Baranyi model with coupled secondary models

library(tidyverse)
library(biogrowth)

Theoretical background

The Baranyi model has demonstrated its reliability for describing the growth of bacterial populations. When combined with the Ratkowsky secondary model, it can describe how the microbial population grows for different temperatures.

The parameters of the Baranyi-Ratkowsky model must be estimated from empirical data. The most common approach is to perform a set of experiments at different temperatures. Then, the model is fitted in a two-steps fashion. First, the primary Baranyi model is fitted to the data obtained at each temperature. Next, based on the values of \(\mu\) obtained for each temperature, the parameters of the Ratkowsky model are estimated. An independent secondary model for the lag phase is defined in a similar way, based on the estimates of \(\lambda\) from the primary model.

\[ \sqrt \mu = b\left(T-T_{min} \right)^2 \]

A recent study revisited the Baranyi-Ratkowsky model, reaching two main conclusions (Garre et al. 2025). The first one is that an inverse-squared relation is the only secondary model for \(\lambda\) compatible with the Baranyi-Ratkowsky model:

\[ \lambda = B\frac{1}{(T-T_{min})^2} \]

The second conclusion is that there is a coupling between the secondary models for \(\mu\) and \(\lambda\). Namely, parameter \(B\) is defined as:

\[ B = \frac{\ln \left(1 + 1/C_0 \right)}{b} \]

Please note that the slope of the Ratkowsky model (\(b\)) also appears in this equation. Then, the secondary model for the lag phase only introduces an additional parameter (\(C_0\)), which is related to the hypotheses of the Baranyi model related to the lag phase.

This opens new possibilities for fitting the Baranyi-Ratkowsky model, exploiting the coupling between both parameters to improve model robustness.

The fit_coupled_growth() function

The biogrowth package includes the fit_coupled_growth() function to fit the Baranyi-Ratkowsky model considering the coupling between both parameters. This function has the following arguments:

Two-steps fitting of the Baranyi-Ratkowsky model

The two-steps fitting mode considers that the primary models have already been fitted to the dataset. Therefore, a table of values of \(\mu\) and \(\lambda\) estimated for each temperature are already available. It must be presented as a tibble (or data.frame) with three columns:

The package includes the example_coupled_twosteps as an example of the data format.

data("example_coupled_twosteps")
example_coupled_twosteps
#> # A tibble: 4 × 3
#>    temp    mu lambda
#>   <dbl> <dbl>  <dbl>
#> 1    10  0.03  20.1 
#> 2    15  0.07   9.61
#> 3    25  0.45   5.4 
#> 4    30  0.94   2.68

The fit_coupled_growth() function uses nonlinear regression. Therefore, it requires an initial guess for every model parameter. This must be defined as a named numeric vector with the following elements: Tmin (theoretical minimum temperature for growth), b (slope of the Ratkowsky model), logC0 (decimal logarithm of the initial value of the theoretical substance \(C_0\) of the Baranyi model):

guess <- c(logC0 = -1, b = .1, Tmin = 5)

Once the initial guess has been defined, the fit_coupled_growth function can be called:

fit_twosteps <- fit_coupled_growth(example_coupled_twosteps, 
                   start = guess,
                   mode = "two_steps")

It returns an instance of FitCoupledGrowth with several S3 methods. The print() method shows a summary of the fit:

fit_twosteps
#> Baranyi model fitted from data considering coupling between secondary models
#> 
#> Fitting approach: two_steps 
#> 
#> Estimated parameters:
#>       logC0           b        Tmin 
#> -0.52237218  0.02905222  0.88742143

Detailed statistics of the fit are shown by the summary() method:

summary(fit_twosteps)
#> 
#> Parameters:
#>        Estimate Std. Error t value Pr(>|t|)   
#> logC0 -0.522372   0.303201  -1.723   0.1455   
#> b      0.029052   0.004318   6.728   0.0011 **
#> Tmin   0.887421   2.581341   0.344   0.7450   
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 0.3047 on 5 degrees of freedom
#> 
#> Parameter correlation:
#>        logC0      b   Tmin
#> logC0 1.0000 0.3079 0.8129
#> b     0.3079 1.0000 0.7759
#> Tmin  0.8129 0.7759 1.0000

The plot() method compares the model fitted against the experimental data.

plot(fit_twosteps)

This method returns a ggplot object, so it can be editted with further layers:

plot(fit_twosteps) + scale_y_sqrt()

The fit_coupled_growth() function also includes the known argument that allows fixing any model parameter. For instance, we can fix parameter Tmin to 5ºC:

known <- c(Tmin = 5)

Please note that now we need to update the initial guess, as no parameter can appear both as an initial guess and a known parameter:

guess <- c(logC0 = -1, b = .1)

We can now fit the model by calling fit_coupled_growth() passing also this parameter

fit_twosteps_fixed <- fit_coupled_growth(example_coupled_twosteps, 
                   start = guess,
                   known = known,
                   mode = "two_steps")

Now, the statistical summary only provides information on two parameters, as Tmin was fixed.

