Variance‑based global sensitivity analysis of a multi‑population, single‑chamber microbial fuel cell operating in continuous flow mode at steady state

Microbial fuel cells (MFCs) are environmentally friendly devices which are used to convert chemical energy in organic wastes to electrical energy. MFCs have a strong non-linearity that requires a very sophisticated controlling system. Consequently, this makes optimization and performance study of MFCs a difficult process. For better estimation of the constants used for optimization of MFCs, global sensitivity analysis is performed. The global sensitivity method based on Sobol’s indices coupled with Monte Carlo simulations was applied on multi-population, single-chamber MFC operating in a continuous flow at steady state for the first time. In this paper, first-order and total-order sensitivity indices were used to visualize the impacts associated with six main parameters resulted from the maximization of power density using Matlab. Such parameters are maximum anodophilic-specific growth rate, half-rate constant of anodophilics, curve steepness factor, mediator half-rate constant, number of electrons transferred per mole mediator and decay rate constant of anodophilic bacteria. The results showed that the curve steepness factor has almost no impact on the power density of MFC. While all other studied, factors are sensitive parameters that impact the power density of MFC. It is worth mentioning that maximum anodophilic growth rate and the number of electrons transferred per mole of mediator are the most sensitive parameters that affecting the power density production having total indices of 0.74 and 0.624, respectively. While the half-rate constant of anodophilics, mediator half-rate constant and decay rate constant of anodophilics have almost similar impact by having total-order indices of 0.127, 0.144 and 0.192, respectively. The findings herein are critical in understanding and further model improvement of microbial fuel cells as the most impacting parameters on MFC power density can be optimized further to reduce uncertainty associated with the experimental parameters in the model.


Introduction
Microbial fuel cells (MFCs) are renewable bio-electrochemical cells that use microbial activity to convert chemical energy that is reserved in wastewater to an electrical energy [1,2].Therefore, MFCs are considered as climateneutral technology for wastewater treatment and energy production [3].This cell contains a bio-anode and/or a biocathode which are separated from one another by a proton exchange membrane.Both compartments of the cell are connected with an outer circuit which allows the electrons to flow through [4].The anode compartment oxidation of the organic material takes place and most of the MFCs use the organic substances as electron donors.As a consequence of oxidation reactions carbon dioxide and/or methane, protons and electrons are produced [5].Eventually, the electrons produced from the oxidization process (electrical energy) travel through the outer circuit producing a flux of electrons.Electrons are transferred directly to the cathode or through a mediator [6].A charge balance occurs in MFCs due to the transfer of electrons through the outer circuit and transfer of protons through the membrane from the anode to the cathode.Based on the electron acceptor in the cathode compartment, reduction reaction occurs.If oxygen is the electrons acceptor; then, water is produced from the reduction reaction [7].This technology has been developed as a result of the huge attentions on sustainable processes [1,2].MFCs as a sustainable process have potentials to reduce water pollution, increase the ecological restoration and generate power.However, the industrial applications of MFC are not widespread.The main reason behind this is the lower power density generation compared to other fuel cells [8].Therefore, optimization of MFC mathematical models using computational methods is required to save laboratory efforts and figure out the optimal conditions at which the power density can be maximized [9].Many studies were conducted on MFC optimization throughout the years.Abu-Reesh, I. conducted a study on a MFC for both single-and multi-objective optimization for a dual chamber which operates at steady state [10].In this study, a mathematical model-based optimization was performed to assess the performance of MFC.The performance is assessed by power density (PD), current density production (CD) and substrate removal efficiency (SRE) using Matlab.Optimizing the single objective of the previous functions were 2.04 W/m 2 , 11.08 A/m 2 and 73.6%, respectively.While multi-objective optimization provided pareto-optimal solutions.Another study was accomplished by Almatouq, A. and Babatunde, O. in which they optimized the operational conditions of a mediator-less, dual chamber MFC for both concurrent power production and phosphorus removal [11].They used mathematical model-based optimization coupled with full-factorial design and central composite designs for the design of experiments.Results showed that operating at high chemical oxygen demand and cathode aeration flow rate improved both objectives of the study [11].Boas, J. et al. [12] optimized a single-chamber MFC using a pure culture of Lactobacillus pentosus and dairy industry wastewater to study the impact of the operating factors and the design of the cell on its performance.The performance was evaluated by power output and the COD removal rate.The studied operating factors were batch cycle and yeast extraction concentration, anode area and membrane thickness.The maximum achieved power density was of 5.04 ± 0.39 mW/m 2 given a COD removal rate between 42 and 58%.This was achieved at anode electrode area of 61 cm 2 , a batch cycle of 48 h and 50 mg/L of yeast extract [12].On the same track, the selection of parameters to be used in the computational methods influence the performance of MFCs [13].Sensitivity analysis is an important tool to reveal the factors that have most impact on the output function [14].This assists to evaluate the accuracy and quantify the uncertainty in the estimated model empirical parameters [15].Two main types of sensitivity analysis can be carried out in studying any model.Local sensitivity analysis (LSA) and global sensitivity analysis (GSA) [16].LSA is carrying out the sensitivity over one point in any selected parameter hyperspace while varying the selected parameter and held other parameters fixed [17], whereas GSA scans the entire defined range of all of the parameters of the study [18,19].Sensitivity analysis assists the extent to which each input parameter impacts the uncertainty on the output functions [20].The Sobol-Saltelli method which is a variance-based method is used to evaluate the relative contribution of each of the parameters to the variance of the output function.This method excluding the interrelations with other parameters in the first-order indices and including all the impacts and interrelation in the total order indices while taking into consideration the non-linearity [20].Yankai, Y. et al. [21] expound the global sensitivity analysis based on variance method on a dual chamber, single population under steady-state mode MFC.They studied the impact of eight parameters on the power density by conduct Sobol's indices.The results show that the cathodic charge transfer coefficient is the most crucial and impacting factor on the output, while the electrical conductivity of the aqueous solution is the least effective.By using the same methodology, a comprehensive global sensitivity analysis using variancebased method is carried out over [7] MFC model.In this work, first-and total-order sensitivity indices of six chosen parameters for optimization of a multi-population, singlechamber MFC with biofilm operating at steady-state mode were estimated.In addition, Monte Carlo simulations were used for predicting the behaviour of each of the parameters on the power density output [14].

