Generalized logistic model of bacterial growth

ABSTRACT This work proposes a new mathematical model describing the dynamics of growing bacterial cultures. The model, described by a first order non-linear differential equation, as a generalization of the logistic equation, was compared with the most studied mathematical models. All models were numerically implemented and fitted to the experimental data, collected from the incubation of a bacterial strain of Pseudomonas fluorescens, to obtain the growth parameters. The experimental data showed the lowest fit error for both the Baranyi–Roberts and new models, which turned out to be equivalent. Simulations of the fitting algorithm were also implemented and repeated for a large number of initial guesses of the parameters, chosen in order to test the fitting and convergence performances. The innovative feature that makes the new model easier to use than Baranyi–Roberts model is definitely its simple and manageable analytical form and its good performance in terms of convergence time.


Introduction
The modelling of microbial growth is certainly one of the most important goals of predictive microbiology [1].In fact, as a result of growth, the risk associated with the presence of a microbial danger (a pathogenic microorganism or a toxic product of its metabolism) can increase significantly, both in terms of probability for the consumer of contracting a disease and in terms of probability as a food business operator to exceed the microbial load limits established by national or international regulations [2,3].The fundamental tools for food industries to predict and control the stages of growth, survival and death of bacteria are mathematical models.
Bacteria in nature can be present both in free form (planktonic form) and organized in structures called biofilms, which are aggregations of microorganisms, microbial cells, and their extracellular material that form thin films adhering to surfaces.These surfaces, called support surfaces, can be made up of biological material (such as a leaf or an animal tissue) or abiotic material (plastic, glass, etc.).The formation of the biofilm helps the bacterial community to protect itself from environmental stresses (such as temperature variations and dehydration), favours the retention of nutrients and creates a real communication between the bacteria [4].The control of bacterial growth and therefore of bacterial biofilm formation is of fundamental importance both in the health sector, for the treatment of pathologies (cystic fibrosis, tuberculosis, etc.), the improvement of antibiotics, and the safety of medical devices [5], and in the industrial sector, to prevent the obstruction and deterioration of the pipelines of the plants (e.g.aqueducts) [6].
The growth of a bacterial culture is characterized by the succession of four phases: lag, exponential, stationary and death phases.The lag phase, which can have a variable length depending on the type of bacteria and on the changes in the physico-chemical environment such as temperature, pH, water activity and nutrient availability, is characterized by cellular activity but not growth.In fact, before the bacteria start to duplicate, when they are transferred in growth conditions different from the starting ones, it is necessary that they synthesize all the enzymes and coenzymes necessary for growth; the lag phase is therefore an adaptation phase in which there is no cell duplication.After the lag phase, bacterial cells enter the exponential phase.This is the time when cells divide by binary fission and double in number after each generation time.Metabolic activity is high as DNA, RNA, cell wall components and other substances required for growth are generated for division.The progressive exhaustion of available resources and the accumulation of waste products then determine the end of the exponential phase and the entry into the stationary phase of growth, in which there is a balance between dividing cells and cells that begin to die, hence the number of cells in culture remains constant.As nutrients become less available and waste products increase, the number of dying cells continues to rise.Finally, in the death phase, the number of living cells decreases exponentially and population growth experiences a sharp decline.As dying cells lyse or break open, they spill their contents into the environment making these nutrients available to other bacteria.This helps the spore-producing bacteria survive long enough for spore production.Spores are able to survive the harsh conditions of the death phase and become growing bacteria when placed in an environment that supports life [4,7].
It is therefore important to develop mathematical models that predict microbial growth and therefore give a quantitative estimate of the changes in a bacterial population under certain experimental conditions.Predictive models have been studied in the literature, including the logistic model, Baranyi-Roberts, Reparameterized Gompertz and Huang models, which describe bacterial growth as a function of time, and this paper in particular aims to develop an alternative model capable of predicting bacterial growth as a function of time, which is a generalization of the logistic model.In particular, the new model improves the logistic model because it manages to overcome the limits of the latter and at the same time has a simpler analytical form than the other three most studied models.Furthermore, this work aims to elaborate a computational analysis of the mathematical models most used to describe the microbial growth and the new model, proving that the latter is equivalent to the Baranyi-Roberts model in terms of accuracy.Moreover, it has been demonstrated that the analytical simplicity of the model does not represent a limit for its efficiency and in particular does not affect its performance in terms of convergence of the fitting algorithm.
Hence, the objective of this research was to develop an alternative and more costeffective non-linear mathematical model than those already known in literature, which describes the dynamics over time, of a sample of Pseudomonas fluorescens (P.F.) bacteria.The P.F.sample was incubated at room temperature in Tryptic soy broth (TSB) culture medium for about 20 h.Since a medium containing a bacterial inoculum becomes more and more turbid with the progress of growth, and this increase in turbidity is proportional to the number of cells present, optical density (OD) measurements were made to monitor the quantitative variation of the biofilm and therefore of the biomass [8].

