Use of Stochastic Models To Assess the Effect of Environmental Factors on Microbial Growth

We present a novel application of a stochastic ecological modelto the study and analysis of microbial growth dynamics as influencedby environmental conditions in an extensive experimental dataset. The model proved to be useful in bridging the gap betweentheoretical ideas in ecology and an applied problem in microbiology.The data consisted of recorded growth curves of Escherichiacoli grown in triplicate in a base medium with all 32 possiblecombinations of five supplements: glucose, NH₄Cl, HCl, EDTA,and NaCl. The potential complexity of 2⁵ experimental treatmentsand their effects was reduced to 2² as just the metal chelatorEDTA, the presumed osmotic pressure imposed by NaCl, and theinteraction between these two factors were enough to explainthe variability seen in the data. The statistical analysis showedthat the positive and negative effects of the five chemicalsupplements and their combinations were directly translatedinto an increase or decrease in time required to attain stationaryphase and the population size at which the stationary phasestarted. The stochastic ecological model proved to be useful,as it effectively explained and summarized the uncertainty seenin the recorded growth curves. Our findings have broad implicationsfor both basic and applied research and illustrate how stochasticmathematical modeling coupled with rigorous statistical methodscan be of great assistance in understanding basic processesin microbial ecology.

INTRODUCTION

Top

Introduction

Mathematical modeling coupled with rigorous statistical methodscan be of great assistance in understanding the interactionof organisms with their physical and biological environment(8, 47, 48, 50, 52, 56). Studies in the field of predictivemicrobiology have shown that successful modeling requires bothadequate models and thorough data sets (57; for extensive reviews,see references 10, 42, and 55). In predictive microbiology atwo-step modeling approach is used (reference 55 and citationstherein). Primary models describe the basic rules of how microbialnumbers change over time (52). Next, these simple models areused to derive secondary models that account for the effectof a set of factors in microbial growth.

The forces of the environment and the events of reproductionand growth are themselves stochastic in nature (3, 21, 24, 32,36, 39, 44, 45, 51, 53), yet simple ecological models, suchas the Verhulst logistic equation, result in deterministic predictions.However, smooth convergence to asymptotic results is not whatis usually seen, even in rigorous experimental settings (17,31). Hence, a more realistic alternative to population growthmodeling is to confront stochastic equations with the data athand (references 5, 14, and 17 and citations therein and references41 and 47).

In this paper, we describe a novel application of a stochasticpopulation model to analyze how environmental conditions influencemicrobial growth dynamics using an extensive experimental dataset. The primary model used was the stochastic Ricker (SR) equation(27, 51). Our secondary model was a novel use and applicationof the SR model in an analysis of variance (ANOVA)-like formatto account for differences between and within experimental conditions.We use rigorous statistical methods to characterize the effectsof chemical supplements on the following aspects of bacterialgrowth: the rate of attaining stationary phase, the populationdensity at which stationary phase occurs, and the variabilityassociated with the growth process. The experimental data weretime series measures of Escherichia coli growth in a basal mediumamended with all 32 possible combinations of five supplements,namely, glucose, NH₄Cl, HCl, EDTA, and NaCl (Fig. 1). All supplementcombinations were tested in triplicate. The methods and resultspresented here have broad implications and can be applied toany situation in which the experimentalist wishes to determinethe factors that affect growth responses of microbial isolates.

View larger version (54K):
[in this window]
[in a new window]
FIG. 1. Plotted are the three replicated growth curves (OD x 100 as a function of time [min]) of each of the 32 design types (or treatment combinations). Each replicated growth curve is marked with a different line type. In the titles above each plot, the factors present or absent are indicated by 0 or 1, respectively, with the sequence of treatments being glucose, NH₄Cl, HCl, EDTA, and NaCl.

The organization of this paper is as follows. In the first section,we present the experimental procedures used. Next, we presenta theoretical background section that is divided into threeparts. In the first part, we provide a brief description ofthe SR model and maximum-likelihood parameter estimation methods.In the second part, we present details of how we derived varioussecondary models to explain the data. Among those models, thebest one was chosen using information theoretic criteria (1,15, 54). This model selection procedure is described in thethird part of the theoretical background section. After presentingthe results, we discuss how the SR model compares to other stochasticmodels and how it can be used to approach the following importantproblems in predictive microbiology: (i) building a stochasticlag phase model and accounting for the three main phases ofmicrobial growth, (ii) accounting for time-varying environmentalfactors, (iii) accounting for microbial interactions in themodel predictions, and (iv) using an adequate model selectiontool that accounts for overparameterization and considers allthe models simultaneously instead of stepwise.