MFC model description
The mathematical model of MFC used in this study was developed by Pinto, R. P. et al. [7] which can be implemented for optimizing productivity.The model of Pinto, R. P. et al. [7] is a simple, one dimensional, unified model which can be applied on bioelectorchemical cells to optimize either the power generation of MFC or hydrogen generation of MEC.The model combines biochemical reactions in the anode compartment, ordinary differential equations of mole balances of substrate, anodophilic bacteria and methanogenic bacteria, Monod equation for growth kinetics and Butler-Volmer equation.The model considered the co-existence of microbial population of anodophilic (electricigenic) and methanogenic bacteria in a single-anode chamber, air cathode membrane-less that is operating in a continuous flow mode.Acetate was used as the only substrate in the system because of its complete degradation compared to other substrates and both microorganisms are competing for it [1].A uniform distribution of both carbon source and microbial population is proposed in the anode compartment with ideal mixing.The anodophilic microorganisms are the only electron producer in the model [5].The anodophilic are placed attached to the biofilm in the anode compartment, while the methanogens are both attached to the biofilm and suspended in the anode compartment.Intracellular mediators (NAD + /NADH) are assumed to be involved in the charge transfer process from the carbon source (acetate) [4].The proposed model focuses on the bio-electrochemical reactions in the anode compartment and its associated kinetics are considered as the rate-limiting of MFC [7].A schematic diagram of the used MFC is shown in Fig. 1.A combination of electrochemical and biochemical reactions is taking place in MFC.Electrochemical reactions deal with electron production and transfer from the oxidation of organic matter, while biochemical reactions deal with the microorganisms' metabolism and substrate reduction [5].The main model assumptions and considerations are listed as follows [7]: 1. Operating conditions are constant in this study following a temperature of 25 °C and 1 atm.2. The biofilm effect is considered in the design equations and the substrate gradient through the biofilm is neglected.3. Direct and indirect (NAD + /NADH) mechanisms for electron transfer from the carbon source (acetate) to the anode electrode are assumed.4. A balance between the microbial poppulation in the biofilm is assumed.5. Stationary phase is assumed; therefore, equilibrium will be achieved between biofilm growth and washout so that the biofilm will reach its steady state thickness and washout will be equal to the net biofilm growth.Washout is the state at which the dilution rate of the system which is equal to the volumetric flow rate over the volume of the reactor becomes greater than the maximum specific growth rate.Then, the rate of removing cells from the system will be greater than the specific growth rate of the cells.This causes cells to be removed from the reactor.
Substrate degradation and electron generation reactions by anodophilic bacteria are illustrated in Eq. (1). Figure 2 provides a schematic diagram for the anode compartment presenting all the associated reactions [9].
In Fig. 2, anodophilic bacteria are represented in blue and they are attached only to the biofilm, they oxidize the substrate (acetate) which releases electrons that are transferred directly or indirectly by mediators as shown in Eq. (2).Equation (3) shows the reduction reaction of electron acceptor (oxygen) in the cathode compartment.While methanogenic bacteria are represented in green in Fig. 2 and they are both attached to the biofilm and suspended in the bulk.The substrate conversion by methanogenic bacteria is represented by Eq. ( 4).The conceptual reactions for oxidizing the substrate are as follows: If the substrate is acetate, then the reaction over the anode and cathode for anodophilic microorganisms is shown as follows [7]: Anode: Air cathode: While for methanogenic bacteria, the anodic reaction is as follows:

