Application of model fitting technique to enhance bacterial regrowth potential (BRP) measurement for drinking water supply monitoring

The bacterial regrowth potential (BRP) method was utilised to indirectly measure the assimilable organic carbon (AOC) as an indicator for the assessment of the microbial regrowth potential in drinking water distribution systems. A model using various microbial growth parameters was developed in order to standardise the experimental interpretation for BRP measurement. This study used 82 experimental BRP data sets of water samples collected from the water treatment plant to locations (customer taps) in the distribution system. The data were used to model the BRP process (growth curve) by a data fitting procedure and to obtain a best-fitted equation. Statistical assessment and validation of the model obtained equation by fitting these 82 sets of data were conducted, and the results show average R values were 0.987 for treated water samples (collected at the plant prior to chlorination) and 0.983 for tap water (collected at the customer taps). The F values obtained from the F-test are all exceeded their corresponding F critical values, and the results from the t-test also showed a good outcome. These results indicate this model would be successfully applied in modelling BRP in drinking water supply systems.

. These parameters are able to give analytical information to describe bacterial growth characteristics and hence the BRP of water (Withers & Drikas 1998a).
Growth factor (f) ¼ Turbidity (final) Turbidity (initial) Growth rate (m) ¼ Turbidity (final) À Turbidity (initial) DTime Normally, a microbiological growth curve has a lag phase, maximum growth rate and asymptotic phases to describe its behaviour (Zwietering et al. 1990). Hambsch & Werner (1993) determined the growth rate in the BRP method using Equation (2). However, the exponential phase is estimated with uncertainty, due to the difficulty of not being able to consistently determine the initial and final turbidity values on the curve and it can lead to relatively inaccurate determination of the growth rate. To develop a more consistent approach and improve accuracy, both the lag phase and stationary phase are required to be considered as part of the maximum growth rate calculation. Therefore, using the statistical modelling approach to develop a mathematical curve to describe the behaviour of BRP growth curve should improve the consistency of measurements of the maximum growth rate, lag phase and AOC values and enable more reliable comparisons between water samples (van der Kooij et al. 2017).
The Weibull regression has been used to model indeterminate growth for plants, fish and molluscs. López et al. (2004) and Swintek et al. (2019) indicated the Weibull model achieved a better correlation coefficient for microbial growth experimental data than the commonly used Gomperts model/sigmoid curves (Tjørve & Tjørve 2017).
The aim of this study was to develop a mathematical model using BRP data to improve understanding of the behaviour of bacterial regrowth in drinking water distribution systems. The objectives of this study are (1) to use experimental data to determine best-fitted models using the curve-fitting software, (2) to understand bacterial growth curve by developing/extracting microbial growth parameters from the determined model and (3) to validate and evaluate the developed model and compare it with the Weibull model.

MATERIALS AND METHODS
A series of BRP experiments (82 sets) were conducted using water samples obtained from six different distribution systems to generate a database and use this to develop a mathematical model to further study the BRP technique.

