Dynamic Modelling of Tetrazolium Based Microbial Toxicity Assay – A Parametric Proxy of Traditional Dose-Response Relationship

19 Microbial toxicity of test substances in tetrazolium assay is often quantified while referring to 20 IC 50 values. However, the implication of such estimates is very limited and can differ across 21 studies depending on prevailing test conditions. In this work, a factorial design-based end-point 22 microbial assay was adopted, which suggests significant interaction ( P = 0.041) between 23 inoculum and tetrazolium dose on formazan production. A dynamic model framework was 24 proposed to be incorporated in toxicity assay, that not only captures the nonlinear 25 interdependency between biomass, substrate, formazan content but also measure the toxicity 26 in terms of inhibition parameter. Microbial growth, glucose uptake, formazan production in 27 presence and absence of heavy metal toxicant (Cu 2+ ) in designed batch studies were utilized 28 for sequential estimation of model parameters and their bootstrap confidence intervals. A 29 logistic growth model (R 2 >0.96) with multiplicative inhibition terms fit the experimental data 30 reasonably well. Dynamic relative sensitivity analysis revealed that both microbial growth and 31 formazan production profiles were sensitive to toxicant inhibition parameter. The adoption of 32 a dynamic model framework as a stable index for the toxic potential of test substances can be 33 extended to design a versatile, robust in-vitro assay system. 34 35


Introduction
Evaluation of toxicity for new molecules, formulation, materials, toxicants or xenobiotics using appropriate organism is an integral part of risk characterization and assessment (Gabrielson et al. 2003;Wadhia et al. 2007).Biological assay systems with a microorganism, invertebrates, plants, or higher organisms have been used to measure the change in test organism or its response in presence of the toxicant in a quantitative manner (Iqbal 2016;Bartlett et al. 2017;Özgür et al. 2018;Luan et al. 2020).Toxicity assay using microbes as bioindicators remained simple, fast, versatile, reproducible, and can be extended to environmental, food, or biomedical applications (Polo et al. 2011;Hassan et al. 2016).Depending on the microorganism, environment and biological response, appropriate analytical techniques such as respirometric (Vasiliadou et al. 2018), fluorescent and luminescent bioassays (Roslev et al. 2015;Abbas et al. 2018), microbial fuel cell-based sensor (Dávila et al. 2011), methanogenic activity assay (Gonzalez-Estrella et al. 2013), microbial growth (Baek and An 2011) can be adopted.Microbial reduction of tetrazolium has been widely used to measure microbial metabolic activity in toxicity assay (Wang et al. 2010;Lupu and Popescu 2013;van Tonder et al. 2015).
The prototype 2,3,5-triphenyl-tetrazolium chloride (TTC) is one of widely studied redox indicator, which gets reduced into water-insoluble red formazan by microbial dehydrogenase in metabolically active cell (Junillon et al. 2014;Sabaeifard et al. 2014).The amount formazan extracted from the microbial cells grown in presence of toxicants would be lower than that of the control (Sydow et al. 2018).The inhibition of tetrazolium reduction with increasing dose of toxicant has also been used to assign a numerical indicator of the pollutant e.g.IC50 valuethe effective concentration causing 50% reduction in measured response (Ma et al. 2018).
Despite being rapid and simple, selection of time for final reading in endpoint assay require optimization of the whole experiment (Stefanowicz-Hajduk and Ochocka 2020).In some cases, IC50 for a 24-hr incubation period can be different from 48-h incubation.A toxicant cannot be tagged with a unique IC50 value as the outcome of tetrazolium assay depends on a multitude of interacting variables, such as the initial concentration of microorganism, substrate, TTC, and the incubation time (Ghaly and Mahmoud 2007;Sylvester 2011).The formazan yield has been known to increase with TTC concentration but beyond a critical point the yield declines (Ghaly and Mahmoud 2007).Assumption of a proportional relation between the amount of formazan deposit to the number or activity of living cells has already been questioned (Etxeberria et al. 2011).Intracellular accumulation of formazan can compromise the integrity of the plasma membrane (Bernas and Dobrucki 2004;Lü et al. 2012) and cell viability.
Accumulated formazan is a result of non-linear dynamics of substrate consumption, microbial growth, and toxicity of formazan as well as an external toxicant, integrated over the incubation period.Constant rate of formazan production can only be presumed when kinetics is known in a priori (Januszek et al. 2015;Małachowska-Jutsz and Matyja 2019).Substrate inhibition models based on modified Michaelis-Menten kinetics or empirical variants have also been adopted to explain TTC reduction data (Ghaly and Mahmoud 2007).The inhibition constant or parameters estimated for toxicants while modelling of temporal change in biomass, substrate, formazan concentration can be an invaluable static, toxicity indicator, resilience to variation incurred in endpoint toxicity assay.
In this study, a comprehensive dynamic model was developed to simulate the temporal change in biomass, substrate, formazan concentration during a toxicity assay.Model parameters were sequentially estimated using independent experimental data sets, designed to explore incremental subspace of complex model hierarchy.Appropriate kinetics (and parameters) for growth and substrate consumption were selected from a set of representative models that corroborate experimental data.The model fitness was ascertained using cross-validation, and the confidence intervals (CIs) of all estimated parameters, their sensitivities were also evaluated.

