# Simulation of anaerobic digestion processes using stochastic algorithm

- Jegathambal Palanichamy
^{1}Email author and - Sundarambal Palani
^{2}

**12**:121

https://doi.org/10.1186/s40201-014-0121-7

© Palanichamy and Palani; licensee BioMed Central Ltd. 2014

**Received: **10 August 2013

**Accepted: **24 August 2014

**Published: **4 September 2014

## Abstract

### Background

The Anaerobic Digestion (AD) processes involve numerous complex biological and chemical reactions occurring simultaneously. Appropriate and efficient models are to be developed for simulation of anaerobic digestion systems. Although several models have been developed, mostly they suffer from lack of knowledge on constants, complexity and weak generalization. The basis of the deterministic approach for modelling the physico and bio-chemical reactions occurring in the AD system is the law of mass action, which gives the simple relationship between the reaction rates and the species concentrations. The assumptions made in the deterministic models are not hold true for the reactions involving chemical species of low concentration. The stochastic behaviour of the physicochemical processes can be modeled at mesoscopic level by application of the stochastic algorithms.

### Method

In this paper a stochastic algorithm (Gillespie Tau Leap Method) developed in MATLAB was applied to predict the concentration of glucose, acids and methane formation at different time intervals. By this the performance of the digester system can be controlled. The processes given by ADM1 (Anaerobic Digestion Model 1) were taken for verification of the model.

### Results

The proposed model was verified by comparing the results of Gillespie's algorithms with the deterministic solution for conversion of glucose into methane through degraders. At higher value of `τ` (timestep), the computational time required for reaching the steady state is more since the number of chosen reactions is less. When the simulation time step is reduced, the results are similar to ODE solver.

### Conclusion

It was concluded that the stochastic algorithm is a suitable approach for the simulation of complex anaerobic digestion processes. The accuracy of the results depends on the optimum selection of tau value.

## Keywords

## Background

Anaerobic Digestion (AD) is the process by which the complex form of organic matter such as carbohydrates, fats and proteins are converted into simpler form by the cells of microorganisms in the absence of oxygen. Energy production, high organic loading and low sludge production are major advantages of AD process. The energy produced can replace fossil fuel use, and also has positive effect on reduction of global warming. Modeling is a powerful tool which can be applied to simulate various processes occurring in the digester. Models are applied for parameter estimation also. Using the simulation results it is easy to predict and avoid digester failure. The modeling results give useful guidelines for the design of the digester also. The reaction system in an anaerobic digester is complex with many sequential and parallel steps. The reactions can be Biochemical or Physico – chemical in nature which involves species of higher and lower concentrations. A stochastic approach can be applied to simulate these reactions in exact manner. The complex organic matter which is called substrate is converted into simpler form through various steps by living cells called biomass. These cells grow at suitable environmental conditions of pH, temperature etc. They interact with the environment and substrates in a complicated way. Generally the biochemical processes include acidogenesis, acteogenesis, and anaerobic oxidation of Volatile Fatty acids, methanogenesis and extracellular hydrolysis step. The reaction kinetics of growth and decay of biomass and conversion of substrates from one form to another are detailed in the following sections.

### Modeling of anaerobic digestion

Various models have been developed based on five major categories: models considering: (a) non-ionised Volatile Fatty Acids and total VFA inhibition; (b) H_{2} as regulator for VFA production; (c) ammonia inhibition. The growth of methanogenic population is greatly affected by un-dissociated acetic acid, un-ionised VFA and total VFA which cause a drop in pH. Several models have been developed taking the substrate (un-dissociated acetate and VFA) as inhibition factor [7],[8]. The factor which regulates the amount of fatty acids generation is the liquid phase redox potential which is expressed as a function of H_{2} partial pressure. Due to sudden increase in organic loading, the accumulation of VFA takes place, since acetogens grow at a slower rate than the acidogens. This will increase pH which in turn the H_{2} partial pressure is increased. This will cause further accumulation of acids and thus methane generation is reduced. Few models have been developed based on the H_{2} partial pressure as inhibition factor [9]–[12] and manure as substrate and generated ammonia as inhibition factor [13]–[15]. Ammonia inhibits the methanogenesis process, thus acetic acid is accumulated. This in turn inhibits acetogenesis process and thus the total VFA accumulates. The reduction in pH causes decrease in ammonia concentration and the inhibition of methanogenesis process. Thus the ammonia inhibition is a self-regulatory.

In all the above models, many species are involved in more than one reaction. The reaction kinetics is solved by formulated Ordinary Differential Equations (ODE). To determine the rate of concentration of a particular species in time, all the reactions in which the species is involved are to be included during the formulation of ODEs. The complexity increases with increase in number of reactions and species. To avoid this, Stochastic Algorithm has been applied where the state of the system is updated based on the current state of the system and the transition probability.

