ABSTRACT
A few bacterial cells may be sufficient to produce a food-borne illness outbreak, provided that they are capable of adapting and proliferating on a food matrix. This is why any quantitative health risk assessment policy must incorporate methods to accurately predict the growth of bacterial populations from a small number of pathogens. In this aim, mathematical models have become a powerful tool. Unfortunately, at low cell concentrations, standard deterministic models fail to predict the fate of the population, essentially because the heterogeneity between individuals becomes relevant. In this work, a stochastic differential equation (SDE) model is proposed to describe variability within single-cell growth and division and to simulate population growth from a given initial number of individuals. We provide evidence of the model ability to explain the observed distributions of times to division, including the lag time produced by the adaptation to the environment, by comparing model predictions with experiments from the literature for Escherichia coli, Listeria innocua, and Salmonella enterica. The model is shown to accurately predict experimental growth population dynamics for both small and large microbial populations. The use of stochastic models for the estimation of parameters to successfully fit experimental data is a particularly challenging problem. For instance, if Monte Carlo methods are employed to model the required distributions of times to division, the parameter estimation problem can become numerically intractable. We overcame this limitation by converting the stochastic description to a partial differential equation (backward Kolmogorov) instead, which relates to the distribution of division times. Contrary to previous stochastic formulations based on random parameters, the present model is capable of explaining the variability observed in populations that result from the growth of a small number of initial cells as well as the lack of it compared to populations initiated by a larger number of individuals, where the random effects become negligible.
INTRODUCTION
Often bacterial contamination of foods starts with a small number of bacteria that are capable of adapting and proliferating by repeated divisions on a given food matrix. At low cell concentrations, standard deterministic models fail to predict the variability of the bacterial population. This is so because at low initial cell numbers, heterogeneity between individuals and its influence on the division times become relevant and have a net influence on the population. Consequently, the behavior of individual bacteria cannot be neglected when assessing possible health risks along the food chain, either during storage or during distribution.
Recently, attention has been drawn to the need for modeling and simulation methods to observe and describe the variability of single-cell behavior and small populations (1, 2) in order to produce realistic estimations of safety risks along the food chain, for instance, during storage and distribution or during food processing.
In this paper, connections are established between individual bacterial growth and division, the corresponding distributions of times to division, and population growth curves that in the long term can aid the quantification, on a probabilistic basis, of microbial risk and product shelf life.
A number of modeling approaches for bacterial population dynamics have been built around the concept of times to division and particularly the first lag time to division of cell populations (3–6). The development of analytical techniques capable of measuring single-cell parameters such as cell length renewed interest in modeling single-cell growth. Experimental techniques for single-cell studies include turbidimetry (see for example, references 7 and 8), lithographic techniques (9), and flow cytometry (10). Recently, time-lapse microscopy has been successfully applied to obtain quantitative information on colonial growth dynamics originated from single cells (2). Based on the flow chamber microscopy technique proposed in reference 10, models for individual cell growth have been developed as described in references 4 and 11. They are essentially adaptations of the now classical model proposed in reference 12 to describe the growth of bacterial populations before attaining the stationary phase. It consists of two ordinary differential equations that can be written as
(1)
(2) where y in equation 1 represents the natural logarithm of the population size and μ the maximum specific growth rate. Variable a(t) in equations 1 and 2 is known as the adjustment function and relates to a certain physiological state of the population. This variable has been introduced to describe the gradual adaptation of the cells to the new environment. Initially, it takes a small value which increases at a rate proportional to ν (in equation 2) up to a maximum of one. This variable induces a delay in the growth of the population size which depends on the inverse of v (the larger the rate, the smaller the delay becomes) and is employed to model the observed initial lag time.
In the context of a single-cell growth, the state variable y, which relates to the size of the population in reference 12, is identified as a critical variable that determines cell growth and characterizes division when it reaches a particular threshold. In previous works, it has been interpreted as cell length (see, e.g., references 4 and 11) or cell DNA content (13). Similarly, the adjustment function in reference 12 is used in references 4 and 11 to model the adaptation of a given cell to the environment.
Stochastic fluctuations in gene transcription and translation within a cell or in response to cell-to-cell disturbances are considered the main sources of heterogeneity among individual cells within a population (14, 15). In order to capture such “behavioral noise” (2), previous approaches combined deterministic equations for single-cell or population dynamics of the form of equations 1 and 2, with random parameters and/or initial random conditions (both described by appropriate probability distribution functions).
In reference 11, the random nature of cell division is modeled by imposing a uniformly distributed random length threshold at which each cell divides. This seems to be also the case in reference 4, although parameter estimation for cells subject to different heat shock treatments is based on a deterministic model. However, the connections between such models and the observed distributions of division times have not been clearly discussed yet.
Inspired by reference 3, a model of cell population dynamics is suggested in reference 16 that explicitly includes a lag time distribution function. This is in agreement with other well-known approximations such as the Weibull or gamma probability distribution functions (8, 17). Recently, a model similar to that proposed in reference 12 with parameters following random distributions has been employed as suggested in reference 2 to simulate population growth rates.
In reference 5, a population dynamics stochastic model has been proposed for Escherichia coli that makes use of time-to-division distribution functions to compute the time (the stochastic variable) at which cell division occurs. This information is provided to the algorithm either directly from experimental data (e.g., flow chamber microscopy or turbidimetry) or from a gamma distribution, which has been previously fitted to experimental data.
In this work, connections between individual bacterial growth variables and distribution of times to division are established by a stochastic version of the models proposed in reference 4 and in references 11 and 18. The underlying premise is that cell growth is the result of a large number of biochemical reactions taking place on a microscopic domain (thus involving a relatively small number of molecular species). Standard assumptions can then be invoked to relate a chemical (or biochemical) microscopic master equation to its mesoscopic (chemical Langevin) counterpart (19), which is nothing but a stochastic differential equation (SDE) (20).
Hence, heterogeneity between individuals or fluctuations within each individual (e.g., due to behavioral noise) is represented by an SDE that will result from adding a stochastic component to the specific growth rate. The introduction of stochasticity on the growth rate has a strong biological interpretation: SDEs are usually employed in a systems biology context to collect the random effects of fluctuations on the system. Gillespie (19) gives convincing arguments to show how the accumulation of stochastic events can be cast into SDEs, stating that the aggregated effect of many events at the cell level can be captured and described by SDE models. Models based on SDE systems have been used previously to describe cell population growth and division for plankton (13) and bacteria (21). This approach has also been adopted in the study of bacterial systems under the action of bacteriophage (22) or antibiotics (23).
In the context of cell growth and division, such representation, which seeks the aggregation of the undergoing biochemical processes during the cell cycle, is shown to reproduce reasonably well the time to division distributions observed and reported in the literature.
It is well known from stochastic systems theory (24) that the collective effect of an SDE system, namely, the evolution of the probability distribution associated with the random state variables, can be computed as the solution of two partial differential equations (PDE): the so-called Kolmogorov equations, forward or backwards in time.
We made use of one such equation, the backward Kolmogorov, to characterize time to division distributions (TTD)—including first time to division—and to efficiently estimate parameters of the underlying stochastic dynamics from the experimental distributions. It is important to emphasize here that from a computational point of view, this way of approaching the model calibration problem is particularly efficient. This is so because it only requires the solution of a partial differential equation, as opposed to obtaining a complete ensemble of realizations by repeatedly solving the SDE to reconstruct the density function associated to the distribution of times to division (e.g., by Monte Carlo methods).
Similar approaches have been employed previously, although in the context of gene expression networks to extract kinetic parameters associated with individual cells from protein distributions obtained by cell population measurements (25). The use of stochastic methods related to Kolmogorov equations in predictive microbiology has been suggested previously (26), albeit in the context of population growth dynamics, as a means to describe variability of the environment as well as uncertainty due to limitations of the measurement equipment. However, to the best of our knowledge, these methods have not been employed so far either to model single-cell growth kinetics or to estimate parameters based on experimental TTD. Our model assumes that bacterial division is subject to stochastic fluctuations which integrate the effect of transcription at the level of each individual. From that point of view our approach follows the assumptions implicit in reference 25, in which a master equation is employed to connect parameters associated with a stochastic process (a signaling network) to observations obtained at the cell population level. In this work, we used observations obtained at the population level (time-to-division distributions) to estimate parameters of the stochastic process that describes division.
MATERIALS AND METHODS
Stochastic model for single-cell growth and division.In order to characterize the variability of single-cell kinetics within a population, a stochastic (SDE-based) version of the single-cell growth model discussed in reference 4 is proposed. The model is formally similar to the one represented by equations 1 and 2, although y in our case relates with size (length) of a single cell instead of number of individuals, and a with its corresponding adjustment function reflecting the physiological state of the cell. Growth starts at a given initial length, x(0) (equal to x0), and a(0) is equal to a0,with a0 being a small quantity provided that the cell undergoes a first division. The adjustment function induces a delay in x which mimics the initial adaptation of the cell to the new environment (the lag phase). After this period, growth proceeds exponentially up to a threshold value, x(T), equal to xdiv, which determines the division of the cell in two daughter cells. The time, T, at which such value is reached defines the time to division, which when it occurs for the first time also includes the lag phase period.
The model we propose assumes that the specific cell growth rate, μ, is subject to a stochastic fluctuation δW characterized by a Wiener process (see, for example, reference 20). For the sake of completeness, the main characteristics of the Wiener process and its role in the statistical properties of the corresponding random variable Y(t) are discussed in the supplemental material. Accordingly, the dynamics for single-cell growth is written as a linear time-dependent stochastic differential equation:
(3) where {Y(t)∈R; t > 0} denotes the random variable which takes values y (y = ln x) at t with a probability P(y,t | y0,0). Function P(y,t | y0,0) represents the conditional probability of Y(t) being equal to y given Y(0) being equal to y0. Function a(t) is the particular solution of equation 2 for a(0) (equal to a0), being of the form
(4)
The first term at the right-hand side of equation 3 is known as the drift and collects the deterministic size growth dynamics. Stochasticity is added in the second term on the right-hand side of equation 3, where a new parameter, ξ, is postulated that expresses the intensity of the stochastic fluctuation.
Although in practice it is accepted that the size (length) to division (and thus the size [length] of the resulting daughter cells as well) is randomly distributed (11), a minimum critical length seems essential to trigger the process (4). In the formulation we present in this paper, this random effect is aggregated together with other sources of variability within and between individuals, into the stochastic part of the growth dynamics (the second term on the right-hand side of equation 3).
Similarly to what has been proposed in reference 26 in the context of stochastic population dynamics, more elaborated formulations of the stochastic model 3 are possible, which can incorporate distinct sources of stochasticity. These might include, in addition to Y, a random equivalent of the adjustment function, which defined as A(t) would take values a(t) at t with probability PA (a,t | a0,0). Adding this new state, the system would result into the following set of time-independent SDEs:
(5)
(6) where v and ε are the corresponding parameters associated with the stochastic adjustment function.
In a quite similar way, other stochastic variables (states) describing the different sources of biological variability, such as the critical size to initiate division, fluctuations in the initial size of the daughter cells, or cell-to-cell interactions, can be included in the description. Such extensions, however, are not considered in the present work, as they involve extra parameters which call for additional experimental information, difficult to collect or unavailable in the existing literature.
In this paper, cell growth and division are modeled by equations 3 and 4. Note, however, that under appropriate experimental design conditions, the method we propose in this work can be extended in a straightforward manner to estimate those parameters.
Model calibration.Single-cell parameters μ and ξ in equation 3 and the initial value of the adjustment function a0 in equation 4 can be estimated from experimental time-to-division (TTD) distributions [which we refer to hereafter as Λ̂(T)] as those reported in references 4 and 5), for instance, by a least-squares method. Essentially it aims at the minimization of the differences between Λ̂(T) and a theoretical TTD Λ(T;θ), defined as a function of the parameters. Formally the problem can be stated as follows:
(7) where J(θ) is the objective function to be minimized (i.e., the integral over time to division of the square errors between data and model) and θ represents the set of parameters to be estimated.
Assuming that growth and division accept a description based on equations 3 and 4, the theoretical TTD distribution can be computed by the formula
(8) where P(d,T | y0,0) corresponds to the probability of the cell length reaching for the first time a value exp(d) at t equal to T, given the cell initial size being x0 [so that y0 is equal to ln(x0)]. Such probability is obtained from the solution of a partial differential equation (PDE) (equation A-16 in the supplemental material), known in stochastic calculus as the backwards Kolmogorov equation, with appropriate boundary and final conditions for every T in the interval (0,∞).
In this work, the partial differential equation A-16 with the corresponding boundary conditions has been solved with a finite difference discretization scheme (see http://www.matmol.org/) to approximate the original partial differential equation (PDE) (27). For all cases, a mesh consisting of 501 elements was considered enough to accurately approximate the equation. Time integration of the resulting set of ordinary differential equations has been performed in Matlab with a standard ODE solver (ode15s).
During first division, cells undergo adaptation, so a0 is a positive number (smaller than 1) to be estimated. After the second, third, and fourth division times, cells are assumed to be adapted to the environment so that a0 is equal to 1 and the estimation reduces to the computation of the growth rate (μ) and the intensity of the stochastic fluctuation (ξ).
Optimizers fminsearch and fmincon from the Matlab optimization toolbox have been employed to solve the least-squares minimization problem (equation 7), leading both methods to the same results for the cases considered.
In order to test the performance of the proposed model to reproduce time to division distributions, two sources of experimental data have been employed. One is taken from reference 4 and corresponds to the distributions of times to first division for Listeria innocua under different heat shock durations ranging from no shock to a 5-min shock duration. The other source comes from reference 5 and corresponds to the distributions of times to first division and successive division times (up to the fourth division) for Escherichia coli at 25 and 32°C.
Simulation of bacterial population growth.Population growth from a given number of colony formation individuals over a given time horizon has been simulated by assigning to each bacterium (including its corresponding offspring) equations 3 and 4 to be solved from its initial size to its size of division. The process is then repeated for each new offspring over the time horizon. A number of numerical solution methods for solving the stochastic differential equations are available (see reference 28). In this work, the Euler-Maruyama algorithm has been selected for its simplicity. For convenience, simulation of population dynamics from the proposed single-cell stochastic model has been performed on a cluster composed of 12 processing nodes (openSUSE 11.0 Linux with 23.5 GB of RAM) and 160 processors in total, using the SGE task manager to distribute the calculations between them.
RESULTS AND DISCUSSION
Theoretical TTD distributions versus SDE realizations.A computational experiment has been performed to show that equation 8 provides a precise representation of the TTD distributions obtained by a number of individual bacteria growing according to equations 3 and 4 up to a critical length or size. Theoretical arguments that support such equivalence can be found in the framework of Ito calculus. For the interested reader, these are summarized in the supplemental material.
Cells with initial length (x0) of 8 were assumed to divide at double their length, i.e., to attain a division length (xdiv) of 16. The time to division of each realization (cell) is calculated as the time the cell, which grows according to equation 3, first reaches the xdiv. To simulate cell growth, equation 3 is solved numerically with the Euler-Maruyama method (28) until the first time to division, namely, the time variable Y attains a value (d) that is identical to ln(xdiv). The corresponding distribution of times to division (TTD) is constructed by repeating the simulation for a sufficient number of cells (each cell constitutes a realization) in order to produce a representative ensemble. In the present study, ensembles comprised around 10,000 realizations.
Parameters employed in the simulation (μ = 0.6060 h−1 and ξ = 0.0821) are on the order of those obtained from data taken from the literature for E. coli (see next subsection and Table 2 for further details). In order to evaluate the effect of the adjustment function on the resulting TTD, two values were considered: an a0 of 0.1, which would correspond to a cell on its first division time, and an a0 of 1 for a cell completely adapted that is growing at its maximum specific rate.
Figure 1a and b present the resulting TTD distributions obtained from the SDE (bars) and from the proposed theoretical distribution 8 (lines) for a cell population undergoing adaptation and completely adapted, respectively. Parity plots in Fig. 1c and d comparing cumulative distributions obtained from the SDE (QSDE) and the proposed theoretical distribution (QΛ) prove that a perfect match exists between the two approaches. This is in agreement with the Kolmogorov-Smirnov tests applied with a 0.05 significance level to check coincidence between the samples of times to division and the cumulative distribution associated with Λ (equation 8).
Distribution of times to division for an ensemble of 10,000 realizations obtained from the SDE (equations 3 and 4) (bars) and from equation 8 (lines). Model parameters are μ of 0.6060 h−1 and ξ of 0.0821. Panel a represents distributions for cells under adaptation (a0 = 0.1000). Panel b represents cells completely adapted (a0 = 1.0000). The corresponding parity plots comparing the cumulative distributions computed from the SDE (QSDE) and the Λ distributions (equation 8) (QΛ) are presented in panels c and d for both scenarios.
It must be noted here that apart from the fact that the backward Kolmogorov PDE (equation A-16 in the supplemental material) offers direct information of the statistical properties of the systems, its solution is by far much more efficient from a computational point of view than is computing a significant number of realizations by means of the SDE, which usually demands a sufficiently populated ensemble to be constructed beforehand in order to be representative. For the scenarios considered in Fig. 1, the ensemble must be on the order of thousands of realizations. As an indication, solving the PDE in Matlab on a standard personal computer takes on the order of seconds and never more than a few minutes. On the other hand, computing just one representative ensemble for any of the two cases considered requires around 1 h of computing time. For this reason, the use of the backward Kolmogorov equation is preferred to Monte Carlo methods for stochastic model calibration purposes, which usually require repeated evaluations of objective function 7.
Model performance to describe experimental TTD.We investigated the capability of our stochastic model of single-cell growth to describe the distribution of times to division observed in experiments. The proposed model consists of SDE 3 together with adjustment function 4. Model parameters include the growth rates, (μ and ν), the intensity of the stochastic fluctuation (ξ), the initial length (x0 [thus b = ln x0]), the length at division (xdiv) (so that d = ln xdiv), and the initial value of the adjustment function (a0). Following reference 4, it will be assumed that ν is equal to μ. As we have shown above for the computational experiment, the TTD distribution produced by our model can be precisely computed by means of equation 8.
Based on cell length data reported in reference 4, an average value for cell length at an xdiv of 16 (in the units of pixels as the authors report) was selected in all simulations. Division is assumed to result in two daughter cells of similar length, so that x0 is equal to 8 pixels. Cell viability is taken into account in the model by setting a minimum cell length below which it is assumed that the cell will die.
In order to compute the theoretical distribution of times to division (equation 8), we set y0 as equal to ln x0 and solve the PDE A-16 discussed in the supplemental material for different time horizons (T). In this study, a minimum length of 4 pixels was selected, so that the domain for variable x lay in the interval between 4 and 16 pixels. Hence, boundary conditions A-14 with c being equal to ln(4) and d being equal to ln(16) must be imposed (see the supplemental material).
A comparison of model predictions with experimental data is presented in Fig. 2 for the first division times of Listeria innocua under no heat shock (Fig. 2a) and for a 5-min shock duration (Fig. 2b) (4). The comparison of the predictions of our model and the corresponding fittings provided by gamma functions for Escherichia coli data at different division times and temperatures (5) is presented in Fig. 3 and 4. Table 1 summarizes the gamma function parameters given by reference 5, in the form of mean time and standard deviation of the distributions.
Comparison between the experimental distribution (bars) of times to first division for Listeria innocua (4) and the theoretical distribution of times to division (equation 8) (lines). Estimated model parameters are presented in Table 2. (a) No heat shock; (b) 5-min heat shock.
TTD distributions (1st to 4th times to division) for Escherichia coli at 25°C (5). (a) 1st division; (b) 2nd division; (c) 3rd division; (d) 4th division. In each plot, experimental distributions (bars) and fittings to gamma distributions (dashed lines) reported in reference 5 are compared with the corresponding model estimations (solid lines).
TTD distributions (1st to 3rd times to division) for Escherichia coli at 32°C (5). (a) 1st division; (b) 2nd division; (c) 3rd division. In each plot, experimental distributions (bars) and fittings to gamma distributions (dashed lines) reported in reference 5 are compared with the corresponding model estimations (solid lines).
Gamma function parameters given by reference 5a
Plots in Fig. 3 and 4 show the experimental data (bars) together with the corresponding fittings to gamma functions provided in reference 5 (dashed lines), as well as model estimations (solid lines). Parameter values obtained from model calibration are presented in Table 2, including the resulting summation of square errors (SSqED—model), which corresponds with the final value attained by the objective function (denoted as J in equation 7). A quantitative measure of the agreement between gamma distributions and the experimental data is given in the last column of Table 2, through the corresponding summations of square errors (SSqED—gamma).
As can be seen in the figures, there is quite good agreement (corroborated by Kolmogorov-Smirnov tests) between our model and the experimental data for all division times at both temperatures. In fact, our model is capable of fitting experimental data for all divisions at 25 and 32°C, better than gamma distributions: as shown in Table 2, the values of the summation of square errors (SSqED—model) for 2nd to 4th division times are 1 order of magnitude smaller than those corresponding to the gamma distributions (SSqED—gamma). For the 1st division time, model fittings to the corresponding experimental data are only slightly better than those provided by gamma distributions (values for SSqED—model and SSqED—gamma are of the same order of magnitude), probably due to larger errors or uncertainties in the experimental measurements.
A sensitivity analysis test was performed on the computed parameters to calculate how their changes affect the objective function. A typical shape of the objective function J for different parameter sets is presented in Fig. 5, showing in all cases a clear minimum. This implies that just one optimal set of parameters complies with a given TTD curve, which suggests that the model parameters are identifiable.
Objective function J in expression 7 is plotted on a three-dimensional plot for different values of μ and ξ to show the presence of a minimum. Experimental TTD data employed to compute the objective function correspond with the second division times at 32°C for Escherichia coli (5).
Finally, it is worth noting how similar the values are for the growth rates for the second to fourth times to division (μ ≈ 0.6 to 0.7 h−1 at 25°C or μ ≈ 0.8 to 1.0 h−1 at 32°C) to the values corresponding to the first time to division (μ ≈ 0.24 h−1 at 25°C or μ ≈ 0.60 to 0.70 h−1 at 32°C). The intensities of the stochastic fluctuations, on the other hand, remain quite uniform in the range of 0.1 to 0.2 for ξ at 25°C and 0.2 to 0.3 for ξ at 32°C for 1st to 4th times to division.
Stochastic modeling of population dynamics.To demonstrate the capability of the proposed model to describe population growth, we followed reference 5 and made use of the parameters obtained from the TTD at 32°C for the first, second, and third divisions to simulate growth curves for Escherichia coli. As in reference 5, the results of the simulation are compared with experimental data measured by viable counts and different inoculum sizes, showing excellent agreement, as can be seen in Fig. 6. The proposed stochastic model can therefore be used to provide an effective link between single cell growth and cell population dynamics.
Growth curves at 32°C for Escherichia coli measured by viable counts (dots) for different inoculum size versus those generated by simulation (lines). Experimental data are taken from reference 5.
It can be argued that the model described by equation 3 is not particularly superior to other constructions fitted to the experimental data, although it has been shown already to produce a more consistent fit (see Fig. 3 and 4). In this regard, it must be said that the proposed model has a clear biological interpretation, whereas a gamma function approximation is a mere empirical fitting. The model makes use of a critical variable which might be related to length or size (e.g., via the total amount of DNA) of single cells, evolves during their cell cycle, and determines division in a way that incorporates cell variability within a bacterial population.
Moreover, the use of a differential formulation instead of a particular solution makes the model predictive even under changes in environmental variables (such as temperature or pH, for instance), provided that experimental data are available to calibrate secondary models that relate stress variables to parameters of the single-cell growth model. From this point of view, the approach presented can be particularly useful to predict population growth under variable environmental conditions.
Finally, it must be remarked that although alternative stochastic methods, and particularly the one proposed in reference 5 to simulate population growth from a given TTD, constitute legitimate approaches, the computational efficiency is inferior to the one based on the use of equation 3. The direct use of the TTD distribution (as in reference 5) requires one to obtain and keep track of the set of random division times, which becomes quite inefficient as the size of the population increases. Handling division by means of stochastic differential equations as in our approach becomes more efficient both from a mathematical and from a computational point of view.
Effect of initial number of cells.In order to test the model's ability to reproduce the behavior for small populations and to predict population growth as a function of the initial population size, we made use of the experimental data reported in reference 2 for Salmonella enterica at 25°C. Data were obtained by the authors using time-lapse microscopy and comprise time-to-division distributions for first, second, and third division times (presented in the plots in Fig. 7 as bars) as well as plots of colonial growth evolution starting from single individuals.
Comparison between experimental TTD distributions for Salmonella enterica taken from reference 2 (bars) and the model for the 1st, 2nd, and 3rd division times (panels a, b, and c, respectively).
Stochastic model 3 was calibrated according to the procedure described in Materials and Methods. Model parameters for the different division times are presented in Table 3. Comparison between experimental and estimated TTD distributions are plotted in Fig. 7, showing a quite reasonable agreement, especially when one takes into account the relatively scarce experimental information available.
Model parameters estimated for the data provided by reference 2 for Salmonella enterica at 25°C
Figure 8 presents a few simulations of population growth curves initiated from one cell. As it can be seen in the figure, the resemblance with the experimental observations reported in reference 2 (Fig. 4 in that article) is remarkable.
Detailed representation of a few realizations for population growth starting from one single bacterium.
Our model is also able to explain the effect of the initial number of cells on the variability of population growth. Figure 9 presents simulations starting from 1 to 100 individuals, showing that variability in population growth is reduced as the number of cells initiating the population increases.
Ten thousand simulated realizations for population growth curves for Salmonella enterica starting from different initial cell numbers: 1 single bacterium (a), 2 single bacteria (b), 10 individuals (c), and 100 individuals (d).
From this point of view, the proposed model is in agreement with deterministic growth descriptions involving large population sizes and with the fact that random effects become negligible as the size of the population increases. Note, however, that this is not the case when models with random parameters are employed instead to simulate population growth curves.
To illustrate this fact, a model similar to that proposed in reference 12, with lag time and specific growth rate parameters following random distributions, has been employed (as suggested in reference 2) to simulate population growth rates. Simulations are depicted in Fig. 10 for different numbers of initial cells using the random distributions for model parameters reported in reference 2. This model predicts larger variabilities as the initial number of cells decreases (Fig. 10), due to the random nature of its parameters; however, given sufficient time, variability spreads out no matter the number of initial cells, which is in contradiction to the observed experimental behavior of large microbial populations (variability becomes negligible through the law of large numbers). In contrast, the model we propose achieves a constant (stationary) variability range which reduces with population size and becomes negligible for populations initiated with large number of individuals, which is in agreement with classical bacterial growth kinetics (in this example, for numbers above 100 individuals, as can be seen in Fig. 9d). This effect can also be seen in Fig. 11a and b, which depict probability distribution functions at different times for populations starting from 1 and 100 individuals, respectively. As observed in Fig. 11a and b, distribution width shrinks as the initial number of individuals increases. Note that in both situations, distributions evolve to a constant shape function which travels right in the x axis (representing log of cell numbers).
Ten thousand simulated realizations for population growth curves for Salmonella enterica using an ordinary differential equation with random parameters for different initial cell numbers. Parameter distribution is taken from reference 2. (a) One single bacterium; (b) 2 single bacteria; (c) 10 individuals; (d) 100 individuals.
Probability and cumulative probability distribution functions for Salmonella enterica populations computed with the proposed model (parameters are in Table 3) starting from one individual (a and c) and from 100 individuals (b and d). Panels c and d depict the cumulative probability for the population to attain 1,000 individuals as a function of time for 1 and 100 bacteria, respectively.
Finally, it must be remarked that our model allows the computation of probability (and cumulative probability) distributions associated to the population as a function of time in a straightforward manner (Fig. 11). As an example, let us suppose that for the case considered (Salmonella enterica at 25°C), 1,000 individuals define a risk threshold. Figure 11c and d indicate that the probability of a bacterial population being larger than 1,000 individuals is negligible before 10 days if growth is initiated with 1 cell, while such time is reduced to more than half (less than 5 days) if the initial number of cells is 10. From this perspective, the plots of the resulting cumulative probability for scenarios involving different initial bacteria can give indications on expected shelf life, which may be of help for risk assessment and shelf life evaluations.
Conclusions.A stochastic version of an individual bacterial growth model has been proposed to describe cell heterogeneity for a given population and to predict time-to-division distributions. Distributions are reconstructed by solving the backward Kolmogorov equation, which is a partial differential equation establishing the functional relationship of probability density with bacterial length (size) and time.
The method is computationally efficient since it only requires the solution of a partial differential equation (which takes only a few seconds on a standard personal computer) to reconstruct the density function of the distribution and thus the distribution of times to division while still maintaining information about the stochastic nature of the process. Evidence of the model's ability to cope with the observed behavior is given for the distributions of division times of Escherichia coli and Listeria innocua reported in the literature.
In addition, the capability of the model to describe population growth has been demonstrated, explaining population variability increase as the initial number of cells is reduced, thus providing a bridge that links individual bacterial growth and the growth of bacterial populations and opens the door for applications in risk assessment and prediction of shelf life.
ACKNOWLEDGMENTS
This work has been funded by the Spanish Ministry of Science and Innovation through mobility grant Salvador de Madariaga (PR2011-0363) and by project ISFORQUALITY (AGL2012-39951-C02-01).
FOOTNOTES
- Received 6 May 2014.
- Accepted 8 June 2014.
- Accepted manuscript posted online 13 June 2014.
- Address correspondence to Antonio A. Alonso, antonio{at}iim.csic.es.
Supplemental material for this article may be found at http://dx.doi.org/10.1128/AEM.01423-14.
REFERENCES
- Copyright © 2014, American Society for Microbiology. All Rights Reserved.