MATERIALS AND METHODS

Top

Materials and Methods

Growth curve experiments.
The organism used for the experiments was Escherichia coli HB101.As a basal medium, a diluted and modified Luria-Bertani (LB)broth was used (1 g/liter tryptone, 0.5 g/liter yeast extract,and 10 g/liter NaCl). The basal medium was supplemented withall 32 possible combinations of five chemicals in the followingconcentrations: glucose, 8 g/liter; NH₄Cl, 10 g/liter; HCl,0.0005 M; Na₂EDTA · 2H₂O, 0.05 g/liter; and NaCl, 20g/liter. Each growth medium was inoculated with a 10^–3dilution of an overnight (22-h) culture of the organism grownin undiluted LB broth at 37°C and shaken at 200 rpm. A 2⁵factorial design with 32 experimental treatments was achievedby transferring triplicate 200-µl aliquots of each mediumto a 96-well microtiter plate. This plate was sealed with abreathable membrane (Breathe-Easy; Diversified Biotech, Boston,Mass.) and incubated in an automated plate reader (PowerWaveHT; Bio-Tek Instruments, Inc., Winooski, Vt.). The plate wascontinuously shaken between readings at level 4, and opticaldensity (OD) values of the cultures were measured at 15-minintervals. Thus, 96 growth curves, each consisting of 193 measurements,were obtained. The experimental design is summarized in theplot titles in Fig. 1. Since the maximum recorded OD value was0.362, we assume a linear relationship between OD and populationsize throughout this work.

Theoretical background. (i) Primary model: the SR model.
We modeled the growth dynamics of each treatment combinationusing a stochastic logistic growth model based on the Rickermodel (51). The Ricker model has a long history in populationecology modeling (35) and has been used as a discrete versionof the well-known Verhulst logistic differential equation (40).The Ricker equation expresses the one-step-ahead populationsize as a function of the current population size and includesa density-dependent effect:

(1)
LetN_t be the population size at time t. The parameter a can bethought of as the speed at which equilibrium is reached (mathematically,it is the eigenvalue of the Jacobian matrix, found after linearizationaround the equilibrium [39, 40]). The parameter b representsthe effect of the current population size on the rate of growth.For example, if b is <0, then the current population sizehas a negative effect on growth, thus expressing density dependence.In that case, the model has a nontrivial stable equilibriumwhose position is given by N = –a/b. The dynamic behaviorof this model up to and including chaos is well known (17, 27,40).

A stochastic version of this model bridges the gap between theoreticalideas and available data sets. The stochastic formulation ofthe deterministic model leads to a hypothesis that explainshow the departures from the deterministic predictions occurin a given time series of population growth (17). The sourcesof variability may come from observational error or from processerror. Observational error relates to the fact that at a giventime, the total number of individuals cannot be known exactly,and the experimenter's estimate is based on a sample. Processerror comes from uncertainty inherent in the process of growthitself. Population dynamics theory further subdivides the processerror into "demographic" and "environmental" stochasticity.Demographic stochasticity represents the variability due torandom contributions of births, deaths, and migrations of individualsin the population (17, 33). Environmental stochasticity representsthe effect of external factors on the individuals of the population(23). For example, in Baranyi's stochastic models (5, 6, 7),the equations' coefficients are random variables that representthe biological variability between the individual cells in thepopulation. In the ecological theory context, this is a modelfor demographic stochasticity (23). Cushing et al. (17) haveshown that considering both demographic and environmental stochasticityis essential to adequate population dynamics modeling but thatintegrating both sources of uncertainty in a single model isnot an easy task.

Because we were dealing with an experiment in which a givenenvironment was imposed on bacterial populations, we chose tomodel environmental stochasticity alone. The stochastic Rickermodel is written as follows (27):

(2)
where ${sigma}$ E_t is a random shock to the population growth rate at time t,and the E_t are independent and normally distributed with meanzero and variance one. Because the current population size dependson the previous observation, this model has the Markov propertyand in the log scale it becomes a first-order nonlinear autoregressivemodel:

(3)
where X_t = ln N_t.

