# 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:

- `fit_data` defines the data for the fitting. Its format will depend on the type of fit (see below).
- `start` initial guess for the model parameters.
- `known` vector of fixed model parameters.
- `mode` the function defines two types of fitting: one-step or two-steps. This argument allows switching between both.
- `weight` type of weight to apply for the fitting (only for two-steps fitting).
- `...` ignored
- `logbase_mu` logarithmic base for the definition of `mu` (only used for two-steps fitting). By default, `exp(1)` (natural logarithm).
- `logbase_logN` logarithmic base for the definition of `logN` (only used for one-step fitting). By default, `10` (decimal logarithm).

### 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:

- temp: the temperature of the experiment
- mu: the value of \(\mu\) estimated at that temperature
- lambda: the value of \(\lambda\) estimated at that temperature

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:

- `temp`: the storage temperature
- `time`: the elapsed time of each sample
- `logN`: the log-microbial concentration. By default, decimal logarithm, although this can be modified using the `logbase_logN` argument.

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:

- `logN0`: log of the initial microbial concentration
- `logNmax`: log of the maximum microbial concentration
- `b`: slope of the Ratkowsky model
- `logC0`: decimal logarithm of the initial value of theoretical variable \(C\) from the Baranyi model
- `Tmin`: theoretical minimum growth temperature.

    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