summary(fit_twosteps_fixed)
#> 
#> Parameters:
#>       Estimate Std. Error t value Pr(>|t|)    
#> logC0 0.003418   0.152971   0.022 0.982896    
#> b     0.035797   0.003965   9.029 0.000103 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 0.3635 on 6 degrees of freedom
#> 
#> Parameter correlation:
#>         logC0       b
#> logC0  1.0000 -0.8703
#> b     -0.8703  1.0000

We can check the impact on the fitted curve using the plot() method.

plot(fit_twosteps_fixed)

One-step fitting of the Baranyi-Ratkowsky model

The one-step fitting mode estimates secondary models directly from the log-microbial concentrations. Hence, the data must be defined as a tibble (or data.frame) with three columns:

The package includes the example_coupled_onestep as an example dataset:

data("example_coupled_onestep")
head(example_coupled_onestep)
#> # A tibble: 6 × 3
#>     time  temp  logN
#>    <dbl> <dbl> <dbl>
#> 1     0      6  2.22
#> 2  2159.     6  1.73
#> 3  4317.     6  2.39
#> 4  6476.     6  2.46
#> 5  8635.     6  4.02
#> 6 10793.     6  5.20

The function uses nonlinear regression for parameter estimation. Therefore, it requires initial guesses for every model parameter. They must be defined as a named vector including:

guess <- c(logN0 = 2, logNmax = 8, b = 0.04, logC0 = -4, Tmin = 5)

Once the data and guess have been defined, the model can be fitted by passing mode = "one_step" to fit_coupled_growth():

fit_onestep <- fit_coupled_growth(example_coupled_onestep,
                             start = guess,
                             mode = "one_step")

It returns an instance of FitCoupledGrowth with several S3 methods. The print() method shows a summary of the fit:

fit_onestep
#> Baranyi model fitted from data considering coupling between secondary models
#> 
#> Fitting approach: one_step 
#> 
#> Estimated parameters:
#>       logN0     logNmax           b       logC0        Tmin 
#>  2.01474226  7.95150807  0.04095209 -4.31665354  5.00139993

Detailed statistics of the fit are shown by the summary() method:

summary(fit_onestep)
#> 
#> Parameters:
#>           Estimate Std. Error t value Pr(>|t|)    
#> logN0    2.0147423  0.0444055   45.37   <2e-16 ***
#> logNmax  7.9515081  0.0586472  135.58   <2e-16 ***
#> b        0.0409521  0.0005485   74.66   <2e-16 ***
#> logC0   -4.3166535  0.2097269  -20.58   <2e-16 ***
#> Tmin     5.0013999  0.0079633  628.06   <2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 0.3186 on 157 degrees of freedom
#> 
#> Parameter correlation:
#>              logN0    logNmax       b      logC0       Tmin
#> logN0    1.000e+00 -0.0221747  0.1570 -0.3750199 -4.268e-06
#> logNmax -2.217e-02  1.0000000 -0.1961  0.1659043  4.431e-04
#> b        1.570e-01 -0.1961348  1.0000 -0.9448999  1.318e-01
#> logC0   -3.750e-01  0.1659043 -0.9449  1.0000000  6.326e-04
#> Tmin    -4.268e-06  0.0004431  0.1318  0.0006326  1.000e+00

The plot() method compares the model fitted against the experimental data.

plot(fit_onestep)

Any model parameter can be fixed using the known argument, which must be a numeric vector:

known <- c(logN0 = 3)

Please note that the initial guess must be updated, so no parameter appears both as fixed as an initial guess:

guess <- c(logNmax = 8, b = 0.04, logC0 = -4, Tmin = 5)

Then, the function can be called again, including the additional argument.

fit_onestep <- fit_coupled_growth(example_coupled_onestep,
                             start = guess,
                             known = known,
                             mode = "one_step")

Please note how the statistical summary now excludes parameter logN0

summary(fit_onestep)
#> 
#> Parameters:
#>          Estimate Std. Error t value Pr(>|t|)    
#> logNmax  7.932401   0.120174  66.008  < 2e-16 ***
#> b        0.042519   0.001638  25.964  < 2e-16 ***
#> logC0   -5.953190   0.652192  -9.128 3.14e-16 ***
#> Tmin     5.001981   0.016893 296.102  < 2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 0.6743 on 158 degrees of freedom
#> 
#> Parameter correlation:
#>            logNmax        b    logC0      Tmin
#> logNmax  1.0000000 -0.18978  0.17204 0.0002007
#> b       -0.1897806  1.00000 -0.98304 0.0960106
#> logC0    0.1720404 -0.98304  1.00000 0.0018395
#> Tmin     0.0002007  0.09601  0.00184 1.0000000

The impact of the overall model fitting can be evaluated using the plot() method.

plot(fit_onestep)

References

Garre, A., Valdramidis, V., Guillén, S., 2025. Revisiting secondary model features for describing the shoulder and lag parameters of microbial inactivation and growth models. International Journal of Food Microbiology 111078. https://doi.org/10.1016/j.ijfoodmicro.2025.111078