This article was downloaded by: [University of Chicago Library] On: 26 December 2014, At: 14:05 Publisher: Taylor & Francis Informa Ltd Registered in England and Wales Registered Number: 1072954 Registered office: Mortimer House, 37-41 Mortimer Street, London W1T 3JH, UK

Environmental Technology Publication details, including instructions for authors and subscription information: http://www.tandfonline.com/loi/tent20

Estimation of parameter uncertainty for an activated sludge model using Bayesian inference: a comparison with the frequentist method ab

Živko J. Zonta , Xavier Flotats a

ab

ac

& Albert Magrí

GIRO Technological Centre, Mollet del Vallès, Barcelona, Spain

b

GIRO Joint Research Unit IRTA-UPC, Department of Agrifood Engineering and Biotechnology, Universitat Politècnica de Catalunya, UPC-BarcelonaTECH, Parc Mediterrani de la Tecnologia, Castelldefels, Barcelona, Spain c

IRSTEA, UR GERE, Rennes, France Accepted author version posted online: 23 Dec 2013.Published online: 22 Jan 2014.

Click for updates To cite this article: Živko J. Zonta, Xavier Flotats & Albert Magrí (2014) Estimation of parameter uncertainty for an activated sludge model using Bayesian inference: a comparison with the frequentist method, Environmental Technology, 35:13, 1618-1629, DOI: 10.1080/09593330.2013.876450 To link to this article: http://dx.doi.org/10.1080/09593330.2013.876450

PLEASE SCROLL DOWN FOR ARTICLE Taylor & Francis makes every effort to ensure the accuracy of all the information (the “Content”) contained in the publications on our platform. However, Taylor & Francis, our agents, and our licensors make no representations or warranties whatsoever as to the accuracy, completeness, or suitability for any purpose of the Content. Any opinions and views expressed in this publication are the opinions and views of the authors, and are not the views of or endorsed by Taylor & Francis. The accuracy of the Content should not be relied upon and should be independently verified with primary sources of information. Taylor and Francis shall not be liable for any losses, actions, claims, proceedings, demands, costs, expenses, damages, and other liabilities whatsoever or howsoever caused arising directly or indirectly in connection with, in relation to or arising out of the use of the Content. This article may be used for research, teaching, and private study purposes. Any substantial or systematic reproduction, redistribution, reselling, loan, sub-licensing, systematic supply, or distribution in any form to anyone is expressly forbidden. Terms & Conditions of access and use can be found at http:// www.tandfonline.com/page/terms-and-conditions

Environmental Technology, 2014 Vol. 35, No. 13, 1618–1629, http://dx.doi.org/10.1080/09593330.2013.876450

Estimation of parameter uncertainty for an activated sludge model using Bayesian inference: a comparison with the frequentist method Živko J. Zontaa,b , Xavier Flotatsa,b and Albert Magría,c∗ Technological Centre, Mollet del Vallès, Barcelona, Spain; b GIRO Joint Research Unit IRTA-UPC, Department of Agrifood Engineering and Biotechnology, Universitat Politècnica de Catalunya, UPC-BarcelonaTECH, Parc Mediterrani de la Tecnologia, Castelldefels, Barcelona, Spain; c IRSTEA, UR GERE, Rennes, France

a GIRO

Downloaded by [University of Chicago Library] at 14:05 26 December 2014

(Received 8 August 2013; accepted 12 December 2013 ) The procedure commonly used for the assessment of the parameters included in activated sludge models (ASMs) relies on the estimation of their optimal value within a confidence region (i.e. frequentist inference). Once optimal values are estimated, parameter uncertainty is computed through the covariance matrix. However, alternative approaches based on the consideration of the model parameters as probability distributions (i.e. Bayesian inference), may be of interest. The aim of this work is to apply (and compare) both Bayesian and frequentist inference methods when assessing uncertainty for an ASM-type model, which considers intracellular storage and biomass growth, simultaneously. Practical identifiability was addressed exclusively considering respirometric profiles based on the oxygen uptake rate and with the aid of probabilistic global sensitivity analysis. Parameter uncertainty was thus estimated according to both the Bayesian and frequentist inferential procedures. Results were compared in order to evidence the strengths and weaknesses of both approaches. Since it was demonstrated that Bayesian inference could be reduced to a frequentist approach under particular hypotheses, the former can be considered as a more generalist methodology. Hence, the use of Bayesian inference is encouraged for tackling inferential issues in ASM environments. Keywords: activated sludge model; parameter uncertainty; Markov chain Monte Carlo; Bayesian; Hessian

Nomenclature Activated sludge model: parameters bH endogenous decay coefficient of biomass, 1/d bSTO endogenous decay coefficient of storage products, 1/d fXI production of inert organic materials in endogenous respiration, mg COD-XI /mg COD-XH K1 regulation constant of biomass as function of the storage products-to-biomass ratio, mg COD-XSTO /mg COD-XH K2 lumped parameter related to the affinity of biomass towards the storage products-to-biomass ratio (mg COD-XSTO /mg COD-XH )2 KS substrate affinity constant, mg COD-SS /L kSTO maximum storage rate of biomass, 1/d YH,S yield coefficient of biomass for growth on substrate, mg COD-XH /mg COD-SS YH,STO yield coefficient of biomass for growth on storage products, mg COD-XH /mg COD-XSTO YSTO yield coefficient of biomass for storage on substrate, mg COD-XSTO /mg COD-SS δ efficiency of the oxidative phosphorylation, mol ATP/mol NADH2

∗ Corresponding

author. Email: [email protected]

© 2014 Taylor & Francis

μMAX,S μMAX,STO τ

maximum growth rate of biomass on substrate, 1/d maximum growth rate of biomass on storage products, 1/d first-order time constant, min

Activated sludge model: variables SS readily biodegradable organic substrate, mg COD/L t time, min XH heterotrophic biomass, mg COD/L XI inert, non-biodegradable organics, mg COD/L XSTO intracellular storage products, mg COD/L Inference notation C covariance matrix D data vector E expectation operator f (x) function F(X ) cumulative distribution function H Hessian matrix K −S Kolmogorov–Smirnov statistic

Downloaded by [University of Chicago Library] at 14:05 26 December 2014

Environmental Technology L(θ ; D) Norm p p(D|θ ) p(θ ) p(θ |D) Si Si,j STi SS U (ai , bi ) V x X Y χ2 ε θ σ2

likelihood function Normal distribution probability likelihood prior posterior main effect index first-order interaction effect index total effect index sum-of-squares uniform distribution variance variable stochastic variable model chi-square statistic error factor/parameter vector error variance parameter