Setting a equal to 0 and b equal to 0 defines a discrete-timeBrownian motion process with zero drift where there is no populationdensity feedback typical of ecological processes (11, 23, 27).When a is not equal to 0 and b is equal to 0, the model is alsoa discrete-time Brownian motion with added drift. If a is notequal to 0 and b is not equal to 0, the model includes (positiveor negative) density dependence. In particular, the model inwhich a is not equal to 0 and b is <0 represents a stochasticlogistic growth. Under this model, the population no longerattains a single deterministic equilibrium, as in the Rickerequation, but instead, it approaches a "cloud of points" (60),a stationary distribution which can be approximated by a gammaprobability density function (20, 24) whose mean is –a/b.The point –a/b represents a center for return tendencies:it is the population abundance at which the average change isN_t, conditional on N_t_–₁ being zero, thus accounting forthe stationary phase. Note that the stochastic Ricker modeldoes not include a term to account for the lag and death phasesof bacterial growth in a batch culture, but as will be shown,the model serves as a very good approximation of these growthphases during the monitored period, despite the fact that thedeath phase is often evident. We later propose and briefly exploresimple modifications to the stochastic Ricker model and othermodels to account for the processes occurring during these twophases of bacterial growth in batch cultures.

The problem of connecting a time series of observations withthe proposed model in equation 2 requires the specificationof a likelihood function. For each of the treatments, we estimatedthe corresponding parameters of the SR model as follows: supposepopulation abundances of a single culture are observed fromtime 0 to q. The likelihood function is a function of the unknownmodel parameters ${theta}$ = [a, b, ${sigma}$ ²]'. It is the joint probabilitydensity function for the vector of random variables X = [X₁,X₂, ..., X_q]' conditional on X₀ = x₀ and evaluated at the recordedvector of values x = [ln n₀, ln n₁, ..., ln n_q]' = [x₀, x₁,..., x_q]' (25, 27). Because of the Markov property, the jointprobability density function of the observed data, given ourproposed model, is just the product of the individual transitiondensity functions. The maximum-likelihood parameter estimates(MLEs) are the parameter values that make the observed data"most probable" or "most likely," i.e., they maximize the likelihoodfunction for time series data, which is given by (27):

(4)
Let Y_t = X_t – X_t_–₁ denote the one-stepdifferences of the logarithmic population size. Dennis and Taper(27) show that the MLEs of a, b, and ${sigma}$ ² under this model areidentical to the least-squares estimates obtained by a linearregression of the observed y_t on n_t_–₁, where t = 1,2,..., q. This makes parameter estimation a straightforward task,and many commercial statistical packages can be used. We notethat the confidence intervals returned by those packages arenot correct in this case, because the observations at each timestep are not independent of each other. Furthermore, the valueb = 0 is at the edge of the set of values b < 0 for whichthe stochastic process N_t is ergodic, thus violating one ofthe regularity conditions under which a ${chi}$ ² approximation to thelikelihood ratio test is valid. Using the ${chi}$ ² approximation thusyields inflated type I error rates (27). Instead, confidenceintervals for the parameter estimates have to be found usingthe parametric bootstrap (PB) (29, 38). The PB of a stochasticmodel of population dynamics involves the following steps (26,27). From a given data set, the ML estimates of the parametersare calculated for the chosen stochastic model and used to simulatemany time series data sets (e.g., 2,000) of the same lengthas the original data set, using the same model. The ML estimatesare then found for each of these data sets. Their histogram(or kernel density estimate) represents an estimate of theirsampling distribution, from which various descriptive statisticscan be obtained (e.g., the 2.5 and 97.5 percentiles).

(ii) Secondary model: an ANOVA variant of the SR model.
In predictive microbiology, secondary models describe the dependenceof primary-model parameters on environmental factors. Here,the primary-model parameters are contained in the stochasticgrowth rate function (i.e., the terms in brackets in equation2). Hence, the natural way to rigorously determine whether acertain combination of environmental factors affects the growthrate is to adopt an approach that is conceptually "ANOVA-like"and test whether there are differences in the estimated modelparameters within and between experimental scenarios. Inferenceusing this approach was greatly facilitated in this study becauseour data set was a perfect 2⁵ factorial design with three replicatesof each experimental treatment. If the general hypothesis thateach of the 32 treatments has a distinct effect on the populationdynamics of bacteria is supported by the data, then the speedof attaining stationary phase, the population size around whichthe cell counts vary in stationary phase, and the variabilityof the process (i.e., each of the SR model parameters) shoulddiffer between the 32 experimental treatments. To estimate themodel parameters for a single treatment, the joint likelihoodof the three replicated time series has to be considered. Becausea treatment's replicates were independent, the likelihood functionfor a given experimental treatment was taken as the productof the individual likelihoods. Note that although the idealsituation is having a full factorial design, just like in thetraditional ANOVA methods, the overall procedure could be modifiedto account for unbalanced designs.