Media preparation, culture growth and toxicity assay conditions
Overnight grown fresh Escherichia coli (ATCC 25922) in sterile nutrient broth medium centrifuged at 3500 × g for 15 min, re-dispersed in distilled water, was used as the test microorganism in all toxicity studies.Autoclaved BH media (pH 7) supplemented with filtersterilized glucose stock (10 g l -1 ) was used for further studies.TTC and CuSO4, if required in experimental sets, were added from the aqueous stock solution (5 g l -1 , 50 mg l -1 , respectively).

Influence of initial biomass, glucose, TTC on tetrazolium reduction
Influence of biomass concentration (X1), glucose (X2), and TTC dosage (X3) on tetrazolium reduction was evaluated by adopting a 2 3 full factorial experimental design as in Table 1.
Different volume freshly grown (and resuspended) E. Coli culture, glucose and TTC solution was transferred into a set of 15 ml-sterile centrifuge tubes.The tubes were wrapped with aluminium foil, incubated in a shaker incubator for 2 h under dark condition.Following centrifugation at 3500 g for 10 min, the cell palette was processed for formazan extraction and analysis.Analysis of variance was carried out to determine the individual effects as well as the interaction effects of the factors on the responses.Data processing and plotting were carried out using the software Minitab 18.

Experimental design for the discretization of dynamic model
The purpose of the design is to explore temporal change in microbial growth, glucose consumption, and TTC production (in presence or absence of toxic chemical), which can be used to model the dynamic behaviour of the process.The experiment was independently conducted in sets of 100-ml Erlenmeyer flasks with 50 ml of BH media (Table 2).Glucose was added at final concertation of 1 g l -1 , while 2 ml TTC (5 g l -1 stock) and 250 µl CuSO4 (50 mg l -1 ) were added if included in the design.All flasks were incubated at 37ºC in a shaker incubator for 16 h.On every 2 h interval, 2 ml sample was collected, optical density was measured at 600 nm, and kept in the refrigerator for further analysis.After completion of the experiment, samples were centrifuged at 3500 g for 15 min, the supernatant was collected for residual glucose analysis and cell palette was processed for formazan extraction.

Analytical Methods
The optical densities of the samples were measured at 600 nm (path length 1 cm) using a UVvisible spectrophotometer.The measurements were converted to dry-weight based biomass concentration (in g l -1 ) using a standard prepared in the identical range of concentrations.
calibration Residual glucose concentration in the culture supernatant was measured using dinitro-salicylic acid (DNSA) method (Kona et al. 2001) using a standard curve prepared with known glucose concentration (0.2-1 g l -1 ).Concentrated HCl (20 µl) was added to cell palette to release intracellular formazan, which was dissolved in 2 ml methanol (Lopez Del Egido et al. 2017).The sample extract was centrifuged at 3500 g for 10 min and absorbance of supernatants were recorded 480 nm against reagent blank.

Formulation of model structure
Microbial growth on glucose is assumed to be substrate limited process, and a flexible and generalized rate expression is preferred.High formazan concentration (F) and/or presence of toxicant (M) is assumed to inhibit microbial growth, where the increasing concentration of formazan crystal inside a cell can reduce viability.That phenomenon can be modelled as product inhibition perspective (Lin et al. 2008) as shown in generalized specific growth rate expression as in Eq. 1.
Where µg is a mathematical expression for glucose-limited microbial growth (explained later).
The second and third terms on the right-hand side represent growth inhibition due to formazan, and toxicant (e.g., heavy metal), respectively.The / FX r is the relative concentration of formazan (F) to biomass (X) and f K is the inhibition constant.The inhibition effect would be negligible at an early phase of incubation when biomass accumulated formazan concentration is very low (i.e., rF/X ≈ 0).However, at a very high level of intracellular formazan (i.e., rF/X >> Kf), the growth rate becomes inversely proportional to formazan content.The Km is inhibition constant for the toxicanta dose that would cause a 50% reduction in microbial growth rate than in the toxicant free environment.The spiked toxicant in the present study is heavy-metal, that is not degraded metabolized, and hence its concentration (M) is assumed to remain constant, timeinvariant.
The dynamic change in biomass, glucose, and formazan concentration is given by Eq. 2-4.
Glucose consumption and formazan production are (from tetrazolium reduction) modelled as growth associated or non-associated process as in Eq. 3 and 4, respectively.The term YX/G on the right-hand side in Eq. 3 represent true biomass yield coefficient, while the second term (β) represents maintenance requirement.The rate of formazan production ceases once formazan concentration attains its theoretical asymptotic maximum (Fmax) assuming 1:1 molar reduction of added tetrazolium in the culture broth.

