In Kinbiont is possible to simulate and fit the bacterial growth with any Ordinary Differential Equations System. We can broadly divide the classes of possible mathematical models into the following:
Pages = ["index.md"]
Depth = 2
In this section we show some of the harcoded models of Kinbiont.jl but note that you can input any custom model both has analytic functio both as ODE.
In general, NL fitting is preferred to ODE fitting in the following cases:
- Since it is faster when you have to analyze a large dataset
- When you do not trust initial conditions (e.g., initial inoculum under the detection limit of the instrument). ODE fit needs to fix the initial condition of data on the first (or an average of the first) time point of data, and this could lead to errors.
In this case, we are supposed to know the analytic formula of microbial growth; in particular, we have implemented some models from "Statistical evaluation of mathematical models for microbial growth" and added some piecewise models. They are:
-
Exponential
$$N(t) = N_0 \cdot e^{\mu \cdot t}$$ where
$\mu$ is the growth rate, and$N_0$ is the starting condition. -
Gompertz
$$N(t) = N_{\text{max}} \cdot e^{-e^{-\mu \cdot (t - t_{\text{L}})}}$$ where
$\mu$ is the growth rate,$N_{\text{max}}$ is the total growth, and$t_{\text{L}}$ is the lag time. -
Logistic
$$N(t) = \frac{N_{\text{max}}}{1 + \left( \frac{N_{\text{max}}}{N_0} - 1 \right) \exp\left( - \mu \cdot t \right)}$$ where
$\mu$ is the growth rate,$N_0$ is the starting condition, and$N_{\text{max}}$ is the total growth. -
Richards model
$$N(t) = \frac{N_{\text{max}}}{[1 + \nu \cdot e^{-\mu \cdot (t - t^{\text{L}})}]^{\frac{1}{\nu}}}$$ where
$\mu$ is the growth rate,$N_{\text{max}}$ is the total growth,$t_{\text{L}}$ is the lag time, and$\nu$ is a shape constant. -
Weibull
$$N(t) = N_{\text{max}} - (N_{\text{max}} - N_0) \cdot e^{-(\mu \cdot t)^{\nu}}$$ where
$\mu$ is the growth rate,$N_0$ is the starting condition,$N_{\text{max}}$ is the total growth, and$\nu$ is a shape constant. -
Morgan
$$N(t) = \frac{N_0 \cdot K^{\nu} + N_{\text{max}} \cdot t^{\nu}}{K^{\nu} + t^{\nu}}$$ where
$N_0$ is the starting condition,$N_{\text{max}}$ is the total growth, and$K$ and$\nu$ are shape constants. -
Bertalanffy
$$N(t) = N_0 + (N_{\text{max}} - N_0) \cdot (1 - e^{-\mu \cdot t})^{\frac{1}{\nu}}$$ where
$N_0$ is the starting condition,$N_{\text{max}}$ is the total growth,$\mu$ is the growth rate, and$\nu$ is a shape constant. -
Piece-Wise Linear-Logistic
$$N(t) = \begin{cases} N_0, & t < t_{\text{L}} \ N(t) = \frac{N_{\text{max}}}{1 + \left( \frac{N_{\text{max}}}{N_0} - 1 \right) \exp\left( - \mu \cdot (t - t_{\text{L}}) \right)}, & t_{\text{L}} \leq t \end{cases}$$
where
$N_0$ is the starting condition,$N_{\text{max}}$ is the total growth,$\mu$ is the growth rate, and$t_{\text{L}}$ is the lag time. -
Piece-wise Exponential-Logistic
$$N(t) = \begin{cases} N_0 \exp{(\mu_0 \cdot t)}, & t < t_{\text{L}} \ \frac{N_{\text{max}}}{1 + \left( \frac{N_{\text{max}}}{N_0 \exp{(\mu_0 \cdot t_{\text{L}})}} - 1 \right) \exp\left( - \mu \cdot (t - t_{\text{L}}) \right)}, & t_{\text{L}} \leq t \end{cases}$$
where
$N_0$ is the starting condition,$N_{\text{max}}$ is the total growth,$\mu$ is the growth rate,$t_{\text{L}}$ is the lag time, and$\mu_0$ is the growth rate during the lag phase.
To call these models use the string present in this table, the parameters will be returned in the same order of this table.
| Model Name | Parameters List | String to call |
|---|---|---|
| Exponential | "NL_exponential" |
|
| Gompertz | "NL_Gompertz" |
|
| Logistic | "NL_logistic" |
|
| Richards model | "NL_Richards" |
|
| Weibull | "NL_Weibull" |
|
| Morgan | "NL_Morgan" |
|
| Bertalanffy | "NL_Bertalanffy" |
|
| piece-wise linear-logistic | "NL_piecewise_lin_logistic" |
|
| piece-wise exponential-logistic | "NL_piecewise_exp_logistic" |
For a general idea of the properties of models, consult the following table:
| Model Name | Has Lag? | Is Piecewise? | Has Stationary Phase? |
|---|---|---|---|
| Exponential | No | No | No |
| Gompertz | Yes | No | Yes |
| Logistic | No | No | Yes |
| Richards model | Yes | No | Yes |
| Weibull | No | No | Yes |
| Morgan | No | No | Yes |
| Bertalanffy | No | No | Yes |
| Piece-wise linear-logistic | Yes | Yes | Yes |
| Piece-wise exponential-logistic | Yes | Yes | Yes |
If undecided between different models, please use the model selection function.
The models implemented in Kinbiont are the following:
- Exponential:
- Hyper Gompertz:
- Hyper Logistic:
-
Bertalanffy-Richards:
$$\frac{d N(t)}{dt} = \frac{\mu}{N_\text{max}^n } \cdot \left ( N_\text{max}^n - N^n(t) \right)$$
where
- Logistic:
- Adjusted Logistic:
- Gompertz:
- Baranyi Richards:
- Baranyi Roberts:
- Piece-wise Adjusted Logistic:
$$\frac{d N(t)}{dt} =
\begin{cases}
\text{const.} , N(t) & t < t_{\text{L}} \
\mu \left( 1 - \left( \frac{N(t)}{N_{\text{max}}} \right)^m \right) , N(t) & t \geq t_{\text{L}}
\end{cases}$$
where
- Triple Piece-wise Adjusted Logistic:
where
- Triple Piece-wise:
where
- Triple Piece-wise Exponential:
where
- Four Piece-wise Exponential:
where
- Heterogeneous Population Model (HPM): $$\begin{cases} N(t) = N_1(t) + N_2(t), \ \frac{d N_1(t)}{dt} = - r_{\text{L}} \cdot N_1(t), \ \frac{d N_2(t)}{dt} = r_{\text{L}} \cdot N_1(t) + \mu \cdot N_2(t) \cdot \left(1 - \frac{N_1(t) + N_2(t)}{N_{\text{max}}}\right), \end{cases}$$
where
Note that these models assume that the cells are in two states:
- Exponential Heterogeneous Population Model:
$$\begin{cases} N(t) = N_1(t) + N_2(t) \ \frac{d N_1(t)}{dt} = - \text{r}{\text{L}} , N_1(t) \ \frac{d N_2(t)}{dt} = \text{r}{\text{L}} , N_1(t) + \mu , N_2(t) \end{cases}$$
where similarly to the HPM model,
- Adjusted Heterogeneous Population Model:
where
- Heterogeneous Population Model with Inhibition:
where
Note that these models assume that the cells are in two states:
- Heterogeneous Population Model with Inhibition and Death:
where
Note that these models assume that the cells are in three states:
- Heterogeneous Population Model with Inhibition, Death and Resistance:
where
To call these models use the string present in this table, the parameters will be returned in the same order of this table.
| Model Name | Parameters | String to call |
|---|---|---|
| Exponential ODE | label_exp, well, model, gr, th_max_gr, emp_max_gr, loss |
"exponential" |
| Hyper Gompertz | label_exp, well, model, gr, N_max, shape, th_max_gr, emp_max_gr, loss |
"hyper_gompertz" |
| Hyper Logistic | label_exp, well, model, doubling_time, gr, N_max, shape, th_max_gr, emp_max_gr, loss |
"hyper_logistic" |
| Von Bertalanffy ODE | label_exp, well, model, alpha, beta, a, b, th_max_gr, emp_max_gr, loss |
"ode_von_bertalanffy" |
| Bertalanffy-Richards | label_exp, well, model, gr, N_max, shape, th_max_gr, emp_max_gr, loss |
"bertalanffy_richards" |
| Logistic | label_exp, well, model, gr, N_max, th_max_gr, emp_max_gr, loss |
"logistic" |
| Adjusted Logistic | label_exp, well, model, gr, N_max, shape, th_max_gr, emp_max_gr, loss |
"alogistic" |
| Gompertz | label_exp, well, model, gr, N_max, th_max_gr, emp_max_gr, loss |
"gompertz" |
| Baranyi Richards | label_exp, well, model, gr, N_max, lag_time, shape, th_max_gr, emp_max_gr, loss |
"baranyi_richards" |
| Baranyi Roberts | label_exp, well, model, gr, N_max, lag_time, shape_1, shape_2, th_max_gr, emp_max_gr, loss |
"baranyi_roberts" |
| Piece-wise Adjusted Logistic | label_exp, well, model, gr, N_max, lag, shape, linear_const, th_max_gr, emp_max_gr, loss |
"piecewise_adjusted_logistic" |
| Triple Piece-wise Adjusted Logistic | label_exp, well, model, gr, N_max, lag, shape, linear_const, t_stationary, linear_lag, th_max_gr, emp_max_gr, loss |
"triple_piecewise_adjusted_logistic" |
| Triple Piece-wise | label_exp, well, model, gr, gr_2, gr_3, lag, t_stationary, th_max_gr, emp_max_gr, loss |
"ODE_triple_piecewise" |
| Triple Piece-wise Exponential | label_exp, well, model, gr, gr_2, gr_3, lag, t_stationary, th_max_gr, emp_max_gr, loss |
"ODE_triple_piecewise_exponential" |
| Four Piece-wise Exponential | label_exp, well, model, gr, gr_2, gr_3, gr_4, lag, t_decay_gr, t_stationary, th_max_gr, emp_max_gr, loss |
"ODE_four_piecewise" |
| Diauxic Piecewise Adjusted Logistic | label_exp, well, model, gr_1, N_max, shape_1, lag, linear_const, t_shift, gr_2, N_max_2, shape_2, end_second_lag, lag_2_gr, th_max_gr, emp_max_gr, loss |
"Diauxic_piecewise_adjusted_logistic" |
| Heterogeneous Population Model (HPM McKellar) | label_exp, well, model, gr, exit_lag_rate, N_max, th_max_gr, emp_max_gr, loss |
"HPM" |
| Exponential Heterogeneous Population Model (HPM McKellar) | label_exp, well, model, gr, exit_lag_rate, th_max_gr, emp_max_gr, loss |
"HPM_exp" |
| Adjusted Heterogeneous Population Model | label_exp, well, model, gr, exit_lag_rate, N_max, shape, th_max_gr, emp_max_gr, loss |
"aHPM" |
| Heterogeneous Population Model with Inhibition | label_exp, well, model, gr, exit_lag_rate, inactivation_rate, th_max_gr, emp_max_gr, loss |
"HPM_3_inhibition" |
| Heterogeneous Population Model with Inhibition and Death | label_exp, well, model, gr, exit_lag_rate, inactivation_rate, death_rate, th_max_gr, emp_max_gr, loss |
"HPM_3_death" |
| Heterogeneous Population Model with Inhibition, Death and Resistance | label_exp, well, model, gr, exit_lag_rate, inactivation_rate, death_rate, n_res, shape, th_max_gr, emp_max_gr, loss |
"aHPM_3_death_resistance" |
In the following table, you can find a general description of the properties of the hardcoded ODE models of Kinbiont:
| Model Name | Has Lag? | Is Piecewise? | Has Stationary Phase? | Has Inhibition? | Is Monotonic? | Supposes Multiple States? |
|---|---|---|---|---|---|---|
| Exponential ODE | No | No | No | No | Yes | No |
| Hyper Gompertz | No | No | Yes | No | Yes | No |
| Hyper Logistic | No | No | Yes | No | Yes | No |
| Von Bertalanffy ODE | No | No | Yes | No | Yes | No |
| Bertalanffy-Richards | No | No | Yes | No | Yes | No |
| Logistic | No | No | Yes | No | Yes | No |
| Adjusted Logistic | No | No | Yes | No | Yes | No |
| Gompertz | No | No | Yes | No | Yes | No |
| Baranyi Richards | Yes | No | Yes | No | Yes | No |
| Baranyi Roberts | Yes | No | Yes | No | Yes | No |
| Piece-wise Adjusted Logistic | Yes | Yes | Yes | No | No | No |
| Triple Piece-wise Adjusted Logistic | Yes | Yes | Yes | No | No | No |
| Triple Piece-wise | Yes | Yes | Yes | No | No | No |
| Triple Piece-wise Exponential | Yes | Yes | Yes | No | No | No |
| Four Piece-wise Exponential | Yes | Yes | Yes | No | No | No |
| Diauxic Piecewise Adjusted Logistic | Yes | Yes | Yes | No | No | No |
| Heterogeneous Population Model (HPM McKellar) | Yes | No | Yes | No | Yes | Yes |
| Exponential Heterogeneous Population Model (HPM McKellar) | Yes | No | No | No | Yes | Yes |
| Adjusted Heterogeneous Population Model | Yes | No | Yes | No | Yes | Yes |
| Heterogeneous Population Model with Inhibition | Yes | No | Yes | Yes | Yes | Yes |
| Heterogeneous Population Model with Inhibition and Death | Yes | No | Yes | Yes | No | Yes |
| Heterogeneous Population Model with Inhibition, Death and Resistance | Yes | No | Yes | Yes | No | Yes |
If undecided between different models, the model selection function should be used.
In the stochastic version of the growth models, the growth rate of each population component (denoted as
- Monod Model:
- Haldane Model:
- Blackman Model:
- Tessier Model:
- Moser Model:
- Aiba-Edwards Model:
- Verhulst Model:
In this section, we present some examples of ODEs multidimensional system hardcoded in Kinbiont. Note that these are just examples since you can define custom models:
-
SIR Model (Susceptible-Infected-Recovered)
$$\begin{cases} \frac{dS}{dt} = -\beta S I \ \frac{dI}{dt} = \beta S I - \gamma I \ \frac{dR}{dt} = \gamma I \end{cases}$$ Parameters: Infection rate ($\beta$ ), Recovery rate ($\gamma$ ). -
SIR with Birth and Death (SIR_BD)
$$\begin{cases} \frac{dS}{dt} = -\beta S I + b S - d S \ \frac{dI}{dt} = \beta S I - \gamma I - d I \ \frac{dR}{dt} = \gamma I - d R \end{cases}$$ Parameters: Infection rate ($\beta$ ), Recovery rate ($\gamma$ ), Birth rate ($b$ ), Death rate ($d$ ). -
SIS Model (Susceptible-Infected-Susceptible)
$$\begin{cases} \frac{dS}{dt} = -\beta S I + \gamma I \ \frac{dI}{dt} = \beta S I - \gamma I \end{cases}$$ Parameters: Infection rate ($\beta$ ), Recovery rate ($\gamma$ ). -
Lotka-Volterra Predator-Prey Model
$$\begin{cases} \frac{dP}{dt} = \alpha P - \beta P C \ \frac{dC}{dt} = -\delta C + \gamma C P \end{cases}$$ Parameters: Prey birth rate ($\alpha$ ), Predation rate ($\beta$ ), Predator death rate ($\delta$ ), Predator reproduction rate ($\gamma$ ). -
Lotka-Volterra with Substrate Limitation
$$\begin{cases} \frac{dP}{dt} = \alpha P \frac{S}{S + K} - \beta P C \ \frac{dC}{dt} = -\delta C + \gamma C P \ \frac{dS}{dt} = -\alpha P \frac{S}{S + K} \end{cases}$$ Parameters: Growth rate ($\alpha$ ), Half-saturation ($K$ ), Predation rate ($\beta$ ), Predator mortality ($\delta$ ), Predator efficiency ($\gamma$ ). -
Monod Chemostat Model (Microbial Growth in a Chemostat)
$$\begin{cases} \frac{dX}{dt} = \mu X - D X \ \frac{dS}{dt} = D (S_{\text{in}} - S) - \frac{\mu X}{Y} - m X \end{cases}$$ where
$$\mu = \mu_m \frac{S}{K_s + S}$$ Parameters: Substrate affinity ($K_s$ ), Maintenance coefficient ($m$ ), Yield coefficient ($Y$ ), Max growth rate ($\mu_m$ ), Dilution rate ($D$ ), Substrate inflow ($S_{\text{in}}$ ). -
Droop Model (Nutrient Quota Model)
$$\begin{cases} \frac{dX}{dt} = \mu X - D X \ \frac{dS}{dt} = \rho X - D S + D S_{\text{in}} \ \frac{dQ}{dt} = \rho - \mu Q \end{cases}$$ where
$$\mu = \mu_m \left(1 - \frac{Q_0}{Q}\right)$$ and
$$\rho = \rho_m \frac{S}{K_s + S}$$ Parameters: Growth rate ($\mu_m$ ), Nutrient uptake rate ($\rho_m$ ), Half-saturation ($K_s$ ), Dilution rate ($D$ ), Minimum quota ($Q_0$ ), Substrate inflow ($S_{\text{in}}$ ). -
Synthetic Chemostat Model (Including Biological Inertia)
$$\begin{cases} \frac{dx}{dt} = Y q_s - a_0 r x - D x \ \frac{ds}{dt} = D (s_r - s) - q_s x \ \frac{dr}{dt} = (Y q_s - a_0 r) \left(\frac{s}{K_r + s} - r\right) \end{cases}$$ where
$$q_s = r \frac{Q_s K_s}{K_s + s} + (1 - r) \frac{Q_s' K_s'}{K_s' + s}$$ Parameters: Yield ($Y$ ), Biological inertia ($a_0$ ), Dilution rate ($D$ ), Nutrient uptake coefficients ($Q_s, Q_s'$ ), Saturation constants ($K_s, K_s'$ ), Half-saturation constant for$r$ ($K_r$ ). -
Monod-Ierusalimsky This model describes microbial growth, substrate consumption, and product formation using Monod-Ierusalimsky kinetics.
The specific growth rate
- The first fraction represents substrate-limited growth (Monod equation).
- The second fraction accounts for product inhibition (Ierusalimsky modification).
The effective biomass yield is given by:
Finally, the System Dynamics is specified by:
where
where
The system state variables are:
- Biomass concentration:
$x = u_1$ - Substrate concentration:
$s = u_2$ - Product concentration:
$p = u_3$
The system dynamics are governed by the following parameters:
| Parameter | Description |
|---|---|
| Substrate affinity constant | |
| Product inhibition constant | |
| Maximum specific growth rate | |
| Maximum yield coefficient | |
| Product yield coefficient | |
| Dilution rate | |
| Inflow substrate concentration | |
| Maintenance coefficient |
| Model Name | Parameters List | String to Call |
|---|---|---|
| SIR Model | SIR |
|
| SIR with Birth/Death | SIR_BD |
|
| SIS Model | SIS |
|
| Lotka-Volterra | Lotka_Volterra |
|
| Lotka-Volterra with Substrate | Lotka_Volterra_with_substrate |
|
| Monod Chemostat | Monod_Chemostat |
|
| Droop Model | Droop |
|
| Synthetic Chemostat | Synthetic_Chemostat |
|
| Monod-Ierusalimsky | Monod_Ierusalimsky |
In this case, Kinbiont does not take as input an ODE function but a data structure, and it proceeds to construct the ODE system. Note that the same results in theory could be obtained by declaring a custom ODEs system that behaves in the same way. The data struct is composed in the following way:
model = Kinbiont_Cybernetic_Model(
Bio_mass_conc = 1.0, # Initial biomass concentration
Substrate_concentrations = [3.0, 3.0], # Concentrations of 2 substrates
Protein_concentrations = [0.0, 0.0], # Initial protein concentrations
allocation_rule = threshold_switching_rule, # Dynamic resource allocation rule
reaction = nothing, # No specific reaction function provided
cost = nothing, # No cost function
protein_thresholds = 0.01, # Protein activation threshold
a = [0.1, 0.4], # Synthesis rate constants for proteins
b = [0.00001, 0.000001], # Degradation constants for proteins
V_S = [0.7, 0.1], # Substrate utilization rates
k_S = [0.1, 0.11], # Saturation constants for substrates
Y_S = [0.07, 0.11] # Yield coefficients for biomass per substrate
)Then the ODEs system that will be constructed is the following:
-
Substrate Consumption:
$$\frac{dS_i}{dt} = -\frac{V_{S_i} P_i S_i}{k_{S_i} + S_i} \cdot u_1, \quad \forall i \in {1, \dots, n}$$ -
Protein Synthesis:
$$\frac{dP_i}{dt} = a_i \cdot \text{alloc}i \cdot k{S_i} - b_i P_i u_1, \quad \forall i \in {1, \dots, n}$$ -
Biomass Growth:
$$\frac{dN}{dt} = -\sum_{i=1}^{n} Y_{S_i} \frac{dS_i}{dt} \cdot N$$
Where
The allocation_rule in the cybernetic model defines how resources are allocated to different substrates, proteins, or other components of the system. You can create custom rules based on various factors such as substrate concentration, protein levels, or cost-benefit analysis. Below are some examples of how to define a custom allocation rule:
# Function for proportional allocation based on substrate concentrations
function proportional_allocation_rule(a, b, V_S, k_S, Y_S, P, S, cost, protein_thresholds)
# Normalize substrate concentrations to create allocation vector
alloc = S ./ sum(S)
return alloc
endThis rule allocates resources to substrates proportionally based on their concentrations.
For the moment all substrate follow a Monod-like concentration effect on growth rate, so please leave the option reaction = nothing.
In the case the user want to input the system as a network of reaction, Kinbiont.jl relies on Catalyst.jl to generate the ODE problem to be fitted and solved. A network and its parameters should be declared as the following:
u0 = [:S => 301, :E => 100, :SE => 0, :P => 0]
ps = [:kB => 0.00166, :kD => 0.0001, :kP => 0.1]
model_Michaelis_Menten = @reaction_network begin
kB, S + E --> SE
kD, SE --> S + E
kP, SE --> P + E
end For other examples on how declare a reaction network please consult:Catalyst.jl documentation.