Water samples
Six drinking water distribution systems, Anstey Hill, Barossa, Myponga, Happy Valley, Hope Valley and Little Para, in the South Australian metropolitan area were selected to conduct this study. The drinking water in these systems was supplied by the corresponding treatment plant using conventional treatment processes (coagulation/flocculation, sedimentation or floatation and then sand filtration). Water samples were collected at the water treatment plant from the treated water tank (after filtration and prior to chlorination) and from the end of distribution system (consumers' taps).

Instrument
A Monitek Turbidimeter (Monitec mAOC Analyser Model 251-4 mAOC, Germany) was used for the acquisition of the BRP data. There were four units of this instrument, each one had a four-channel turbidimeter (total of 16 channels) connected with an external computer and automatically registered readings and analysing growth curve through an accompanying computer program. The accompanying program exported the growth curve data (raw data) to a text format file, as time and turbidity (two columns of data).

Experimental procedure
2.2.2.1. Sample and inoculum preparations. Water samples were collected in glass bottles with a quenching agent (4% of sodium thiosulphate) added and transported back to the laboratory with cold pack. The samples were then filtered through pre-rinsed glass fibre membrane filters. The sterilised (autoclaved) inorganic nutrient salt solution prepared using analytical grade chemicals, ammonium chloride, calcium nitrate, calcium chloride, magnesium sulphate, ferrous sulphate, sodium silicate nanohydrate, aluminium sulphate and potassium dihydrogen orthophosphate, with Milli Q water was added to the water samples to ensure AOC was the limiting nutrient (Withers & Drikas 1998b). Subsequently, pre-rinsed polycarbonate membranes were used to filter the sample, and then these membranes were placed in isotonic saline solution to obtain natural consortia as inoculum (Withers & Drikas 1998b).

Data processing
The water industry used the carbon equivalent value as an indicator for water quality, because the carbon equivalent value is more practical in monitoring the water treatment processes (Withers et al. 1996). Hence, an investigation (Withers & Drikas 1998b) to establish a link between an acetate carbon equivalent (ACE) and growth factor had been carried out, and a calibration curve has been established by utilising an acetate standard spiked into different types of water that relates acetate to biomass (Withers & Drikas 1998b). The ACEs were calculated by turbidity change in a bacterial regrowth curve of the water sample, dividing average slope (1.5). It has been determined that an ACE of 80 μg/L is equal to a growth factor of 5, which are thresholds to determine if the water is biologically stable or not (Withers & Drikas 1998b). These researchers carried out further calibrated experiments and then determined that an average slope of 3 (Equation (3)) and an acetate equivalent of 40 μg/L are equal to a growth factor of 5. The raw data, two columns of data (time and turbidity), were exported by the BRP instrument's accompanying program to a text format file. The AOC value was calculated by using the formula function developed in-house in Microsoft Excel™.

Develop a mathematical model
The curve-fitting software TableCurve 2D version 5.01 (Systat Software, Inc., Richmond, CA, USA) was selected to analyse data and determine ideal empirical equations representing the bacterial growth curve for different types of water. TableCurve 2D was applied to find the best models by fitting the experimental data to its database of a wide range of equations in its equation pool. To simplify the process, the automated fitting feature was selected, and 3,665 built-in equations were used to provide the ranks of R 2 , degrees of freedom-adjusted R 2 , maximum error and F-statistic.
The best-fitted equation (Equation (4)) was selected after 10 sets of data (five of them were treated water and another five of them were tap water) which was fitted by utilising TableCurve 2D. Subsequently, this equation was used for fitting all data. Finally, this equation was determined as a relatively ideal model to represent the BRP of waters based on the good R 2 values. Some preliminary data were selected from database for using TableCurve 2D fitting data, and the resulting model was restricted to fit non-linear equations. The logistic dose-response equation (No. 8013) in TableCurve 2D's library was determined the best-fitted model for describing the BRP of waters. This equation is a four-parameter logistic non-linear model; it is a form of Hill equation (Gadagkar & Call 2015). The shape of the curve of this equation visually fitted the bacterial growth curve (Gadagkar & Call 2015). Hence, this equation was determined as the preliminary equation to fit the other sets of data to evaluate if it was the best equation. The equation as follows, was determined to be the best-fitted model to describe the BRP of waters in the drinking water distribution system, where y represents the turbidity (microbial level, measurement time interval is 30 min), and t represents the incubation time. Parameter a is the lower asymptotic value, and parameter b is the difference between the upper asymptotic value and the lower asymptotic value. The coefficient of parameter c is near to x-value of inflection point of the growth curve but is larger than this value. When parameter d is smaller than À1, a sigmoidal shape curve is being able to describe the microbial growth curve. When parameter d is greater than or equal to À1, it is a saturation shape curve that has no lag time. If bacterial population size increases exponentially after a very short time, the growth curve's shape is more likely to be a saturation shape visually, so the growth curve can be fitted by adjusting the parameter d (,À 1) to make the fitted curve be a sigmoidal curve. Due to the AOC value is determined by Equation (3), the model-determined AOC value is calculated using the following equation: where b (turbidity difference) from Equation (4).
Two-column data namely time and turbidity were processed using TableCurve 2D. Every set of experimental data was processed by TableCurve 2D and fitted to the selected Equation (4) (named as BRP-Hill model hereafter) to obtain the coefficient of the equation for the experimental data. For a few sets of data, a rational coefficient for parameter b from TableCurve 2D was not able to be determined, because the software output for the highest R 2 value leads to coefficients of parameter b being extremely high. Microsoft Excel Solver was applied to adjust coefficients for these sets of data. Although the R 2 values after adjustment were slightly lower than the TableCurve 2D fitted, the coefficients of parameters provided for accurate fitting. (4)) The goodness of fit of the BRP-Hill model (Equation (4)) was statistically assessed by the coefficient of determination (R 2 ) values.