### Kinetics of anaerobic digestion

The complex organic matter which is called substrate is converted into simpler form through various steps by living cells called biomass. These cells grow at suitable environmental conditions of pH, temperature etc. They interact with the environment and substrates in a complicated way. The reaction kinetics of growth and decay of biomass and conversion of substrates from one form to another are detailed in the following sections.

### Kinetics of biomass growth and decay

Once the inoculum is introduced into anaerobic digester, the cells pass through lag phase where they adjust to the new environment. Once they get accustomed to the new environment, they start growing in exponential phase called log phase. During this phase, the specific growth rate remains constant. The specific growth rate is given by dx/dt = μx. When the growth limiting substrate is exhausted, the growth remains stationary and reaches the maintenance mode in the stationary phase. When there is no substrate, the cell population slowly starts decreasing in the death phase. In each growing culture, there is a maximum rate of growth per unit biomass with unlimited substrate in the given environment (μ_{max}).

### Kinetic model of biomass growth

A kinetic model represents the cell population kinetics. The models can be unstructured, structured unsegregated and segregated. In the unstructured model a single substrate is considered as growth limiting one. Multiple substrates are considered in structured model. In unsegregated model, the average properties of the cell population are considered. In segregated model the discrete and heterogeneity of cell populations are considered. Also various kinetic models such as Maltha's Law, Slater model and Monod model are used in modeling AD processes. In Maltha's law the rate of increase of biomass is a function of microbial population only (F(X) = μx). In this model, lag or death phase is not considered. It is assumed that there is unrestricted growth of biomass. In Slater model, final population of biomass is also included. This is represented by logistic equation of biomass growth which relates the specific growth rate μ, biomass concentration X, maximum specific growth rate μ_{max} and final population X_{f} . μ = μ_{max}(1 - X/X_{f}). This is an empirical formula which can be applied in batch studies. Generally, Monod kinetic model is applied in most of the biological wastewater treatment processes. In this model, when one of the substrate concentrations (S) is Limiting, the biomass growth is represented by μ = μ_{max}S/(K_{s} + S). K_{s} is the value of the limiting nutrient concentration at which the specific growth rate is half its maximum value.

### Advantage and limitation of Monod equation

In this model, the kinetic parameters (μ_{max ,}K_{s}) which describe the microbial processes are able to predict the conditions of maximum growth and when the activity will cease. The main disadvantage of this model is that since the kinetic parameters (μ_{max ,}K_{s}) vary with substrate, one set of parameters cannot describe biological process at short and longer retention time. To overcome this limitation, first - order models are used [16].

### Application of Stochastic algorithm for simulation of kinetic reactions in the anaerobic digestor

In this work, four anaerobic microbial groups are considered for degradation: (i) glucose fermenting acidogens; (ii) propionic acid degrading acetogens; (iii) butyric acid degrading actetogens, acetoclatic methanogens and (iv) hydrogenotrohic methanogens.

### Step I: Acidogenesis

After hydrolysis of complex substrate into simpler organic compounds such as glucose, short chain fatty organic matter takes place, acidogens degrade glucose into acetic, propionic and butyric acids. Hydrogen is considered as an inhibitor for acidogens according to Embden- Myerhoof pathway since the correlation of NAD^{+}/NADH in the cells of biomass depend on the concentration of hydrogen.

### Step II: Acetogenesis

In this slowly growing and pH sensitive acteogens oxidize propionic and butyric acids to acetate. At high partial pressure of hydrogen, it acts as an inhibitor in acetogenesis phase. Also acetate inhibition of the propionic and butyric acid degradation step has been considered in numerous studies. So this can be represented by non-competitive type inhibition model.

### Step III - Acetoclastic methanogenic stage

In this step, pH -sensitive and slowly growing acetocalstic methanogens reduce acetate to methane. Here, free ammonia is the inhibitor for the growth of methanogens. The stoichiometric equation and specific growth kinetic reactions are given below:

### Step IV - Hydrogenotrohic methanogenesis

There may be growth limitations due to deficiency of CO_{2} in the reaction system due to digestion of propionic and butyric acids by acetogens. So dual substrate form of monod equation can be applied to represent the specific growth rate of Hydrogenotrohic methanogens.

### Biochemical reactions and their kinetics in the anaerobic digestion system