Anode compartment material balances
The model balances equations used in this paper are represented as follows: where S and S 0 are the substrate concentration and feed substrate concentration respectively (mgS/L); X m and X a are methanogenic and electricigenic microorganism concentrations respectively (mgX/L); t is the time (d); D is the dilution rate (1/d) which is defined as (D = F in /V A ); V A is the volume of anode (L); F in inlet flow rate; max is the maximum specific growth rate (1/d); q max is the maximum substrate consumption rate (mgS/mgX.d)and k d,a and k d,m are decay rates constants for anodophilics and methanogens respectively (1/d).
In this model, the biofilm impacts in the anode compartment are taken into consideration based on two-phase biofilm growth model.For the effects of biomass formation and retention in the layer, (1/d) which is the biofilm retention constant that is calculated below, assuming the cell is at stationary phase and the biofilm reached its steady-state thickness.The maximum attainable biomass concentration (X max ) is hold in the biofilm at this phase.

Intracellular mediators material balance
Intracellular mediators are used for electrons transport; therefore, they exist in both oxidized and reduced forms.
The material balance of the intracellular mediators per anodophilic microorganisms is as follows: Given that: where M ox is the oxidized mediator fraction per anodophilic microorganism (mgM/mgX a ); M red is the reduced mediator fraction per anodophilic microorganism (mgM/mgX a ); M T is the total mediator fraction per microorganism (mgM/mgX a ) which is a constant known value tabulated in the nominal values table.Y M is the mediator yield (mgM/mgS); m is the number of electrons transferred per mol of mediator and is mediator molar mass (mgM/molM).

Kinetics equations
where μ m andμ a are the specific growth rates of metha- nogens and anodophilics microorganisms, respectively (1/d); max,m and max,a are the maximum specific growth rates for both methanogens and anodophilic microorganisms, respectively (1/d); q m andq a are substrate consump- tion rate by methanogens and anodophilic microorganisms, respectively (mgS/mgXd); q max,m andq max,a are the maximum substrate consumption rate by methanogens and anodophilics microorganisms, respectively (mgS/ mgXd).

Electrochemical equations
Theoretically, the MFC output voltage (E MFC ) can be cal- culated by open circuit potential ( E OCV ) with removing the losses which are mainly concentration losses ( η conc ) , ohmic losses (η ohm ) and activation losses (η act ).
Activation losses can generally be neglected, as from the optimization the internal resistance is aimed to be equal to the external resistance.This simplifies the equation of MFC voltage to: where the open circuit potential, concentration losses, internal resistance and MFC current are calculated theoretically as follows: where is a constant and it is approximately equated to zero; then, the current of MFC is reduced to: The power density of MFC can be calculated as: where E OCV is the open circuit potential (V); R max andR min are the maximum and minimum internal resistance ( Ω ); is the mediator molar mass (mgM/molM); E MFC is the MFC voltage (V); conc is the concentration losses (V); R is the universal gas constant (J/K.mol); K r is the curve steepness factor (L/mgX); P MFC is the power density of MFC (mW/L).Table 1 below shows the nominal values of all constants used in the optimization model.