Various hypotheses were contained in the general hypothesisdescribed above, and modified multivariate techniques can beused to identify patterns in the SR model parameter space. Ifthere is no difference between the model parameters found in,e.g., two experimental treatments, the data for these two treatmentscan be pooled and the SR model parameters can be estimated anew.In that fashion, the total number of possible treatment effects,combinations, and interactions can be reduced and the main factorsthat drive the system dynamics can be identified. The patternsof variation in the 96 data sets were summarized in order topropose a set of competing hypotheses consistent with the observationsand to group the data accordingly. To do so, the model parametersfor each of the 96 data sets were calculated as described before.These values were used as coordinates in a three-dimensionalspace. Canonical-variates analysis was used to identify patternsof variation among treatments. The design types were plottedin a two-dimensional canonical-variates space. A careful inspectionof the canonical-variates plot allowed the formulation of sevenhypotheses that were consistent with the data. Each hypothesisgrouped all of the 32 experimental treatments in a particularway and proposed that distinct population dynamics occurredwithin each group. We then tested and compared the degree towhich each of the seven hypotheses was in agreement with theobserved patterns of variation. Confidence intervals using aPB were calculated for the model parameters under each hypothesis(see reference 27 for details).

(iii) Model selection.
Statistical theory and evidence from many biological disciplines(15) show that traditional stepwise regression methods basedon a series of likelihood ratio tests may miss the best modelor hypothesis consistent with the data. Also, a series of pairwisecomparisons can lead to erroneous P values in likelihood ratiotests and inflated type I errors. Therefore, we relied hereon the Akaike information criterion (AIC) (1) and on the Bayesianinformation criterion (BIC) (54) to simultaneously assess thequality of each of those hypotheses to explain the data.

The use of the AIC and BIC for hypothesis selection has a strongtheoretical rooting in information theory. For a given dataset, the AIC gives an estimate of the expected, relative, directeddistance between the fitted model and the unknown true mechanismthat generated the data (15). Thus, the decision rule for modelselection using those statistics is to choose the model withthe lowest AIC or BIC. For a fixed data set, adding more parametersto the model reduces that distance but further increases uncertaintyin the estimation process. That trade-off between underfittingand overfitting is directly expressed in both the AIC and BICas a term that penalizes their scores as a function of the numberof estimated parameters in the model. We note that Schwarz (54)showed that if the true model is within the suite of evaluatedmodels, then the BIC is guaranteed to find it. Hence, we usedboth statistics to derive the conclusions of the hypothesisselection. For the stochastic Ricker model, the AIC is equalto:

(5)
where ln()is the likelihood function evaluated at the MLEs and p is thenumber of model parameters. The BIC is calculated with:

(6)
where q is the number of data points used inthe parameter estimation process. A disagreement between thetwo statistics would indicate that there is not enough evidencein the data to support the best model, and a decision wouldhave to be taken after investigating the type I error ratesof each model using extensive simulations (26, 33a).

Finally, an evaluation of the quality of the best-fitted modelwas done via a residual analysis. Under the proposed model,the residuals should be normally distributed, centered aroundzero, and nonautocorrelated. A strong deviation from normality,if it appears, is an indicator that the current model mechanismis insufficient to explain the available data. All the calculationswere done using MATLAB 6.5.1, release 13 (The MathWorks, Inc.,Natick, Mass.).

RESULTS

Top

Results

Growth curve experiments.
A biological data set was obtained by growing E. coli in a microtiterplate in a base medium that contained all 32 possible combinationsof five supplements (glucose, NH₄Cl, HCl, EDTA, and NaCl) intriplicate. The growth that occurred with each supplement combinationwas recorded using an automated plate reader. Figure 1 showshow the recorded population size (OD x 100) varies for eachof the 32 treatments.

