Stochastic Individual-Based Modeling of Bacterial Growth and Division Using Flow Cytometry.

A realistic description of the variability in bacterial growth and division is critical to produce reliable predictions of safety risks along the food chain. Individual-based modeling of bacteria provides the theoretical framework to deal with this variability, but it requires information about the individual behavior of bacteria inside populations. In this work, we overcome this problem by estimating the individual behavior of bacteria from population statistics obtained with flow cytometry. For this objective, a stochastic individual-based modeling framework is defined based on standard assumptions during division and exponential growth. The unknown single-cell parameters required for running the individual-based modeling simulations, such as cell size growth rate, are estimated from the flow cytometry data. Instead of using directly the individual-based model, we make use of a modified Fokker-Plank equation. This only equation simulates the population statistics in function of the unknown single-cell parameters. We test the validity of the approach by modeling the growth and division of Pediococcus acidilactici within the exponential phase. Estimations reveal the statistics of cell growth and division using only data from flow cytometry at a given time. From the relationship between the mother and daughter volumes, we also predict that P. acidilactici divide into two successive parallel planes.

1. Introduction

Population- and individual-based modeling are usually presented as incompatible approaches, although both describe the same system at different levels (Fahse et al., 1998; Wilson, 1998).

Traditionally, deterministic population-based models have been the underlying method behind predictive microbiology (Baranyi and Roberts, 1995). These models have been successfully applied to, for example, monitoring of food spoilage and microbial safety (Koutsoumanis and Nychas, 2000; Ross et al., 2000), smart sensing of food quality (García et al., 2015, 2017), Quantitative Microbial Risk Assessment (Cassin et al., 1998; Membré and Lambert, 2008), and design and control of food processes (Simpson et al., 1993; Alonso et al., 2013).

Over the last 15 years, stochastic individual-based modeling emerged as a promising tool to produce realistic estimations of safety risks along the food chain by describing the variability of single-cell behavior and small populations (Ferrer et al., 2009; Augustin et al., 2015; Koutsoumanis and Aspridou, 2017). Often, contamination of food starts with a small number of bacteria that adapt and proliferate on a given food matrix. At low cell concentrations, standard deterministic population 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 cells cannot be neglected when assessing possible health risks along the food chain, either during storage or during distribution.

There are still many challenges in individual-based modeling, including the lack of information about single-cell behavior inside a population. The emergence of individual-based modeling was possible thanks to two main factors: (1) the increase of computer processing power and (2) the availability of single-cell measurements using new techniques such as the “mother machine” microfluidic device (Wang et al., 2010). However, the information from single-cell measurements is limited in those techniques where cells have to be isolated from the population. That was illustrated for example by Gangan and Athale (2017), who showed the difference in single-cell growth in a “mother machine” or in a population.

In this work, we hypothesize that population statistics of cell volume encode information about single-cell growth and division that can be used for individual-based modeling. For this purpose, we derive a modified Fokker-Planck Equation (forward Kolmogorov equation) describing the population volume distribution. The underlying idea is similar to that in Alonso et al. (2014) who derived a backward Kolmogorov equation to estimate single-cell growth using time-to-division distributions. We should remark that not only our approach is different, but it extends the previous theory to consider not only single-cell growth but also single-cell division. This allows us to simulate the individual-based modeling of bacterial growth and division.

In the first part of the work, we will test theoretically how population statistics obey a modified Fokker-Planck equation that encodes single-cell information. The equivalent individual-based modeling approach is derived in parallel to check consistency. Both models simulate single-cell growth and division assuming that cell volumes grow exponentially and cells divide following the sizer principle, i.e., division occurs at a critical volume (Métris et al., 2005; Alonso et al., 2014; Robert et al., 2014). Whereas, exponential growth of cell volume is a standard principle in bacterial physiology (Fishov et al., 1995), the main trigger of bacterial division is still a matter of controversy (Taheri-Araghi et al., 2015). There are three major paradigms: the sizer, timer, and adder principle depending on whether division is triggered by a certain volume, time, or after growing a given volume. As in most predictive microbiology studies (Métris et al., 2005; Alonso et al., 2014), we focus on fully adapted cells (medium growth is kept constant and measurements are within the exponential phase) and the sizer principle remains the reasonable assumption.