Alternative algorithm to cover the issue caused by insufficient experimental running time
When modelling some BRP data sets by the BRP-Hill model (Equation (4)), an issue was noted where differences between the model fitted AOC value and the instrument built-in determined value were relatively high. This is because the lag phases of some BRP experiments were relatively long and they did not reach the stationary phase during test. In other words, the experimental operator stopped the experiments too soon at the predetermined time, which followed the routine procedure. Hence, a mathematical approach (second derivative) was applied.

Model validation
To validate this model, this study compared the instrument built-in method for AOC values with model-determined AOC values using regression by plotting instrument built-in vs. model-determined and both the slope and R 2 were used, and errors were enumerated for two types of water, treated water and tap water.

Microbial growth parameters development
To describe a bacterial growth curve, the characteristics of the curve need to be studied. Four phases of the growth curve need to be evaluated: (1) the lag time (T Lag ) (the intercept of lower asymptote line and tangent line that passes through the inflection point), (2) the maximum growth rate (μm) (the slope of the tangent line through the inflection point), (3) lower asymptote (start value of the curve) and (4) upper asymptote (maximum values that the curve is able to reach) (Zwietering et al. 1990). In this model, those four phases are able to be calculated (as shown in Figure 1(c)), and this algorithm is according to Zwietering et al.'s (1990) method.
The second derivative of the BRP-Hill model (Equation (4)) with respect to t is obtained as follows: At the inflection point, where t ¼ t i , the second derivative is equal to 0 as follows: The maximum growth rate is the slope of the tangent line through the inflection point, which can be gained by calculating the first derivative at the inflection point: The tangent line through the inflection point can be described as follows: Therefore, the lag time can be calculated: When t ¼ 0, the lower asymptotic value can be obtained from the following equation: When t approaches infinity, the upper asymptotic value can be calculated from the BRP-Hill model (Equation (4)): Therefore, the four phases of the growth curve can be expressed by the parameters of BRP-Hill model (Equation (4)): Maximum growth rate: Lower asymptotic value (initial turbidity): a Upper asymptotic value (maximum turbidity): a þ b (17)

Models comparison
The known growth curve (Weibull model) was selected to compare with the BRP-Hill model (Equation (4)), and Excel Solver Add-in was applied for fitting the Weibull model. Table 1 shows equations of those two models and the four phases of growth curve expressed by parameters of those two equations. Six sets of data are treated water selected from drinking water distribution system. The R 2 values, microbial growth parameters and model-determined AOC values that were obtained from two models for these six sets of data were compared, and the Akaike information criterion (AIC) (López et al. 2004) method has been used to compare the Weibull model and the BRP-Hill model (Equation (4)). The AIC method can be used for a comparison of either nested models or non-nested models (Motulsky & Christopoulos 2004). The AIC is computed by the where SS is the sum of square, n is the number of data points, p is the number of parameters of the fitting model. This equation is the corrected equation of original one (López et al. 2004;Motulsky & Christopoulos 2004). Individual AICc values were obtained for comparison, and the model that has a lower value is more likely to be better fitted (López et al. 2004). The evidence ratio is calculated by comparing two models' AICc values (Equation (19)), for example, when comparing two models, and the evidence ratio is 2, it means the model with the lower AICc value is more than two times likely to be correct than the model with the higher AICc value.