Biochemical reactions and their kinetics in the anaerobic digestion system were assumed to follow first order reactions (hydrolysis), monod type kinetic reactions and inhibition reaction. The complex particulate waste from industries or household is first disintegrated into carbohydrate, protein and lip (both particulate and soluble inert material). During hydrolysis by extra cellular enzymes (hydrolases), monossaccharides, amino acids and Long Chain Fatty Acids (LCFA) are formed. All these bio-chemical extracellular steps were assumed as first order [6]. The first order kinetic model is an empirical relation, which assumes that the hydrolysis rate is a linear function of the available biodegradable substrate at a certain pH and temperature. The acidogenic bacteria turn the products of hydrolysis into simple organic compounds such as short chain Volatile Acids (VA), e.g. propionic, formic, lactic, and butyric and alcohols such as ethanol, methanol, glycerol and acetone. Then two types of acetogenic mechanism can occur [5] (a) acetogenic hydrogenations and (b) acetogenic dehydrogenations. In acetogenic hydrogenations, the organic acids formed are subsequently converted by acetogenic bacteria to acetate as the main product. Acetogenic dehydrogenations include the anaerobic oxidation of volatile long and short chain fatty acids. In this reaction, acetate is formed from the separated carbon atoms. During this process, due to high hydrogen partial pressure, oxidation process can be inhibited. So the hydrogen produced by these organisms is consumed by a hydrogen-utilizing methanogenic group and acetate by an aceticlastic methanogenic group. Almost the 64–70% of methane production is from acetate. Methanogenic bacteria are very sensitive to pH, temperature, loading rate and other compounds. All the substrate uptake reactions are intracellular reactions and they are modeled using Monod kinetics reaction (single Monod, double Monod and also with competitive and non competitive reactions) Methane is considered to be water insoluble, whereas the carbondioxide is partially soluble and partly escapes into gas phase. When temporary accumulation of Volatile fatty acids occurs, the pH of the digester is reduced. This will increase the concentration of unionized VFA in the system. This will inhibit methanogenic activity. Inhibition function includes pH, hydrogen and free ammonia. Hydrogen and free ammonia inhibition can be represented by non-competitive reaction whereas pH inhibition can be represented by empirical equations. The inorganic nitrogen uptake is represented by competitive secondary Monod-kinetics, where the prevention of growth due to limitation of nitrogen and competition for uptake of butyrate and valerate occur.

### Deterministic and Stochastic approaches for reaction simulation

In the deterministic approach for reaction simulation, the time evolution of the system is considered as continuous which is governed by a set of coupled ODEs. But the stochastic approach regards the time evolution similar to a random walk process which is governed by a single differential equation which is called the master equation. In the deterministic approach, given the initial concentration of various species in a free homogeneous macroscopic environment, the concentration at all future time intervals is determined by averaging random fluctuations to produce a predictable deterministic behaviour [17]. All the elementary reactions (first, second order) follow the reaction rate law where the rate of the reaction is always proportional to the concentration of each reactant involved in the reaction. The ODE solvers propagate the system's state using finite time steps. For non-linear reactions, extremely small time steps need to be adopted to keep the numerical exactness. In that case adaptive step size or implicit method is recommendable. The ODE approach is empirically accurate for reaction systems where large concentrations occur and may not be adequate for systems with small concentrations. In the deterministic approach, the species population is described by a continuous state although a chemical reaction involves random collisions between individual species. Also a predictable system is assumed for the reaction rate or velocity and time evolution. These two assumptions are not appropriate for a system with low concentration of species where high relative fluctuations due to stochastic effect occur. Stochasticity in the state change occurs as intrinsic and extrinsic stochasticity. The intrinsic stochasticity is inherent to the system and mainly arises due to low concentration of species and extrinsic stochasticity arises due to the random variation of environmental factors.

### Concept of master equation

*k*

_{ j }(τ;

*X*(

*t*),

*t*) , that the reaction event

*R*

_{μ}occurs in time interval [t,t + τ ) is given by a Poisson random variable, P(

*a*

_{ j }(

*X*(

*t*)),τ ), where

*j*= 1,..

*M*. The flowchart of Gillespie's Tau Leap Method is outlined in Figure 2. The key point of Gillespie's Tau Leap Method lies in the selection of time step. 'τ' will be approaching to 1/

*a*

_{0}(

*X*(

*t*)) when the number of reactions or the concentration of the reactant species is low. In that case it will be equivalent to Gillespie's Direct Method where only one reaction is chosen in a time step. In these cases, the algorithm becomes inefficient.

### Application of a stochastic algorithm for the simulation of anaerobic digestion processes

_{max},

_{glu}= 1.25, μ

_{max},

_{bu}= 0.833, μ

_{max,pro}= 0.542, μ

_{max,ace}= 0.333, μ