Once the theory is established, we combine the modified Fokker-Planck equation with flow cytometry data to find the single-cell behavior of Pediococcus acidilactici within the exponential phase. The food industry is interested in this species for several reasons, including its probiotic attributes (Planas et al., 2004; Standen et al., 2015), its ability to valorize food wastes (Vázquez et al., 2011; Banwo et al., 2013; Scatassa et al., 2015), and its ability to produce a very potent and broad-spectrum bacteriocin (pediocin SA-1) with high capacities as food biopreservative (Ray, 1992; Anastasiadou et al., 2008; Vázquez and Murado, 2008). P. acidilactici has been also selected for being coccoid cells with interest in the food industry. This shape makes easier to find correlations between side scatter and volume, and differ from the well-studied Escherichia coli (model of rod-shaped cells). We should stress that cell volume, membrane area and diameter scale similarly when the cell is rod-shaped, although that is not the case for round cells. For such reason, along this work we consider the term size as equivalent to volume, but not bacterial diameter or membrane area.

2. Materials and methods

This work combines theory with experimental data and requires three types of methodologies to (1) develop models at the single-cell and population level, (2) acquire data with flow cytometry and optical density, and (3) determine the best parameters to reconcile the theory with the experimental results.

2.1. Modeling at single-cell and population levels

2.1.1. Individual-based modeling of single-cell bacterial growth and division

We tested different alternatives with stochastic or deterministic division and growth, that are specific cases of the general individual-based modeling approach we describe in this section.

The model assumes that the growth of the logarithm of the single-cell volume is subject to a stochastic fluctuation δW characterized by a Wiener process (Alonso et al., 2014):

where X represents the volume of cell i (V) in a logarithmic scale, μ represents the growth rate within the exponential phase and ξ is the intensity of the stochastic fluctuation. For the case of deterministic division ξ = 0 and cell volumes grow exponentially.

The division was modeled by adding a new cell to the population and resizing mother and newborn cell to the daughter size. The division event is triggered when the size of one or more cells is greater than a continuous random variable X with statistics defined by the probability density function (pdf) of mother sizes (f(x)):

where n is the number of cells in the population, i runs from 1 to n and the daughter volume is half the mother volume (υ = υ/2).

We tested different probability density functions to describe the statistics of cell division, i.e., for the probability density function of mother sizes f(x). The probability showing the best agreement with the data suggests that the volume of the mothers V is a random variable following a log-normal distribution (Koch, 1966; Amir, 2014):

For simulating the deterministic division, σ was set to zero so that the normal distribution turns into Dirac delta function centered at the logarithm of the mother size x.

Simulations were initialized with a single cell () and run for a given time horizon where a given cell and its offspring grow following (1a) and divide according to the rule in (1b). We selected the Euler-Maruyama algorithm to solve the stochastic differential equations for its simplicity as compared with other numerical methods (Higham, 2001). The bins of all the histograms to represent the population statistics were determined using the Freedman-Diaconis rule (Freedman and Diaconis, 1981). For convenience, simulations of population dynamics from the proposed single-cell stochastic model were 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 among them.

2.1.2. Population modeling using the modified fokker-planck equation

The statistic of the cell sizes in the population is formally described by a probability density function (pdf) p(t, x) that depends on size x and time t. For the sake of clarity, we keep previous subsection notation: υ and x denote size in terms of volume and natural logarithm of the volume, respectively. Subindexes d and m denote daughter and mother respectively, whereas f and f are the corresponding pdfs for sizes. As before, growth rate and fluctuation intensity are denoted by μ and ξ, respectively. The function p(t, x) is the solution of the following modified Fokker-Planck equation:

being

where is the pdf of daughter sizes, the pdf of mother sizes and F its cumulative distribution function. The model is valid only for large domains where no cell sizes are close to their (minimum and maximum) boundaries.

Without the terms of division and normalization, Equation (2a) is the classical Fokker-Planck equation of the stochastic differential Equation (1a). It explains how the change in the distribution of volumes depends on a diffusion term which is proportional to the square of the fluctuation plus a convective term proportional to the cell growth rate (Gardiner, 2004; Alonso et al., 2014). To account for division, we have added two terms proportional to the pdfs of daughter and mother volumes. Normalization is required for p(t, x) to be a pdf. This is performed via the last term in the right hand side.