Primary model: the SR model.
The SR model predictions found through a process error fit weresuperior to those generated by alternative models using an observationerror fit, as is clearly illustrated in Fig. 2. In fact, evenimproved observation error models including a death phase donot fit as well as the SR model. This figure also illustratesthe fact that a process error model fit is conceptually differentfrom an observation error fit in that it naturally yields astochastic one-step-ahead prediction of future population sizes.The advantage of this becomes evident while plotting the predictionsfrom both types of statistical fit (Fig. 2). The residuals obtainedby an observation error fit of the Ricker model were at least1 order of magnitude larger than the process error fit residuals.A typical deterministic prediction using an observation errorfit differs dramatically from the one-step-ahead predictionsgenerated using the stochastic Ricker model parameter estimates.

View larger version (19K):
[in this window]
[in a new window]
FIG. 2. Stochastic (Stoch.) Ricker model predictions compared to the observed data, an observation error fit of the Ricker equation, and an observation error fit of an ordinary differential equation (ODE) model that accounts for the rate of substrate consumption. In the ODE model, n is the population abundance, r is the resource variable, and ${phi}$ _max, k, µ, and ${alpha}$ are constants:

This model assumes that a population is growing according to the well-known Monod function and that individuals are dying at a rate µn.

The estimates of the parameters a, b, and ${sigma}$ ² of the stochasticlogistic growth models fitted to each of the experimental growthcurves are reported in Fig. 3a to c. In Fig. 3d to f, the parametersfor supplement combinations 26 through 32 were omitted for clarity.

View larger version (47K):
[in this window]
[in a new window]
FIG. 3. Box plots of the stochastic Ricker model parameter estimates as a function of the design type for treatments 1 to 32 (a, b, and c) and for designs 1 to 25 (d, e, and f). The bar inside the box indicates the median, the boxes delimit the interquartile range, and the vertical lines extend to the maximum and minimum.

Secondary model: an ANOVA variant of the stochastic Ricker model.
After estimating the parameters using the population dynamicsmodeling approach for each of the 96 growth curves, we groupedthe data in various ways to establish how the data could bestbe explained. We used a graphical approach to develop sevenhypotheses to explain the treatment groupings. Further dataanalysis showed that not all of the 32 different chemical supplementcombinations had a significant effect on the population dynamicsof the microbial cultures. The data were best explained withjust four distinct dynamics, corresponding to the cases in whicheither NaCl or EDTA was present, both were present, and bothwere absent. Furthermore, we found a positive interaction effectbetween these two supplements.

The MLEs of the stochastic Ricker model for each of the 96 datasets were used to run a canonical-variate analysis. Fisher'sdiscriminant function based on the estimated parameter valuesof the stochastic Ricker models identified the vectors in themultidimensional space on which, when the data points were orthogonallyprojected, the experimental designs were maximally separated(34). The test for the relative contributions of the three selectedeigenvalues and the total canonical structure extracted by Fisher'sdiscriminant function indicated that only the first two eigenvalueswere significant (P values of < 0.001, < 0.001, and 0.0684,respectively). Furthermore, they cumulatively explained 87.26%of the variation in the three-dimensional space of [a,b, ${sigma}$ ²].A closer examination of the canonical structure indicated thatthe extracted canonical variate number 3 was not very usefulto explain the variation in the three-dimensional space. Theplot of the designs in the discriminant space (Fig. 4) was usedto identify possible design groupings and to formulate the followinghypotheses.

View larger version (10K):
[in this window]
[in a new window]
FIG. 4. (Left) Experimental treatments plotted in a two-dimensional space according to their mean coordinate values in canonical variates 1 and 2. (Right) Same as before, but treatment 26 is excluded for clarity. The two canonical variates represent the two axes in the three-dimensional space of a, b, and ${sigma}$ ² along which the highest separation of the data points is obtained, i.e., these are the vectors along which most of the variability seen in the values of the model parameters occurs. Hence, two treatments for which the biochemical environment leads to similar growth patterns will appear closer when plotted in the two-canonical-variates space.

(i) Hypothesis 0.
Only two groups (A and B) can be identified: group A consistsof designs 26 to 32, and group B consists of designs 1 to 25.The first group corresponds to the case where EDTA and NaClare both present at the same time regardless of any other supplementaddition. The second group corresponds to the case where EDTAand NaCl are not present at the same time (Fig. 1). Group Adramatically increased the value of a, the rate at which stationaryphase is reached (Fig. 3a), and lowered the mean level (–a/b)at which this occurred by decreasing the value of b (Fig. 3b).