_{max,hyd}= 0.35, K

_{s}= 500, K

_{b}= 200, K

_{p}= 100, K

_{a}= 150, K

_{h}= 150, where μ

_{maxS}are maximum uptake rate of various species, K's are the half saturation value [23]. The yield of product on the substrate (f) and the yield of biomass on substrate (Y) of the various reactions involved in biodegradation are given in the Table 1. Monod reactions are used as rate reactions to determine the next state of the system. The propensity function which is the rate of each reaction a

_{μ}is determined. Then based on the total reaction rate of the system, a reaction is chosen by Poisson random variable. The time step is assumed in such a way that more than one reaction is chosen in a single simulation step. If higher value is chosen, more reactions are chosen where there will be great variation in the propensity function of the chosen reactions. Based on the chosen reactions in a given time step and the stoichiometric constants of products and substrates of the chosen reactions, the number of particles assigned for the reactant and product species involved in the chosen reactions will be modified. So the accuracy of the algorithm depends on the selection of the time step. The simulation results obtained from SBTOOLBOX in Matlab and stochastic algorithm s are shown in Figures 4 and 5.

**Stoichiometric constants for substrates and products in the model**

Sugar/Glucose | Constant | Constant | Constant | Constant | Constant |
---|---|---|---|---|---|

| 10 | 0 | 0 | 0 | 0 |

| 0 | 1 | 0 | 0 | 0 |

| 0 | 0 | 1 | 0 | 0 |

| 0 | 0 | 0 | 1 | 0 |

| 0 | 0 | 0 | 0 | 1 |

| 0 | 0 | 0 | 0 | 0 |

| 0 | 0 | 0 | 0 | 0 |

| 0.117 | 0 | 0 | 0 | 0 |

| 0.243 | 0 | 0 | 0 | 0 |

| 0.369 | 0.752 | 0.5472 | 0 | 0 |

| 0.171 | 0.188 | 0.4128 | 0 | 0 |

| 0 | 0 | 0 | 0.95 | 0.94 |

When simulation time step is more than 1, the number of selections of the chosen reaction is insignificant and the formation of methane is less. Update of sugar takes longer time to reach the steady state. The simulation of formation of acids and formation of methane are rapid. When the tau value is reduced, the simulation results approaches the results of simulation using ODE solver. This is due to the selection of reactions where the propensity function remains constant. At the time step of 0.4, the simulation results get closer to deterministic values. Thus the accuracy of simulation depends on the selection of the optimal value of tau.

## Conclusion

When more number of reactions are involved in a biochemical processes, the rate of change of concentration of species is determined by evolving a rate equation for each species. As the number of species increases, the number of rate of reactions also increases to evolve the change in the concentration of species. But it is a complex process. The difficulty in formulation of Ordinary Differential Equation (ODE) for reaction simulation can be overcome using Gillespie's algorithms where the reactions can be represented in the form of matrix which involves the stoichiometric constants and the reaction rate constants. Stochastic Algorithm (SA) has been identified as a better solution where multiscale concentration species are involved in reactions such as sorption, precipitation, degradation, sequential and parallel decay reactions etc. In this paper, the feasibility of applying SA (Gillespie Tauleap Method) in simulation of various bio chemical processes occurring during anaerobic degradation of complex organic matter has been studied.

## Declarations

### Acknowledgements

The authors acknowledge Water Institute, Karunya University for providing laboratory facilities and support.

## Authors’ Affiliations

## References

