Baranyi and Roberts Model#

The Baranyi and Roberts growth model (1994) is introduced by a two-dimensional vector of state variables \(\mathbf{x}=(x_1, x_2)\), where \(x_1(t)\) denotes the cell concentration of a bacterial population at the time \(t\), and \(x_2(t)\) defines a physiological state of the cells, the process of adjustment (lag-phase):

(1)#\[\begin{split}\begin{alignat}{3} &\dot x_1(t) &&= \frac{x_2(t)}{x_2(t) + 1} \mu^\text{max} \bigg(1 - \frac{x_1(t)}{x_1^\text{max}}\bigg) x(t)\\ &\dot x_2(t) &&= \mu^\text{max} x_2(t) \end{alignat}\end{split}\]

Here \(\mu^\text{max}\) is the maximum growth rate, and \(x_1^\text{max}\) is bacteria concentration at the saturation. To account for the influence of the temperature on the activity of the model, we will use the ‘square root’ or Ratkowsky-type model for the maximum growth rate

(2)#\[\begin{alignat}{3} \sqrt{\mu^\text{max}} = b (T - T_\text{min}), \end{alignat}\]

where \(b\) is the regression coefficient, and \(T_\text{min}\) is the minimum temperature at which the growth can occur. Here \(x_1^\text{max}, b, T_\text{min}\) are parameters that we estimate. And temperature \(T\) is an input of the system.

First of all, import al the necessary libraries.

import numpy as np
from eDPM import *

Define the system of ODEs.

 8def baranyi_roberts_ode(t, x, u, p, ode_args):
 9    # Unpack the input vectors x, u, p for easy access
10    # of their components
11    x1, x2 = x
12    (Temp, ) = u
13    (x_max, b, Temp_min) = p
14
15    # Calculate the maximum growth rate
16    mu_max = b**2 * (Temp - Temp_min)**2
17
18    # Calculate the right hand side of the ODE, store it
19    # in a list [...] and return it.
20    return [
21        mu_max * (x2/(x2 + 1)) * (1 - x1/x_max) * x1,           # f1
22        mu_max * x2                                             # f2
23    ]
24
25def ode_dfdp(t, x, u, p, ode_args):
26    x1, x2 = x
27    (Temp, ) = u
28    (x_max, b, Temp_min) = p
29    mu_max = b**2 * (Temp - Temp_min)**2
30    return [
31        [
32            mu_max * (x2/(x2 + 1)) * (x1/x_max)**2,             # df1/dx_max
33             2 * b * (Temp - Temp_min)**2 * (x2/(x2 + 1))
34                * (1 - x1/x_max)*x1,                            # df1/db
35            -2 * b**2 * (Temp - Temp_min) * (x2/(x2 + 1))
36                * (1 - x1/x_max)*x1                             # df1/dTemp_min
37        ],
38        [
39            0,                                                  # df2/dx_max
40            2 * b * (Temp - Temp_min)**2 * x2,                  # df2/db
41            -2 * b**2 * (Temp - Temp_min) * x2                  # df2/dTemp_min
42        ]
43    ]
44
45def ode_dfdx(t, x, u, p, ode_args):
46    x1, x2 = x
47    (Temp, ) = u
48    (x_max, b, Temp_min) = p
49    mu_max = b**2 * (Temp - Temp_min)**2
50    return [
51        [
52            mu_max * (x2/(x2 + 1)) * (1 - 2*x1/x_max),          # df1/dx1
53            mu_max * 1/(x2 + 1)**2 * (1 - x1/x_max)*x1          # df1/dx2
54        ],
55        [
56            0,                                                  # df2/dx1
57            mu_max                                              # df2/dx2
58        ]
59    ]

Define the parameters of the system p and initial conditions x0.

58    p = (2e4, 0.02, -5.5) # (x_max, b, T_min)
59    x0 = np.array([50., 1.])

Define optimization of 6 time points with lower bound 0.0, upper bound 100.0.

62    times = {"lb": 0., "ub": 100., "n": 6}

Define optimization of one input value (temperature) with lower bound 2.0, upper bound 12.0.

65    inputs = [{"lb": 2.0, "ub": 12.0, "n": 1}]

As an observable, it is pretty common to measure the bacteria count \(x_1\) or the logarithm of this value. For simplicity, we would consider the prior case where the observable is the null-component of the state variable vector \(y(t_i) = x_1(t_i)\).

obs_fun = 0

The resulting Fisher Model:

69    fsm = FisherModel(
70        ode_x0=x0,
71        ode_t0=0.0,
72        ode_fun=baranyi_roberts_ode,
73        ode_dfdx=ode_dfdx,
74        ode_dfdp=ode_dfdp,
75        ode_initial=x0,
76        times=times,
77        inputs=inputs,
78        parameters=p,
79        obs_fun=0,
80        covariance={"abs": 0.3, "rel": 0.1}
81    )

The optimization is then held using relative sensitivities and D-optimality criterion (determinant).

83    fsr = find_optimal(
84        fsm,
85        relative_sensitivities=True,
86        recombination=0.7,
87        mutation=(0.1, 0.8),
88        workers=-1,
89        popsize=10,
90        polish=False,
91    )

Check the structural indentifiability of the Design.

95    check_if_identifiable(fsr)

Save and plot the results of optimization.

95    plot_all_solutions(fsr)
96    json_dump(fsr, "baranyi_roberts.json")

The resulting Optimal Experimental Design:

../../_images/Observable_Results_baranyi_roberts_ode_fisher_determinant_rel_sensit_cont_6times_1temps_000_x_00.svg

The output of the Experimental Design optimization procedure. Line plot: the model solution for the observable, scatter plot: the design time points.#