Sensitivity analysis
Measuring the quality of the parameters which can significantly impact the performance of MFC is crucial.This can be carried out by the means of performing sensitivity analysis on the system.The analysis method is based on mathematical rules by generating values of the factor; then, calculating the value of the output corresponds that mathematical model.If a small change in a certain input parameter has resulted in a relatively large change in the output, it is said that the output of the model is sensitive to this input parameter [22,23].Ideally, both uncertainty and sensitivity analysis of a model should be run tandemly [24,25].Figure 3 shows the interrelationship between uncertainty, sensitivity analysis and modelling.
Figure 3 shows that a model which has input data (errorfree for simplicity) from which the parameters are estimated using the specified model, then after the estimation; the best parameter values are considered.At this stage, the model is considered as a true model and it is proceeded to the sensitivity analysis [18].Sobol (1998) defined sensitivity analysis as the method of assigning variation of model outputs to different changes in the input parameters [26].Global sensitivity analysis is computed in this paper.

Global sensitivity analysis (GSA)
Global sensitivity analysis (GSA) is a variance-based analysis in which all the input factors are varied simultaneously over the whole range of interest allocating the uncertainty in the model output to the uncertainties of each of the input parameters.GSA is an important tool for decision-making in complicated systems such as systems in environmental sciences in which many input factors have a strong interaction with one another [21].Evaluating the individual impact of each parameter while the other inputs are constant produces nonaccurate results for the model output.In GSA, while evaluating the impacts of inputs on the model output, it can calculate both single impacts of individual factors called first-order indices and multiple impacts due to factors interactions with one another and it is called total-order indices [21].MFC is one of the convoluted systems in which many factors simultaneously affect the performance of the cell, which can cause multiple sources of uncertainty.To figure out the effect of these uncertain factors on the cell and on the population dynamics, GSA is performed.Therefore, GSA is becoming an essential tool for the assessment of environmental systems [8,23].