(ii) Hypothesis 1.
Because the variability in the data coming from design 26 wasmuch higher than the rest (Fig. 3c), this design can be separatedfrom the rest and constitute another treatment by itself, besidesthe two proposed by hypothesis 0. Hence, this hypothesis proposesto divide the designs into three groups, corresponding to threeseparate population dynamics types.

(iii) Hypothesis 2.
Four groups explain the data: designs 27 to 32, in which EDTAand NaCl are both present and at least one of NH₄Cl and HClis present; design 26; designs 16 and 24; and all other designs.

(iv) Hypothesis 3.
Four experimental groups defined by just two supplements explainthe data. The treatments are EDTA and NaCl, and the experimentalgroups are as follows: EDTA and NaCl absent, EDTA present andNaCl absent, EDTA absent and NaCl present, and both present.This hypothesis suggests that there is a significant interactionbetween EDTA and NaCl and that the effects of that interactiondrive the dynamics of the system.

(v) Hypothesis 4.
Hypothesis 4 considers the same experimental groups as in hypothesis3 but also includes the presence or absence of glucose as afactor. Hence, a total of eight experimental groups are proposedby this hypothesis.

(vi) Hypothesis 5.
All of the supplements except glucose have an effect, and interactionsexist between EDTA plus NaCl and NH₄Cl plus HCl. There are atotal of 16 experimental groups.

(vii) Hypothesis 6: full model.
The data at hand can be best explained by proposing 32 distincttypes of population dynamics types corresponding to each ofthe 32 experimental treatments with three replicates each.

Model selection.
To evaluate the consistency of each hypothesis with the data,the BIC and AIC statistics were calculated from the likelihoodof each hypothesis evaluated at the MLEs. To estimate the parametersunder each hypothesis, we pooled the replicates of designs containedwithin each group proposed by the seven hypotheses. The jointlikelihood function for each group was computed as the productof the individual time series likelihoods.

Both the AIC and the BIC clearly favored hypothesis 3 (Table1), implying that a good explanation of the changes in a, b,and ${sigma}$ ² was obtained while considering just the presence or absenceof EDTA and NaCl. With the parameter estimates and confidenceintervals obtained using PB (Table 2), interaction plots weredrawn (Fig. 5). The interaction plots show that the order ofthe effect of NaCl changed as the state of the EDTA treatmentchanged. Thus, the effects of these two factors were not simplyadditive. When the data coming from the treatments with theabove-mentioned interaction were excluded, then the effectsof the other factors became visible (Fig. 3d to f). The BICvalues of hypotheses 0, 1, 2, 4, and 5 show that even when othertentative explanations may be valid and better than the fullmodel, the patterns of variation in the time series cannot beexplained as well as with hypothesis 3. The confidence intervalsobtained with a PB of hypothesis 5 showed that the interactionsbetween NaCl plus EDTA and NH₄Cl plus HCl were not significant(Fig. 5). Thus, in hypotheses 4 to 6, the treatment effectswere just additive and could not generate as much variationas the one generated by the interaction between EDTA and NaCl.The second-best hypothesis was the reduced hypothesis (hypothesis0), which considered just two groups. This hypothesis had thebenefit of a reduced number of parameters to estimate. However,with just two groups it is not possible to estimate the interactioneffects.

View this table:
[in this window]
[in a new window]
TABLE 1. Hypothesis selection^a

View this table:
[in this window]
[in a new window]
TABLE 2. Stochastic Ricker model parameter estimates and PB confidence intervals for the model specified by hypothesis 3^a

View larger version (24K):
[in this window]
[in a new window]
FIG. 5. Interaction plots for hypothesis 3 and hypothesis 5. (Left) Closed circles represent absence of NaCl, open circles, presence of NaCl. (Right) Open circles represent absence of both NH₄Cl and HCl, closed circles, presence of both chemicals. Parametric bootstrap confidence limits are marked with crosses.