Simulations were performed using the finite difference discretization scheme in http://www.matmol.org/ (Vande Wouwer et al., 2014) for the x domain. For all cases, the discretization scheme consisted of 501 elements. That was considered enough to approximate the equation since further refinements resulted in negligible improvements in the accuracy of the results. First derivatives were calculated using an upwind 5 points in the stencil and second derivatives with centred 5 points in the stencil. Due to the hard non-linearity at division, not only is a refined mesh in x required, but also a stiff time integrator. Ode15s in Matlab (Shampine and Reichelt, 1997) was selected for time integration of the resulted set of ordinary differential equations after the spatial discretization.

2.2. Data acquisition and analysis

2.2.1. Microbiological methods

Pediococcus acidilactici NRRL B-5627 was kindly provided by the Northern Regional Research Laboratory (Peoria, IL, USA). Stock cultures of bacterium were stored at -80 °C in MRS commercial medium (Pronadisa, Hispanlab S.A., Spain) with 25% glycerol. The inoculum to study the growth dynamics of P. acidilactici was prepared as follows:

One hundred and fifty milliliters of cellular suspension from the cryotube was transferred to 5 mL of MRS fresh medium and then incubated at 30 °C in an orbital shaker at 200 rpm for 16 h.

From the obtained culture, 1 mL was added to an Erlenmeyer flask with 150 mL of MRS fresh medium and fermented at 30 °C/200 rpm for 22 h.

From the previous cultivation, serial 10-fold dilutions were prepared in peptone-buffered solutions, and 0.1 mL samples were plated (MRS agar medium) in triplicate and incubated at 30 °C for 48 h.

Five individual colonies from plates were isolated and transferred to 5 Erlenmeyer flasks with 200 mL of MRS fresh medium and cultivated at 30 °C/200 rpm. Samples from flasks were taken each hour up to 17 h (except at 12 and 16 h). All samples were separated in two aliquots, one of them was prepared for cytometer evaluation following the indications described in the next section. The other aliquot was centrifuged at 4,000 g for 15 min and the sediment washed twice and re-suspended in distilled water at an appropriate dilution to measure the optical density at 700 nm. The dry weight was estimated from a calibration curve ().

The percentage of viable cells smaller than certain diameters was calculated during the exponential phase. The cultures at 8 h were filtered, under sterile conditions, through 1.2 and 1 μm glass microfiber filters (Filter-Lab, Filtros Anoia S.A., Barcelona, Spain). Thus, in these filtered solutions and in the final unfiltered culture (control), viable cells (colony forming units per mL) were quantified by count on MRS-agar plate as it was mentioned in the previous paragraph.

2.2.2. Flow cytometer data acquisition

The abundance and size of P. acidilactici were determined with a BD FACSCalibur flow cytometer (BD Biosciences, San José, CA, USA) equipped with a laser emitting at 488 nm. Bacteria samples were fixed with a P+G solution (1 % paraformaldehyde + 0.05 % glutaraldehyde) at 10 % final concentration for 15 min in the dark. Then, samples were quickly frozen in liquid nitrogen and stored at -80 °C. Prior to analysis, bacteria were stained with SybrGreen I DNA dye (5 mM final concentration) and diluted adequately. Bacteria were detected in the flow cytometer by their signature in a plot of Side Scatter (SSC) vs. FL1 (green fluorescence). All the reagents and chemicals were purchased from Sigma-Aldrich S.A. (St. Louis, MO, USA).

2.2.3. Estimation of P. acidilactici population growth

The logistic equation was used to fit the population growth data (Zwietering et al., 1990; Peleg and Shetty, 1997). Biomass and number of cells were obtained respectively by dry weight and cytometry of P. acidilactici:

where G is the P. acidilactici growth as biomass or cells (g L−1 or cells mL−1). G represents the maximum growth or plateau phase (g L−1 or cells mL−1), λ is the lag phase (h), t denotes the time of culture (h) and is the specific growth rate of the population (h−1) with being the maximum specific growth rate (g L−1 h−1 or cells mL−1 h−1).