Subscripts/superscripts true value k data size (1, . . ., N ) m mode ∼p non-problematic p problematic s sample size (1, . . ., S) θ∼i all factors but θi θi i-th factor (1, . . . , P) ∗

Abbreviations ATP: adenosine triphosphate COD: chemical oxygen demand, mg O2 /L NADH2 : reduced nicotinamide adenine dinucleotide VFA: volatile fatty acid

1. Introduction Respirometry is a fast and relatively inexpensive technique applied in the field of biological wastewater treatment for the characterization of activated sludge. The technique is based on the measurement of the oxygen uptake rate (OUR) due to microbial activity and its use is common for the compositional typification of wastewaters and the calibration of biokinetic models.[1–7] Working with activated sludge models (ASMs), the assessment of the parameter accuracy (or uncertainty) within a frequentist confidence region is quite usual. [2–9] Briefly, the frequentist procedure is random while a parameter is assumed as a fixed unknown quantity. Once the optimal parameter value is estimated, its uncertainty is assessed from the covariance matrix, which can be estimated from a local linearization in parameters of the model outcome (i.e. involving the Fisher information matrix (FIM)) or from a quadratic expansion of the goodness-of-fit function (i.e. through the Hessian).

1619

While frequentist methods assume parameters as single unknown values, the Bayesian approach treat them as random variables defined by a probability distribution. In the Bayesian framework, the initial degree of belief over the parameter uncertainty is expressed by a prior distribution (before experimental data are collected), while the likelihood distribution indicates how likely it would be to observe data, given a particular parameter value. The updated degree of belief over the parameter is represented by the posterior distribution, which is an aggregation of the information contained in the prior and the likelihood distributions. If the prior is non-informative (i.e. flat and uniform distribution) and the model structure is fixed, then the relative posterior may be assumed as only dependent on the available data through the likelihood. In this case, the posterior mode is known as the maximum likelihood estimator (MLE). In case of considering a linear model in parameters and a large enough data set, the frequentist and Bayesian procedures will yield the same result.[10] From a Bayesian perspective this fact justifies the common scientific practice of interpreting frequentist inference (point estimates and confidence regions) in a Bayesian fashion (i.e. probability statements about confidence regions). However, for nonlinear models (including ASM-type models) used with sparse data the frequentist approach may underestimate parameter uncertainty.[11,12] Since respirometric data used for calibration are often sparse in relation to model complexity, their use do not guarantee the practical identifiability of all the parameters included in an ASM-type model.[13] Once the evidence provided by data is established, the onset of nonidentifiability may be considered as gradual because there are parameters: (i) informative only for a given model regime but data do not force the model in this regime, and/or (ii) informative in groups that cannot be resolved into individual components (i.e. correlated parameters).[14] In this regard, the higher the practical identifiability of a given parameter, the higher the accuracy of its estimate. The covariance matrix (which is estimated by the inversion of the FIM or the Hessian) fully describes the practical identifiability of parameters when dealing with linear models. Thus, if it contains zero eigenvalues some model parameters are non-identifiable.[14] Alternatively, if the FIM or the Hessian matrices are badly conditioned (i.e. high ratio between the highest and lowest eigenvalue of the matrix), it may become interesting to re-think the model and re-run the numerical analysis or, in other cases, to get more data.[13] Thus, complementary data to the OUR profiles can be considered in order to enhance the inferential procedure. This is the case for compositional variables (apart from dissolved oxygen (DO)), storage products or the pH (titrimetry).[15,16] In order to simplify the analysis, rather than re-thinking the model or collecting more data, the problematic (non-identifiable) parameters are sometimes fixed at a nominal value before re-running the uncertainty estimation.

Downloaded by [University of Chicago Library] at 14:05 26 December 2014

1620

Ž.J. Zonta et al.

In the wastewater treatment community, local trajectory sensitivity analysis (i.e. based on the use of the FIM) is commonly used to assess practical identifiability of the model parameters [13–17] although for nonlinear models nearzero eigenvalues of the covariance matrix remain only as indicative of non-identifiability.[14] In this context, a local derivative-based approach may be inconclusive.[18] Alternatively to the local trajectory sensitivity analysis, global sensitivity analysis (GSA) can be applied.[19] Then, when the GSA is performed over the goodness-of-fit function,[20] it is possible to assess the parameters (alone or under interactions) leading the model fitting to the experimental data. The frequentist and Bayesian inference methods were already compared by Omlin and Reichert using simple synthetic data sets of a Monod model for microbial substrate conversion.[11] They concluded that the frequentist technique is superior in the case of identifiable parameters (computational efficiency). However, in the case of poor parameter identifiability, the conceptual advantage of the estimation of parameter distributions and the use of prior knowledge make the Bayesian approach more recommendable. The aim of this paper is to apply (and compare) both Bayesian and frequentist inference methods to an ASM which considers microbial substrate conversion under aerobic conditions according to simultaneous storage and growth processes.[2] This model has a more complex structure than the Monod model, including parameters which cannot be perfectly classified as identifiable or non-identifiable. The inference exercise is carried out exclusively considering OUR profiles (three data sets) and with the aid of GSA. 2.

Materials and methods

2.1. Global sensitivity analysis Sensitivity analysis deals with the study of how uncertainty in the model output can be apportioned to different sources of uncertainty in the model input factors.[21] Particularly, GSA focuses on the output uncertainty over the entire range of possible values of the input factors. When considering nonlinear models, the variance-based methods are well suited to account for interactions between factors. In this regard, given a model of the form Y = f (θ ), with Y a scalar, a variance-based main effect for a generic factor θi (i = 1, . . ., P) can be written as [21] Vθi (Eθ∼i {Y |θi }),

(1)

where θi is the i-th factor and θ∼i denotes the vector of all factors but θi . The meaning of the inner expectation operator, E, is that the mean of Y is taken over all possible values of θ∼i while keeping θi fixed. The outer variance, V , is taken over all possible values of θi . When Equation (1) is normalized by the unconditional total variance

V (Y ) = Vθi (Eθ∼i {Y |θi }) + Eθi (Vθ∼i {Y |θi }), the associated sensitivity measure (main effect index) is obtained, which can be written as Si =

Vθi (Eθ∼i {Y |θi }) . V (Y )

(2)

In a similar way, the first-order interaction effect index can be written as   Vθi,j (Eθ∼i,j Y |θi,j ) Si,j = . (3) V (Y ) Another important sensitivity measure is the total effect index, defined as STi =

Eθ∼i (Vθi {Y |θ∼i }) Vθ (Eθ {Y |θ∼i }) = 1 − ∼i i , V (Y ) V (Y )

(4)

which measures the main, first- and higher- order interaction effects of factor θi . In probabilistic GSA, the factor θ is a stochastic variable characterized by a distribution g(θ ) that describes the prior knowledge over θ . In this case, the input factor θ is the model parameter to be inferred while the function f (θ ) is the sum-of-squares defined in the next section ‘Bayesian Inference’. Intuitively, an individual parameter θi can be estimated with precision if it leads the model fitting to the available data (i.e. Si is high). Those parameters involved in interactions (i.e. STi >> Si ), but which cannot be estimated precisely (i.e. Si is low), should be still considered as important to get a good fitting. 2.2. Bayesian inference In the Bayesian framework, the concept of probability, p, is defined as the degree of belief or the plausibility that a proposition is true and is quantified as a real positive number in the range [0, 1]. For the determination of the probability of a continuous parameter, θ , once given data D, with the k-th element Dk (k = 1, . . ., N ), and considering the prior information on the parameter θ, inference concerning θ is based on its posterior distribution [10] p(θ |D) ∝ p(θ )p(D|θ),

(5)

which depends both on the ‘subjective’ belief over the parameter θ through the prior p(θ ) and on data D through the likelihood p(D|θ). The prior p(θ ) is formulated ‘before’ data are observed while the likelihood p(D|θ) is a conditional probability depending on θ , with D held fixed. Note that commonly the likelihood p(D|θ) = L(θ ; D) is seen as a mathematical function of θ , with parameters D. At this point, the uncertainty in the future model predictions y = f (θ ) can be inferred as   p(y|D) = p(y, θ |D)dθ = p(y|θ)p(θ|D)dθ, (6) where y and D are assumed to be conditionally independent given the value of θ .

Environmental Technology A typically used measurement model for the observations D is Dk = fk (θ ) + ε, ε ∼ Norm(0, σ 2 ),

(7)

where the error term ε is an independent and identically distributed (IID) random variable with a homoscedastic error variance parameter, σ 2 . Once at that point, the corresponding likelihood for this measurement model takes the form   1 p(D|θ , σ 2 ) ∝ exp − 2 SS(θ ; D) , 2σ

Downloaded by [University of Chicago Library] at 14:05 26 December 2014

N 

(Dk − fk (θ ))2

the most probable value of θ (i.e. the MLE) is the leastsquare (LS) estimate, θm . Moreover, if f (θ ) is a linear model in parameters the shape of χ 2 is parabolic, and thus the posterior is multi-variate Gaussian (MVG). In this context (since the posterior is completely defined by the covariance matrix, C), is possible to estimate the uncertainty of the LSestimate θm from the Hessian matrix, H,[24] according to χ 2 (θm ) (13) 2H(θm )−1 . N −P Summarizing, the frequentist uncertainty estimation procedure for θ was derived from a more general framework (i.e. Bayesian framework), under the following clearly stated hypotheses: C=

(8)

where SS(θ ; D) =

(9)

k=1

is the sum-of-squares. For a general likelihood function, SS(θ ; D) corresponds to twice the negative log-likelihood, 2ln(p(D|θ )).[22] In practice, the measurement model which is assumed should be tested against the relative distribution of residuals. Finally, in Equation (5) the posterior p(θ |D) is only referred as proportional to the product of a given prior and likelihood. The reason is that, for the commonly used numerical methods, the posterior is not provided in a closedform solution. For example, those algorithms based on Markov chain Monte Carlo (MCMC) produce correlated samples {θ (s) , s = 1, . . ., S} (i.e. a Markov chain) which have in the posterior p(θ|D) equilibrium distribution. Based on the sample θ (s) , the posterior is ‘recovered’ within a histogram or kernel-density approximation. Thus, the quality of such approximation will improve with the number of samples from the posterior.[23] 2.3.

Frequentist inference

Under specific assumptions, the frequentist approach is a particular case of the Bayesian inference. Assuming that a uniform distribution is a practical choice for the prior, Equation (5) becomes [10] p(θ |D) ∝ p(D|θ ) = L(θ; D).

(10)

The set of θ which is most likely that which maximizes L(θ; D); a result that fits with the maximum likelihood principle. Thus, if the error ε is assumed as a normal distribution according to Equation (7), the likelihood-dominated result becomes  p(θ |D) ∝ exp − 12 χN2 , (11) where χ 2 = χN2 =

N 1  (Dk − fk (θ ))2 σ2

1621

(1) The prior probability p(θ ) is uniform. (2) The measurement error term ε is Norm(0, σ 2 IN ), with N -by-N identity matrix IN . (3) The LS-estimate θm is asymptotically θ ∗ (the true value of θ ). (4) f (θ ) is a linear model in parameters (at least locally around θm ). In routine applications involving nonlinear models in parameters, the linear approximation of the uncertainty of θ may still fit because the above hypotheses are only slightly violated or because the posterior exhibits an approximate MVG shape. 2.4.

ASM structure

The applied model was first proposed by Sin et al. [2] as a modified version of the ASM3.[25] Basically, it assumes that the growth process for the heterotrophic biomass (under aerobic conditions) is simultaneous to the intracellular storage process. The model also includes a time-depending term (i.e. concerning processes of the formation of intracellular storage products and aerobic growth on external substrate) to simulate the transient response of the observed OUR.[26] The growth of biomass on XSTO is assumed to occur under strictly famine conditions.[2] Thus, the corresponding outcome in terms of OUR is given by Equation (14) 

SS 1 − YSTO · (1 − e−t/τ ) · kSTO · · XH OUR = YSTO KS + S S 

1 − YH,S SS · (1 − e−t/τ ) · μMAX,S · + · XH YH,S KS + SS 

1 − YH,STO KS · μMAX,STO · + YH,STO KS + S S

(12)

k=1

is the chi-square statistic with N degrees of freedom. Maximizing the likelihood is equivalent to minimizing χ 2 , and

·

(XSTO /XH )2 · XH K2 + K1 · (XSTO /XH ) · XH

+ (1 − fXI ) · bH · XH + bSTO · XSTO ,

(14)

1622

Ž.J. Zonta et al.

According to the assumptions made by Sin et al.,[2] the three yield coefficients (YSTO , YH,S , and YH,STO ) were considered to depend exclusively on the efficiency of the oxidative phosphorylation (δ) 4.2 4δ − 2 4δ − 2 4.5 · , YSTO = · , and 4.2δ + 4.32 4 4.5δ 4 4.5δ − 0.5 4.2 YH,STO = · , (15) 4.2δ + 4.32 4.5

YH,S =

Downloaded by [University of Chicago Library] at 14:05 26 December 2014

and the model parameters μMAX,S and μMAX,STO were assumed to be equals. Furthermore, it was considered that the concentrations of DO and ammonium were high enough not to affect kinetics.

2.5.

Model implementation and computational analyses The model was implemented in Matlab (Mathworks Inc., USA) as a Simulink S-function Cmex code block. A free-Matlab toolbox called ‘Adaptive Robust Numerical Differentiation’ [27] was used for the estimation of H. The estimation routine is based on a finite-difference, fourthorder Romberg-extrapolation method with an adaptive routine for the determination of the step-size perturbation parameters. Strictly positive uniform priors were assumed in order to preserve the physical meaning of the parameters. The measurement model for the experimental data was defined according to Equation (7). The error variance σ 2 was considered as an unknown stochastic parameter, which was characterized by the inverse gamma distribution prior [22,23] with mean of 0.01 and an accuracy of 4. It was verified that the computed posterior was insensible to the relative choice of the gamma prior parameter. The posterior p(θ |D) was approximated with a sample size S of 15,000 (burn-in sample of 5000) and obtained after the convergence of the delayed rejection adaptive metropolis (DRAM) sampler.[22–28] The Markov chain was considered stationary when the value of the Geweke’s convergence diagnostic [29] was higher than 0.9. For a one-dimensional case, the K − S statistic is defined as [30] K − S = supx |FS1 (x1 ) − FS2 (x1 )|,

(16)

where K − S is used to quantify the distance between the empirical and theoretical cumulative distribution functions FS1 and FS2 with S IID observations from the stochastic variable X1 . In a two-dimensional case (with a new stochastic variable X2 ), the Peacock’s algorithm considers the four quadrants (x1 < X1 , x2 < X2 ), (x1 < X1 , x2 > X2 ), (x1 > X1 , x2 < X2 ) and (x1 > X1 , x2 > X2 ), adopting as the final K − S statistic the largest of the four differences between the two empirical cumulative distributions.

2.6. Data sources The experimental data, which were used for the numerical analysis, consisted of three OUR profiles (data sets A–C). In all cases, the respirometric assay was carried out under low substrate-to-biomass ratio (SS (0)/XH (0) = 0.05 − 0.2) and with feast plateaus lasting about 15 min. The feast periods after substrate dosage were followed by famine periods; that is when the OUR suddenly drops from its maximum level to a level that still is higher than the endogenous OUR, and thereafter gradually decreases to the endogenous level. (1) Data set A was provided by the authors Sin et al. [2] It consisted of two pulses of acetate, each one equivalent to 40 mg COD/L, which were added according to an optimal experimental design to activate sludge collected from a municipal wastewater treatment plant (WWTP) aiming at nitrogen removal. XH (0) was calculated as 800 mg COD/L and XSTO (0) was measured as 6.8 mg COD/L. (2) Data set B was digitalized from Hoque et al. [5]. It consisted of a single pulse of acetate equivalent to 50 mg COD/L (pH 7.8, 20◦ C), which was added to the activated sludge collected from a WWTP. In this case, XH (0) was calculated as 900 mg COD/L. (3) Data set C was obtained in the laboratory and consisted of a single pulse of acetate equivalent to 40 mg COD/L (pH 8.0, 20◦ C), which was added to the activated sludge purged from a lab-scale sequencing batch reactor (SBR). The SBR was fed with raw leachate under a schedule based on the dosage of 4 pulses per day, a loading rate of 1 g COD/L-d and intermittent aeration. The composition of the leachate was equivalent to 9.81 g CODVFA /L (25% acetate, 9% propionate, 52% butyrate and 14% valerate), 48% CODVFA /COD and 1.01 g N/L. The respirometric test was carried out in a 2.5-L liquid phase principle, flowing gas, static liquid (LFS) respirometer.[31] The OUR values were estimated offline from the DO measurements (Inolab 740 – CellOx 325, WTW, Germany) by applying an optimal local polynomial filtration paradigm called ‘Lazy Learning’.[32] The initial content of storage products in biomass, XSTO (0), was assumed as 7.6 mg COD/L. Similar to the previous works (data sets A and B), once fixed the parameters bSTO (0.2), bH (0.2), and fXI (0.2) according to the values given in the ASM3, the initial concentration of biomass, XH (0), was calculated from the endogenous OUR as 214 mg COD/L. The ASM used in this research was first described by Sin et al. [2] (data set A) and later by Hoque et al. [5] (data set B). In both researches, parameter estimation was carried out exclusively considering OUR measurements. Although such approximation may be constrained by the lack of

Downloaded by [University of Chicago Library] at 14:05 26 December 2014

Environmental Technology

Figure 1. Schema of the uncertainty analysis procedure followed in this research including Bayesian and frequentist inferences.

data,[33] the same approach is considered here since the main target of this work remains the application of Bayesian and frequentist inferential procedures using the same data. For all the aforementioned respirometric tests, a sub-group of measurements was considered for the inference (from 1.4 to 2.8 min/sample). Final length of data sets A, B and C was 70, 49 and 72 measurements, respectively. 2.7.

Procedure to assess and compare parameter uncertainty The parameter vector to be assessed was defined as θ = [τ , KS , kSTO ,μMAX,S , δ, K1 , K2 ]. Moreover, θ will be split in two sub-clusters during the analysis according to θ = [θ ∼p , θ p ], where θ ∼p includes ‘non-problematic’ model parameters and θ p includes potentially ‘problematic’ model parameters (non-identifiable). Once at this point, the uncertainty analysis procedure was applied as follows (Figure 1): Step 1: Assessment of the influence of the model parameters in the outcome by GSA. The evaluation was performed on the basis of the analysis of the three OUR profiles by means of Bayesian GSA [34] using a free software tool,[35] which is able to compute the main effect index (Si ), the first-order effect index (Si,j ) and the total effect index (STi ). The distribution g(θi ) was assumed uniform – i.e. U (ai , bi ) – for all the parameters. Particularly, τ ∼ U (0.1, 5), KS ∼ U (0.1, 10), kSTO ∼ U (0.1, 12), μMAX,S ∼ U (1, 20), δ ∼ U (1, 8), K1 ∼ U (10−2 , 1) andK2 ∼ U (10−4 , 10−2 ). Limits of the uniform distribution were chosen wide enough to include the parameter values reported in the literature.[2,5] The parameter sub-clusters θ ∼p and θ p will be defined after the analysis. Step 2: Calculation of θm by means of Bayesian inference. A sample from the posterior p(θ|D) was obtained after running the DRAM routine and θm was calculated from the kernel-density approximation of the MCMC sample. Residuals were subsequently analysed by means of histograms. Auto-correlation plots [36] were performed in order to

1623

check the appropriateness of the measurement model (i.e. including the error ε as a normal distribution). ∼p,∗ Step 3: Calculation of θm by means of Bayesian inferp p ence (once fixed θm ). The subset θ p was fixed to θm and p the posterior p(θ ∼p |D, θm ) was approximated through the ∼p,∗ DRAM routine. The mode estimate θm of the posterior was subsequently calculated. ∼p ∼p,∗ Step 4: Calculation of θm and θm by means of frequentist inference. The covariance matrix C for the ∼p parameter point-estimation θm and the covariance matrix ∼p,∗ C∗ relative to θm were estimated according to Equation ∼p (13), while θ p was fixed at θm . The MVG approximations ∼p ∼p Norm(θm , C) and Norm(θm , C∗ ) were computed. Step 5: Contrast of Bayesian inference vs. frequentist inference (K − S statistic). The posterior distribution ∼p p(θ ∼p |D) was compared with its approximation Norm(θm , C). This case was named ‘full-case’. Alternatively, the p posterior p(θ ∼p |D, θm ) was compared with its approxi∼p mation Norm(θm , C∗ ). This case was named ‘reducedcase’. Differences between posteriors and relative approximations were evaluated for both case studies using the two-dimensional K − S statistic.

3. Results and discussion 3.1. Global sensitivity analysis The main effect index (Si ) and the total effect index (STi ) for the three assessed data sets A–C are presented in Table 1. The parameters δ and μMAX,S lead the model fitting to the experimental data, since they have the highest indexes Si (δ: 53.25–58.75%; μMAX,S : 4.78–14.56%) and STi (δ: 79.82–89.91%; μMAX,S : 29.32–35.71%). The Table 1. GSA indexes for the parameters of the model. Si : main effect index, STi : total effect index (units: %). Data set Aa Parameter

Si

STi

Data set B

Data set C

Si

Si

STi

STi

τ 2.97 6.0 3.25 8.96 0.38 2.66 KS 1.24 4.0 1.23 3.45 1.21 8.12 kSTO 0.01 0.01 0.08 0.08 0.06 0.70 μMAX,S 14.56 32.83 13.42 29.32 4.78 35.71 δ 58.44 79.23 58.75 79.82 53.25 89.91 K1 0.09 2.80 0.10 4.10 0.23 7.05 K2 0.05 2.87 0.06 3.65 0.28 5.98 Total main 77.4 76.9 60.2 effectb [a] Total first-order 19.0 18.2 32.1 interaction effectc [b] Sum [a+b] 96.3 95.1 92.3 a Data

set A: Sin et al. [2]; data set B: Hoque et al. [5]; data set C: own data. b Total main effect: S . i c Total first-order interaction effect: S . Individual first-order i,j interaction effect indexes are not shown.

Downloaded by [University of Chicago Library] at 14:05 26 December 2014

1624

Ž.J. Zonta et al.

majority of the interaction effects in which these parameters are involved is attributed to first-order interactions (i.e. Si,j indexes are high, with total values of 18.2–32.1%). Since these two parameters have relatively high main effects, both are expected to be well estimated during the subsequent inference exercise. On the other hand, other parameters with lower indexes are assumed to be assessed with less accuracy. According to this analysis, and also taking into account the results reported by Sin et al. [2] and Hoque et al. [5] (who found high correlation between parameters K1 and K2 ), the parameter θ was split in θ ∼p = [τ , KS , kSTO , μMAX,S , δ] and θ p = [K1 , K2 ] for the subsequent inferential exercise. In accordance with this assumption, the condition number [37] for H relative to θ ∼p was low, and thus the estimation of the covariance matrix could be considered as ‘non-problematic’ still including the parameter kSTO . In this research, GSA was introduced in order to provide an insight into what to expect from the inferential

procedures, and the analysis was performed only once. However, the methodology applied here could be still enhanced by repeating the sensitivity analysis and the parameter inference iteratively. Moreover, the use of informative priors, or additional experimental data (i.e. concentration of storage products, pH, etc.) could also help in the analysis. 3.2. Bayesian inference The marginal posteriors of θ for data sets A-C (Figure 2) are depicted in Figure 3. All the relative modes θm (95% credibility) for data sets A and B are similar to those values initially reported by the respective authors (Table 2). The only exemption is the first-order time constant τ for data set A. Here, the mode is 1.67 ± 0.36 min while the value estimated by the original authors was 0.51 ± 0.07 min. Such difference may be caused by the pre-processing

Figure 2. OUR profiles (data sets A–C) and predictive envelope for the model outcomes relative to the measurement error (95% credibility). The graph on the left corresponds to data set A, the middle graph corresponds to data set B and the graph on the right corresponds to data set C.

Figure 3. Marginal posteriors of parameter for data set A [-], data set B [–] and data set C [-·-]. The Y -axes qualitatively represent the probability density, which ranges from 0 to 1.

Environmental Technology Table 2.

Comparison of parameter uncertainty estimation including Bayesian inference and literature results. Data set A

Parameter

Downloaded by [University of Chicago Library] at 14:05 26 December 2014

τ KS kSTO μMAX,S δ K1 K2

1625

Data set B

Data set C

Bayesian inference

Sin et al. [2]

Bayesian inference

Hoque et al. [5]

Bayesian inference

1.67 (0.36) 0.86 (0.37) 1.95 (0.71) 1.32 (0.68) 2.91 (0.59) 0.103 (0.088) 7.1 × 10−4 (3.7 × 10−3 )

0.51 (0.07) 0.70 (0.10) 2.02 0.97 2.56 (0.08) 0.102 (0.012) 1.2 × 10−3 (3 × 10−4 )

1.13 (0.33) 0.81 (0.92) 2.25 (0.65) 2.14 (0.61) 4.58 (0.60) 0.026 (0.038) 9.0 × 10−4 (1.2 × 10−3 )

1.50 (0.02) 1.23 (0.06) 2.53 1.99 4.36 (0.13) 0.051 (0.003) 1.7 × 10−7 (4.8 × 10−7 )

1.87 (0.22) 6.42 (1.83) 7.68 (2.17) 13.0 (1.51) 5.00 (0.20) 0.140 (0.094) 3.1 × 10−3 (5.4 × 10−3 )

sub-sampling, which was applied in order to reduce the auto-correlation of the residuals. The modes obtained for the efficiency of the oxidative phosphorylation (δ) when analysing data sets B (4.58 ± 0.60) and C (5.00 ± 0.20) are higher than the theoretically expectable values, which are of between 1 and 3.[2–39] However, Hoque et al. [5] found the same limitation for the parameter δ and reported a final value of 4.36 ± 0.132. These high values call into question the mechanistic significance of the parameter according to the original conception of the mathematical model. Such values of the parameter δ also lead to high yield coefficients. For all cases, YSTO (0.80–0.90) is found to be the highest yield coefficient. Direct biomass growth on acetate (YH,S ) (0.57–0.75) is estimated as very similar to the indirect biomass growth via intracellular storage products (YSTO · YH,STO ). Thus, a narrow ratio of 0.96–0.98 is calculated between both pathways, which are in agreement with the findings of Beun et al.[38] The parameters KS , kSTO and μMAX,S reach higher modes for data set C than for data sets A and B. Furthermore, the inference results for data set C are the only providing kSTO lower than μMAX,S . Thus, different numerical values are obtained during the assessment for a given parameter depending on the OUR profile being analysed. This is because parameter values are expected to be highly influenced by many factors, such as wastewater composition, bioreactor operating conditions, dominant microbial communities, etc.[16] In this regard, the activated sludge used in data set C (obtained from a lab-scale SBR treating high-loaded leachate) was very different with respect to sludges used in data sets A and B (obtained from WWTPs). On the other hand, the boundary conditions considered during the inference exercise will decisively affect the results. This may be the case of XSTO (0), which was measured/assumed (depending on the case), although it was reported as a potential source of problems for the identification of parameters when using OUR profiles alone.[2] The case of XH (0) is similar, since it was calculated taking into account the endogenous OUR and assuming the values for fXI and bH as proposed in ASM3.[25] Any mismatch in these boundary conditions will be propagated during the inference exercise and consequently may result in non-realistic results.

Predictive envelopes for the OUR model outcomes relative to the measurement error are also depicted in Figure 2. The wideness of the envelope is particularly significant for data set A. Furthermore, the high value of σ (0.052 ± 0.009) obtained for such data set (Figure 3) indicates that the model have some problems in reproducing the dynamic phase of the OUR response. Indeed, Sin et al. [2] already reported that the model was unable to perfectly fit in the second peak in the OUR profile. In this study, the model was unable to perfectly fit in the first peak. Although this is a contradiction in the results, both cases seem to be possible according to the predictive envelope. Lower values of σ were obtained for data sets B (0.035 ± 0.007) and C (0.018 ± 0.003) due to a more simple experimental design, among other possible reasons. 3.3. Frequentist inference Parameter assessment is also conducted assuming the posterior distribution as MVG. It could be observed that the condition number for H relative to θ was higher than the machine precision, which implies inappropriateness in the estimation of the covariance matrix. On the other hand, the condition number for H relative to θ ∼p takes a value around 103 , and thus the estimation of the covariance matrix (once fixed θ p ) can be considered as ‘non-problematic’. This fact evidences that the model is overparametrized with respect to OUR data. Yet, all parameters are still necessary in order to maintain the structure of the model (and to reproduce data correctly). 3.4. Residuals analysis The validity of the above statistical inference is dependent on the assumptions taken for the measurement model defined in Equation (7). In Figure 4, the residuals analyses are presented on basis of histograms and auto-correlation plots. The histogram of the residuals for each data set is then compared with respect to the corresponding measurement model (Figure 4, top). However, from the comparison between the histograms and the measurement model it may be still difficult to assess whether the normality

Ž.J. Zonta et al. 30

60

10

20

40

5

0 −0.1 −0.05 0 0.05 Residuals

−0.2 −0.1 0 0.1 0.2 Residuals

Downloaded by [University of Chicago Library] at 14:05 26 December 2014

0

5 Lags

10

0.1

0.5

0

−0.5 0

20 0 −0.1 −0.05 0 0.05 Residuals

0.1

0.5

ACF

ACF

0.5

−0.5 0

10

ACF

0

Density

15 Density

Density

1626

5 Lags

10

0

−0.5

0

5 Lags

10

Figure 4. Residuals analysis. Top: histograms of the residuals and the relative most probable measurement model, Norm(0, σm2 ), in solid lines. The measurement model is given within its 95% prediction intervals in dashed lines, since its parameter σ 2 is estimated during the inferential procedure. Bottom: ACF plots with 95% significance intervals. The graphs on the left correspond to data set A, the middle graphs correspond to data set B and the graphs on the right correspond to data set C.

assumption over ε is appropriate, or instead should be rejected. Kleinbaum et al. [40] reported that the confidence intervals used in a regression analysis are robust in the sense that only extreme departures of the distribution of the residuals from normality yield spurious results. On the other hand, the auto-correlation function (ACF) of the residuals (Figure 4, bottom) is within the significance levels, and thus the independence requirement of the measurements is fulfilled. It seems interesting to remark that, if all the samples from a respirometry device with a high frequency of sampling were used raw and without further processing, the residuals would be highly auto-correlated, which would lead to the underestimation of parameter uncertainty. In order to avoid this, aside from sub-sampling, they can be applied whitening techniques as well as the correction of the residuals by means of auto-regressive models which permit their proper weighting.[41] 3.5. Comparison of Bayesian vs. frequentist inference In order to contrast the parameter uncertainty assessment under both inferential procedures, the results obtained through the Bayesian posterior are compared with respect to those obtained through the MVG approximation. The ‘full-case’ comparison takes into account the parameter θ = [θ ∼p , θ p ], whereas the ‘reduced-case’ comparison only accounts for θ ∼p (while θ p is assumed as perfectly known according to the values provided by the original authors). The comparison is performed over a two-dimensional parameter space. Thus, considering data set B as illustrative example, and for the ‘reduced-case’ comparison, Figure 5 shows the MCMC sample, the kernel-density p approximation of the posterior p(θ ∼p |D, θm ) and its linear

∼p