Mathematical modelling of bacterial growth in waters
Examples of growth curves of two water samples showing details of the experimental data and the 95% confidence interval of the fitted growth curves are shown in Figure 2. The model-determined AOC values, maximum growth rates and lag times are shown in Figure 2. These water samples were collected from the Myponga drinking water distribution system (treated water from the plant and tap water from customer taps, respectively) in winter. This model provided good fitting to the data of these waters, and the R 2 values are all greater than 0.99.

Coefficient of determination (R 2 ) for the BRP-Hill model (Equation (4))
The BRP-Hill model (Equation (4)) provided overall good fitting to the data for plant treated water (after filtration and prior to chlorination) and tap water samples, shown in Table 2. Although tap water's minimum R 2 value was smaller than 0.900, only one set of tap water's R 2 values are smaller than 0.900. These data indicate the BRP-Hill model (Equation (4)) fits suitably the experimental data for the BRP of the waters.

)) and Weibull model's equations and tree phases of growth curve and model-determined AOC value that expressed by their parameters
Model BRP-Hill model (Equation (4) AQUA -Water Infrastructure, Ecosystems and Society Vol 00 No 0, 8 Corrected Proof

Alternative algorithm to calculate model-determined AOC
With pre-set analysis times, some cases of BRP analyses with relatively long lag phases were stopped before the plateau phase was reached. Hence, the original method for AOC value determination and model-determined values have unacceptable differences. For instance, Figure 3(a) shows an example of a large difference between the method built-in AOC value and the model-determined AOC value of 233.8 and 293.1 μg/L ACE for a data set acquired using Equations (5) and (3), respectively. These two AOC values were calculated from the model-determined turbidity difference (b) and experimental turbidity difference. It can be seen that the end point of predicted curve (Pred-Tursolid curve)'s vertical axis value is 0.76, and the upper asymptotic value is 0.93 (coefficient of parameter a plus coefficient of parameter b). Hence, the model predicted curves have not reached the plateau phase, and that was the reason for the relatively large difference between the AOC values. It can be inferred that for such experiments the full monitoring period to plateau stage was not completed, and that made experimental AOC values lower than model-determined values, which was based on a plateau being attained.

Corrected Proof
When parameter d is smaller than À1, the BRP-Hill model (Equation (4)) is a sigmoidal curve with an upper asymptote (Panik 2014), and it has a bell-shaped first derivative (Basterretxea et al. 2004). As seen in Figure 4, for Example 2, that when the BRP experiment is run over a longer time to reach the stationary phase, the circular dot on x-axis (the cross point of tangent line 2 and x-axis) is a distinct point where the curve did not reach the plateau (when growth curve (top solid line)'s x-value was smaller than circular dot's x-value. Hence, if the experiment stops before this point (cross point of second tangent line and x-axis), the full monitoring period is not completed.
For Example 1, the triangular dot on the curve shown in Figure 3(b) is the second inflection point of first derivative curve and from the BRP-Hill model (Equation (4)), x″ ¼ 87.4. Therefore, the x-value of the second inflection point of curve of first derivative equalled 87.4.
where x″ is the x-value of second inflection point of first derivative curve.  (4)) for Example 1, where the dotted curve is the first derivative curve and the dotted curve is the second derivative curve, circular dot is the second inflection point of second derivative curve and the triangular dot is the corresponding point of first derivative curve for circular dot (x 0 ¼ x 00 ¼ 87.4).
AQUA -Water Infrastructure, Ecosystems and Society Vol 00 No 0, 10

Corrected Proof
The second inflection point of first derivative curve was selected as the first point to calculate the average AOC value. Some sets of data were investigated, and the first derivative values of this point are not very small, which means there was modest gradient of this point on the model curve. Subsequently, calculating the average turbidity value from x″ to the end point, and obtain the difference between this average value and parameter a, then dividing by 3, gives a result that equals 198.0 μg/L ACE. The instrument built-in AOC value equals 212.0 μg/L ACE, the model-determined AOC value equals 234.7 μg/L ACE, and hence the result from this approach was the closest value to the instrument built-in AOC value.

Model validation (AOC values)
The BRP-Hill model (Equation (4)) was investigated to determine the accuracy of how it represented the experimental data. The 82 sets of data were analysed by comparing the instrument built-in values and model-determined values in terms of AOC. There were seven sets of data's model-determined AOC values being calculated by applying the alternative algorithm. The model's parameter values can be utilised for determining the AOC value on the condition where experiments were stopped after arriving at the plateau phase, and the alternative algorithm was used to determine AOC values for experiments that were stopped before reaching the plateau phase. The linear relationships between the AOC value obtained from instrument built-in and the AOC value determined from the model's parameter values for the two types of water are presented in Figure 5. The slope values of the two regression lines are close to 1, and the R 2 values of the plant treated water and tap water were 0.97 and 0.98, respectively. These values indicate that the results determined by the instrument built-in method and model-determined results are in good agreement.  In this study, the Weibull model frequently used in predictive microbiology (Cunha et al. 1998) was selected for comparison for the proposed BRP-Hill model (Equation (4)). In this study, Gomperts model and another frequently used logistic model (Fujikawa et al. 2004) have been trialled to fit the experimental data, but every R 2 value obtained for those experimental data was smaller than 0.9. Hence, those two models were not been considered further in this study. Table 3 provides details of six sets of data of plant treated waters and tap waters selected from six drinking water distribution systems (Anstey Hill, Barossa, Myponga, Happy Valley, Hope Valley and Little Para). It can be seen that both the BRP-Hill model (Equation (4)) and the Weibull model described well these growth curves. Half of the BRP-Hill model (Equation (4))'s determined maximum growth rates were closer to the observed growth rates, and half of the Weibull model's determined maximum growth rates were closer to observed growth rates (Table 3). For the lag time, all of the BRP-Hill model (Equation (4))'s determined lag times were closer to observed lag time than the Weibull model. There are four sets of data where the determined AOC values by the BRP-Hill model (Equation (4)) are closer to the observed AOC values than the Weibull model (Table 3).

Comparison of the BRP model with the Weibull model
There were 24 sets of data for comparing AICc values between the BRP-Hill model (Equation (4)) and the Weibull model selected randomly from all drinking water distribution system data for two types of water. Seventeen out of the 24 sets of data's BRP-Hill model (Equation (4)) had smaller AICc values, and the mean evidence ratio of BRP-Hill model and Weibull model for these 17 sets of data was 10.79 (data not shown). When some sets of data were fitted by Weibull model, their coefficient d of equation was smaller than 1. The curve is saturation shape when d 1, and the curve is sigmoid curve when d . 1 (Swintek et al. 2019). However, the microbiological growth curve generally is a sigmoidal curve (Fujikawa et al. 2004).
In summary, the BRP-Hill model (Equation (4)) was found to be better for describing bacterial regrowth curve in water than the Weibull model. Although both the Weibull model and the BRP-Hill model (Equation (4)) have good (high) R 2 values, the Weibull model cannot provide a good fit when experimental data plotted curve has a more saturation shape, and as well, the BRP-Hill model had better AICc scores.

CONCLUSIONS
A BRP-Hill model (Equation (4)) that provides microbial growth parameters to describe the bacterial growth curve has been developed with some alternative algorithms to adapt the large difference between the method built-in AOC value and the Table 3 | BRP-Hill model (Equation (4)), Weibull model-determined maximum growth rate, lag time, lower asymptotic value and upper asymptotic value, and observed maximum growth rate, lag time and AOC value for these six sets of data  (4)) fitted data sets (n ¼ 82) were close to 1. The average R 2 values were 0.987 for plant treated water and 0.983 for tap water (as shown in Table 2). Furthermore, the BRP method built-in AOC values and model-determined AOC values for the 82 sets of data were compared, and the result show this model fitted these data well. The BRP-Hill model (Equation (4)) parameters can help improving understanding of BRP and also of the AOC in waters. This suggests that the developed model can be a helpful tool to understand the behaviour of bacterial regrowth of drinking water and enable prediction where some data is not available. However, the BRP-Hill model (Equation (4)) was developed based on experiments carried out at a constant temperature (22°C); thus, all the BRP experiments were conducted at this temperature.
The BRP-Hill model (Equation (4)) may not be able to fit data from the BRP experiment conducted in a different environment such as at another constant temperature (30°C).