Selection of Growth model
In absence of tetrazolium/formazan or heavy metals, the glucose-limited growth of microorganism G  was modelled using Logistic, Monod, Contois, and Tessier kinetics as Where Xmax is the asymptotic biomass concentration in logistic growth.The Kd, Kc, Kt are kinetic constant relevant to Monod, Contois and Tessier models.

Formulation of the dynamic model and parameter estimation strategy
A complete set of coupled ordinary differential equations (ODEs) reflecting the rate of change for all state variables is shown in Eqs. ( 1)-( 5).Simultaneous estimation of all eight parameters from a single experimental data set using biomass, glucose and formazan concentration profiles could be inherently inaccurate.A discretized approach, where only a subset of parameters is required to be estimated from the hierarchical model subspace, utilizing the data from a designed experiment (Table 2), was adopted as outlined in following steps.
Step 1. Experimental set containing only biomass, glucose data were used to select the ideal growth model (i.e., Logistic, Monod, Contois, or Tessier) as well as their relevant parameters.
The inhibition terms for formazan and toxicant as in Eq. 1 can be dropped.
Step 2. Experimental set containing biomass, glucose, formazan data was used to estimate growth inhibition term (KF) in Eq. 1, and formazan production parameters (γ, δ).Growth related parameters were adopted from an earlier step, the biomass yield coefficient or maintenance requirement of glucose ( / XG Y , β) was re-calibrated.
Step 3. Experimental set containing biomass, glucose, formazan, and toxicant was used to estimate growth inhibition term (Km) in Eq. 1 while adopting all other parameters from step 2.
The numerical solution of the dynamic model was approached from an initial value state variable and reasonable guess of parameters (within suitable space boundaries) in Matlab® ode45 solver.Subsequently, a normalized residual sum of squares (NSSE) between experimental and modelled profiles was minimized (Eq.6) while changing relevant parameter(s) in a constrained minimization algorithm using fmincon function in Matlab®.
( ) where u is the number of measured state variables (i.e., two or three in this study), n is the number of data points in each measured ( , ij x ) or modelled ( , ˆij x ) profiles.The x is the average measured value of the ith profile.

Cross-validation, confidence interval estimation and sensitivity
As a check against overfitting, and to assess the model's ability to predict data not used in model fitting (i.e., its capacity for out-of-sample prediction), a five-fold cross-validation analysis was carried out and coefficient of determination (R 2 cv) was estimated.The confidence interval (CI) for parameter estimates in the dynamic model data was processed using a residualbased, non-parametric bootstrap resampling approach (Dogan 2007).The dynamic relative sensitivity, i.e. relative sensitivity at different time intervals, for the model parameters was calculated using the finite difference approximation (Yue et al. 2006) with a repeated solution (twice) of the model for each parameter.The sensitivity analysis represented here is essentially a univariate, where a change in modelled responses (e.g., biomass, substrate, product) is mapped from the perturbation of single parameter value while keeping other parameters on hold at their original estimate.Local relative sensitivity as a function of time allows a direct comparison of responses at different states or across different parameters (Eq.7).
,, ,, where yi,t is the i th response at time t and θj is the j th parameter of which sensitivity is being measured.

Impact of initial biomass, substrate, and tetrazolium on formazan concentration
A factorial analysis was performed on the data to study the influence of inoculum size, glucose concentration, TTC dose and their interactions on formazan production and the results are presented in Table 3. Analysis of variance performed on the data set (Table 4) revealed that inoculum size and TTC dose and their 2-way interactions did significantly affect formazan production (P<0.05).This is also supported by the Pareto chart (Fig. 1) showing order of influencing parameters.The plots for individual factor effects on formazan productivity showed that increasing amounts of inoculum was favourable while increasing amount of glucose or TTC can inhibit actively growing microbial culture, resulting in lower TTC reduction into formazan.Similar toxic effect of tetrazolium salts has been made in previous studies (Junillon and Flandrois 2014).A reversible increase in the duration of lag phase has also been observed for gram-positive Listeria monocytogenes culture growing in presence of different tetrazolium salts (Junillon and Flandrois 2014).

Dynamic model parameters in toxicity assay
The microbial biomass and glucose concentration profile from experimental sets without TTC or toxicant (i.e., cupper salt) were modelled using Logistic growth, Monod, Contois, and Tessier models.The fitness of the chosen models against experimental data has been shown in The experimental and modelled profiles (for biomass, substrates and formazan) are also shown in Fig. 2, depicting the suitability of the selected kinetics for modelling observations in toxicity assay.The modelled profiles are in good agreement with experimental data, except for substrate concentration which is somehow over-predicted at the late phase in a batch experiment involving both TTC and toxicants.
The least-square estimates of substrate consumption, formazan production and microbial growth inhibition (due to formazan and toxicant) parameters along with their boot-strap CIs are presented in Table 6.The high cross-validation R 2 for both biomass and substrate concentration profiles suggests no over-prediction in the proposed logistic growth model.
Certainly, the estimates for µmax (0.3674 h -1 ) and Xmax (1.2658 g l -1 ) can be extended to other experimental sets containing TTC and/or cupper, as long as the growth-limiting substrate remains identical.However, estimates biomass yield and maintenance coefficient can change in presence of formazan and/or toxicant (Şengör et al. 2009).Those parameters were recalibrated from additional experimental sets containing inhibitor /toxicants.It should be noted that all parameters were not estimated from a single batch experiment, rather sequentially estimated that improves the accuracy (Wang et al. 2018).
The biomass growth yield coefficient (0.747 g g -1 ) and formazan inhibition constant KF, (30.963 mg g -1 ) were estimated from experimental sets containing glucose, and TTC.The yield coefficient estimates seem to be in the same order of growth-associated stoichiometric true yield coefficient (0.57 g g -1 ) as reported for E. coli in glucose-limited continuous culture (Kayser et al. 2005).It should be noted that the biomass yield coefficient reported in the present study does not consider the corrections for the maintenance requirement.Though literature advocates bacterial growth inhibition in presence of tetrazolium salts (Villegas-Mendoza et al. 2015), no attempts (either mechanistic or empirical) has been made to model the kinetics of formazan toxicity.No comparison or discussions on the estimate of KF (30.963 mg g -1 ) can be made in the light of existing literature.
The estimates formazan production parameters (γ and δ respectively) suggests the process be primarily growth associated one.The non-growth associated formazan production rate (1.257e - 05 mg g -1 h -1 ) is very trivial and accumulated intracellular formazan within the incubation period for in-vitro toxicity assay is likely to be far below than the KF.In other words, non-growth associated formazan production alone cannot exert any toxicity.A significant positive relationship between specific growth rate and formazan production rate has also been reported in other studies (Villegas-Mendoza et al. 2019).Though formazan production in batch cultures increases with time, the formazan production rate decreases more rapidly.
The inhibition parameter pertinent to metal toxicity has been estimated to be 152.39mg l -1 from the third batch study containing TTC as well as Cu +2 .The regression coefficients for biomass, substrate and formazan concentration profiles (0.856, 0.782, 0.964) in this least square estimation exercise were inferior in comparison to the other two batches.This could be attributed to constrained parameter subspace for Km, as all other parameters were already prefixed, leaving a lower degree of freedom.Utgikar et al modelled toxicity of copper on sulfatereducing mixed microbial cultures at two perspectives, i.e., a toxic effect causing the death of cells or inhibitory effect using an inverse exponential function (Utgikar et al. 2003).The author reported inhibition constant 1140 mg l −1 for copper, which is about 10-fold higher than the current estimate.The reported half inhibitory concentrations (IC50) for Cu(II) in enzymatic or respirometric assays vary widely (<0.2 to 200 mg l -1 ) in the literature (Esquivel-Rios et al.

2014
).This difference or discrepancy is very logical considering the variance in microbial culture, test conditions and the mathematical model adopted.