Growth models in literature
Among the most studied primary models in literature, it is worth mentioning Baranyi-Roberts, Reparameterized Gompertz and Huang models.Primary models describe how population density changes with time in a specified environment and are depicted as microbial growth or death curves.
The Huang model is expressed as or equivalently if the quantity y t ð Þ taken into consideration is expressed as the natural logarithm of the cell concentration: where t is the time h ð Þ, λ is the lag phase duration h ð Þ, y t ð Þ represents the natural logarithm of microorganism count ln CFU=g ð Þ in eq. 3 ð Þ or the optical density OD ð Þ of the bacteria in eq. 1 ð Þ, which is related to the cell number, depending on the type of the considered measured quantity, α is a constant which defines the transition from the lag phase to the exponential phase of the growth curve (α ¼ 4Þ, y 0 is the initial microorganism count, y max is the bacterial count at the stationary phase and μ max is the maximal specific growth rate h À 1 ð Þ [9].The Baranyi-Roberts model [10] is expressed as or equivalently if the quantity y t ð Þ taken into consideration is expressed as the natural logarithm of the cell concentration: where h 0 is a dimensionless parameter which represents the initial physiological state of the cells and it is equals to μ max λ, and the other variables t; y t ð Þ; y 0 ; y max ; and μ max are defined as those used in the Huang model [11].
The Reparameterized Gompertz model [1,12] is expressed as μ max e ymaxÀ y 0 λÀ t ð Þþ1 (7) where all parameters λ; t; y t ð Þ; y 0 ; y max ;and μ max are identical to those used in the Huang model.

New model derivation
The new proposed model describes the evolution over time of a bacteria population, on the basis of the following assumptions, which allow to easily use it [13,14].
• We disregard the possible spatial distribution of the bacterial cells assuming that the bacterial population is homogeneous.• All cells are equally reproductive.
• The reproduction process is continuous.
• Each new cell is immediately fertile.
We consider the following model variables: y � t ð Þ is cell concentration, and S the concentration of the nutrient.Assuming that the absence of nutrient has a limiting effect according to the Monod equation [15][16][17] the instantaneous specific growth rate μ, which can be explained as the instantaneous measure of the birth rate per cell per unit time as a function of substrate concentration, is given by: where μ max is the maximum specific rate of potential growth h À 1 ð Þ, K s is the 'half-velocity constant' that is the substrate concentration corresponding to 1  2 μ max , and l ¼ S K s þS is a limiting factor.
Since the culture conditions are adjusted so the substrate concentration is always in excess, nutrient limitation can be excluded, hence also to simplify and use the model easily, we can neglect the nutrient modelling and we can introduce a simple limiting L t ð Þ which describes the transition from the exponential to the stationary phase, depending on the parameter y max called carrying capacity which stands for the maximum cell density [10]: It should be noted that the maximal specific growth rate μ max is typical of the species as it is related to the reproductive mechanism and to the incubation temperature of that species, it increases with temperature in the temperature range considered and is accompanied by a decrease in the duration of time for which the growth rate is near this maximum, and it depends also on the nature of the substrate, instead the carrying capacity y max depends mainly on the environmental conditions but not on temperature, it remains almost constant to the temperature range considered [18,19].Moreover, we introduced the parameter y s that has the function of regulating the growth rate of the model while preserving and respecting its mathematical properties (such as existence of the inflection point and convergence to the carrying capacity).In particular, the parameter y s represents a concentration, and, as shown in Figure 1, it is a shift parameter that regulates the duration of the lag phase but also the position of the point of maximum growth and therefore of the time and the value of the bacterial concentration at which the bacterial growth rate begins to decrease and enters the stationary phase.Hence, the model under consideration reads as follows: and its analytical solution is written as: where y m ¼ y max : The Equation ( 11)can be expressed in the following form: where are constant.We notice that the proposed model presents the same number of parameters as the existing ones, thus it does not introduce any additional complexity from a biological or computational point of view.
With simple calculations, it can be shown that the Baranyi-Roberts and the new proposed models are both a generalization of the logistic equation, in fact, we can derive the logistic equation from both when we assume the parameter λ ¼ 0 in the first and y s ¼ 0 in the second [14,20].It has been verified that the logistic equation has some limitations and often fails to model data satisfactorily, so for this reason it has been subject to some revisions and refinements over the years [21].
Furthermore, it is possible to verify analytically that they are equivalent, that is, even if the parameters do not represent the same quantities, indeed the unit measure of λ and y s are time and concentration, respectively, these two parameters are actually related to each other.As Figure 2(b) shows, while λ regulates the lag phase, y s models the lag phase together with the bacterial growth rate.In fact, according to the Baranyi-Roberts model  the cell concentration dynamics of a bacterial culture is described by the following nonlinear differential equation [14]: where and after simple computations we get the following relationship which correlates the Equations (10,14,15): where y 0 is the initial microorganism count, μ max BR and μ max NM are the maximum specific rate of potential growth of the Baranyi-Roberts model and the new model, respectively.From Equation ( 16) we can verify that the parameters λ and y s of the two models 4 ð Þ À 5 ð Þ and 11 ð Þ have some properties in common.In fact, assuming that the lag phase is quite long, from Equations (4, 5) we get lim Moreover, from Equation ( 16) if λ !þ1, it follows that y s !y 0 , and so lim Similarly, if we assume the lag phase is very small, we have that lim , that is the logistic equation, and from Equation ( 16) if λ !0, it follows, assuming that y s !0, therefore lim In conclusion, the two models describe similarly long and short lag phases of bacterial growth, which can then be considered equivalent.

Analysis of the model
The function y t ð Þ, as it is defined in Equation (14, 15) must be positive and must also have a positive first derivative, hence it must be monotonically increasing.Moreover, y t ð Þ must have an inflection point that marks the change in the direction of the growth rate of the model dy t ð Þ dt defined by Equation (10) [7].We then determine the definition interval of the parameter y s that satisfies the above mentioned properties of y t ð Þ.First of all, for y t ð Þ to be defined, should be y s �y 0 ; y max .Analysing the first derivative of y t ð Þ: where c; b are defined as 13 ð Þ, we have that _ y t ð Þ � 0 $ y s < y 0 < y m Furthermore, studying the second derivative of y t ð Þ we observe that: and y :: where Hence, has the characteristics of an inflection point for y t ð Þ.Now let us verify what value y s must assume in order that y t ð Þ maintains the inflection point.We have already verified that y t ð Þ is monotonically increasing if and only if y s < y 0 , and now we want to compute the lower bound of the defining interval of y s .
is defined if and only if 2y 0 À y m � y s < y 0 .However, we need to analyse two possible cases: • If 2y 0 À y m � 0, the initial condition is quite close to the maximum concentration y max , therefore we can state that y s 2 2y 0 À y m ; y 0 ½ ½. • If 2y 0 À y m < 0, since y s is a concentration, we can state that y s 2 0; y 0 ½ ½.
The equilibrium points of the autonomous differential equation of the model 10 ð Þ are: In order to study the stability of the equilibrium points, let us apply a small perturbation to the equilibrium point y � , that is y t ð Þ ¼ y � þ εỹ t ð Þ, with ε arbitrarily small.Linearizing around the equilibrium point, the reproduction function can be written as follows: Therefore, in the neighbourhood of y � it holds that: which has the following solution: where ỹ0 is the initial perturbation.
follows that y s is an unstable equilibrium point, moreover y m is an asymptotically stable equilibrium point and it is also an attractor, and its basin of attraction is the interval y 0; y m � � : Another information that allows us to understand the dynamics of the response to perturbations around the equilibrium y max is the characteristic return time t R .In fact, it was studied that any enclosed population experiencing density-dependent mortality or fecundity will tend to return towards an equilibrium level, and in particular, the characteristic return time t R gives us a measure of the rate at which this equlibrium is reached [22].
From the equation 23 ð Þ with y � ¼ y m we have: that is equivalent to Hence, according to [22], Equation (26) gives us a measure of the rate at which this equilibrium point y m is reached.Larger values of μ max b lead to a more rapid approach to y m .

Trend of growth and rate of growth
From Equation (10) it is easy to deduce that the growth rate dy t ð Þ dt is positive, and it assumes its maximum value at the inflection point of y t ð Þ, which in this case is obtained from the following expression: Figure 2(a) shows how by varying the parameter μ max the population growth rate and the duration of the lag phase also vary, whereas the ordinate of the inflection point remains constant [23].
Instead, as we can see in Figure 2(b), by varying the parameter y s , both the population growth rate and the location of the inflection point vary.In particular, the inflection point represents the point where the population growth rate starts to decrease resulting in the convergence of the total biomass to the carrying capacity y max .As y s increases, the growth becomes slower, indeed the abscissa and the ordinate of the inflection point increase.It also affects the lag phase, in fact as y s increases, the lag phase increases.

Bacteria sample
The analysed models were numerically implemented by comparing them with the experimental data obtained using Pseudomonas fluorescens bacteria samples purchased from the collection of the Leibniz Institute DSMZ -German Collection of Microorganisms and Cell Cultures (Braunschweig, Germany).These are rodshaped gram-negative strains, which secrete a yellow-green fluorescent pigment which can be applied to study the influence of antibacterial agents on bacterial proliferation in static conditions [24,25].For growth curve determination, some bacterial suspensions were collected and transferred to Petri dishes containing nutrient agar to permit the growth of bacterial colonies.P.F. was incubated at room condition (21-22°C) under a laminar chamber till the evident formation of colonies on the agar plate and stored at 4°C until used.A loopful of bacterial biomass was transferred from the nutrient agar plate to a volume of 30 g/L Tryptic soy broth (TSB) as a bacterial medium for 24 h at room temperature (R. T.).Bacteria stocks were centrifuged for 10 min at 4000 rpm.The supernatant was discarded, and desired media such as TSB was added to the test sample.Three replicates of each sample were prepared for this study.

Bacterial growth curve measurements
The growth curve measurements were performed by measuring the OD of the bacteria sample, which is related to bacterial concentrations in cultures [26].Before analysis, their initial optical density (OD) was adjusted identically for all replications.Then, average values of OD for Pseudomonas fluorescens strain at 570 nm were recorded every 30 min for the first 8 h and then every 60 min, with a gap of 7 h, for the last 6 h (plateau reached) to monitor bacterial growth.Finally, the obtained control and sample curves were compared.Absorbance growth results from bacteria growth.

Fitting algorithm and model parameter estimation
The developed model was numerically implemented in Matlab and used to fit the experimental data in order to estimate the model parameters y s ; y max and μ max .A nonlinear least square fitting, using the optimization function lsqnonlin, based on the Levenberg-Marquardt algorithm, was used to adjust the model's unknown parameters aiming at minimizing the error function, which is the mean squared value of the relative error between the observed growth data and the results of the numerical analysis [27].For # ¼ y s ; y max ; μ max À � the error function is expressed as In Equation ( 29) y i is the experimental measurement of bacterial concentration at time t i , ŷi # ð Þ represents the bacterial concentration predicted by the model at same time t i and N is the number of measurements collected during the experiments.
Figure 3 shows the result of the model fitting, in which the experimental data concern the variation of the optical density of a P.F.sample, with initial OD about 0.182, during an incubation interval of about 20 h.In particular, as can be seen from Figure 3, the experimental measurement has a gap where the data is missing, due to the instrument and the manual measurement technique used to collect the experimental data which therefore lead to an interruption of the measurement process, nevertheless using the parameters estimated by the fitting algorithm, the model was used to simulate the trend of the experimental data in that data gap.It can be seen from Figure 3, despite the lack of a part of the data, the model, like the others studied, is able to predict satisfactorily the experimental data.
Figure 4 shows that the growth data of Pseudomonas fluorescens on TSB culture, at storage room temperature of about 22°C, were also fitted to the other three known models: Baranyi-Roberts, Reparameterized Gompertz and Huang models, to compare the estimated model parameters and the statistical parameters for model validation with the proposed new model.
Table 1 shows the modelled growth parameters and some statistical parameters calculated for the validation of the model: root mean square error (RMSE), sum of squared errors (SSE), accuracy factor (A f ) and the bias factor (B f ).The root mean square error is given by: RMSE ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi P N i¼1 y i À ŷi it measures the average magnitude of the error, that is, it provides information on how close the predicted data are to the observed values.The RMSE of the new model is 0.0072, and being very small, in particular almost close to zero, indicates that the fit is good.The comparison between the measured and predicted data was based also on the accuracy and bias factors.The accuracy and bias factors provide an indication of the average deviation between the model predictions and observed results [28] and their closeness to a value of 1 is an effective and practical measure of predictive model validity.These statistical characteristic indices are expressed as follows: where y i is the value of the data experimentally measured, ŷi is the value predicted by the model at the same time t i as the data were measured and N is the number of measurements.Table 1 shows the A f and B f values for the new model which are 1:012 and 0:9967 respectively, indicating how the model fits satisfactorily with the experimental data.
As shown in Table 1 all the statistical parameters analysed (RMSE, MSEP, A f , B f ) show larger values for the Gompertz and Huang models, moreover they assume the same values for our new proposed model and Baranyi-Roberts model.Hence, the latter two

Model validation
The developed new model was used to simulate the microbial growth of other bacterial strains in response to different environmental conditions, i.e. the culture environment and the incubation temperature.A comparison of the experimental data, obtained at two different temperatures, 20°C and 25°C, respectively, from the bacterial growth of Pseudomonas spp. in which culture four bacterial strains were combined, P. fluorescens, P. miguale, P. tolaasii and P. agarici [29], and of the predictions of the new model are shown in Figures 5, 6.The Pseudomonas spp.strains were incubated in Agaricus bisporus fruiting bodies.
In particular, as shown in Figures 5, 6, the experimental data of the bacterial growth are represented by the natural logarithm of microorganism count ln CFU=g ð Þ as a function of time.
In both cases the model can fit satisfactorily the experimental data with a root mean square error of 0.0039 and 0.0105 for data obtained at 20°C and 25°C respectively and can be used to analyse the bacterial growth curves of various bacterial strains at different incubation temperatures.
It has been verified that even in the absence of a part of the experimental data, such as those in Figure 3, the model satisfactorily predicts the experimental data, managing to simulate the missing data.

Performance comparison of Baranyi-Roberts and new models
We have verified that the new proposed model and the Baranyi-Roberts model provide equivalent fitting results, both in terms of estimated model parameters and statistical parameters for measuring the goodness of the model.Furthermore, by comparing the expressions of the two models, Equations (4-6) and Equations (14,15) one immediately realizes that the new proposed model is much simpler, more intuitive, and easier to use than the Baranyi-Roberts model and therefore the probability of making computational errors will be less.This is definitely an aspect that makes us prefer the new model over the others already known.
Another aspect to be analysed in order to assess the efficiency of the models concerns the evaluation of their performance in terms of convergence of the fitting algorithm.The convergence speed of the fitting algorithm according to the two analytical models was evaluated by implementing fitting algorithms with initial guesses of the parameters uniformly distributed in intervals like I μ max ; ε μ max À � where ε μ max is a fixed percentage of variation of the parameter μ max and the same for the other parameters.The goal was to calculate the number of iterations required by the fitting algorithm, according to each model, to converge to the solution starting from initial guesses chosen randomly and so trying to simulate the greatest number of cases that can actually occur when making a choice on the initial parameters.
Figure 7 shows the results obtained by 10.200 simulation runs, considering each implementation of the fitting algorithm with different initial guesses, each chosen assuming a maximum error of 40% with respect to the estimated parameters.The distributions corresponding to the number of iterations necessary for the convergence of the fitting algorithm for each of the two models are also shown.As can be seen from Figure 7, the distribution of the data relative to the new model is slightly shifted to the left with respect to the distribution of Baranyi-Roberts model, it has an average of about 11 iterations and the largest number of analysed cases converges with a number of iterations less than or equal to this average, however both the two models converge with a rather small number of iterations.Hence, following this analysis we have verified that the simplicity of the new model does not even affect the performance in terms of convergence of the fitting algorithm, therefore both models show satisfactory results and therefore a rapid convergence.

Conclusions
In this work the authors proposed a new mathematical model to predict the growth of Pseudomonas fluorescens.The new proposed model was studied and compared with the models known and studied in literature.The developed model was numerically implemented in Matlab and used to fit the experimental data obtained by incubating a bacterial strain of P.F. with TSB culture for about 20 h.By comparing the error of the fitting algorithm and the estimated parameter values of the models, we verified that the new model and the Baranyi-Roberts model fit satisfactorily to the initial data with a lower root mean square error than the other two considered models.In particular, it has been verified that the new model and the Baranyi-Roberts model are equivalent both from a numerical and theoretical point of view, and although the two parameters y s for the new proposed model and λ for Baranyi-Roberts model represent different quantities we have shown that there is a relationship between the new parameter introduced by the proposed model and the lag phase parameter λ used in the other three models, being able to conclude that the shift parameter y s is able to model the lag phase together with the bacterial growth rate in a very satisfactory way.We have also demonstrated that the new model and the Baranyi-Roberts model both show satisfactory performances in terms of convergence speed of the fitting algorithm, therefore we can conclude that the new developed model is a generalization of the logistic model and therefore manages to overcome the limitations of this last but at the same time its simple and intuitive analytical expression does not affect its performance in terms of accuracy and convergence of the fitting algorithm which are almost equivalent to those of the Baranyi-Roberts model.

Figure 1 .
Figure 1. a ð Þ Growth rate over cell concentration diagram of a population of Pseudomonas fluorescens.b ð Þ Cell concentration (blue line) and growth rate (red line) over time diagram of a population of Pseudomonas fluorescens.It is assumed for both that y 0 ¼ 0:182OD, y max ¼ 1:095OD, y s ¼ 0:128 and μ max ¼ 0:56h À 1 .

Figure 2 .
Figure 2. Cell concentrations (solid-blue lines) and growth rates (dashed-red lines) over time a ð Þ varying the maximal specific growth rate μ max , b ð Þ varying y s (chosen in the range 2y 0 À y m ; y 0 ½ ½).

Figure 3 .
Figure 3. Results of the model fitting of the growth of Pseudomonas fluorescens inoculated with TSB culture at room temperature (about 22°C).The x-markers are the measured data during the growth of P.F.sample, o-markers are the data predicted by the model.

Figure 4 .
Figure 4. Comparison of results of the three models fitting of the growth of Pseudomonas fluorescens inoculated with TSB culture at room temperature (about 22°C).The x-markers are the measured data during the growth of P.F.sample, o-markers are the data predicted by the model.

Figure 5 .
Figure 5.Comparison of experimental data and of the new model prediction for the bacterial growth of Pseudomonas spp. on Agaricus bisporus fruiting body incubated at temperature of 20°C.

Figure 6 .
Figure 6.Comparison of experimental data and of the new model prediction for the bacterial growth of Pseudomonas spp. on Agaricus bisporus fruiting body incubated at temperature of 25°C.

Figure 7 .
Figure 7. Distribution of the number of iterations of the fitting algorithm to converge to the solution for the new model (blu line) and the Baranyi-Roberts model (red line).

Table 1 .
Estimated parameters of the four different growth models for P.F.models were more suitable for analysing the growth curve than the Reparameterized Gompertz and Huang models, as shown in Table1and Figures3, 4. In particular, by comparing the fitting results and the statistical parameters in Table1, we can conclude that the new model and the Baranyi-Roberts model have the same accuracy in terms of prediction of the experimental data and therefore the analytical simplicity of the new model does not affect its performance.