Notebook source code:
examples/substrate/depletion.ipynb
Run the notebook yourself on binder
Biomass maximization with substrate depletion¶
by S. Psalmon and B. Schall (Polytech Nice Sophia, Applied math. dep)¶
Microorganisms must assign the resources available in their environment to different cellular functions. In nature, it is assumed that bacteria have evolved in a way that their resource allocation strategies maximize their biomass. In our case, the substrate contained in the environment is finite. We define a self-replicator model describing the dynamics of a microbial population growing inside a closed bioreactor of fixed volume. At the beginning of the experience, there is an initial mass of substrate \(S\) inside the bioreactor, that is gradually consumed by the bacterial population. transforming it into precursor metabolites \(P\) building blocks for the creation of proteins and other large molecules. Precursors are transformed into components of the gene expression machinery \(R\) (ribosomes, RNA polymerase etc …), into enzymes that make up the metabolic machinery \(M\), and into a metabolite of interest \(X\) that is excreted from the cell.
While the production of proteins \(M\) and \(R\) is catalyzed by the enzymes \(R\), the absorption of \(S\) and synthesis of X are both catalyzed by \(M\), following these dynamics.
Where \(s\), \(p\), \(r\), \(x\) are the mass fractions of the total bacteria mass of \(S\), \(P\), \(R\) and \(X\), and \(V\) the bacteria’s volume. Here, \(u\) is the allocation control, regulating the production of \(R\) and \(M\). We suppose that compared to \(m\) and \(r\), \(p\) is really small, so we can write \(V = r + m\), that is why \(m\) is not explicit in the system above.
\(w_M\), \(w_R\) and \(w_X\) are function following the Michaelis-Menten kinetics (one of the best-known models of enzyme kinetics), such as :
Where \(K_1,K_2,k_1,k_1\) are constants depending on the bacteria’s environment.
Test of the model with a constant control¶
[1]:
using DifferentialEquations
using Plots
using JuMP, Ipopt
using Statistics
[2]:
# Giordano et al. (2015)
e_M = 3.6 # 1/h
k_R = 3.6 # 1/h
K_R = 1 # g/L
beta = 0.003 # L/g
# Yegorov et al. (2018)
k_X = 1 # 1/h
K_X = 1 # g/L
k_M = 4.32 #1.2*k_R
K_M = 33.33 #0.1/beta
K = beta*K_R
K_1 = beta*K_X
K_2 = beta*K_M
k_1 = k_X/k_R
k_2 = k_M/k_R
t0 = 0.0 # initial time
tf = 40 # final time
p0 = 0.001 # p0
r0 = 0.8 # r0
V0 = 0.003 # V0
s0= 0.3
Eps = s0+(p0+1)*V0
rmin = r0*(V0/Eps)
rmax = 1-(1-r0)*(V0/Eps)
a=0.5
u0 = [s0, p0, r0, V0]
function derivatives!(du, u, p, t)
s, p, r, V = u
du[1] = (-k_2*s/(K_2+s))*(1-r)*V
du[2] = (k_2*s/(K_2+s))*(1-r) - (p/(K+p))*r*(p+1)
du[3] = (a-r)*(p/(K+p))*r
du[4] = (p/(K+p))*r*V
end
tspan = (0.0, 40.0)
prob5 = ODEProblem(derivatives!,u0,tspan)
alg = RK4()
sol5 = DifferentialEquations.solve(prob5, alg, saveat = 0.01, abstol = 1e-9, reltol = 1e-9)
a = 0.6
prob6 = ODEProblem(derivatives!,u0,tspan)
alg = RK4()
sol6 = DifferentialEquations.solve(prob6, alg, saveat = 0.01, abstol = 1e-9, reltol = 1e-9)
a = 0.7
prob7 = ODEProblem(derivatives!,u0,tspan)
alg = RK4()
sol7 = DifferentialEquations.solve(prob7, alg, saveat = 0.01, abstol = 1e-9, reltol = 1e-9);
[3]:
s = plot(ylabel = "s", fmt = :png)
s = plot!(sol5, vars = 1, label = "u = 0.5")
s = plot!(sol6, vars = 1, label = "u = 0.6")
s = plot!(sol7, vars = 1, label = "u = 0.7")
p = plot(ylabel = "p", fmt = :png)
p = plot!(sol5, vars = 2, label = "u = 0.5")
p = plot!(sol6, vars = 2, label = "u = 0.6")
p = plot!(sol7, vars = 2, label = "u = 0.7")
r = plot(ylabel = "r", fmt = :png)
r = plot!(sol5, vars = 3, label = "u = 0.5")
r = plot!(sol6, vars = 3, label = "u = 0.6")
r = plot!(sol7, vars = 3, label = "u = 0.7")
V = plot(ylabel = "V", fmt = :png)
V = plot!(sol5, vars = 4, label = "u = 0.5")
V = plot!(sol6, vars = 4, label = "u = 0.6")
V = plot!(sol7, vars = 4, label = "u = 0.7")
display(plot(s, p, r, V, layout = (2, 2)))
┌ Warning: To maintain consistency with solution indexing, keyword argument vars will be removed in a future version. Please use keyword argument idxs instead.
│ caller = ip:0x0
└ @ Core :-1
Volume maximisation (no metabolite production)¶
[4]:
# JuMP model, Ipopt solver
sys = Model(optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 3))
set_optimizer_attribute(sys, "tol", 1e-9)
set_optimizer_attribute(sys, "max_iter", 100000)
# Parameters
# Giordano et al. (2015)
e_M = 3.6 # 1/h
k_R = 3.6 # 1/h
K_R = 1 # g/L
beta = 0.003 # L/g
# Yegorov et al. (2018)
k_X = 1 # 1/h
K_X = 1 # g/L
k_M = 4.32 #1.2*k_R
K_M = 33.33 #0.1/beta
K = beta*K_R
K_1 = beta*K_X
K_2 = beta*K_M
k_1 = k_X/k_R
k_2 = k_M/k_R
t0 = 0.0 # initial time
tf = 10 # final time
p0 = 0.001 # p0
r0 = 0.8 # r0
V0 = 0.01 # V0
s0=0.1
Eps = s0+(p0+1)*V0
rmin = r0*(V0/Eps)
rmax = 1-(1-r0)*(V0/Eps)
N = 500 # grid size
Δt = (tf-t0)/N # Time step
JuMP.@variables(sys, begin
0. ≤ s[1:N] ≤ 1. # s
0. ≤ p[1:N] ≤ 1. # p
rmin ≤ r[1:N] ≤ rmax # r
0. ≤ V[1:N] ≤ Eps # V
0. ≤ u[1:N] ≤ 1. # allocation (control)
end)
# Objective
@objective(sys, Max, V[N])
# Initial conditions
@constraints(sys, begin
s[1] == s0
p[1] == p0
r[1] == r0
V[1] == V0
end)
# Dynamics, Crank-Nicolson scheme
@NLexpression(sys, wR[j = 1:N], p[j]/(K+p[j]))
@NLexpression(sys, wM[j = 1:N], k_2*s[j]/(K_2+s[j]))
for j in 1:N-1
@NLconstraint(sys, # s' = -wM*(1-r)*V
s[j+1] == s[j] + 0.5 * Δt * (-wM[j]*(1-r[j])*V[j] - wM[j+1]*(1-r[j+1])*V[j+1]))
@NLconstraint(sys, # p' = wM*(1-r) - wR*r*(p+1)
p[j+1] == p[j] + 0.5 * Δt * (wM[j]*(1-r[j]) - wR[j]*r[j]*(p[j]+1) + wM[j+1]*(1-r[j+1]) - wR[j+1]*r[j+1]*(p[j+1]+1)))
@NLconstraint(sys, # r' = (u-r)*wR*r
r[j+1] == r[j] + 0.5 * Δt * ((u[j]-r[j])*wR[j]*r[j] + (u[j+1]-r[j+1])*wR[j+1]*r[j+1]))
@NLconstraint(sys, # v' = wR*r*V
V[j+1] == V[j] + 0.5 * Δt * (wR[j]*r[j]*V[j] + wR[j+1]*r[j+1]*V[j+1]))
end
# Solve for the control and state
println("Solving...")
status = optimize!(sys)
println("Solver status: ", status)
p1 = value.(p)
r1 = value.(r)
s1 = value.(s)
V1 = value.(V)
u1 = value.(u)
println("Cost: ", objective_value(sys))
Solving...
******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
Ipopt is released as open source code under the Eclipse Public License (EPL).
For more information visit https://github.com/coin-or/Ipopt
******************************************************************************
Total number of variables............................: 2500
variables with only lower bounds: 0
variables with lower and upper bounds: 2500
variables with only upper bounds: 0
Total number of equality constraints.................: 2000
Total number of inequality constraints...............: 0
inequality constraints with only lower bounds: 0
inequality constraints with lower and upper bounds: 0
inequality constraints with only upper bounds: 0
Number of Iterations....: 65
(scaled) (unscaled)
Objective...............: -1.0522150088843175e-01 1.0522150088843175e-01
Dual infeasibility......: 1.2925054411097729e-15 1.2925054411097729e-15
Constraint violation....: 5.9158233867151466e-13 5.9158233867151466e-13
Variable bound violation: 0.0000000000000000e+00 0.0000000000000000e+00
Complementarity.........: 9.0909095584056773e-11 9.0909095584056773e-11
Overall NLP error.......: 9.0909095584056773e-11 9.0909095584056773e-11
Number of objective function evaluations = 79
Number of objective gradient evaluations = 66
Number of equality constraint evaluations = 79
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 66
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 65
Total seconds in IPOPT = 5.886
EXIT: Optimal Solution Found.
Solver status: nothing
Cost: 0.10522150088843175
[5]:
t = (1:N) * Δt
s_plot = plot(t, s1, xlabel = "t", ylabel = "s", legend = false, fmt = :png)
p_plot = plot(t, p1, xlabel = "t", ylabel = "p", legend = false, fmt = :png)
r_plot = plot(t, r1, xlabel = "t", ylabel = "r", legend = false, fmt = :png)
V_plot = plot(t, V1, xlabel = "t", ylabel = "V", legend = false, fmt = :png)
display(plot(s_plot, p_plot, r_plot, V_plot, layout = (2, 2)))
u_plot = plot(t, u1, xlabel = "t", ylabel = "u", fmt = :png)
[5]:
Proportions as a function of time (1/2)¶
[6]:
y2 = ones(length(s1))
y3 = zeros(length(s1))
plot(t,[s1/Eps (r1.*V1)/Eps+s1/Eps V1/Eps+s1/Eps y2 ] ,xlabel = "t", ylabel = "Proportion of Sigma", fillrange = [y3 s1/Eps s1/Eps+(r1.*V1)/Eps V1/Eps+s1/Eps], fillalpha = 0.6, c = [1 7 3 4] , label = ["s" "rv" "mV" "Vp"], legend = :topleft, fmt = :png)
[6]:
Proportions as a function of time (2/2)¶
[7]:
@gif for i ∈ 1:(length(t))
pie(["s", "rV", "mV","Vp"], [(s1/Eps)[i], ((r1.*V1)/Eps)[i], ((V1 - (r1.*V1))/Eps)[i],(V1.*p1/Eps)[i]], title = "")
end every 2
┌ Info: Saved animation to
│ fn = /Users/ocots/Boulot/recherche/Logiciels/dev/ct/gallery/examples/substrate/tmp.gif
└ @ Plots /Users/ocots/.julia/packages/Plots/gl4j3/src/animation.jl:137
[7]: