Global Sensitivity Analysis of the Biological Part of the Integrated BTEX Bioremediation Model

Abstract

The greatest challenge when working with mathematical models describing biological processes is to choose appropriate values of kinetic parameters. This applies especially to the models of environmental biological processes, such as bioremediation. Parameter values are usually obtained from various sources and under specific experimental conditions and constitute sources of uncertainties. A useful tool for dealing with model parameter uncertainties is sensitivity analysis. It studies how a given model depends on its input factors. In this work, the nonstationary Extended Fourier Amplitude Sensitivity analysis was used to analyze the effect of 28 kinetic parameters on 18 metabolite concentrations and 16 reaction rates of the biological part of the integrated BTEX bioremediation model. First-order sensitivities and total effect indices were estimated for three time points of the process to identify the places of the metabolism control shift. Analyses revealed that the model is nonadditive and there are no strong interactions between the parameters. Results of this study indicate the importance of using the global sensitivity analysis before bioremediation process planning.

Introduction

Biological systems consist of complex processes that can be described with different mathematical models. Mathematical models of biological processes (Reed et al., 2008; Peskov et al., 2012; Wayman and Varner, 2013) are mainly based on the simplification of the complex control and regulatory mechanism, but ensure good insight into processes taking place in biological systems. Mathematical models are very useful for hypothesis testing and selection among different possible regulatory structures (Gombert and Nielsen, 2000).

According to Wiechert et al. (2002), mathematical models of biological systems have several aims as follows: (i) structural understanding (Geng et al., 2014), (ii) exploratory simulation (Essaid et al., 1995), (iii) interpretation and evaluation of measured data (Schirmer et al., 2000), (iv) systems analysis (Geng et al., 2013), (v) prediction and design (Schirmer et al., 1999), and (vi) optimization (Nicol et al., 1994). There are several phases in the development of a mathematical model of a biological process. First, the purpose of the model must be defined. Thereupon, the network structure and diagram of the model, kinetic rate expressions (mostly enzyme-catalyzed biochemical reactions), and model parameters should be established to describe the time-dependent change of model variables (Kent et al., 2013). The final step of model development includes model validation (Almquist et al., 2014).

Due to the enormous interest in environmental protection, mathematical models dealing with biological degradation of pollutants are being used more often. Activated sludge models are widely applied in describing wastewater treatments (Henze et al., 1987, 1995, 1999; Gujer et al., 1999). Bioremediation, a technology for degradation of harmful chemicals, has gained popularity lately due to its economic advantages, environmental friendliness, and ease of application. Some methodical models describing bioremediation processes are given in literature. Yaghmaei (2002) developed a model for the in situ bioremediation of contaminated soil; Li et al. (1994), and Surajudeen (2012) presented mathematical models describing the bioremediation of oil-contaminated soil. The most detailed bioremediation model available in literature was proposed by Maurer and Rittmann (2004).

Mathematical analysis of complex reactions and intraregulations of those processes allows easier process planning and optimization. It is important to mention that all these models are complex and overparameterized. Therefore, currently the greatest challenge in the mathematical modeling of a biological process is to choose the values of the kinetic parameters (Jia et al., 2012). Parameter values are usually obtained from various sources and are estimated under specific experimental conditions and consequently constitute sources of uncertainties (Kent et al., 2013; Vanlier et al., 2013).

Sensitivity analysis is a useful tool for dealing with these uncertainties. According to Saltelli et al. (2010), sensitivity analysis is the study of how uncertainty in the model output can be apportioned to different sources of uncertainty in the model input factor, such as model parameters. Sensitivity analysis methods can be categorized according to the outcome of the sensitivity measures: qualitative or quantitative methods, local or global methods, and methods that depend on or are independent of the model characteristics (Homma and Saltelli, 1996; Schwieger, 2004).

Global sensitivity analyses are especially useful in dealing with complex models with a large set of model parameters focusing on the output uncertainties over the entire range of values of the input parameters simultaneously (Cosenza et al., 2014). According to Iooss and Lemaître (2015), global sensitivity analysis methods can be divided into (i) screening methods (based on a discretization of the inputs in levels, allowing fast exploration of the code behavior), (ii) importance measure methods (methods based on the analysis of the linear models and method based on functional decomposition of variance), and (iii) deep exploration of the model. Both Sobol' and Fourier Amplitude Sensitivity Test (FAST) methods belong to a group of variance-based global sensitivity analysis methods (Chan et al., 1997).

In variance-based methods, the variability or uncertainty associated with an important input parameter is propagated throughout the model, resulting in a large contribution to the overall output variability (Chan et al., 1997). Comparing the Sobol' and FAST methods, the advantage of FAST is that the evaluation of sensitivity coefficients can be carried out independently for each parameter by only one simulation because all the terms in Fourier expansion are mutually orthogonal (Saltelli and Bolado, 1998), while Sobol' uses the Monte Carlo integration procedure. FAST is also computationally more efficient than Sobol' because of it robustness, especially at small sample size (Chan et al., 1997; Saltelli et al., 1999).

FAST was originally developed for a functional model whose input parameters were uniformly distributed, but over the years, the procedure to transform the nonuniform distributions into uniform was developed (Fang et al., 2003), allowing the application of the FAST methodology in a wide range of mathematical models (from engineering to physics to chemistry).

There are several studies in literature that apply sensitivity analysis methods to bioremediation process models. Benekos et al. (2007) applied a sampling-based sensitivity analysis method on a complex model that couples chemical degradation and transformation with movement of four chlorinated ethenes undergoing bioremediation. The sampling-based approach is based on Monte Carlo simulations and involves the generation and exploration of a mapping from uncertain analysis input to uncertain analysis results (Helton et al., 2006). Zhao et al. (2010) used the multiparametric sensitivity analysis (MPSA) to identify the most important parameters for the uranium concentration in the electron load and nonelectron load models of Geobacter-mediated bioremediation. According to the basic concept of the MPSA, for each parameter of the analyzed model, a comparison of the distributions of the parameter values associated with the acceptable and unacceptable is performed. If the two distributions are not statistically different, the parameter is classified as insensitive; otherwise, the parameter is classified as sensitive (Nassar, 2011). The extension of the model describing the Geobacter-mediated bioremediation, proposed by Zhao et al. (2011), integrated the model on the electron capacitance of Geobacter with a comprehensive simulator of coupled fluid flow, hydrologic transport, heat transfer, and biogeochemical reactions and analyzed it by applying the Extended FAST method combined with the Morris screening method. Based on available information on the use of the sensitivity analysis method in the bioremediation process, it can be concluded that although suitable, the Extended FAST is not widely used because of the lack of tools addressing the computational challenge (there is no commercially available software package performing the Extended FAST analysis).

In this work, the Extended FAST was used for analysis of the kinetic parameters of the biological part of the integrated BTEX bioremediation model. The motivation for this study was the lack of available information in literature on the use of the nonstationary Extended FAST method on mathematical models describing a bioremediation process. The aim of this study is to determine the first-order sensitivities and total effect indices of model metabolites and reactions at three time points to identify the possible shift in metabolism control.

Materials and Methods

Simulation of integrated BTEX bioremediation model

For analysis and simulation of the BTEX bioremediation processes, only the biological part of the Coupled Biological and Chemical (CBC) mathematical model, proposed by Maurer and Rittmann (2004), was used in this work. The proposed model describes the basic microbiological and chemical reactions that affect the presence and concentrations of BTEX intermediates. The biological part of the model is based on modified Monod kinetics and describes the BTEX degradation by bacteria, methanogenesis, fermentation, oxidation of fermentative products, and biomass decay. The specific reactions described with the model are (i) aerobic growth, (ii) denitrification, (iii) iron reduction, (iv) sulfate reduction, (v) fermentation, (vi) lysis of fast-growing heterotrophic bacteria that respire oxygen and nitrate, (vii) lysis of iron-reducing bacteria, (viii) lysis of sulfate-reducing bacteria, (ix) lysis of acetoclastic methanogens, (x) lysis of hydrogen-oxidizing methanogens, (xi) cleavage of acetate, (xii) methanogenesis, (xiii) aerobic growth on organic fermentation products, (xiv) anoxic growth on organic fermentation products, (xv) iron reduction on organic fermentation products, and (xvi) sulfate reduction on organic fermentation products. The rate terms of the microbiological reactions based on modified Monod kinetics describe the biomass synthesis and electron donor utilization, taking into account the effects of multiple substrates, nutrients, and inhibitors (Maurer and Rittmann, 2004). The chemical part of the model describes the dissolution and precipitation of calcium carbonate, conversion of iron(II) monosulfide into thermodynamically stable pyrite, and the reductive dissolution of solid iron(III) by oxidation of sulfide. The model includes 16 dissolved components (BTEX, acetate, oxygen, nitrate, ammonium, sulfate, methane, hydrogen, calcium, iron(II), inorganic carbon, and inorganic sulfide), distinguishes five biomass types (fast-growing heterotrophic bacteria, sulfate-reducing bacteria, iron-reducing bacteria, acetoclastic methanogens, and hydrogen-oxidizing methanogens), and has five inorganic solid phases (calcium carbonate, iron(III) oxide, iron(II) monosulfide, pyrite, and elementary sulfur). The biological part of the CBC model used for analysis includes 18 balances in the form of ordinary differential equations. The time change of each metabolite concentration is the difference between the rate of formation and the consumption of individual molecules. The list of model balances derived by multiplying the stoichiometric matrix with the reaction rate vector, proposed by Maurer and Rittmann (2004), is given in Table 1, whereas the numerical values and units of the analyzed kinetic parameters are shown in Table 2.

Table 1.

Table of Model Balances Based on the Stoichiometric Matrix and Reaction Rate Vector Given by Maurer and Rittmann, 2004

The list of balances is derived by multiplying the stoichiometric matrix with the reaction rate vector proposed by Maurer and Rittmann, 2004.

Table 2.

List of Kinetic Parameters of the BTEX Bioremediation Mathematical Model Analyzed in This Work (Maurer and Rittmann, 2004)

	Symbol	Description	Value ^a	Unit
1.	K _O	Half-saturation coefficient for oxygen	3.00 × 10⁻⁶	mol_O2/L
2.	K _HC	Half-saturation coefficient for BTEX	0.25 × 10⁻⁶	mol_C/L
3.	K _F	Half-saturation coefficient for fermentation products	1.00 × 10⁻⁵	mol_O2/L
4.	K _N	Half-saturation coefficient for nitrogen and biomass growth	5.00 × 10⁻⁶	mol_N/L
5.	μ _O	Heterotrophic growth rate constant growing on oxygen	1.60 × 10⁻¹	1/day
6.	μ _O,F	Aerobic heterotrophic growth rate constant with acetate	6.00 × 10⁰	1/day
7.	K _NO	Half-saturation coefficient for nitrate	3.00 × 10⁻⁶	mol_N/L
8.	η _NO	Reduction coefficient for the anoxic maximum growth rate	8.00 × 10⁻¹	—
9.	K _Fe2O3	Half-saturation coefficient for iron(III)	1.00 × 10⁻³	mol_C/mol_Fe
10.	μ _Fe	Heterotrophic growth rate constant growing on iron(III)	8.00 × 10⁻²	1/day
11.	μ _Fe,F	Heterotrophic growth rate constant with acetate for iron(III)	1.00 × 10⁻¹	1/day
12.	K _SO	Half-saturation coefficient for sulfate	1.00 × 10⁻⁶	mol_S/L
13.	K _IS1	Inhibition coefficient for sulfide for acetoclastic methanogens	3.40 × 10⁻³	mol_S/L
14.	K _IS3	Inhibition coefficient for sulfide of acetate SBR	1.10 × 10⁻³	mol_S/L
15.	K _IS4	Inhibition coefficient for sulfide of SBR	21.20 × 10⁻³	mol_S/L
16.	K _IS5	Inhibition coefficient for sulfide of fermentation	1.00 × 10⁻³	mol_S/L
17.	μ _SO	Aerobic heterotrophic growth rate constant growing on sulfate	8.00 × 10⁻²	1/day
18.	μ _SO,F	Heterotrophic growth rate constant with acetate for iron(III)	1.00 × 10⁻¹	1/day
19.	K _H2,an	Inhibition coefficient for hydrogen and fermentation	5.00 × 10⁻⁶	mol_H2/L
20.	η _an	Reduction coefficient for the fermentative maximum growth rate	1.00 × 10⁻¹	—
21.	b _BO	Lysis rate constant for X_BO	1.60 × 10⁻²	1/day
22.	b _BF	Lysis rate constant for X_BF	6.00 × 10⁻¹	1/day
23.	b _BS	Lysis rate constant for X_BSO	8.00 × 10⁻³	1/day
24.	b _BAc	Lysis rate constant for X_BAc	2.14 × 10⁻⁴	1/day
25.	b _BH2	Lysis rate constant for X_BH2	5.00 × 10⁻²	1/day
26.	μ _Ac	Growth rate constant of acetoclastic organisms	9.00 × 10⁻²	1/day
27.	μ _H2	Growth rate constant of methanogens	8.00 × 10⁻¹	1/day
28.	K _H2	Inhibition coefficient for hydrogen utilization of methanogens	1.50 × 10⁻⁷	mol_H2/L

The information on origin of initial value of kinetic parameters is given by Maurer and Rittmann, 2004, in detail.

Simulations of the biological part of the CBC model were performed using W.R. Mathematica 9.0 software with initial concentrations presented by Kao and Wang (2000). They presented the results of intrinsic bioremediation at a gasoline spill in Dublin, Bladen County, North Carolina, in 1995. Due to the appearance of nonaqueous phase, the liquid hydrocarbons underneath the spill site were continuously leaching into the groundwater. Twenty-one monitoring wells were installed to observe the BTEX distribution. Each BTEX plume was divided into four zones (source zone, mid-plume area, denitrification area, and down-gradient oxidized area) for easier interpretation of the results. Based on the results presented by Kao and Wang (2000), it can be seen that other electron acceptors (e.g., nitrate and ferric ions) were used to oxidize dissolved contaminants after oxygen depletion. In addition, by reducing the reduction potential in the source zone and mid-plume zone, it reflects the change for oxidizing to reducing conditions. Groundwater samples were collected from four representative wells. To evaluate the trend of the BTEX degradation process, Kao and Wang (2000) analyzed the distribution of microorganisms within plume. Analysis revealed the high total anaerobe population, denitrifiers, iron reducers, sulfate reducers, methanogens, and small aerobe population (Kao and Wang, 2000). For CBC model simulations, the groundwater sample composition presented as MW1 was used.

Global sensitivity analysis Extended FAST

The importance of individual parameters of the biological part of the CBC model was analyzed numerically by the Extended FAST method (Saltelli et al., 1999, 2000, 2010) using W. R. Mathematica 9.0 software. The Extended FAST method allows calculation of two sensitivity indices (first-order effect index and total effect index). The system analysis is based on evaluation of the first moments (joint expected value E and variance σ²) from the multidimensional joint probability distribution p(y₁, y₂…. y_n). The effect of each individual process variable, x_i, on each of the system outputs, y_k, is evaluated by the ratio of the ensemble variance of the conditional expected value relative to the total ensemble dispersion.

The Extended FAST method is based on the transformation of the multidimensional parameter space into one-dimensional space of a single parameter, s, by the nonlinear transformation of each parameter according to Equation (1). In the Extended FAST method, the sensitivity coefficients are evaluated by a search curve that scans the space of the input factors so that each factor is explored with a selected integer frequency (Saltelli et al., 2000): \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} { x_ { \rm { i } } } = \frac { 1 } { 2 } + \frac { 1 } { { \rm { \pi } } } \cdot { \rm { arcsin } } \left( { { \rm { sin } } \left( { { \rm { \pi } } \cdot { \omega _ { \rm { i } } } \cdot s + { \varphi _ { \rm { i } } } } \right) } \right) \tag { 1 } \end{align*} \end{document}

where s is the sampling parameter in the range of s ∈[–1, 1], ω_i are selected frequencies, and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$ \varphi_{\rm i}$$ \end{document} are randomly selected phase angles. The phase angles are random numbers from the range \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$ \varphi_{\rm i}\,$$ \end{document} ∈[–π, π], while the frequencies are integers and form an incommensurate set. A sequence of the scan parameter s generates a large set of random and uncorrelated sequences of the input variables. The input variables (transformed kinetic parameters) were sampled with a frequency range from 8 to 64 Hz, while 200 Hz was selected as the highest frequency in the Fourier analysis.

Responses of the output variable are expanded into Fourier series by which the overall variance, D, of the output function is decomposed into summands of squares of Fourier coefficients, A_j and B_j [Eqs. (2 )–(4)]: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}Var \left( Y \right) = 2 \mathop \sum \limits_{j = 1}^{ \infty} { \left( {A_{ \rm{j}}^{ \rm{2}} + B_{ \rm{j}}^{ \rm{2}}} \right) } \tag{2}\end{align*} \end{document} \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} { A_ { \rm { j } } } = \frac { 1 } { { 2 \pi } } \int \limits_ { - \pi } ^ { \pi } { y \left( s \right) } \cdot \cos \left( \, { j \cdot s } \right) \cdot { \rm { d } } s \tag { 3 } \end{align*} \end{document} \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} { B_ { \rm { j } } } = \frac { 1 } { { 2 \pi } } \int \limits_ { - \pi } ^ { \pi } { y \left( s \right) } \cdot \sin \left( \, { j \cdot s } \right) \cdot { \rm { d } } s \tag { 4 } \end{align*} \end{document}

Individual parameter contributions in the total dispersion are calculated by the corresponding harmonics [(Eq. (5)]. The summation is meant to include all the harmonics related to the frequency associated to the input variable considered (Saltelli and Bolado, 1998), from minimum frequency for transformation of the selected parameter to maximum frequency defined for the Fourier analysis. \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}Var_{ \rm{i}} \left( Y \right) = 2 \mathop \sum \limits_{j = 1}^{ \rm m} { \left( {{ A^{2}_{{ \rm{j}}{ \omega }_{ \rm{i}}}}} + B_{{ \rm{j}}{ \omega }_{ \rm{i}}}^{2} \right) } \tag{5}\end{align*} \end{document}

The relative sensitivity coefficient or first-order effect index, S_i, is a measure of the global influence of parameter x_i on the output and is given as the ratio of partial dispersion and total dispersion [(Eq. (6)]: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} { S } _ { \rm i } = { \frac { { Var } _ { \rm i } \left( Y \right) } { Var \left( Y \right) } } \tag { 6 } \end{align*} \end{document}

On the other hand, the total effect index is defined as follows [(Eq. (7)]: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} { S_T } _ { \rm { i } } = 1 - { \frac { { Var_ { \left( { - i } \right) } } \left( Y \right) } { Var \left( Y \right) } } \tag { 7 } \end{align*} \end{document}

and describes the amount of the variance that stays undescribed when only variable x_i is varied.

Results and Discussion

Nonstationary Extended FAST global sensitivity analysis was applied for identification of the most important kinetic parameters of the biological part of the integrated BTEX bioremediation model. The first-order indices and total effect indices were calculated for all 18 metabolites and 16 reactions included into the model at three time points. The first-order effect index measures the relative importance of an individual kinetic parameter in total uncertainty and can be seen as indicating where to direct efforts to reduce the uncertainty. On the other hand, the total effect index is linked to the amount of the output variance that would remain unexplained if an individual kinetic parameter was left free to vary over its uncertainty range while all other parameters have been learnt. In practice, total sensitivity indices are used in model building to identify unessential variables and reduce the number of uncertain model inputs. By definition, total sensitivity indices are usually larger than first-order indices for nonadditive models and equal for additive models (Chan et al., 1997; Sallteli et al., 2005). The biological part of BTEX bioremediation was simulated using the initial concentrations measured at the gasoline spill site presented by Kao and Wang (2000) to define the model steady state. Two time points, shorter than steady-state time (t_f = 0.25 day), were chosen for the Extended FAST analysis to analyze the possible control shift. The first-order effect indices for the most important parameters at three time points (t_f = 0.01 day, t_f = 0.1 day, t_f = 0.25 day) for metabolites included in the biological part of CBC model are given in Fig. 1, while corresponding total effect indices are given in Fig. 2. The subfigures in Figs. 1 and 2 are as follows: (a) presents carbon sources, (b) nitrogen source, (c) dissolved gases, (d) ions, and (e) particulate components.

FIG. 1.

First-order effect indices for metabolites included in the biological part of the CBC model: (a) carbon sources, (b) nitrogen sources, (c) dissolved gases, (d) ions, and (e) particulate components. CBC, Coupled Biological and Chemical.

FIG. 2.

Total effect indices for metabolites included in the biological part of the CBC model: (a) carbon sources, (b) nitrogen sources, (c) dissolved gases, (d) ions, and (e) particulate components.

Model metabolites

Carbon sources

For better insight into the effect of parameters on the concentration of metabolites, the results are grouped into five clusters (carbon sources, nitrogen sources, dissolved gases, ions, and particulate components). The first group to be analyzed was the carbon sources group (hydrocarbons, S_HC; dissolved methane, S_CH4; organic fermentation products, S_F). The first-order sensitivity indices are presented in Fig. 1a and total effect indices in Fig. 2a. By analyzing the first-order sensitivity indices, it can be observed that the most important parameters and their variance contribution change over time. In the case of hydrocarbons (S_HC) at t_f = 0.01 day and t_f = 0.1 day, the three most important parameters are the half-saturation coefficient for oxygen (K_O), half-saturation coefficient for BTEX (K_HC), and the half-saturation coefficient for fermentation products (K_F) included into the reaction rates of microbial BTEX degradation by sulfate reduction, cleavage of acetate in methanogenesis, and oxidation of the organic fermentation products. Based on the obtained results, it can be concluded that in case of oxygen depletion, sulfate reduction and methanogenesis are the most dominant mechanisms affecting the hydrocarbon concentration. At t_f = 0.25 day (at steady state), a slight control shift can be observed; the reduction coefficient for the anoxic growth rate contribution (η_NO) in denitrification reactions becomes the most important kinetic parameter, followed by the heterotrophic growth rate constant growing on oxygen (μ_O). Thirty-two percent of variance in hydrocarbon output at t_f = 0.25 day is explained by η_NO and μ_O. The sums of the first-order sensitivities for the three points were 0.890, 0.894, and 0.924, indicating the nonadditive model. No significant variation in terms of the total effect indices (Fig. 2a) for hydrocarbon molecules had been detected during the BTEX bioremediation process. At the t_f = 0.01 day and t_f = 0.1 day, there were no influential parameters (total effect index larger than 0.01). By prolonging the simulation until steady state, only ten (η_NO, μ_O, K_O, K_HC, K_F, μ_O,F, K_N, K_NO, K_IS5, and K_Fe2O3) of 29 parameters were influential.

In the case of dissolved methane, the same trend as for hydrocarbons was observed. At t_f = 0.01 day and t_f = 0.1 day, the three most important parameters were K_O, K_HC, and K_F, explaining the 37% variance in the dissolved methane output. Interestingly, it was observed that all 28 parameters contributed to the variance of methane output. The smallest contribution to the variance was observed for the parameters included into the biomass decay reactions. At steady state, the dissolved methane concentration mostly depends (40%) on the aerobic heterotrophic growth rate constant with acetate (μ_O,F) and the half-saturation coefficient for nitrate (K_NO), both included into reactions of acetate consumption (aerobic growth on acetate and methanogenesis reactions). By analyzing the total effect index factor, it was found that at t_f = 0.01 day and t_f = 0.1 day, 13 parameters (included into fermentation, methanogenesis, and biomass decay reactions) were noninfluential, while at t_f = 0.25 day, 19 parameters were noninfluential. This result could be attributed to the fact that there is no high interaction among parameters.

Organic fermentation products (S_F) were included into all reactions of BTEX bioremediation, except fermentation and biomass decay. Analysis of the first-order sensitivities showed that S_F concentration depended mostly on K_SO at t_f = 0.01 day and on K_H2 at t_f = 0.1 day and t_f = 0.25 day (not visible on the presented charts due to ranking according to the first-order effect indices of the first analyzed component). The half-saturation coefficient, K_SO, was included into the kinetics of sulfate reduction on acetate, and by increasing K_S_O, the reactions had slowed down, whereas S_F output concentration increased. On the other hand, the half-saturation coefficient for hydrogen utilization of methanogens (K_H2) was part of the kinetics of the methanogenesis reaction and had the same effect on S_F output concentration as K_SO. At t_f = 0.01 day, K_SO explained 85% of S_F output variability (24 parameters were noninfluential), and at t_f = 0.1 day and t_f = 0.25 day, K_H2 explained 84% (25 parameters were noninfluential) and 74% (22 parameters were noninfluential), respectively. In terms of the total effect index, at t_f = 0.01 day, all parameters except K_SO were considered noninfluential. At t_f = 0.1 day, 25 parameters were noninfluential, but on t_f = 0.25 day, all parameters were influential, indicating the shift in metabolism control.

Nitrogen sources

The biological module of the CBC model of BTEX bioremediation proposes two nitrogen sources; nitrite and nitrate equivalents (S_NO) and ammonium (S_NH). The fast-growing bacteria utilizing oxygen and nitrogen metabolize the BTEX and acetate as a carbon source, but in conditions of oxygen and nitrogen limitation, they ferment BTEX to acetate (Maurer and Rittmann, 2004). By analyzing the results of the Extended FAST on nitrogen sources (Figs. 1b and 2b) in terms of the first-order sensitivities, it could be observed that η_NO (reduction coefficient for anoxic maximum growth rate) and K_H2 (half-saturation coefficient for hydrogen utilization of methanogens) were the most important parameters for both S_NO and S_NH at all three analyzed time points. By increasing the numerical value of η_NO, the denitrification reactions were faster and larger amounts of the nitrogen source were consumed. In contrast, by increasing the numerical value of K_H2, the reaction rate of the methanogenesis had reduced. In the case of S_NH at t_f = 0.01 day, K_H2 contributed to 69% of the output variability, and there were 26 noninfluential parameters. The number of noninfluential parameters was the same at t_f = 0.1 day and t_f = 0.25 day. The sum of the first-order sensitivities for both metabolites at all three analyzed time points was less than one. By analyzing the total effect indices, it was observed that in the case of S_NO, there was no significant variability and all parameters were influential at all three analyzed time points. In the case of S_NH, K_H2 contributed to approximately 60% of S_NH output concentration variability and other parameters could be considered noninfluential (total effect index less than 0.01).

Dissolved gases

BTEX-metabolizing microorganisms can utilize only the components in dissolved form. The amount of available dissolved oxygen determines the metabolic processes; therefore, it is very interesting to define the effect of kinetic parameters on the oxygen concentration. The results of the first-order sensitivities and total effect indices are given in Figs. 1c and 2c. Regarding the first-order sensitivities, it can be noticed that the dissolved oxygen concentration depended mostly on μ_O (contributing to 98% of output variability) at all three analyzed time points. This result was expected because the higher value of μ_O means that the dissolved oxygen is consumed faster and vice versa. It was also observed that there were 26 noninfluential parameters with the possibility of model reductions. In the case of total effect indices, the μ_O was observed to be the most important at t_f = 0.01 day (68%). By prolonging the simulation, at t_f = 0.1 day and t_f = 0.25 day, all other parameters, except K_HC, K_NO, and K_F, became influential (data not shown) on dissolved oxygen concentration. The high value of the total effect index for μ_O at t_f = 0.01 day, and μ_O and μ_Fe at t_f = 0.1 day and t_f = 0.25 day, indicates their indirect influence.

The second dissolved gas included into the model was hydrogen. According to the first-order sensitivities, the dissolved hydrogen output concentration mostly depended on K_O and K_HC at t_f = 0.01 day, on K_HC and K_SO (data not shown) at t_f = 0.1 day, and on K_HC and K_iS4 at t_f = 0.25 day. According to the model balances, dissolved hydrogen was a fermentation product included into microbial BTEX degradation and consumed in methanogenesis. The global sensitivity analysis revealed the effect of the parameters included into the sulfate reduction reaction (K_SO and K_iS3) on the hydrogen concentration at t_v = 0.1 day and t_f = 0.25 day. The first-order sensitivities were around 0.187 (data not shown). This interconnection is probably regulated through hydrogen sulfide, which is included into methanogenesis and the sulfate reduction reaction. Hydrogen sulfide is the inhibitor in both reactions. With regard to the total effect index, there were some changes over time. At t_f = 0.01 day, there were only two influential parameters, K_O and K_HC; at t_f = 0.1 day, all parameters were influential with the maximum influence of K_HC; and at t_f = 0.25 day, K_HC was the most influential, and there were 22 noninfluential parameters.

Ions

Inorganic carbon (H₂CO₃, HCO₃⁻), sulfate, sulfide (SO₄²⁻, HS⁻), and iron (Fe²⁺) were consumed by the microorganisms in ionic form. Sulfate (SO₄²⁻) and sulfide (HS⁻) concentrations change through aerobic and anaerobic sulfate reduction, iron (Fe²⁺) concentration through aerobic and anaerobic iron reduction, while inorganic carbon takes part in BTEX biodegradation, methanogenesis, and oxidation of the organic fermentation product. According to the Extended FAST analysis (first-order sensitivity), the sulfate and iron ions, as well as inorganic carbon concentrations, mostly depended on K_O and K_HC, while the sulfide concentration mostly depended on μ_O at t_f = 0.01 day and t_f = 0.1 day. These results can be explained by the fact that in aerobic conditions, microbial BTEX degradation is a dominating mechanism and the half-saturation constants for oxygen and BTEX affect how fast the oxygen is consumed. A slight change in control was noticed at t_f = 0.25 day, where K_HC and μ_O,F were observed to be the most influential parameters for sulfate, iron ions and sulfide, and K_O and K_HC for inorganic carbon. At steady state (t_f = 0.25 day), the available oxygen is consumed and fermentative metabolism is dominated. Analogous results were obtained for the total effect indices with the largest variability for sulfide, indicating that sulfide concentration mostly varied during the bioremediation process. The results also indicate possible model reduction regarding iron ions due to the fact that at the steady state, there were 23 noninfluential parameters; K_HC and μ_O,F contributed to 57% of variability.

Particulate components

The analyzed model of BTEX bioremediation distinguishes five biomass types: (i) X_BO, fast-growing heterotrophic bacteria; (ii) X_BS, sulfate-reducing bacteria; (iii) X_BF, iron-reducing bacteria; (iv) X_BAc, acetoclastic methanogens; and (v) X_BH2, hydrogen-oxidizing methanogens. In the biological part of the CBC model, iron(III) oxide was also included as a particulate component (X_Fe). Regarding the first-order sensitivities, it could be observed that at t_f = 0.01 day, all parameters were considered influential, but with the small contribution for all five particulate components. K_O and K_HC contributed to approximately 30% of output variability for all five particulate components. At t_f = 0.1 day, there was a slight change in model control, μ_O,F and K_iS3 were shown to be the most important parameters for X_Fe; μ_O,F and K_F for X_BO; and K_HC and μ_O,F for X_BF, X_BS, X_BAc, and X_BH2. In the case of X_Fe, there were 22 noninfluential parameters, while in the case of X_BO, there were 14. The model control, again, slightly changed with the prolongation of the simulation to t_f = 0.25 day. At t_f = 0.25 day, K_HC and μ_O,F were shown to be the most influential on X_Fe; K_HC and K_NO on X_BO; and K_HC and K_O on X_BF, X_BS, X_BAc, and X_BH2. The significant influence of K_HC on almost all biomass types indicates the importance of the microbial BTEX degradation processes. Regarding the total effect indices, it could be observed that at t_f = 0.01 day, there was no significant variability in the total effect indices. It was also observed that there were 23 noninfluential parameters for all five particulate components. The variations of the first-order indices and total effect indices at t_f = 0.1 day and t_f = 0.25 day indicate that particulate component concentrations are not homogeneous during the bioremediation process.

Comparison of groups

Results of the global sensitivity analysis revealed that carbon source concentrations over time mostly depended on the half-saturation coefficient for oxygen (K_O), half-saturation coefficient for BTEX (K_HC), the half-saturation coefficient for fermentation products (K_F), and on the half-saturation coefficient for sulfate (K_SO). In the case of nitrogen sources, η_NO (reduction coefficient for the anoxic maximum growth rate) and K_H2 (the half-saturation coefficient for hydrogen utilization of methanogens) were the most important parameters, while in the case of dissolved gases, μ_O (growth rate constant growing on oxygen) was shown to be the most important.

By comparing the sums of the first-order sensitivities for all analyzed metabolite groups and for all analyzed time points, it could be concluded that the sum of the first-order sensitivities is always less than one. The results indicate that the model is nonadditive and that interactions among model factors do occur. The presence of interactions can have significant implications on the interpretation of the model. In practice, it makes it more difficult to predict the consequences of changing the value of a variable (kinetic parameter value), particularly if the variable interacts with those that are hard to measure or difficult to control. According to Cosenza et al. (2014), the high interactions between factors can be explained with a wide factor variation range. By broadening the variation range of each model factor, the interactions among factors and nonlinear behavior of the model become more apparent. This explanation can be applied to the result obtained in this work because 50% variation of kinetic parameter values around nominal values was tested. The sums of the total effect indices, in all cases, were greater than the corresponding sums of the first-order sensitivities. There is no significant variability between the mentioned sums for specific metabolite at analyzed time points. The sums of the first-order sensitivities were in the range 0.823–0.986, while sums of total effect indices varied in the range 1.896–14.181. The maximum total effect indices sum was obtained for dissolved oxygen, indicating the influence of dissolved oxygen on the process. In the work presented by Yaghmaei (2002), sensitivity analysis was applied on the model that consisted of three nonlinear differential equations describing in situ bioremediation of contaminated soil. The results of applied sensitivity analysis revealed the importance of initial contaminate concentration on the time and mechanism of remediation, which is in consistence with results obtained in this study where the importance of the oxygen concentration on remediation mechanism was demonstrated.

Model reaction

For better insight into the complex inter-regulation and cumulative effects of kinetic parameters on the bioremediation process, the Extended FAST analysis was applied on model reactions. The reactions of the biological module of the BTEX bioremediation model were divided into four groups (microbial BTEX degradation, biomass decay, methanogenesis, and oxidation of organic fermentation products), and the results are presented accordingly in Figs. 3 and 4.

FIG. 3.

First-order effect indices for reactions included in the biological part of the CBC model: (a) microbial BTEX degradation, (b) biomass decay, (c) methanogenesis, and (d) oxidation of organic fermentation products.

FIG. 4.

Total effect indices for reactions included in the biological part of the CBC model: (a) microbial BTEX degradation, (b) biomass decay, (c) methanogenesis, and (d) oxidation of organic fermentation products.

Microbial BTEX degradation

The microbial BTEX degradation group included aerobic growth, denitrification, iron reduction, sulfate reduction, and fermentation. In the case of first-order sensitivity, the aerobic growth reaction mostly depended on K_O, K_HC, and K_F (half-saturation constants for oxygen, BTEX, and fermentation products) at all three analyzed time points. These three parameters contributed to approximately 35% of output variability. The mentioned constants were included into the kinetic expression describing the aerobic growth rate. All three parameters had the same effect; the aerobic growth rate decreased with the increase in the value of an individual parameter. An analogy in the total effect indices was observed. There was no significant variance in the total effect indices and all parameters were influential.

In conditions of oxygen absence, microorganisms use available electron acceptors and perform metabolic reactions. As a result, the anoxic growth reaction was mostly influenced by μ_O, μ_O,F, and K_O at t_f = 0.01 day. The faster the oxygen consumption (larger μ_O value), the sooner the anoxic metabolism is activated. By prolonging the simulation to t_f = 0.1 day and t_f = 0.25 day, the Extended FAST analysis revealed that denitrification depended mainly on K_N, the half-saturation coefficient for nitrogen, and biomass growth (approximately 80% of output variability). Regarding the total effect indices, the highest variability was observed for the denitrification reaction over time.

Sulfate-reducing bacteria and iron-reducing bacteria are able to reduce sulfate and solid iron, respectively, and use them as an electron acceptor. These microorganisms use BTEX and acetate as electron donors and carbon sources. The results of the global sensitivity analysis suggested that the most important parameters for the iron reduction reaction rate change over time. At t_f = 0.01 day, the most important parameter was K_NO; at t_f = 0.1 day, μ_O,F; and at t_f = 0.25 day, K_O. The aforementioned half-saturation constants were included into the kinetic expression describing iron reduction, both having an inhibitory effect. The importance of μ_O,F could be explained by the competition between aerobic growth on acetate and iron reduction. Perhaps aerobic growth on acetate is a preferable metabolism. On the other hand, μ_O was the most influential parameter for the sulfate reduction reaction rate at t_f = 0.01 day and t_f = 0.1 day (∼70% contribution), while at t_f = 0.25 day, K_O was the most influential. The sum of the first-order sensitivities for the mentioned three time points varies from 0.821 to 0.882, indicating a nonadditive model. No relevant variations of total effect indices had been detected during the simulation time for both iron reduction and sulfate reduction reaction rates.

In the fermentation reaction, BTEX is transformed to acetate and hydrogen. According to the first-order sensitivity, fermentation was mostly influenced by K_O and μ_O,F at t_f = 0.01 day, by K_O and K_N at t_f = 0.1 day, and by K_O and K_HC at t_f = 0.25 day. At t_f = 0.01 day, K_O and μ_O,F contributed to 57% of output variability, and there were 22 noninfluential parameters. By increasing the values of K_N and K_O, the fermentation reaction had slowed down. Total effect indices showed the same trend with no significant variability.

Biomass decay

The BTEX bioremediation model also describes microorganism decay through reactions of biomass decay. Each reaction describing five types of biomass decay was analyzed separately. Lysis of X_BFe, X_BAc, and X_BH2 was influenced (∼99%) only with one kinetic parameter K_NO, K_F, and η_NO, respectively, at all three analyzed time points. All other 27 kinetic parameters were considered noninfluential on the mentioned reactions. The sum of the first-order sensitivity for the mentioned reactions at the analyzed time points was close to 1 (0.99), indicating the unimportant effect related to interactions. K_NO is considered an inhibition constant in the iron reduction reaction performed by X_BFe biomass. By increasing the K_NO value, the mentioned reaction had slowed down and the biomass lysis activated. K_F had the same effect on the reaction of acetate cleavage performed by X_BAc and η_NO on the denitrification reaction. The reaction describing X_BO and X_BSO biomass decay was mostly sensitive to Ko and μ_Fe. The total effect indices indicate the same trend as first-order sensitivity. Regarding the total effect indices, there was possible model reduction with respect to lysis of X_BFe, X_BAc, and X_BH2.

Methanogenesis

Methanogenesis includes two reactions, cleavage of acetate and methanogenesis by oxidation of hydrogen and reduction of inorganic carbon. During the bioremediation process, there was a change in the most influential parameters that affected methanogenesis reactions. At the beginning of the process, both reaction rates were mostly influenced by μ_O (sensitivity coefficient = 0.469 and 0.242), and there were 20 noninfluential parameters. μ_O affects the oxygen consumption rate and indirectly acetate production. The faster the oxygen is consumed, the faster the fermentation is activated and, respectively, the faster is the acetate cleavage reaction. At t_f = 0.1 day, μ_O,F, K_NO, and K_HC were shown to be the most influential, indicating the competition between methanogenesis and acetate oxidation. Interestingly, at the steady state, the FAST analysis revealed that all kinetic parameters had approximately the same influence on the analyzed reactions. The same effect was noticed for the total effect indices. At steady state, there were minimum variations in the total effect index values for both reactions. Some variations of the total effect indices were observed at t_f = 0.1 day and t_f = 0.25 day for reaction of acetate cleavage.

Oxidation of organic fermentation products

During BTEX bioremediation, the microorganisms oxidize acetate using the same electron acceptors as for the BTEX through four reactions: (i) aerobic growth on acetate, (ii) anoxic growth on acetate, (iii) iron reduction on acetate, and (iv) sulfate reduction on acetate. Global sensitivity analysis regarding the first-order sensitivities revealed that there was a slight change in process control. At t_f = 0.01 day and t_f = 0.1 day, K_NO was shown to be the most important parameter for all four reactions, contributing to approximately 30% of output variability. As mentioned previously, it had an inhibitory effect on the mentioned reactions. At steady state (t_f = 0.25 day), the individual parameter contributions to variance were low, and the most influential parameter was observed to be K_O. The half-saturation coefficient for oxygen had an inhibitory effect on anoxic growth on acetate, iron reduction on acetate, and sulfate reduction on acetate. This indicates that aerobic growth on acetate is a preferable process.

In the case of total effect indices, significant availabilities could be observed at t_f = 0.01 day. There were 26 noninfluential parameters for aerobic and anoxic growth on acetate, and 24 noninfluential parameters for iron and sulfate reduction on acetate. Interestingly, at t_f = 0.1 day and t_f = 0.25 day, all parameters were influential, with no significant variability.

Comparison of groups

By analyzing the results of the global sensitivity analysis for model reactions, it can be seen that aerobic growth reaction mostly depended on K_O, K_HC, and K_F (the half-saturation constants for oxygen, BTEX, and fermentation products). It was also observed that reactions included into methanogenesis mostly depend on μ_O (growth rate constant growing on oxygen), while in the case of the oxidation of organic fermentation products, K_O was the most influential parameter.

Like for the models of metabolites, the sums of the first-order sensitivities for all analyzed reactions at all analyzed time points were less than one. The sums of the first-order sensitivities ranged from 0.63 for the reaction of the lysis of X_BO at t_f = 0.01 day to 0.997 for the reaction of the lysis of X_BAc at all three analyzed time points. Correspondingly, the sums of the total effect indices ranged from 2.162 for the reaction of fermentation at t_f = 0.01 day to 11.124 for the reaction of acetate cleavage at t_f = 0.01 day.

Conclusions

In this work, the Extended Fourier Amplitude Sensitivity Test was applied for analysis of the importance of kinetic parameters on 18 metabolites and 16 reaction rates of the biological part of the integrated BTEX bioremediation model on three time points. The Extended FAST results of two sensitivity indices (first-order effect and total effect indices) were analyzed. According to the sum of the first-order effect indices for both metabolites and reactions at three analyzed time points, it can be concluded that the model is nonadditive, but the relationships between model variables are not strong (sums of the first-order sensitivities for metabolites are in the range 0.823–0.986 and, for reactions, are in the range 0.63–0.997). Comparing all calculated first-order sensitivities, it can be observed that the most influential parameters on both metabolites and reactions are K_HC and K_O. The obtained results indicate that a BTEX bioremediation process described with the set of equations given in Table 1 is regulated depending on the dissolved oxygen concentration. Although a different set of equations describing the BTEX bioremediation process would demand a new global sensitivity analysis, the obtained results offer a good starting point for model analysis planning.

Footnotes

Author Disclosure Statement

No competing financial interests exist.

References

Almquist

, Cvijovic

, Hatzumanikatis

, Nielsen

, and Jirstrand

(2014). Kinetic models in industrial biotechnology-Improving cell factory performance. Metab. Eng. 24, 38.

Benekos

I.D.

, Shoemaker

C.A.

, and Stedinger

J.R.

(2007). Probabilistic risk and uncertainty analysis for bioremediation of four chlorinated ethenes in groundwater. Stoch. Environ. Res. Risk. A.21, 375.

Chan

, Saltelli

, and Tarantola

(1997). Sensitivity analysis of model output: Variance-based methods make the difference. Presented at the 1997 Winter simulation conference, Atlanta, Georgia, USA, December 7–10.

Cosenza

, Mannina

, Vanrolleghem

P.A.

, and Neumann

M.B.

(2014). Variance-based sensitivity analysis for wastewater treatment plant modelling. Sci. Total Environ. 470, 1068.

Essaid

H.I.

, Bekins

B.A.

, Godsy

E.M.

, Warren

, Baedecker

M.J.

, and Cozzarelli

I.M.

(1995). Simulation of aerobic and anaerobic biodegradation processes at a crude oil spill site. Water Resour. Res. 31, 3309.

Fang

, Gertner

G.Z.

, Shinkareva

, Wang

, and Anderson

(2003). Improved generalized Fourier amplitude sensetivity test (FAST) for model assessment. Stat. Comp. 13, 221.

Geng

, Boufadel

M.C.

, and Wrenn

(2013). Mathematical modeling of the biodegradation of residual hydrocarbon in a variably-saturated sand column. Biodegradation. 24, 153.

Geng

, Boufadel

, Personna

, Lee

, Tsao

, and Demicco

(2014). BioB: A mathematical model for the biodegradation of low solubility hydrocarbons. Mar. Pollut. Bull. 83, 138.

Gombert

A.K.

, and Nielsen

(2000). Mathematical modelling of metabolism. Curr. Opin. Biotech. 11, 180.

10.

Gujer

, Henze

, Mino

, and Van Loosdrecht

M.C.M.

(1999). Activated sludge model No. 3. Water Sci. Technol. 39, 183.

11.

Helton

J.C.

, Johnson

J.D.

, Sallaberry

C.J.

, and Storlie

C.B.

(2006). Survey of sampling-based methods for uncertainty and sensitivity analysis. Reliab. Eng. Syst. Safe. 91, 1175.

12.

Henze

, Grady

C.P.L.

, Gujer

, Marais

G.V.R.

, and Matsuo

(1987). Activated sludge model No. 1. IAWPRC Scientific and Technical Report No. 1. IAWPRC, London.

13.

Henze

, Gujer

, Mino

, Matsuo

, Wentzel

M.C.

, and Marais

G.V.R.

(1995). Activated sludge model No. 2. IAWO Scientific and Technical Report No. 3. IAWO, London.

14.

Henze

, Gujer

, Mino

, Matsuo

, Wentzel

M.C.

, Marais

G.V.R.

, and Van Loosdrecht

M.C.M.

(1999). Activated sludge model No. 2d, ASM2d. Water Sci. Technol. 39, 165.

15.

Homma

, and Saltelli

(1996). Importance measures in global sensitivity analysis of nonlinear models. Reliab. Eng. Syst. Saf. 52, 1.

16.

Iooss

, and Lemaître

(2015). A review on global sensitivity analysis methods. In Meloni

and Dellino

, Eds., Uncertainty Management in Simulation-Optimization of Complex Systems: Algorithms and Applications. New York: Springer US, p.101.

17.

Jia

, Stephanopoulos

, and Gunawan

(2012). Incremental parameter estimation of kinetic metabolic network models. BMC Syst. Biol. 6, 142.

18.

Kao

C.M.

, and Wang

C.C.

(2000). Control of BTEX migration by intrinsic bioremediation at a gasoline spill site. Water Res. 13, 3413.

19.

Kent

, Neumann

, Kummer

, and Mendes

(2013). What can we learn from global sensitivity analysis of biochemical systems?. PLoS One. 8, 1.

20.

K.Y.

, Xu

, Cawley

W.A.

, Colapret

J.A.

, Bonner

J.S.

, Ernest

, and Verramachaneni

P.B.

(1994). Field test and mathematical modelling of bioremediation of an oil-contaminated soil-Part 2: Mathematical modelling. Waste Manage. 14, 571.

21.

Maurer

, and Rittmann

B.E.

(2004). Formulation of the CBC-model for modelling the contaminants and footprints in natural attenuation of BTEX. Biodegradation. 15, 419.

22.

Nassar

M.A.

(2011) Multi-parametric sensitivity analysis of CCHE2D for channel flow simulations in Nile River. J. Hydro-Env. Eng. 5, 187.

23.

Nicol

J.P.

, Wise

W.R.

, Molz

F.J.

, and Benefield

L.D.

(1994). Modeling biodegradation of residual petroleum in a saturated porous column. Water Resour. Res. 30, 3313.

24.

Peskov

, Moqilevskaya

, and Demin

(2012). Kinetic modelling of central carbon metabolism in Escherichia coli. FEBS J. 18, 3374.

25.

Reed

M.C.

, Thomas

R.L.

, Pavisic

, James

S.J.

, Ulrich

C.M.

, and Nijhout

H.F.

(2008). A mathematical model of glutathione metabolism. Theor. Biol. Med. Model. 5, 8.

26.

Saltelli

, and Bolado

(1998). An alternative way to compute Fourier amplitude sensitivity test (FAST). Comput. Stat. Data Anal. 26, 445.

27.

Saltelli

, Annoni

, Azzini

, Campolongo

, Ratto

, and Tarantola

(2010). Variance based sensitivity analysis of model output. Design and estimator for the total sensitivity index. Comput. Phys. Commun. 181, 259.

28.

Saltelli

, Tarantola

, and Campolongo

(2000). Sensitivity analysis as an ingredient of modelling. Stat. Sci. 15, 377.

29.

Saltelli

, Tarantola

, and Chan

(1999). A quantitative, model independent method for global sensitivity analysis of model output. Technometrics. 41, 39.

30.

Schirmer

, Butler

B.J.

, Roy

J.W.

, Frind

E.O.

, and Barker

J.F.

(1999). A relative-least-squares technique to determine unique Monod kinetic parameters of BTEX compounds using batch experiments. J. Contam. Hydrol. 37, 69.

31.

Schirmer

, Molson

J.W.

, Frind

E.O.

, and Barker

J.F.

(2000). Biodegradation modelling of a dissolved gasoline plume applying independent laboratory and field parameters. J. Contam. Hydrol. 46, 339.

32.

Schwieger

(2004). Variance-based sensitivity analysis for model evaluation in engineering surveys. Presented at INGEO 2004 and FIG Regional Central and Eastern European Conference on Engineering Surveying, Bratislava, Slovakia, November 11.

33.

Surajudeen

(2012). Mathematical modelling of bioremediation of soil contaminated with spent motor oil. J. Emerg. Trends Eng. Appl. Sci. 3, 654.

34.

Vanlier

, Tiemann

C.A.

, Hilbers

P.A.J.

, and van Riel

N.A.W.

(2013). Parameter uncertainty in biochemical models described by ordinary differential equations. Math. Biosci. 246, 305.

35.

Wayman

J.A.

, and Varner

J.D.

(2013). Biological systems modelling of metabolic and signalling networks. Curr. Opin. Chem. Eng. 2, 365.

36.

Wiechert

(2002). Modelling and simulation: Tools for metabolic engineering. J. Biotech. 94, 37.

37.

Yaghmaei

(2002). Mathematical modelling of natural in situ bioremediation to estimate initial contaminate concentration effect. IJE Trans. 2, 105.

38.

Zhao

, Scheibe

T.D.

, and Mahadevan

(2011). Model based analysis of the role of biological, hydrological and geochemical factors affecting uranium bioremediation. Biotechnol. Bioeng. 108, 1537.

39.

Zhao

, Scheibel

T.D.

, Fang

, Lovley

D.R.

, and Mahadevan

(2010). Modeling and sensitivity analysis of electron capacitance for Geobacter in sedimentary environments. J. Contam. Hydrol. 112, 30.