approximation Norm(θm , C∗ ). In general terms, the MVG approximation seems to estimate well the shape of the Bayesian posterior. However, the MVG approximation overestimates the uncertainty with respect to the Bayesian posterior in those cases where it is truncated as strictly positive (see for example, KS vs. kSTO ). This is because the MVG distribution is defined over an unbounded parameter space. On the other hand, if one considers the case where the inference is not affected by the positivity constraint in parameters (see for example, τ vs. μMAX,S ), the MVG approximation may underestimate the parameter uncertainty in agreement with those observations reported by Omlin and Reichert [11] or Vrugt and Bouten.[12] The above results are also confirmed for data sets A and C (data not shown); the MVG approximation is a reasonable estimate for the Bayesian posterior. The quality of the frequentist approximation is subsequently analysed through the K − S statistic for both the ‘reduced-case’ (Figure 6, top) and the ‘full-case’ (Figure 6, bottom), and considering the three OUR data sets. In the ‘reduced-case’ study, the MVG approximation is shown as a reasonable approximation for the posterior. However, in the ‘full-case’ study, the higher values of the K − S statistic suggest that the MVG approximation is not so appropriate. It is worth to remark that in the ‘full-case’ the Bayes’ posterior considers θ p through a uniform prior distribution of finite range, while the corresponding MVG approxima∼p tion Norm(θm , C) is estimated by fixing the problematic p p parameter θ at θm (modal value; derived from the posterior distribution). In other words, during the Bayesian inference θ p is weakly specified, while during the frequentist inference θ p is considered as perfectly known. Such assumption of a perfectly known θ p is questionable (i.e. if

1627

Figure 5. Two-dimensional comparison (95% credibility) between the posterior distribution p(θ ∼p |D, θm ) [-] estimated from the MCMC ∼p samples [ ] and the linear approximation Norm(θm , C∗ ) confidence ellipses [-·-]. This case corresponds to the ‘reduced-case’ for data set B. p

K−S statistic

0.5 0.4 0.3 0.2 0.1 0 0.5 K−S statistic

Downloaded by [University of Chicago Library] at 14:05 26 December 2014

Environmental Technology

0.4 0.3 0.2 0.1 0

(τ,

kS

TO

)

(τ,

K

S

(k )

ST O

,K

S

(τ, )

(k

μ

MA

X,S

)

ST O

(K

,μ

MA

S

X,S

)

(τ,

δ)

,μ

MA

X,S

)

(k

ST O

,δ)

(K

S

,δ)

(μ

MA

X,S

,δ)

Figure 6. The two-dimensional K − S statistic for the ‘reduced-case’ (top) and the ‘full-case’ (bottom) comparisons relative to data set A (black bars), data set B (grey bars) and data set C (white bars). Figure 5 is evaluated at the top (grey bars). The linear approximation is reasonable in the ‘reduced-case’ while it is unsatisfactory in the ‘full-case’.

θ p is problematic because its value is highly uncertain, then how is it possible to assume it as perfectly known?) but necessary to make the computation of the covariance matrix

feasible. Hence, the disparity in the assumptions about θ p between both inferential procedures (i.e. vague knowledge vs. perfect knowledge) may justify the quantitative

1628

Ž.J. Zonta et al.

differences observed through the K − S statistic. Moreover, since the model is nonlinear in parameters, the goodness of the MVG approximation to the Bayesian posterior might change just with the location of the estimation (which depends on the available data).

Downloaded by [University of Chicago Library] at 14:05 26 December 2014

4.

Conclusions

This research introduces the use of Bayesian inference for the assessment of parameter uncertainty when addressing issues of practical identifiability in the field of wastewater treatment modelling. Although the use of Bayesian inference (based on the definition of model parameters as probability distributions) is already described in the scientific literature, it can be still considered as a new approach when using ASM-type models. The Bayesian method was then compared with the commonly applied methodology of the frequentist inference (based on the estimation of optimal parameter values and confidence regions). The inference exercise was carried out considering: (i) a mathematical model accounting for simultaneous intracellular storage and biomass growth under aerobic conditions, (ii) experimental data exclusively consisted of respirometric OUR profiles (three data sets) and (iii) the aid of GSA. Bayesian inference allowed a probabilistic description of the assessed parameters. However, the exclusive use of OUR data led to limitations in the practical identifiability of parameters. Since it was demonstrated that Bayesian inference could be reduced to a frequentist approach under particular hypotheses, the first methodology can be considered as more generalist. Working with ASM-type models, there are many situations where is not possible to precise a priori the suitability of the frequentist methods to deal with the estimation of parameter uncertainty. Use of Bayesian inference in such scenarios may become advisable. Acknowledgements This research was supported by ACC1Ó/Generalitat de Catalunya under programme INNOEMPRESA-2009 (YT-CONIT project IEINN09–1–0180) and the Spanish Ministry of Science and Innovation (ADAMOX project CTM2010–18212). Authors thank Dr A. Guisasola (UAB) for providing the OUR profile used as data set A.

References [1] Vanrolleghem PA, Spanjers H, Petersen B, Ginestet P, Takacs I. Estimating (combinations of) activated sludge model no. 1 parameters and components by respirometry. Water Sci Technol. 1999;39(1):195–214. [2] Sin G, Guisasola A, De Pauw DJW, Baeza JA, Carrera J, Vanrolleghem PA. A new approach for modelling simultaneous storage and growth processes for activated sludge systems under aerobic conditions. Biotechnol Bioeng. 2005;92:600–613. [3] Magrí A, Flotats X. Modelling of biological nitrogen removal from the liquid fraction of pig slurry in a sequencing batch reactor. Biosyst Eng. 2008;101:239–259.