- Andrews JF:
**A dynamic model of the anaerobic digestion process.***J Sanitary Eng Div Proc Am Soc Civil Eng SA*1969,**1:**95–116.Google Scholar - Andrews JF, Graef SP: Dynamic modelling and simulation of the anaerobic digestion process. In Anaerobic Biological Treatment Processes, Advances in Chemistry Series. 105th edition. Edited by Gould RF. Washington DC: 1971:126-162. Amer. Chem. Soc, doi:10.1021/ba-1971–0105.ch008.Google Scholar
- Olafadehan OA, Alabi AT:
**Modelling and simulation of methanogenic phase of an anaerobic digester.***J Eng Res*2009,**13**(2):1–16.Google Scholar - Kalyuzhnyi SV:
**Batch anaerobic digestion of glucose and its mathematical modeling. II: description, verification and application of model.***Biore Technol*1997,**59:**249–258. 10.1016/S0960-8524(96)00125-3View ArticleGoogle Scholar - Gavala HN, Angelidaki I, Ahring BK:
**Kinetics and modeling of anaerobic digestion process.**In*Biomethanation I*. 81st edition. Edited by: Ahring BK. Berlin: Springer, Berlin; 2003:57–93. Available from: 10.1007/3–540–45839–5_3 10.1007/3-540-45839-5_3View ArticleGoogle Scholar - Batstone DJ, Keller J, Angelidaki I, Kalyuzhnyi SV, Pavlostathis SG, Rozzi A, Sanders WT, Siegrist H, Vavilin VA:
**The IWA Anaerobic Digestion Model No 1(ADM 1).***Water Sci Tech*2002,**45**(10):65–73.Google Scholar - Moletta R, Verrier D, Albagnac G:
**Dynamic modelling of anaerobic digestion.***Wat Res*1986,**20:**427–434. 10.1016/0043-1354(86)90189-2View ArticleGoogle Scholar - Hill DT, Barth CL:
**A dynamic model for simulation of animal waste digestion.***J Water Pollut Contr Fed*1977,**10:**2129–2143.Google Scholar - Mosey FE:
**Mathematical modelling of the anaerobic digestion process: regulatory mechanisms for the formation of short-chain volatile acids from glucose.***Wat Sci Technol*1983,**15:**209–232.Google Scholar - Costello DJ, Greenfield PF, Lee PL:
**Dynamic modelling of a single- stage high-rate anaerobic reactor - I: model derivation.***Wat Res*1991,**25:**847–858. 10.1016/0043-1354(91)90166-NView ArticleGoogle Scholar - Costello DJ, Greenfield PF, Lee PL:
**Dynamic modelling of a single-stage high-rate anaerobic reactor - II: model verification.***Wat Res*1991,**25:**859–871. 10.1016/0043-1354(91)90167-OView ArticleGoogle Scholar - Pullammanapallil P, Owens JM, Svoronos SA, Lyberatos G, Chynoweth DP:
**Dynamic model for conventionally mixed anaerobic digestion reactors.***AIChE Annu Meet*1991,**paper 277c:**43–53.Google Scholar - Siegrist H, Renggli D, Gujer W:
**Mathematical modelling of anaerobic mesophilic sewage sludge treatment.***Wat Sci Technol*1993,**27:**25–36. 10.1021/es00038a601View ArticleGoogle Scholar - Angelidaki I:
*Anaerobic thermophilic biogass process: the effect of lipids and ammonia*. The Technical University of Denmark, Copenhagen; 1992.Google Scholar - Angelidaki I, Ellegaard L, Ahring BK:
**A mathematical model for dynamic simulation of anaerobic digestion of complex substrates: focusing on ammonia inhibition.***Biotechnol Bioeng*1993,**42:**159–166. 10.1002/bit.260420203View ArticleGoogle Scholar - Mussati M, Aguirre P, Fuentes M, Scenna N:
**Aspects on methanogenic biofilm reactor modeling.***Latin Am Appl Res*2006,**36:**173–180.Google Scholar - Turner TE, Schnell S, Burrage K:
**Stochastic approaches for modeling in vivo reactions.***Comput Biol Chem*2004,**24:**165–178. 10.1016/j.compbiolchem.2004.05.001View ArticleGoogle Scholar - Gillespie D:
**Concerning the validity of the stochastic approach to chemical kinetics.***J Stat Phy*1977,**16:**311–318. 10.1007/BF01020385View ArticleGoogle Scholar - Gillespie D:
**Approximate accelerated stochastic simulation of chemically reacting systems.***J Chem Phy*2001,**115**(4):1716–1733. 10.1063/1.1378322View ArticleGoogle Scholar - Cao Y, Gillespie D, Petzold L:
**Multiscale stochastic simulation algorithm with stochastic partial equilibrium assumption for chemically reacting systems.***J Comput Phy*2005,**2:**395–411. 10.1016/j.jcp.2004.12.014View ArticleGoogle Scholar - Rathinam M, Petzold L, Cao Y, Gillespie D:
**Stiffness in stochastically reacting systems: the implicit tau-leaping method.***J Chem Phy*2003,**119**(24):12784–12794. 10.1063/1.1627296View ArticleGoogle Scholar - Jegathambal P, Becker T, Spiller M, Kongeter J, Mohan S:
**Multicomponent reaction modeling using a stochastic algorithm.***J Comput Visual Sci*2009,**12**(2):51–61. 10.1007/s00791-007-0080-yView ArticleGoogle Scholar - Jeong HS, Suh CW, Lim JL, Lee SH, Shin HS:
**Analysis and application of ADM1 for anaerobic methane production.***Bioproc Biosyst Eng*2005,**27**(2):81–89. 10.1007/s00449-004-0370-4View ArticleGoogle Scholar

## Copyright

This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly credited. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.