Non-linear least-squares method (quasi-Newton) was applied for growth data modeling. Confidence intervals from the parametric estimates (Student's t-test) and consistency of mathematical models (Fisher's F-test) and residual analysis (Durbin-Watson test) were evaluated by “SolverAid” macro (Levie's Excellaneous website: http://www.bowdoin.edu/~rdelevie/excellaneous)

2.3. Model calibration

Reliable single-cell parameters of P. acidilactici are unknown and were estimated by minimizing the distance between the modified Fokker-Planck Equations (2a–2d) and the data from the flow cytometry. The unknown parameters are the growth rate μ, the fluctuation intensity ξ, the statistics of the mother distribution (x and σ) and the parameter relating mother and daughter sizes ω. The estimated parameters were used to simulate single-cell dynamics based on the individual-based modeling approach (1a–1d).

The method of least squares was used to define the distance between the model and the data. Essentially, it aims at minimizing the differences between the stationary distribution of sizes calculated with (2a–2d) and the stationary distribution estimated from the data. The experimental distribution was obtained using different replicates at one given time using histograms with a number of bins given by the Freedman-Diaconis rule. For noisy data, optimization could lead to a multimodal problem, i.e., it has several sub-optimal solutions (Vilas et al., 2017). To assure convergence to the global solution in a reasonable time, the global optimizer Enhanced Scatter Search (eSS) was employed (Egea et al., 2009).

3. Results and discussion

3.1. Population statistics of single-cell growth and division obey a modified fokker-planck equation

3.1.1. Deterministic growth and division

We first simulate the individual-based modeling of deterministic growth and division. We assume that cell volumes grow exponentially and cells divide following the sizer principle, i.e., division occurs at a critical volume (Métris et al., 2005; Alonso et al., 2014; Robert et al., 2014). The adder model is more realistic while cells are adapting to the growth media (lag phase), but it becomes a sizer when, as in our case, cells are fully adapted within the exponential phase (see Figure 1B in Sauls et al., 2016). It should be noted that considering an adder model would complicate considerably the derivation of the equivalent Fokker-Planck equation without altering the final results.

Figure 1A shows the simulations of the deterministic single-cell dynamics. All the cells have the same volume because they grow and divide at the same velocity and time. For this example only 3 parameters are required: the daughter volume (υ) and the mother volume (υ = 4), depicted in dashed blue and red lines, and the single-cell growth rate μ = 0.7. At time 0 there is only one cell that grows until reaching the critical volume of division (or mother volume). This cell divides into two cells of half their volumes (daughter volume). Therefore, at time 1 there are two cells that cannot be distinguished in the figure because their dynamics overlap. The process is repeated until reaching a population of 32 cells at time 5.

Figure 1

Population statistics of bacterial growth obey a modified Fokker-Planck equation from which we can extract information about single-cell behavior. (A) Simulation of the deterministic individual-based modeling approach shows that cells grow exponentially at the same rate and divide at the same time into two cells with equal volume when reaching a critical size (the mother volume). (B) The statistics of the population volume obey a modified Fokker-Planck equation that oscillates between mother and daughter volumes without reaching a stationary distribution. (C) Individual-based modeling with stochastic growth assumes that the logarithm of the volume is subject to a stochastic fluctuation δW characterized by a Wiener process. (D) Simulation of the equivalent modified Fokker-Plank now shows that the population volumes evolve to a stationary distribution. (E) Individual-based modeling can be used to estimate the histogram of the stationary distribution, but at the expenses of expensive computations that scale exponentially with time and linearly with the number of cells in the population. The distribution is sharp and skewed to the right and encodes single-cell features such as the mother and daughter volumes, fluctuation, and growth ratios. (F) The modified Fokker-Planck simulates the continuous shape of the stationary distribution (red line) in a efficient way that is independent on the number of cells within the population. (G) The individual-based modeling simulates stochastic growth and division. (H) The resulting stationary distribution of the population volumes is smooth and equivalent when calculated using the individual-based modeling and the modified Fokker-Planck equation.

The dynamics of the population volume distribution obey a modified Fokker-Planck equation that, as shown in Figure 1B, consists of a pulse that oscillates between the daughter and mother volumes. Each color line represents the distribution at a different time and are simulated using the partial derivative Equation (2a) with stochastic parameters set to zero (ξ = 0 and σ = 0) and the single-cell parameters in Figure 1A (υ = 2, υ = 4, μ = 0.7).

This model with deterministic growth and division is invalid since the population statistics fluctuates instead of evolving to a stationary distribution. Fishov et al. (1995) explained how balanced exponential growth implies steady-state growth and stationary frequency distribution of the various components that constitute the cell. In words by Painter and Marr (1968) “the distribution of each intensive random variable is time-invariant.”

3.1.2. Stochastic growth and deterministic division

Single-cell measurements of bacteria suggest that growth is a stochastic process (see for example Figure 1A in Deforet et al., 2015 for Pseudomonas aeruginosa). Alonso et al. (2014) proposed an individual-based modeling approach reproducing such single-cell dynamics. They assumed that the logarithm of the volume is subject to a stochastic fluctuation δW characterized by a Wiener process following Equation (1a). We use the same assumption to model stochastic growth. Figure 1C shows simulations of the single-cell dynamics with stochastic growth and deterministic division. After the first division (time = 0.9) the dynamics of the two daughter volumes differ and the same happens with division times.

We confirm how simulations of the population statistics with stochastic growth now evolve to a stationary distribution (Figure 1D) as predicted by (Fishov et al., 1995). At time 0 the distribution is a Dirac delta function that spreads and moves between mother and daughter volumes. The rate of spread depends on the intensity of the fluctuation ξ whereas the velocity of the moving pulse is determined by the growth parameter μ. Single-cell parameters are as in previous section (υ = 2, υ = 4, μ = 0.7) except for the fluctuation that it is now ξ = 0.1.

The stationary distribution can be calculated either using a large number of single-cell simulations (1a–1d) or solving the modified Fokker-Planck equation (2a–2d). We plot in Figure 1E the histogram for a population of 1e6 cells using the individual-based modeling approach. The mode of the histogram coincides with the daughter volume whereas the end of the histogram is the mother volume. The shape depends on the growth rate μ and the fluctuation ξ. Normalizing the area of this histogram we obtain the blue probability density function in Figure 1F. The modified Fokker-Planck equation, on the other hand, calculates directly the probability density function (red line). Both approaches coincide as shown in Figure 1F.

Results with stochastic growth and deterministic division evolve to a stationary distribution, but sharper than observed experimentally. In fact, the distribution of E. coli is commonly approximated for some authors by a smooth log-normal distribution (Kaya and Koser, 2009; Athale and Chaudhari, 2011). In addition, the end of the distribution (the mother volume) is not exactly the double of the mode (the daughter volume) in experiments.

3.1.3. Stochastic growth and division

A great amount of works in the literature assumes stochastic division and focus on the distribution of mother or daughter sizes (Amir, 2014; Taheri-Araghi, 2015; Taheri-Araghi et al., 2015; Sauls et al., 2016). Some works measure symmetric distributions close to a normal distribution (Sauls et al., 2016), whereas others assumed asymmetric distributions. That is the case of Koch (1966) and Amir (2014) who concluded that the daughter volume distribution is log-normal.

We extended our model considering that cells divide stochastically with a certain probability, either normal (Sauls et al., 2016) or log-normal (Amir, 2014). In other words, we moved from a strict sizer model to a sizer model with stochastic division. For the simulations in Figure 1G we assumed log-normal division. Now cells may divide before or after reaching the volume of 4 (red dashed line). As the model was implemented in the logarithm of the volume X, the log-normal stochastic division in the logarithm becomes a normal probability with mean υ = 4 and a standard deviation that we assumed to be σ = 0.1. The remaining parameters are kept as in the previous section (υ = 2, v = 4, μ = 0.7, ξ = 0.1).

As shown in Figure 1H, results from the individual-based modeling coincide with the modified Fokker-Planck also for stochastic growth and division.

We should note how it is critical to assume stochastic growth and division to obtain realistic and smooth stationary distributions where the mode is larger than the daughter volume (υ).

3.1.4. Comparison of individual-based modeling and population modeling with the modified fokker-planck equation

Individual-based modeling is a bottom-up approach providing valuable information at the single-cell level, but it requires parameters that cannot be easily measured (Ferrer et al., 2009; Augustin et al., 2015). The modified Fokker-Planck equation here presented focuses on population statistics that can be measured by flow cytometry. Comparisons between this equation and the experimental data are sufficient to estimate single-cell parameters that can be used for individual-based modeling.

The modified Fokker-Planck equation directly provides the evolution of the volume distribution without the need of predefining a probability density function. When populations are not large enough, the histograms calculated with individual-based modeling are too poor to extract relevant statistics. It is then usually preferred to assume a certain family of probability density functions. The precision of the modified Fokker-Planck equation, however, depends only on the discretization method to solve the partial differential equation, and works for large and small populations whenever the assumptions of the Wiener process are satisfied (Gardiner, 2004).

In addition, individual-based modeling is characterized by requiring long computational times which make its use prohibitive in applications that demand many model evaluations (An et al., 2017), such as parameter estimation. The computation time of individual-based modeling grows exponentially with time, whereas the growth is linear for the equivalent modified Fokker-Planck equation. Note that the individual-based modeling approach requires one equation per cell. As cells grow exponentially, computation time scales linearly with the number of cells and exponentially with time. The modified Fokker-Planck equation, however, is a unique partial differential equation (PDE). Its computation time will depend on the degree of discretization and the simulation time. For the examples in Figure 1, computational times (2–6 s) are similar until time 15 (population of less than 3e4 cells) for both approaches. However from this time the modified Fokker-Plank equation becomes more efficient in orders of magnitude.

Moreover, the partial differential equation for stochastic growth has a diffusion term allowing efficient simulations using the appropriate techniques. Classical discretization methods transform the partial differential equation into a large number of ordinary differential equations. When the original equation is diffusive, a number of methods are at hand to take advantage of this property and significantly reduce the number of ordinary differential equations, thus reducing computational times (Trefethen, 2000; García et al., 2008).

3.2. Flow cytometry allows estimation of volume distributions for P. acidilactici

Flow cytometry is the standard technique for fast acquisition of population statistics. It is commonly employed to estimate different mammalian cell characteristics such as cell size using forward scattered light. This technique is also useful for bacteria, but as their diameters are close to the light wavelength (488 nm or 0.5 μm), side scattered light has better resolution and is preferred. Figure 2A shows how side scattered light (y-axis) discriminates among different bead diameters. Events in red, green, pink and blue correspond with beads of diameters 0.2 , 0.5 , 1 , 2 μm, respectively. Sizes smaller than 0.2 μm were below the detection limit of the device.

Figure 2

Flow cytometry is an efficient technique to extract volume distributions of coccoid bacteria such as Pediococcus acidilactici. (A) Side scatter light (y-axis) discriminates among round beads of different diameters represented in red, green, pink, and blue for diameters 0.2, 0.5, 1, and 2 μm, respectively. (B) Bead diameter correlates with side scatter as a second order polynomial (Julià et al., 2000; Prats et al., 2010) and can be transformed into bacterial diameter following the linear relationship in Chandler et al. (2011). (C) Flow cytometry correlogram and gating (red box) of fluorescent dye Sybrgreen with side scatter light for Pediococcus acidilactici after 8 h of growth. (D) Growth kinetics of P. acidilactici in terms of biomass and cell counting shows that the selected time where we took the sample (t = 8 h) is within the exponential phase (error bars are the confidence of intervals for n = 5 and α = 0.05). (E) Estimated diameter distribution of five replicas of the population of Pediococcus acidilactici at time 8. (F) Population volume distributions of Pediococcus acidilactici at time 8.

Estimation of bacterial diameters from side scattered light requires two steps: (1) to find the correlation between the bead diameter and side scattered light and (2) to transform the bead diameters to bacterial diameters. Figure 2B shows that the bead diameter is a second order polynomial of side scatter (Julià et al., 2000; Prats et al., 2010). However, bacteria and polystyrene beads have different refractive indexes. To correct the differences in the refractive index we make use of the linear relationship in Chandler et al. (2011). The figure shows in the right y-axis the final relationship between bacterial diameter (d) and side scatter (SSC), outlined in the following expression:

where SSC is the side scatter channel, a = 3.2911 and b = −0.2769 are the parameters provided in Chandler et al. (2011) to correct the differences in the refractive index and the estimated parameters of the second-order polynomial.

We acquired and processed side scatter data of P. acidilactici at one sampling time after 8 h of growth. Figure 2C shows the flow cytometry correlogram of Sybrgreen fluorescence and side scatter for one of the replicates. Red points indicate beads used to count the number of events.

The red box defines the gating where viable cells lie. We used two sources of information to define the gating: Sybrgreen fluorescence and experiments counting viable cells at different diameters. Sybrgreen was helpful to determine those events with too low DNA material to be consistent with a viable cell (upper and lower horizontal lines of red box). They probably represent either dead cells from the lag phase that have lost some of their DNA material, or free DNA detected as an event. Only with this gating, most diameters were between 1 and 2.5 μm as reported in the literature (Holt, 1994). However, the first calibrations of the model suggested that smaller cells were not able to divide. We passed cells through a 1.2 μm filter and found that only about 2% of the cells were viable. Consequently, we did not consider cells smaller than 1.2 μm (left vertical line of the red box).

In order to assume that the selected sampling time (t = 8 h) was within the exponential phase of growth, we estimated growth kinetics in terms of biomass and cell counting. Figure 2D shows the experimental curves with the standard sigmoid growth pattern for lactic acid bacteria. Both cases are described by the logistic equation (3) (R2 = 0.993–0.999). p-values from Fisher's F-test show consistency and robustness of the logistic to appropriately describe these profiles (Table 1). It is noted that no autocorrelation was observed in the fittings (data not shown). All parameters were statistically significant (t-Student test). The production of biomass was 25–30% lower as compared to previous cultures (Vázquez and Murado, 2008) which may be due to the minimum inocula employed in the present study. Typically, inocula used for the production of bacteriocins from lactic acid bacteria including P. acidilactici are much more populated, reaching values of 105−107 cfu/mL and longer productive periods (Vázquez et al., 2008).

Table 1

Summary of the parameter values obtained from the fittings of P. acidilactici growths (biomass and cells production) to the logistic Equation (3).

Parameters	Biomass	Cells
Gm	0.990 ± 0.021gL−1	(10.14 ± 0.72) × 108cells mL−1
μp	0.904 ± 0.021h−1	1.163 ± 0.394h−1
λ	8.23 ± 0.18h	6.72 ± 0.61h
μmp	0.224 ± 0.018gL−1h−1	(2.95 ± 0.93) × 108 cells mL−1h−1
R2	0.999	0.993
p-values	<0.0001	<0.0001

Statistical parameters R.

Five data replicates at t = 8 h and the relationship in (4) were combined to estimate the volume distribution of P. acidilactici. Figure 2E shows the bacterial diameter distribution of different replicates at time 8 h calculated from relationship (4). Volume distributions in Figure 2F are calculated by applying the volume equation of a sphere to the diameter histograms. Error bars show the mean and standard deviation of the histogram values for the five replicates. The distribution is smooth and skewed to the right, as expected from the modeling analysis in section 3.1. A similar behavior was observed in volume distributions of E. coli growing inside a population (Gangan and Athale, 2017).

3.3. The modified fokker-planck equation is combined with flow cytometry data to find the single-cell behavior of P. acidilactici

The volume distributions of the P. acidilactici at time 8 h provide enough information to find the single-cell parameters of the modified Fokker-Planck equation. The problem consists in finding the best set of parameters of the modified Fokker-Planck equation that represents the experimental volume distributions by solving a least square problem.

Preliminary computations demonstrated the inability of the model in (2a–2d) to reproduce the data, suggesting that the volume of each daughter is approximately one-fourth of the mother volume, i.e., υ = (1/4)υ. This hypothesis contradicts the principle of binary fission and had to be rejected because, if a mother gives two cells of this size, total volume would be destroyed and not conserved during division. This would imply that the specific growth rate of the population μ differs from the growth rate of the single-cell volumes μ, contradicting common observations in rod-shaped bacteria (Taheri-Araghi et al., 2015; Harris and Theriot, 2016). We could devise that this hypothesis, however, may be plausible in coccoid cells since cell volume, membrane area, and diameter scale differently. In fact, assuming that cells are spheres, simple calculations indicate that membrane area, instead of volume, is conserved if υ = (1/4)υ. However, it is well-known that coccoid bacteria create membrane (septal growth) before division (Pinho et al., 2013; Monteiro et al., 2015), suggesting again, that volume is conserved.

It resulted that the data pointed out to another mechanism to explain the daughter volumes: P. acidilactici, like Pediococcus pentosaceus and other coccoid cells (Zhou et al., 2010; Pinho et al., 2013; Monteiro et al., 2015), have two planes of division (Turner et al., 2010). That means that cells undergo two consecutive divisions and, at the time-scales of interest, one mother divides into four daughters with fourth of the mother volume. Hence we derive the individual-based and population-based modeling approaches considering two planes of division (see Appendix).

Figure 3A compares experimental data with both models at the population and single-cell levels considering two planes of division. Table 2 shows the single-cell parameters used in the simulations. The three elements exhibit the same steady-state volume distribution. The black line represents the experimental data (see also Figure 2F) and the red line and histogram are the solutions of the modified Fokker-Planck distribution and individual-based modeling, respectively.

Figure 3

The modified Fokker-Planck equation allows us to estimate single-cell behavior of Pediococcus acidilactici from acquisition based on cytometry data at one sampling time within the exponential phase. (A) The stationary distribution of Pediococcus acidilactici (black line) coincides with the modified Fokker-Planck equation (red line) and the individual-based modeling (blue histogram). The single-cell parameters to simulate both models were obtained by minimizing the differences between data and the stationary Fokker-Plank equation. (B) Single-cell dynamics of Pediococcus acidilactici with the estimated parameters. During the cycle of a single-cell, the volume growths four times and it divides in four daughters.

Table 2

Bounds and estimations of the best set of parameters of the modified Fokker-Planck equation to reproduce the experimental volume stationary distribution of P. acidilactici.

θ	θmin	θmax
υm = 4.5574	1	9
μ = 1.1619	0.7	1.5
ξ = 0.13439	0.075	0.2
σ = 0.3	0	0.3

We assume that bacteria have 2 division planes and therefore one mother produces four daughters of 1/4 of the mother's volume. Upper bound on the standard deviation of the statistics of division σ was considered 0.3 to avoid overlapping between distribution of mother sizes and daughter sizes.

All the estimated parameters are within the selected bounds (see Table 2) except for the standard deviation of the distribution of mother volumes (σ). Probably, this parameter tries to accommodate the errors for assuming that the mother divides into four perfectly round daughters. For coccoid cells with two planes of division, there is a short transition where cells are not completely round (Pinho et al., 2013; Monteiro et al., 2015). For the time-scales here considered this is a simplification that seems appropriate. Moreover, we have tried to reduce this bound resulting in similar estimations but with worse fit. For all these reasons we have considered an upper bound for the standard deviation of the mother distribution sufficiently small to avoid relevant overlapping between the mother and the daughter distributions. In this way, we force a scenario where the probability of having daughter volumes greater than mother volumes is low. Note that a standard deviation of 0.3 is reasonable attending to other mother distributions in the literature (Amir, 2014).

We also validate the results confirming that the specific growth rate of the population, μ in (3), coincides with the growth rate of the cell volume, μ in (1a) and (2a). Both population and volume should increase by one-fourth in one cell cycle. The specific growth rate of the population was estimated using the growth curves of P. acidilactici in Figure 2D. The calculated specific growth rate (1.1619 h−1) is similar to the rate using cell counting with cytometry (1.163 ± 0.394h−1) and larger than the rate estimated from the biomass growth curve (0.904 ± 0.021h−1).

Individual-based modeling of P. acidilactici can now be implemented using the estimated single-cell parameters in Table 2. Figure 3B depicts the simulations starting with one cell at the mean of the daughter volumes (υ = υ/4 = 1.14). We observe how the fluctuations, dependent on ξ, resembles experimental single-cell dynamics in the literature (Deforet et al., 2015). In order to obtain such simulations, it is critical to consider that during division one mother gives four daughters, volume is conserved and that volumes grow four-fold during a cell-cycle.

4. Concluding remarks

In this work, we have developed a modified Fokker-Planck equation describing the statistics of a population from their single-cell parameters. The model is based on the assumptions that cell volumes grow exponentially and cells divide following the sizer principle. We have tested and numerically compared the modified Fokker-Planck with its equivalent individual-based modeling approach. Simulations resulted critical to understand several observed phenomena during the exponential phase of growth of bacteria, including the necessity of considering stochasticity to obtain a distribution of volumes that is time-invariant and similar to experimental observations.

The modified Fokker-Planck equation is also a powerful tool to estimate the behavior of single-cells inside populations. Instead of requiring single-cell measurements, we make use of flow cytometry to find the volume distribution of a population of P. acidilactici within the exponential phase. The combination of the modified Fokker-Planck equation with data provides information about the growth rate and stochasticity of the single-cell volume, as well as the statistics of the mother and daughter volumes. For a good correspondence between model and data, it is fundamental to assume that the P. acidilactici have two planes of division and its volume and population numbers grow by a fourth during a cell cycle.

The proposed methods allow efficient analysis of flow cytometry data to find single-cell behavior. In fact, they pave the way for studying cell heterogeneity in numerous applications in food microbiology, such as in quantitative risk assessment and prediction of shelf life.

Author contributions

MG and AA: designed the theoretical study, developed the methodology, and drafted the manuscript; JV: designed the experimental study and analyzed the data; IT: collected the data in the flow cytometry; MG: conducted the modeling experiments. All authors approved the final version of the manuscript.

Conflict of interest statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.