Bootstrap parameter distribution and sensitivity
The distribution of bootstrap parameters in the proposed toxicity model is shown in Fig. 3, where a near symmetric normal/gaussian distribution for μmax, Xmax (from glucose consumption and biomass growth data) assures better credibility.The distribution of parameters such as KF and YX/G are somewhat right-skewed, which can be related to the inaccuracy of the estimate, or poor parameter sensitivity.The distribution of γ and Km exhibits a marginal left-skewness, implying the presence of outlier in experimental data (Joshi et al. 2006).Such deviation can be linked to the hierarchical parameter estimation approach and nonlinearity of the dynamic model.
Sensitivity analysis of the toxicity model parameters on biomass, substrate, and formazan content is shown in Fig. 4. The relative sensitivities of the different parameters not only vary between profiles but also with time for a given response.The growth profile appears to be most sensitive toward µmax, Kf, Km.This is very logical as the specific growth rate is directly dependent on those parameters.Moreover, the relative sensitivities have been shown to increase with time.This can be attributed to a higher amount of accumulated biomass and formazan at later phase of microbial growth.As the growth model is not directly coupled with substrate consumption, growth data remains insensitive to β, δ, or YX/G.The negative relative sensitivity of substrate profiles with Xmax reflects faster depletion of the substrate.The substrate consumption profile remains insensitive to β, δ, which are non-growth associated parameters for glucose consumption or formazan production, respectively.Entire sets of parameters that had a significant influence on growth (i.e., µmax, Kf, Km, γ) appears to impact formazan production as well.However, the sensitivity of γ on the formazan profile is positive which is very reasonable from their kinetic relationship.As formazan content increases, the absolute sensitivity of γ declines over time.Conceivably, parameters in the dynamic model are temporarily separated where any parameter(s) relevant to a later stage does not influence the preceding process responses.

External Validation of the Proposed Model
A dynamic model for toxicity assay was validated using an independent TTC free experimental set spiked with 50 mg l -1 Cu 2+ while keeping the other process conditions unchanged.The experimental and modelled data for biomass and residual glucose concentrations are presented in Fig. 5.The experimental and modelled data were in close agreement, particularly at the initial growth phases.However, at the later stage, experimental biomass and glucose profiles drifted toward the lower and higher side respectively as compared to the model.Although an explanation for such discrepancies is not available, such deviations were within reasonable limits of uncertainty due to measurement error.

Conclusions
Mathematical formalization of tetrazolium based microbial toxicity assay and projection of estimated kinetic model parameter as an invariable and reproducible alternative to traditional toxicity assay reporting has been proposed.Growth inhibition was modelled as a function of intracellular formazan content and toxicant concentration.Parameters, sequentially estimated from batch experimental data, were credible with narrow confidence intervals.Though a limited combination of growth, substrate consumption and inhibition models were only adopted in this study, the approach can easily be extrapolated to much larger and multivariate data set incorporating mechanistic insight of the toxicity assay.Parametric representation of result can offer a unified format of reconciliation, reporting and validation for in-vitro toxicity assay.
Table 1.Factorial design to examine the dependence of TTC reduction on biomass, substrate and TTC concentration.

Figure 1 .
Figure 1.Influence of inoculum, glucose and TTC dose on formazan concentration as depicted in (a) main effect plot as a mean response, (b) interaction plot for each combination of two factors, and (c) Pareto chart where bars represent the absolute value of the effect of each factor and the dotted line indicates the limit of significance (α = 0.05).

Figure 2 .
Figure 2. Experimental (symbols) and model-predicted (solid line) profiles for biomass, substrates and formazan concentration (if applicable) in batch studies involving sets with (a) no TTC or toxicant, (b) TTC but no toxicant, and (c) with TTC and toxicant.

Figure 3 .
Figure 3. Distribution of bootstrap parameter estimates (N=5000) relevant to microbial growth, substrate consumption, formazan production in the proposed toxicity model.

Table 5 .
Lindstrom et al. 1998ens et al. 2011)odel to biomass or substrate profiles (respective R 2 of 0.9663, 0.9633) is better than the other three models.It should be noted that the biomass growth is not coupled with substrate concertation in the Logistic model structure, but the substrate depletion is proportion to biomass growth.The logistic model has successfully been used to explain growth and substrate consumption by E. coli in batch studies under substrate limited conditions(Li et al. 2010;Wittgens et al. 2011)or population growth model in microbial culture in toxicity assay (E.Lindstrom et al. 1998).

Table 2 .
Hierarchical experimental design to decipher tetrazolium-based toxicity model

Table 3 .
Factorial design matrix and formazan concentration

Table 4 .
Analysis of variance performed on the formazan production data.514

Table 5 .
Performance of different kinetic models fitting microbial growth and glucose utilization data in absence of tetrazolium or toxicant.

Table 6 .
The bootstrap confidence interval of parameter estimates for the dynamic model structure at different modelling stages and the regression coefficients $ $ four-fold cross-validation R 2 (R 2 CV) with 10 Monte-Carlo repetitions for validation.