Once chosen, the best hypothesis was further scrutinized viaa residual analysis. The model residuals for the best hypothesis,hypothesis 3, were approximately normally distributed and centeredaround zero, suggesting that the stochastic Ricker model wasa good first approximation to differentiate between variouspopulation dynamics (Fig. 6). However, some deviations fromthe assumptions occurred in the extreme data points; at theearly stage of the process, during the lag phase; and as expected,at the late stage of the process, during the death phase. Thosedeviations were accompanied in some cases with an autocorrelationpattern, which suggested a significant effect of an observationalerror.

View larger version (17K):
[in this window]
[in a new window]
FIG. 6. Plotted are the quantile-quantile graphs of the residuals for each of the four experimental design groups in hypothesis 3 (see "Model selection" for details). The points far from the straight lines in the lower right plot represent the individual model residuals that deviated from the model normality assumption.

DISCUSSION

Top

Discussion

Uncertainty in the trajectories of population growth seen inour 96 microbial growth curves was effectively explained andsummarized by the stochastic Ricker model. It was shown thatthe positive and negative effects of the five supplements andcombinations of them directly resulted in an increase or decreasein the time required to attain stationary phase and the populationsize at which the stationary phase started. Furthermore, thepotential complexity of 2⁵ experimental treatments and theireffects was reduced to 2², as just the metal chelator EDTA,the presumed osmotic pressure imposed by NaCl, and the interactionbetween these two factors were enough to explain the variabilityseen in the data. These findings seem to be in agreement withprevious studies (30, 49, 57).

Model selection via the AIC and BIC overcomes two major problemsreported in the literature. The first one is overparameterization.In the search for more mechanistic models, many parameters areintroduced, for example, to model the probability of growthas a long nonlinear function (57). Baranyi et al. (9) have alreadywarned about the effects of overparameterization in predictivemicrobiology. Our results show that by using the AIC or BIC,a good compromise between the number of parameters includedin the model and the quality of the predictions is reached.The second problem is stepwise regression as a model selectiontool. As mentioned before, it has been repeatedly shown thatstepwise regression can lead to an erroneous choice of the bestavailable model. The information criteria used here overcomethat problem by evaluating all the models at the same time,thereby allowing a proper identification of the best availablemodel.

Process error fitting is a more realistic and mechanistic modelingapproach. It is a realistic approach because it aims to explainthe variability seen in time series of population abundanceby proposing stepwise stochastic changes in the growth rateas a function of the experimental medium (Fig. 2). It is a mechanisticapproach because it effectively explains and reproduces theobserved growth patterns through proposed biological processeswell rooted in first principles. Our model basically statesthat exponential growth and density-dependent effects are randomand change as a function of the growth medium. The effects ofthe environmental factors were directly translated into changesin the basic growth characteristics.

Modeling the process error is not a mere fitting technique.Deterministic mathematical models (e.g., Gompertz and logisticequations) are often proposed as the true underlying "perfect"mechanism and fitted to the data via nonlinear fitting techniques(8, 13, 16, 19, 52, 59). Deviations from the deterministic predictionsare accounted for as observational uncertainty. On the otherhand, the SR model asserts that the process of population growthis stochastic by nature and can be modeled with a Markov process.The transition probability density function of this Markov processis used to predict trends and also as a natural tool for parameterestimation via ML. In the presence of substantial observationerror, the process error model ML parameter estimates and predictionsare biased, and the converse is also true (25, 28, 43, 46).However, in the context of this paper, because the microtiterplate reader is highly accurate (coefficient of variation =standard deviation/mean = 1.08%) and large populations are beingsampled, observation error can be neglected. Further evidencesupporting the use of the SR model was given here by the quantileresidual plots (Fig. 6), in which generally the assumptionsof normally distributed errors were met. However, we saw deviationsfrom the normal model at the early stages (lag phase) and towardthe end (death phase). Early deviations can possibly be additionallyaffected by the detectability limitations of the use of OD asa surrogate for population size (18), as well as by a lag phase.

The SR model can be modified to accurately predict more complexbiological phenomena of bacterial batch cultures, such as thelag and death phases, the effect of time-varying environmentalfactors, and multiple-species microbial interactions. Well-knownideas in ecology can be translated into stochastic populationmodels to explain the data. As an example, suppose that at smallinitial population sizes bacteria cannot grow well due to harshenvironmental conditions, e.g., low pH values. This situationmight change at larger population sizes, as metabolic productscould increase the pH. That is, at intermediate population sizes,the population might grow much better than at small populationsizes. Furthermore, a critical initial population size mightexist below which bacteria cannot grow and above which theysucceed. This phenomenon is well known in ecology as the Alleeeffect (2, 20), and stochastic population models exist in theliterature that account for it (20, 21). A stochastic modelingapproach to this problem could, for instance, lead to accurateestimates of that critical initial population size below whichbacteria fail to grow.

