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):
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
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:
The output of the Experimental Design optimization procedure. Line plot: the model solution for the observable, scatter plot: the design time points.#