[4] Coello-Oviedo MD, Barragán-Sánchez J, Aragón-Cruz C, Quiroga-Alonso JM. A new approach to toxicity determination by respirometry. Environ Technol. 2009;30: 1601–1605. [5] Hoque MA, Aravinthan V, Pradhan NM. Assessment on activated sludge models for acetate biodegradation under aerobic conditions. Water Sci Technol. 2009;60:983–994. [6] Ubay Çokgör E, Sözen S, Insel G, Orhon D. Respirometric evaluation of biodegradation characteristics of dairy wastewater for organic carbon removal. Environ Technol. 2009;30:1169–1176. [7] Gillot S, Choubert J-M. Biodegradable organic matter in domestic wastewaters: comparison of selected fractionation techniques. Water Sci Technol. 2010;62:630–639. [8] Checchi N, Marsili-Libelli S. Reliability of parameter estimation in respirometric models. Water Res. 2005;39: 3686–3696. [9] Marsili-Libelli S, Tabani F. Accuracy analysis of a respirometer for activated sludge dynamic modelling. Water Res. 2002;36:1181–1192. [10] D’Agostini G. Bayesian inference in processing experimental data: principles and basic applications. Rep Prog Phys. 2003;66:1383–1419. [11] Omlin M, Reichert P. A comparison of techniques for the estimation of model prediction uncertainty. Ecol Modell. 1999;115:45–59. [12] Vrugt JA, Bouten W. Validity of first-order approximations to describe parameter uncertainty in soil hydrologic models. Soil Sci Soc Am J. 2002;66:1740–1751. [13] Dochain D, Vanrolleghem PA. Dynamical modelling and estimation in wastewater treatment processes. Padstow, UK: IWA Publishing; 2001. [14] Renard B, Kavetski D, Kuczera G, Thyer M, Franks SW. Understanding predictive uncertainty in hydrologic modeling: the challenge of identifying input and structural errors. Water Resour Res. 2010;46:1–22. [15] Petersen B, Gernaey K, Vanrolleghem PA. Practical identifiability of model parameters by combined respirometric-titrimetric measurements. Water Sci Technol. 2001;43(7):347–355. [16] Ni B-J, Yu H-Q. Simulation of heterotrophic storage and growth processes in activated sludge under aerobic conditions. Chem Eng J. 2008;140:101–109. [17] Guisasola A, Sin G, Baeza JA, Carrera J, Vanrolleghem PA. Limitations of ASM1 and ASM3: a comparison based on batch oxygen uptake rate profiles from different fullscale wastewater treatment plants. Water Sci Technol. 2005;52(10–11):69–77. [18] Saltelli A, Ratto M, Tarantola S, Campolongo F. Sensitivity analysis practices: strategies for model-based inference. Reliab Eng Syst Saf. 2006;91:1109–1125. [19] Sin G, Gernaey KV, Neumann MB, van Loosdrecht MCM, Gujer W. Global sensitivity analysis in wastewater treatment plant model applications: prioritizing sources of uncertainty. Water Res. 2011;45:639–651. [20] Ratto M, Tarantola S, Saltelli A. Sensitivity analysis in model calibration: GSA-GLUE approach. Comput Phys Commun. 2001;136:212–224. [21] Saltelli A, Annoni P, Azzini I, Campolongo F, Ratto M, Tarantola S. Variance based sensitivity analysis of model output. Design and estimator for the total sensitivity index. Comput Phys Commun. 2010;181:259–270. [22] Laine M. Adaptive MCMC methods with applications in environmental and geophysical models [PhD thesis]. Lappeenranta, Finland: Department of Mathematics and Physics, Lappeenranta University of Technology; 2008.

Downloaded by [University of Chicago Library] at 14:05 26 December 2014

Environmental Technology [23] Gelman A, Carlin JB, Stern HS, Rubin DB. Bayesian data analysis, texts in statistical science. 2nd ed. Boca Raton (FL): Chapman & Hall/CRC; 2004. [24] Marsili-Libelli S, Guerrizio S, Checchi N. Confidence regions of estimated parameters for ecological systems. Ecol Modell. 2003;165:127–146. [25] Gujer W, Henze M, Mino T, van Loosdrecht M. Activated sludge model no.3. Water Sci Technol. 1999;39(1): 183–193. [26] Vanrolleghem PA, Van Daele M, Dochain D. Practical identifiability of a biokinetic model of activated sludge respiration. Water Res. 1995;29:2561–2570. [27] D’Errico J. Adaptive robust numerical differentiation [Internet]. Matlab File Exchange. 2006 [revised 2013 Oct 22]. Available from: http://www.mathworks.com/mat labcentral/fileexchange/13490-adaptive-robust-numericaldifferentiation [28] Haario H, Laine M, Mira A, Saksman E. DRAM: efficient adaptive MCMC. Stat Comput. 2006;16:339–354. [29] Geweke J. Evaluating the accuracy of sampling-based approaches to the calculation of posterior moments. In: Bernado JM, Berger JO, Dawid AP, Smith AFM, editors. Bayesian statistics 4. Oxford, UK: Oxford University Press; 1992. p. 169–194. [30] Peacock JA. Two-dimensional goodness-of-fit testing in astronomy. Mon Not Roy Astron Soc. 1983;202:615–627. [31] Spanjers H, Vanrolleghem PA, Olsson G, Dold PL. Respirometry in control of the activated sludge process: principles. Scientific and Technical Report 7. London: IAWQ, International Association on Water Quality; 1998. [32] Bontempi G, Birattari M, Bersini H. Lazy learning for local modelling and control design. Int J Control. 1997;72: 643–658.

1629

[33] Dizdaroglu-Risvanoglu G, Karahan O, Ubay-Cokgor E, Orhon D, van Loosdrecht MCM. Substrate storage concepts in modeling activated sludge systems for tannery wastewaters. J Environ Sci Health A Toxic Hazard Subst Environ Eng. 2007;42:2159–2166. [34] Oakley JE, O’Hagan A. Probabilistic sensitivity analysis of complex models: a Bayesian approach. J R Stat Soc Ser B Stat Methodol. 2004;66:751–769. [35] O’Hagan T. The GEM software project [Internet]. Author’s resources. 2009 [revised 2013 Oct 22]. Available from: http://www.tonyohagan.co.uk/academic/GEM/index.html [36] Dodge Y. The concise encyclopedia of statistics. New York (NY): Springer Science and Business Media LLC; 2008. [37] Press WH, Teukolsky SA, Vetterling WT, Flannery BP. Numerical recipes: the art of scientific computing. 3rd ed. New York (NY): Cambridge University Press; 2007. [38] Beun JJ, Paletta F, van Loosdrecht MCM, Heijnen JJ. Stoichiometry and kinetics of poly-β-hydroxybutyrate metabolism in aerobic, slow growing, activated sludge cultures. Biotechnol Bioeng. 2000;67:379–389. [39] Dias JML, Oehmen A, Serafim LS, Lemos PC, Reis MAM, Oliveira R. Metabolic modelling of polyhydroxyalkanoate copolymers production by mixed microbial cultures. BMC Syst Biol. 2008;2:59. Available from: http://www.biomedcentral.com/content/pdf/1752-0509-259.pdf [40] Kleinbaum DG, Kupper LL, Nizam A, Muller KE. Applied regression analysis and multivariable methods. 4th ed. Belmont (CA): Duxbury Press; 2008. [41] Neumann MB, Gujer W. Underestimation of uncertainty in statistical regression of environmental models: influence of model structure uncertainty. Environ Sci Technol. 2008;42:4037–4043.