To account for environmental stochasticity and the fluctuationsover time in resources and in the concentrations of the chemicalsupplements, the stochastic Ricker model could be easily modified.To include the main effects of the chemical factors withoutthe interactions, the growth function in the exponential termof the stochastic Ricker model could be written as follows:

where c₁,c₂,...,c₅ are parameters to be estimatedand S_it is the concentration of the chemical supplement S_i attime t. This model has been successfully used to model longtrends in populations whose growth depends on climatic fluctuations,such as "El Niño" rainfall (26). In the model describedabove, the chemical supplements are included as covariates thatmay enhance or suppress bacterial growth, depending on the valuesof the constants c_i. This model has the advantage that it takesinto account the variation of the concentration of the supplementsover time, and substrate depletion can be easily accounted for.However, to estimate the parameters in this model, the concentrationof each chemical supplement and substrate at each time stepis needed (12). In the context of our experiment, measuringthese concentrations was not experimentally feasible. We optedto sacrifice the number of variables recorded in order to havea large number of data points and a perfectly balanced factorialdesign. From the point of view of improving statistical power,it is preferable to have many replicates of the same univariateprocess rather than to gather data on a large multivariate processwith small sample sizes. Another extension of the model wouldbe to consider several levels of the initial concentration ofthe supplements, and like most ANOVA models, our current approachcan be easily expanded in such a way.

Accounting for the population lag phase (in the sense of Baranyi[5]) is also straightforward using the SR model. During lagphase, bacterial numbers seem to fluctuate randomly around theinitial density. However, the case of zero growth is just aspecial case of the SR model. Indeed, no growth is equivalentto setting a equal to 0 and b equal to 0 in equation 2 (see"Theoretical background" above). Then, accounting for the lagphase amounts to specifying what is known in the statisticsliterature as a stochastic change point problem. While our modeleffectively accounts for the pattern of variability seen inthe data, there are other aspects of growth dynamics that couldbe modeled in future work. For example, delay differential equationmodels, like the one proposed by Baker et al. (4), are usefulin estimating the fraction of cells that are dividing, the degreeof initial synchronization of the cells, and the initial distributionof cells in the cell cycle. It would be interesting to extendBaker's model to a stochastic differential equation and properlyaccount for stochastic effects. However, parameter estimationvia ML for stochastic differential equation models is a wide-openresearch area and presents many challenging statistical problems.These topics constitute materials for future research.

Finally, another of the latest concerns in predictive microbiologyis the consideration of microbial interactions (37, 58). A two-speciesdeterministic discrete Ricker model exists in the literature(40). This model can be easily transformed into a stochasticmodel in the form of equation 2 (22, 28). Preliminary resultsnot presented here show that this competition model can be effectivelyused to analyze the seminal data set of Gause (31).

Stochastic modeling techniques are maturing, and there is noreason why theoretical and applied studies in microbiology shouldbe deprived of such useful mathematical statistics tools andconcepts. The results of this paper have broad implicationsfor both basic and applied research and can be regarded as adifferent starting point to fulfill the urgent need of simplestochastic models of microbial growth.

ACKNOWLEDGMENTS

This research was made possible by National Institutes of Healthgrant number P20 RR 16448 from the COBRE Program of the NationalCenter for Research Resources (to L. J. Forney); a fellowshipfrom the Inland Northwest Research Alliance (INRA) SubsurfaceScience Research Institute, which is funded by the Departmentof Energy under contract DE-FG07-02ID14277 (to F. P. J. Vandecasteele);and the National Science Foundation (NSF DEB-0089756 and NSFDMS-0072198) (to P. Joyce).

We also thank Brian Dennis for reading the manuscript and providingmany insightful comments and ideas.

FOOTNOTES

* Corresponding author. Mailing address: Department of Mathematics and Department of Statistics, University of Idaho, Moscow, ID 83844-1103. Phone: (208) 885-6338. Fax: (208) 885-5843. E-mail: joyce{at}uidaho.edu.

REFERENCES

Top