Sensitivity index estimation
In this study, the objective optimization function (output of the model) is the power density of MFC-based model [7].Based on the proposed model, the variance of some constants defined in Table 1 on power density of MFC was studied.Global sensitivity analysis based on variance decomposition as a screening method which was suggested by Sobol is carried out.The variance-based screening method (Sobol's) is better than the Morris method in terms of precise indices taking into consideration the interrelationships between the input parameters [27].Screening is the process of sieving the factors and identifying the ones which are the least important in their impact on the output function and fix them within their feasible range [28].In this part, the GSA is implemented based on variance-based screening and the sensitivity indices are approximated based on Sobol's method [20].A critical point in such studies is the sample size selection.MFCs are complex systems; therefore, the sample size should not be too small so that the accuracy is lost or not too large and face difficulty to run the simulation.Based on the model input parameters (M) and the number of model evaluations (N) and the sample size (n), N is set as a function of M and n.Therefore, choosing the number of system's total evaluations can lead to the number of the sample size.Some methods use the direct relation between N and n and equate them (i.e.N = n) as it is used in this study [8,29].

Sobol' sequences
Sobol indices is a method of estimation based on the variance calculations of the inputs of a mathematical model.Sobol' technique is basically an example of quasi-random low-discrepancy sequences [30].These sequences are used for GSA extensively for their simplicity and accuracy.This was proposed first by the Russian mathematician Ilya M. Sobol in 1998 [26].The detailed GSA based on Sobol's method is as follows [23,31]: where Y output is the model output; f is the model function and x 1 andx 2 are the chosen parameters to study their influ- ence on the output function.
Based on the variance study, first-order and total-order coefficients are calculated as follows: 1-First-order indices: 2-Total indices can be described as: where ( x i ) is a matrix of all parameters; (x ∼ i ) is the matrix of all parameters but x i ; ( Var x i ) is the variance of argument (•) over x i ; ( E xi ) is the mean of argument (•) over x i ; (E x∼i ) is the mean of argument (•) and (Var x∼ i ) is the vari- ance of argument (•).
Variance and mean value are calculated as follows: where (A and B) are matrices of size N × k.

GSA methodology
The methodology followed for the GSA in this study is as follows: the parameters to be studied were selected based on their expected impact on the power density (output function).After that the probability density function (PDF) of each parameter was predicted for further Monte Carlo analysis.One key of getting good results is the choice of PDF ( 23) type of each parameter.This represents a critical issue which needs an experience or knowledge in the field especially in MFC complicated models.Power density has a normal distribution with most of the factors with time [1]; therefore, normal distribution is used for all of the parameters as a PDF.For global sensitivity analysis, Sobol's method was used to compute first-and total-order indices as the main method of the study.While double loop, Morris, Sobol and FAST methods were applied to compare different methods of calculating total indices.Finally, Monte Carlo simulations were run for each of the parameters to figure out the impact of each parameter on the power density.Input sample generation of size (n = 2000) is created in Matlab for each parameter.The methodology is shown in Fig. 4.

Results and discussion
In this section, the results of GSA which is applied on Pinto's MFC model are analysed and discussed.The impact of input parameters is evaluated under steady-state conditions on power density of the cell.Most representative parameters have been chosen which may have an impact on the performance of MFC represented by the power density.Parameters are shown in Table 2 below.The studied parameters are decay rate constant of anodophilics ( k d,a ), half-rate con- stant of anodophilic microorganism ( K s,a ) , curve steepness factor ( K r ), mediator half-rate constant k M , anodophilic maximum specific growth rate ( max,a ) and number of electrons transferred per mole of mediator (m).The selection of parameters was based on the fact that anodophilics are the only organisms which contribute to electron production and transfer, thus impact the output power density production.The transfer of the electron is via two methods, direct transfer and indirect transfer through the mediator, in turn the half-rate of mediators constant is considered.Also, the number of the electrons transferred per mole of mediator as there are different mediators that can carry different amount of electrons.The initial value used in the simulation initiation is the mean value of each of the parameters which is the proposed value in the model as shown in Table 1, the defined boundaries of both minimum and maximum values in Matlab are shown in Table 2.
First-order and total-order sensitivity indices based on the variance (Sobol's) method are calculated for all parameters in response to their impacts on the power density output using Matlab.This is shown in Fig. 5 below using a bar graph representation.Also, Table 3 shows the values of the indices and the difference between them, as it is expected the total indices are greater than the first order.The indices reflect the degree of impact of the selected input parameters on the output function which is the power density of MFC.First-order indices reflect the level of impact of each parameter stand-alone on the output function, while the total-indices take into consideration the cross-interaction of all the variables with one another and their effect on the output function (power density of MFC).The closer the value of the sensitivity coefficient (index) of the parameter to 1, the more sensitive it is to the output function.
The results of Fig. 5 are shown in Table 3 as stated previously with showing the difference between the first and  total-order indices.From Fig. 5, the sequence of sensitivity of parameters based on their impact on power density can be determined as follows: the maximum specific growth rate of anodophilic bacteria ( max,a ) , number of electrons transferred per mole of mediator (m), decay rate constant of anodophilic bacteria ( k d,a ), half-rate constant of anodo- philic bacteria ( K s,a ), mediator half-rate constant ( k m ) and curve steepness factor ( K r ) .The difference between the first- order and total-order indices indicated that the factor has a larger impact on the power density when interacts with other parameters.Among the parameters, max,a has the greatest impact on the power density.This means that a small change in the value of this parameter will result in a large change in power density compared to other factors, while K r has very weak impact on the power density as a stand-alone factor with a negligible impact on the power density when interacting with other parameters.Decay rate constant of anodophilic bacteria ( k d,a ), half- rate constant of anodophilic bacteria ( K s,a ) and mediator half-rate constant ( k m ) have near but fluctuating impacts on the power density.However, the difference between the first-and total-order impacts for all the three factors is almost the same.Therefore, they have some variety in their impact; the lumped impact is almost the same on the power density function.Based on the studies conducted by Batstone, D. J. et al. [32] and Marcus, A. [33], the value of decay rate constant of the biomass is 2% of the maximum specific growth rate of the biomass.Therefore, the impact of this factor will be less than that of the maximum specific growth rate of the biomass.The net rate of biomass production consists of biomass production and decay.Increased decay rate constant of anodophilic bacteria means more bacteria which is responsible of producing electrons is lost.Thus, a lower power density could be produced.Half-rate constant of anodophilic bacteria ( K s,a ) represents the con- centration of the substrate at which the growth rate equals half of the maximum growth rate in Monod kinetics.Therefore, K s,a is linked with the growth rate anodophilic bacte- ria.Increasing K s,a decreases the growth rate of anodophilic bacteria, hence, decreased the power produced.Therefore, it is expected that both K s,a and k d,a have a decreasing rela- tionship with power density production which is represented using Monte Carlo simulation in the next section.Mediator half-rate constant ( k m ) is the Monod constant when the oxi- dized mediator fraction per anodophilic bacteria is a limiting factor for the growth rate of anodophilic bacteria in addition to substrate concentration, and it is expected to have very similar manner as K s,a , which is also represented by Monte Carlo scatterplots.Number of electrons transferred per mole of mediator (m) appeared to be one of the most impacting factors on the power density.There are many types of artificial mediators which differ in the number of electrons they can carry.Neutral Red (NR) is an artificial mediator that carries two electrons which is shown in Eq. 28, while potassium ferricyanide (Fe(CN 6 ) 3− ) carries one electron, the redox reaction is represented in Eq. (29).Therefore, other mediators can be artificially produced to enhance the power productivity of MFC by carrying more electrons.It is worth mentioning that more than one type of mediator can be used in a MFC [34], if considering them as a group, then the number of electron transferred per one mole of the group would be increased, and this enhances the power production.While curve steepness factor ( K r ), has almost very neg- ligible impact on the power produced from MFC.This is also shown in scatterplot of power density versus the curve steepness factor in Monte Carlo simulation.Finally, the maximum specific growth rate of anodophilics microorganism ( max,a ), this has empirical values that depend on the bacterial species used as well as the environmental conditions.Also, max,a is a result of the true yield of cell production and the maximum anodophilic reaction rate (maximum specific rate of substrate utilization).Therefore, to increase max,a , the yield can be increased by increasing the mass of the cell produced per mass of substrate consumed.Any small increase in the max,a value will increase the value of power density as it is considered a sensitive factor for power density.When the value of max,a increases, the maximum growth rate of anodophilic bacteria increases; then, more bacteria are available to produce electrons if enough substrate is available.Therefore, an increasing relationship of max,a with power density in Monte Carlo simulation is expected.In fact, both K s,a and max,a are Monod kinetics parameters which are estimated using experimental methods.Therefore, a large uncertainty in expected to be associated with their values which are originated from the reaction conditions, number of samples and the measurement errors.For total-order indices, several techniques for sampling are considered to tackle the most accurate results compared to Sobol's method which is shown in Fig. 6.Sampling techniques which are used in this section are Sobol, double loop, Fourier Amplitude Sensitivity Testing (FAST) and Morris.Sobol's sampling technique depends on the direct integral calculations.Which is then developed further to reveal Saltelli's approach by new formula which improved the computational accuracy and called Saltelli-Sobol [20], while the double loop approach is a simple approach and considered as an alternative for the direct techniques after the improvements in the algorithms suggested by Plischke.FAST is a variance-based GSA method which is based on the conditional variances.FAST is considered as one of the most efficient techniques to calculate sensitivities via Monte Carlo integration.However, FAST is limited to the calculations of total effects only [24].Finally, the Morris method can be used to calculate the first-order or total-order indices; however, it is an extension for LSA in the whole parameter range.Therefore, it does not rank the importance of impact of parameters on the model output in GSA and considered as a screening method only [16].Figure 6 shows the different methods in calculating the total indices.It can be noticed that almost all the three methods of double loop, Sobol and FAST give similar results only Morris method underestimates the total indices.This gives an indication that Fig. 6 Total effects index using different data sampling methods the results of GSA based on the Sobol method on this model are correct as the other methods producing similar results.

Monte Carlo simulation
In this paper, Monte Carlo experiments are also performed on MFC model.Monte Carlo is a method of sampling based on the repeated random values of the factors in the study from the distribution that they follow; so that the statistical behaviour of the output function is obtained.In addition, the factors assumed to be independent from one another so that the samples (input values) are extracted from the distribution they follow.In this study, normal distribution is assumed for all of the parameters under the study [28].If Y is assumed the output function (power density), then the input matrix for Y is as follows: Then, calculating Y for each row using the model, then the output matrix is: With this Y matrix and the model inputs, a scatterplot can be produced for each of the model parameters and the output by projecting the N values of the selected output Y vs the N values for each of the input parameters.For MFC model, each of the parameter produce a scatterplot with the model output of power density.These scatterplots are crucial to investigate the behaviour of parameters on the output [16].The scatterplots of the MFC power density versus the three factors of decay rate constant of anodophilic bacteria ( k d,a ), half-rate constant of anodophilic bacteria ( K s,a ) and mediator half-rate constant ( k m ) are represented in Fig. 7.
Figure 7 show K s,a , k d,a and k m relationship with power density using Monte Carlo representation.All the three factors show a discernible linear decreasing relationship with MFC power density.Increasing the decay rate constant of anodophilic microorganisms means increase in the loss of the active biomass, therefore, decrease in the number of biomass which can produce electrons.While the half-rate constant of anodophilics (K s,a ) represents the half saturation concentration in Monod equation which describes the bacterial growth linked to the substrate consumption.
K s,a is the concentration of substrate at which the specific growth rate of the microorganism ( ) becomes half of the maximum specific growth rate ( max ).In addition, k m represents the half-rate constant for the Monod equation representing the oxidized mediators impacts being a limiting factor for the growth rate of anodophilic bacteria as stated before.Therefore, its relationship with the power density is expected to be similar to K s,a .While for the curve steepness factor (K r ), the power density is found to be insensitive to this factor as shown in Fig. 8.It can be noticed that the behaviour of this output function is almost constant at different values of the input and the scatterplot become as a horizontal line towards all the values of this input.K r is a constant which determines the curve steepness factor in the equation of resistance which is identified experimentally in [7] by using the voltage measurements in MFC-1 operation.GSA revealed that power density is insensitive to this factor even though this factor is determined by the voltage values.The change in power density function due to the change in K r is negligible as shown in Fig. 8 which is only 0.04 mW/L for the whole range of the factor.
For maximum specific growth rate ( max,a ), a discernible strong non-linear increasing pattern is shown for the power density.As stated previously, max,a is an important Monod parameter for bacterial growth.As max,a increased, the growth rate of anodophilic bacteria which produce power increase through Monod equation which is shown below.
The relationship between a , max,a and the substrate (S) is non-linear increasing, while increasing max,a increases K s,a which represents the amount of the sub- strate concentration required to reach the half of the maximum growth rate.This impacts the growth rate by decreasing it.Increasing the maximum growth rate and substrate in the system increases the power produced as the nutrient is available for anodophilics to oxidize and produce electrons.The relationship of the power density and max,a is shown in Fig. 9.
For the last parameter which is the number of electrons to be transferred by the mediator, the power density has shown a discernible non-linear increasing relationship in the scatterplot that is shown below in Fig. 10.
The power density can be increased if the number of the electrons transferred per a mediator increase.This mainly depends on the type of the mediator used and partially the substrate.Artificial mediators vary in the number of electrons that they can carry.Therefore, increased the number of electrons transferred per mole of mediator enhances the amount of produced power.Figure 10 shows a discernible non-linear increasing relationship of power density and the number of electrons transferred per mole of mediator.

Conclusions
In conclusion, local and global sensitivity analyses are powerful tools in revealing the impact of different parameters on the output functions.Variance-based global sensitivity analysis using Sobol's method was applied on multi-population, single-chamber MFC in a continuous flow mode at steady state.The conducted study revealed the extent and pattern of influence of some selected parameters on the produced power density.GSA results showed that the maximum specific growth rate and number of electrons transferred per mole of mediator were the most impacting factors on power density and both have a discernible non-linear increasing relationship with power density using Monte Carlo simulation.0.736 was the total indices value for max,a which was greater by 0.072 from the first-order indices for the same parameter.While this parameter was greater than the Number of electrons transferred per mole of mediator by 0.112 in total indices and by 0.12 in first-order indices which had values of 0.624 and 0.544 for total-order indices and first-order indices respectively.Decay rate constant, half-rate constant of anodophilics and mediator half-rate constant had almost the same impact on the output function with values of total-order indices ranging (0.19-0.7), and all showed a decreasing linear relationship with power density upon using Monte Carlo simulation.Finally, curve steepness factor had almost no impact on the results of power density having a value of 1.76 ×10 −3 in total-order indices which is near to zero.Monte Carlo simulation results agreed with the GSA in terms of revealing which factor can affect the power density of MFC using scatterplots.It can be concluded that global sensitivity analysis is very powerful tool and can be successfully applied for parameter analysis and optimization.This study can be utilized further

Fig. 3
Fig. 3 Representation of the sensitivity analysis method

Fig. 7
Fig. 7 Scatterplots of MFC power density vs a half-rate constant of anodophilics (Ks,a).b Decay rate constant of anodophilics (kd,a).c Mediator half-rate constant (km)

Table 2
Selected parameters and their corresponding ranges for sensitivity analysis study

Table 3
First-and total-order indices