In [1]:
import pandas as pd
import numpy as np
import scipy
import seaborn as sns
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from scipy.interpolate import interp1d
from scipy import interpolate
from scipy.optimize import minimize
from scipy import integrate
from matplotlib.pylab import *
from scipy import stats
import lmfit
from scipy import optimize as op
from matplotlib import gridspec
from SALib.sample import saltelli
from SALib.analyze import sobol
import corner
import emcee
from matplotlib.ticker import MaxNLocator
from matplotlib import rcParams
rcParams["savefig.dpi"] = 100
rcParams["figure.dpi"] = 100
rcParams["font.size"] = 14
import matplotlib
matplotlib.rcParams.update({'errorbar.capsize': 3})
import warnings
warnings.filterwarnings("ignore")
%matplotlib inline
Parameter estimation¶
In [2]:
def dy_dt(t, y, gamma, delta, a):
"""
y[0] = dS/dt
y[1] = dx_growing/dt
y[2] = dx_nongrowing/dt
INPUTS:
t: time
y: the time varying concentrations: substrate (S) and biomass (x1 and x2)
"""
F = D*V # L/h
# Change these time values, depending on the question
if t >= 75 and t <= 85:
S_in = 30 ## mg/L
else:
S_in = 30 ## mg/L
S = y[0]
x1 = y[1]
x2 = y[2]
y = np.zeros((3,1))
K = (1-((F/V)/mu_max))*gamma*x1-((F/V)/mu_max)*delta*x2
y[0] = F/V*(S_in - S) - (mu_max*S/(Ks+S))/Y_g*x1-m*x1-m2*x2-((F/V)/mu_max)*delta*x2*a
y[1] = -x1*F/V +(mu_max*S/(Ks+S)) * x1 -K
y[2] = -x2*F/V +K
return y
def residual(p):
"""
Calculation of residual between
model prediction and measurement
"""
v = p.valuesdict()
p_guess = [v['gamma'], v['delta'], v['a']]
y_pred = system_opt(p_guess)
return y_pred - y_true
def system_opt(p_guess):
"""
System definition and optimisation for substrate, growing and non-growing biomass.
y[0] = dS/dt
y[1] = dx_growing/dt
y[2] = dx_nongrowing/dt
INPUTS:
t: time
y: the time varying concentrations: substrate (S) and biomass (x1 and x2)
"""
gamma = p_guess[0]
delta = p_guess[1]
a = p_guess[2]
r = integrate.ode(dy_dt).set_integrator('vode', method='bdf')
# Initial guess for gamma, delta, a parameters
ICs = [0.087, 0.43, 0.05]
t_0 = 0.0
r.set_initial_value(ICs, t_0).set_f_params(gamma, delta, a, )
t_final = 1520.0
dt = 0.2
# Create vectors to store the solutions in;
# add extra space for intial condition
n_steps = np.int((t_final - t_0)/dt) + 1
time = np.zeros(n_steps)
S = np.zeros(n_steps)
x1 = np.zeros(n_steps)
x2 = np.zeros(n_steps)
S[0], x1[0], x2[0] = ICs
k = 1
while r.successful() and r.t < t_final:
r.integrate(r.t + dt)
time[k] = r.t
S[k] = r.y[0]
x1[k] = r.y[1]
x2[k] = r.y[2]
k += 1
t_sim = np.array([0, 12, 17, 18])*24
s_sim = interpolate.PchipInterpolator(time,S, extrapolate=True)(t_sim)
x1_sim = interpolate.PchipInterpolator(time,x1, extrapolate=True)(t_sim)
x2_sim = interpolate.PchipInterpolator(time,x2, extrapolate=True)(t_sim)
arr = np.concatenate((s_sim, x1_sim, x2_sim))
return arr
data = pd.read_excel("dynamics.xlsx", sheet_name='Sheet2')
y_true = data.iloc[:12][['atrazine', 'Growing', 'Non-growing']].dropna().T.values.flatten()*1e-3
In [3]:
V = 2 # L, Reactor volumne
Y_g = 0.028 # g/g, Biomass yield
mu_max = 0.11 # /h, specific growth rate
Ks = 0.237 # mg, Monod affinity
m=0.22 # g/g-h Maintenance demand of growing biomass
m2=0.10 # g/g-h Maintenance demand of non-growing biomass
D = 0.023 # /h , Dilution rate
# Finding global minima using brute-force method
p = lmfit.Parameters()
p.add_many(('gamma', 0.008, True, 0.004, 0.01, None, 0.001),
('delta', 0.004, True, 0.001, 0.008, None, 0.001),
('a', 0.02, True, 0.01, 0.05, None, 0.005))
mi = lmfit.minimize(residual, p, method='brute', nan_policy='omit')
# Polish using local minimizer
mini = lmfit.Minimizer(residual, mi.params, nan_policy='omit')
# first solve with Nelder-Mead
out1 = mini.minimize(method='Nelder')
# then solve with Levenberg-Marquardt using the
# Nelder-Mead solution as a starting point
out2 = mini.minimize(method='leastsq', params=out1.params)
Model simulation using estimated parameters for all conditions¶
In [4]:
def dy_dt(t, y):
"""
y1 = dS/dt
y2 = dx/dt
INPUTS:
t: time
y: the time varying concentrations: substrate (S) and biomass (x)
"""
if t <= 19*24:
F = 0.023*V # ml/hr
elif t >19*24 and t <= 30*24:
F = 0.032*V # ml/hr
elif t>30*24 and t <= 36*24:
F = 0.048*V # ml/hr
elif t>36*24 and t <= 41*24:
F = 0.056*V # ml/hr
else:
F = 0.068*V # ml/hr
S = y[0]
x1 = y[1]
x2 = y[2]
y = np.zeros((3,1))
gamma, delta, a = [ out2.params[i].value for i in out2.params.keys()] #Estimated parameters
K = (1-F/V/mu_max)*gamma*x1-F/V/mu_max*delta*x2
y[0] = F/V*(S_in - S) - mu_max*S/(Ks+S)*(1/Y_g)*x1-m*x1-m2*x2-((F/V)/mu_max)*delta*x2*a
y[1] = -x1*F/V +(mu_max*S/(Ks+S)) * x1 -K
y[2] = -x2*F/V +K
return y
V = 2 # L, Reactor volumne
Y_g = 0.028 # g/g, Biomass yield
mu_max = 0.11 # /h, specific growth rate
Ks = 0.237 # mg, Monod affinity
m=0.22 # g/g-h Maintenance demand of growing biomass
m2=0.10 # g/g-h Maintenance demand of non-growing biomass
D = 0.023 # /h , Dilution rate
S_in = 30 # mg/L
r = integrate.ode(dy_dt).set_integrator('vode', method='bdf')
# Part 1:
ICs = [0.036, 0.43, 0.01]
# Part 2: use the steady-state solution
# ICs = [33.33333333333, 53.3333333333]
t_0 = 0.0
r.set_initial_value(ICs, t_0)
t_final = 1200.0
dt = 2
# Create vectors to store the solutions in;
# add extra space for intial condition
n_steps = np.int((t_final - t_0)/dt) + 1
time = np.zeros(n_steps)
S = np.zeros(n_steps)
x1 = np.zeros(n_steps)
x2 = np.zeros(n_steps)
S[0], x1[0], x2[0] = ICs
k = 1
while r.successful() and r.t < t_final:
r.integrate(r.t + dt)
time[k] = r.t
S[k] = r.y[0]
x1[k] = r.y[1]
x2[k] = r.y[2]
k += 1
for sh in ['Sheet2', 'Sheet3']:
fig, ax = plt.subplots(1,1, figsize=(14,5))
ax.plot(time, S*1e3, 'r', label='Substrate', lw=3)
ax.plot(time, x1*1e3, 'k', label='Growing Biomass', lw=3)
ax.plot(time, x2*1e3, 'b', label='Non-growing Biomass', lw=3)
ax.plot(time, (x1+x2)*1e3, 'y', label='Total biomass', lw=3)
data = pd.read_excel("dynamics.xlsx", sheet_name=sh)
# 'time / d', 'atrazine', 'Biomass'
data['time / d']=data['time / d']*24
data['Biomass']=data['Biomass']*1000
data['Active']=data['Growing']*1000
data['Inactive']=data['Non-growing']*1000
# data = data[data['time / d']> 151*24]
data['time / d'] = data['time / d']-data['time / d'].iloc[0]
# data['time / d']
# data.index= data['time / d']
data[['time / d','atrazine', 'Biomass', 'Active', 'Inactive']].plot.scatter(x= 'time / d',y='atrazine',marker='s', c='r', ax=ax, s=40)
data[['time / d','atrazine', 'Biomass','Active', 'Inactive']].plot.scatter(x= 'time / d',y='Biomass', c='y', ax=ax, s=40)
data[['time / d','atrazine', 'Biomass','Active', 'Inactive']].plot.scatter(x= 'time / d',y='Active', c='k', marker='d', ax=ax, s=40)
data[['time / d','atrazine', 'Biomass','Active', 'Inactive']].plot.scatter(x= 'time / d',y='Inactive', c='b', marker='d', ax=ax, s=40)
ax.set_ylabel(r'Concentration [$\mu$g L$^{-1}$]')
ax.set_xlabel(r'Time [h]')
x=np.array([19, 30, 36, 41])*24
ax.vlines(x=x, ymin= 0, ymax=1070, colors='#d3d3d3', linestyles='solid')
ax.legend(loc=2, ncol=4)
ax.set_xlim(0,t_final);
ax.set_ylim(0,1070)
fig.savefig("dynamics_+str(sh).svg");
Parameter uncertainity estimation by MCMC¶
In [6]:
def Substrate_min1(alpha, ms1, ms2, qs_max, Ks, Yg, mu_max, a, delta):
"""
Calculation of minimum substrate reqirement
"""
mu = 0.023 # /h s
A=mu*alpha/Yg+ms1*alpha+ms2*(1-alpha)+(mu/mu_max)*(1-alpha)*a*delta
Cs_min = Ks*A/(qs_max-A)
return Cs_min
fs=16 ## fontsize for plotting
np.random.seed(123)
# Choose the "true" parameters.
alpha= 0.7
ms1 = 0.22
ms2 = 0.10
qs_max = 4.08
Ks = 237e-6
Yg = 0.028
mu_max = 0.11
a = 0.028
delta= 0.002
N = 50
f_true = 90e-6
alpha1 = np.sort(1*np.random.rand(N))
# Sminerr = 0.001+0.05*np.random.rand(N)
Sminerr = 1.0e-7+0.05e-6*np.random.rand(N)
Smin = Substrate_min1(alpha1, ms1, ms2, qs_max, Ks, Yg, mu_max, a, delta)
Smin += np.abs(f_true*Smin) * np.random.randn(N)
Smin += Sminerr * np.random.randn(N)
# Do the least-squares fit and compute the uncertainties.
A = np.vstack((np.ones_like(alpha1), alpha1)).T
C = np.diag(Sminerr * Sminerr)
#cov = np.linalg.inv(np.dot(A.T, np.linalg.solve(C, A)))
popt, cov = curve_fit(Substrate_min1, alpha1, Smin, maxfev = 500000)
ms1_ls, ms2_ls, qs_max_ls, Ks_ls, Yg_ls, mu_max_ls, a_ls, delta_ls = popt
# Define the probability function as likelihood * prior.
def lnprior(theta):
ms1, ms2, qs_max, Ks, Yg, mu_max, a, delta, lnf = theta
if 0. < ms1 < 0.3 and 0 < ms2 < 0.1 and 1.0 < qs_max < 4.0 and 1.0e-4 < Ks < 3.8e-4 \
and 0.02 < Yg < 0.1 and 0.05 < mu_max < 0.2 and 0.02 < a < 0.1 and 0.001 < delta < 0.01 and -10.0 < lnf < 10.0:
return 0.0
return -np.inf
def lnlike(theta, x, Smin, yerr):
ms1, ms2, qs_max, Ks, Yg, mu_max, a, delta, lnf = theta
# model = m * x + b
model = Substrate_min1(alpha1, ms1, ms2, qs_max, Ks, Yg, mu_max, a, delta)
inv_sigma2 = 1.0/(Sminerr**2 + model**2*np.exp(2*lnf))
return -0.5*(np.sum((Smin-model)**2*inv_sigma2 - np.log(inv_sigma2)))
def lnprob(theta, x, y, yerr):
lp = lnprior(theta)
if not np.isfinite(lp):
return -np.inf
return lp + lnlike(theta, x, y, yerr)
ms1_ls, ms2_ls, qs_max_ls, Ks_ls, Yg_ls, mu_max_ls, a_ls, delta_ls = popt
# Find the maximum likelihood value.
chi2 = lambda *args: -2 * lnlike(*args)
result = op.minimize(chi2, [ms1, ms2, qs_max, Ks, Yg, mu_max, a, delta, np.log(f_true)], args=(alpha1, Smin, Sminerr))
ms1_ml, ms2_ml, qs_max_ml, Ks_ml, Yg_ml, mu_max_ml, a_ml, delta_ml, lnf_ml = result["x"]
print("""Maximum likelihood result:
ms1 = {0} (truth: {1})
ms2 = {2} (truth: {3})
qs_max = {4} (truth: {5})
Ks = {6} (truth: {7})
Yg = {8} (truth: {9})
mu_max = {10} (truth: {11})
a = {12} (truth: {13})
delta = {14} (truth: {15})
f = {16} (truth: {17})
""".format(ms1_ml, ms1, ms2_ml, ms2, qs_max_ml, qs_max, Ks_ml, Ks, Yg_ml, Yg, mu_max_ml,
mu_max, a_ml, a, delta_ml, delta, np.exp(lnf_ml), f_true))
# Set up the sampler.
ndim, nwalkers = 9, 1000
pos = [result["x"] + 1e-4*np.random.randn(ndim) for i in range(nwalkers)]
sampler = emcee.EnsembleSampler(nwalkers, ndim, lnprob, args=(alpha1, Smin, Sminerr))
# Clear and run the production chain.
print("Running MCMC...")
sampler.run_mcmc(pos, 3500, rstate0=np.random.get_state())
samples = sampler.get_chain()
burnin = 50
samples = sampler.chain[:, burnin:, :].reshape((-1, ndim))
fig, axes = plt.subplots(3,3, figsize=(12, 6), sharex=True)
labels = ["$m_s$", "$m\prime_s$", "$q_s^{max}$", "$K_s$", "$Y_G$", "$\mu_{max}$", "$\psi$", "$\delta$", "$\ln\,f$"]
m,n =0,0
for i in range(ndim):
ax = axes[m,n]
ax.plot(samples[:, i], "k")
ax.set_xlim(0, len(samples))
ax.set_ylabel(labels[i])
ax.ticklabel_format(style='sci', axis='x', scilimits=(0,0))
n+=1
if n==3:
n=0
m+=1
fig.tight_layout()
axes[2,1].set_xlabel("step number");
fig.savefig('samples_mcmc.png')
print("Done.")
rcParams["savefig.dpi"] = 300
rcParams["figure.dpi"] = 300
rcParams["font.size"] = 14
rcParams['xtick.major.pad']='18'
rcParams['ytick.major.pad']='18'
# Make the triangle plot.
plt.figure()
fig = corner.corner(samples, labels=["$m_s$", "$m\prime_s$", "$q_s^{max}$", "$K_s$", "$Y_G$", "$\mu_{max}$", "$\psi$", "$\delta$", "$\ln\,f$"],
truths=[ms1, ms2, qs_max, Ks, Yg, mu_max, a, delta, np.log(f_true)], labelpad=280)
Global Senstivity Analysis¶
In [7]:
def ET1(X):
mu = 0.032
alpha, ms1, ms2, qs_max, Ks, Yg, a, delta, mu_max = X[:,0], X[:,1], X[:,2], X[:,3], X[:,4], X[:,5], X[:,6], X[:,7], X[:,8]
A=mu*alpha/Yg+ms1*alpha+ms2*(1-alpha)+(mu/mu_max)*(1-alpha)*a*delta
Cs_min = Ks*A/(qs_max-A)
return Cs_min
alpha= 0.7
factor = 0.3
fl=1-factor
fh=1+factor
ms1 = 0.22
ms2 = 0.10
qs_max = 4.08
Ks = 237e-6
Yg = 0.028
mu_max = 0.11
a = 0.028
delta= 0.002
problem = {'num_vars': 9,
'names': [r'$\alpha$',
r'$m_s$',
r'$m\prime_s$',
r'$q_s^{max}$',
r'$K_s$',
r'$Y_g$',
r'$\psi$',
r'$\delta$',
'$\mu_{max}$' ],
'bounds': [[fl*alpha, fh*alpha],
[fl*ms1, fh*ms1],
[fl*ms2, fh*ms2],
[fl*qs_max, fh*qs_max],
[fl*Ks, fh*Ks],
[fl*Yg, fh*Yg],
[fl*a, fh*a],
[fl*delta, fh*delta],
[fl*mu_max, fh*mu_max]]
}
# Generate samples
param_values = saltelli.sample(problem, 10000, calc_second_order=True)
# Run model (example)
Y = ET1(param_values)
# Perform analysis
Si = sobol.analyze(problem, Y, print_to_console=False, parallel=True, n_processors=4)
res = pd.DataFrame([Si['S1'], Si['S1_conf'], Si['ST'], Si['ST_conf'], Si['S2'], Si['S2_conf']])
res.columns = problem['names']
res.index = ['S1', 'S1_conf', 'ST', 'ST_conf', 'S2', 'S2_conf']
res = res.T
s2= pd.DataFrame(Si['S2']).T
s2.columns = problem['names']
s2.index = problem['names']
sns.set(style="white")
sns.set(font_scale=1.2)
fig, ax = plt.subplots(2, 3, figsize=(12,8))######
ax[0,0].errorbar(res.S1, res.index, xerr=res.S1_conf, fmt='o')
ax[0,1].errorbar(res.ST, res.index, xerr=res.ST_conf, fmt='o')
ax[0,0].set_title('S1');
ax[0,1].set_title('ST');
ax[0,2].set_title('S2');
ax[0,0].set_xlabel('Variance');
ax[0,1].set_xlabel('Variance');
mask = np.zeros_like(s2, dtype=np.bool)
# Generate a custom diverging colormap
cmap = sns.diverging_palette(220, 10, as_cmap=True)
# Draw the heatmap with the mask and correct aspect ratio
sns.heatmap(s2, mask=mask, cmap=cmap, vmax=s2.max().max(), center=0,
square=True, linewidths=.5, cbar_kws={"shrink": .9,'label':'Variance'}, ax=ax[0,2]);
# fig.savefig('sobol1.svg')
def ET2(X):
mu = 0.032
D=0.032
gamma, delta, mu_max, kd1, x2 = X[:,0], X[:,1], X[:,2], X[:,3], X[:,4]
x1=-(mu/mu_max*delta*x2)/(-D+mu-((1-mu/mu_max)*gamma)-kd1)
return x1
problem = {'num_vars': 5,
'names': ['$\gamma$',
'$\delta$',
'$\mu_{max}$',
'$K_{d1}$',
'$x_2$'],
'bounds': [[0.002, 0.1],
[0.001, 0.05],
[0.1, 0.2],
[0.002, 0.006],
[0.001, 0.003]]
}
# Generate samples
param_values = saltelli.sample(problem, 1000, calc_second_order=True)
# Run model (example)
Y = ET2(param_values)
# Perform analysis
Si = sobol.analyze(problem, Y, print_to_console=False, parallel=True, n_processors=4)
res = pd.DataFrame([Si['S1'], Si['S1_conf'], Si['ST'], Si['ST_conf'], Si['S2'], Si['S2_conf']])
res.columns = problem['names']
res.index = ['S1', 'S1_conf', 'ST', 'ST_conf', 'S2', 'S2_conf']
res = res.T
s2= pd.DataFrame(Si['S2']).T
s2.columns = problem['names']
s2.index = problem['names']
# fig, ax = plt.subplots(1, 3, figsize=(15,4))
ax[1,0].errorbar(res.S1, res.index, xerr=res.S1_conf, fmt='o')
ax[1,1].errorbar(res.ST, res.index, xerr=res.ST_conf, fmt='o')
ax[1,0].set_title('S1');
ax[1,1].set_title('ST');
ax[1,2].set_title('S2');
ax[1,0].set_xlabel('Variance');
ax[1,1].set_xlabel('Variance');
mask = np.zeros_like(s2, dtype=np.bool)
# Generate a custom diverging colormap
cmap = sns.diverging_palette(220, 10, as_cmap=True)
# Draw the heatmap with the mask and correct aspect ratio
sns.heatmap(s2, mask=mask, cmap=cmap, vmax=s2.max().max(), center=0,
square=True, linewidths=.5, cbar_kws={"shrink": .9,'label':'Variance'}, ax=ax[1,2]);
fig.tight_layout()
fig.savefig("sobol1+2.svg")
In [8]:
def ET(X):
mu = 0.032
D=0.032
gamma, delta, mu_max, kd2, x1 = X[:,0], X[:,1], X[:,2], X[:,3], X[:,4]
x2=-(1-mu/mu_max)*gamma*x1/(-D-mu/mu_max*delta-kd2)
return x2
problem = {'num_vars': 5,
'names': ['$\gamma$',
'$\delta$',
'$\mu_{max}$',
'$K_{d2}$',
'$x_1$'],
'bounds': [[0.002, 0.1],
[0.001, 0.05],
[0.1, 0.2],
[0.002, 0.006],
[0.001, 0.003]]
}
# Generate samples
param_values = saltelli.sample(problem, 1000, calc_second_order=True)
# Run model (example)
Y = ET(param_values)
# Perform analysis
Si = sobol.analyze(problem, Y, print_to_console=False, parallel=True, n_processors=4)
res = pd.DataFrame([Si['S1'], Si['S1_conf'], Si['ST'], Si['ST_conf'], Si['S2'], Si['S2_conf']])
res.columns = problem['names']
res.index = ['S1', 'S1_conf', 'ST', 'ST_conf', 'S2', 'S2_conf']
res = res.T
s2= pd.DataFrame(Si['S2']).T
s2.columns = problem['names']
s2.index = problem['names']
sns.set(style="white")
sns.set(font_scale=1.2)
fig, ax = plt.subplots(1, 3, figsize=(15,4))
ax[0].errorbar(res.S1, res.index, xerr=res.S1_conf, fmt='o')
ax[1].errorbar(res.ST, res.index, xerr=res.ST_conf, fmt='o')
ax[0].set_title('S1');
ax[1].set_title('ST');
ax[2].set_title('S2');
ax[0].set_xlabel('Variance');
ax[1].set_xlabel('Variance');
mask = np.zeros_like(s2, dtype=np.bool)
# mask[np.triu_indices_from(mask)] = True
# Generate a custom diverging colormap
cmap = sns.diverging_palette(220, 10, as_cmap=True)
# Draw the heatmap with the mask and correct aspect ratio
sns.heatmap(s2, mask=mask, cmap=cmap, vmax=s2.max().max(), center=0,
square=True, linewidths=.5, cbar_kws={"shrink": .9,'label':'Variance'}, ax=ax[2]);
fig.savefig("sobol3.png")
Scenario analysis¶
Constant low substrate availability¶
In [9]:
def Substrate_min2(alpha, ms1, ms2, qs_max, Ks):
mu_max = 0.11
mu=0.023
a = 0.028
delta= 0.002
A=mu*alpha/Yg+ms1*alpha+ms2*(1-alpha)+(mu/mu_max)*(1-alpha)*a*delta
Cs_min = Ks*A/(qs_max-A)
return Cs_min
matplotlib.rcParams.update({'font.size': 11})
qs_max_all = [4, 2]
Ks_all = [220e-6, 100e-6]
ms2_all = [0.02, 0.08, 0.12]
alpha1 = np.linspace(0.05, 1, 14)
ms1 = np.linspace(0.05, 0.3, 14)
ALPHA, MS1 = np.meshgrid(alpha1, ms1) # create the "base grid"
sns.set(font_scale=1.8)
fig, ax = plt.subplots(4,3, figsize=(15,14), sharex=True, sharey=True)
i, j = 0, 0
for Ks in Ks_all:
for qs_max in qs_max_all:
for ms2 in ms2_all:
Cs_MIN= Substrate_min2(ALPHA, MS1, ms2, qs_max, Ks)*1e6
CS = ax[i,j].contourf(ALPHA, MS1, Cs_MIN, levels = np.arange(0, 100, 5), cmap='jet_r')
CS2 = ax[i,j].contour(CS, levels=CS.levels[::2],
colors='black',
origin='lower')
ax[i,j].clabel(CS2, inline=1, fmt='%1.f')
if i == 0: ax[i,j].set_title('Non-growing cell type I' if j==0 else 'Non-growing cell type II' if j==1 else 'Non-growing cell type III')
if j == 0: ax[i,j].set_ylabel('Growing cell type I' if i==0 else 'Growing cell type II' if i==1 else 'Growing cell type III' if i==2 else 'Growing cell type IV')
j+=1
i+=1
j=0
cbar_ax = fig.add_axes([1.05, 0.1, 0.03, 0.8])
plt.colorbar(CS, cax=cbar_ax)
cbar_ax.set_ylabel(r'$Cs_{min}$')
plt.tight_layout()
fig.text(0.5,-0.01, r'Growing fraction of population, $\alpha$', ha="center", va="center")
fig.text(-0.01,0.5, r'Maintenance energy requirement of growing population, $m_s$', ha="center", va="center", rotation=90);
fig.savefig('contour_plot.svg', dpi=300)
Pulse feeding¶
In [10]:
def dy_dt(t, y):
"""
y1 = dS/dt
y2 = dx/dt
INPUTS:
t: time
y: the time varying concentrations: substrate (S) and biomass (x)
"""
sig = t % P < D
S_in = sig*30
F = 0.023*V*sig # ml/hr
S = y[0]
x1 = y[1]
x2 = y[2]
y = np.zeros((3,1))
gamma, delta,a = 0.006,0.002,0.028
qs = m*x1+m2*x2
p = (mu_max*S/(Ks+S)) * x1
mu= mu_max*S/(Ks+S)
K = (1-mu/mu_max)*gamma*x1-mu/mu_max*delta*x2
if F <=0:
y[0]= S - (1/Y_g)*p- qs
y[1] = p -K
y[2] = -x2*F/V +K
if y[0]<=0:
y[0]=S - (1/Y_g)*p
y[1] = p -K-kd*x1
y[2] = K-kd2*x2
else:
y[0] = F/V*(S_in - S) - (1/Y_g)*p -qs-((F/V)/mu_max)*delta*x2*a
y[1] = -x1*F/V +p -K
y[2] = -x2*F/V +K
return y
V = 2000 # ml
Y_g = 0.028 # g/g yield
mu_max = 0.11 #
Ks = 0.237 #mg/L
m=0.22 # g/g-h, maintenance demand of growing biomass
m2=0.10 # g/g-h, maintenance demand of non-growing biomass
kd=m*Y_g # 1/h, endogenous decay rate for growing biomass
kd2=kd/2 # 1/h, endogenous decay rate for non-growing biomass
P = 40
D = 5
r = integrate.ode(dy_dt).set_integrator('vode', method='bdf')
ICs = [0.076, 0.43, 0.1]
t_0 = 0.0
r.set_initial_value(ICs, t_0)
t_final = 800.0
dt = 1
# Create vectors to store the solutions in;
# add extra space for intial condition
n_steps = np.int((t_final - t_0)/dt) + 1
time = np.zeros(n_steps)
S = np.zeros(n_steps)
x1 = np.zeros(n_steps)
x2 = np.zeros(n_steps)
S[0], x1[0], x2[0] = ICs
k = 1
while r.successful() and r.t < t_final:
r.integrate(r.t + dt)
time[k] = r.t
S[k] = r.y[0]
x1[k] = r.y[1]
x2[k] = r.y[2]
k += 1
sns.set(style="white")
sns.set(font_scale=1.2)
fig, ax1 = plt.subplots(1,1, figsize=(11,4))
lns1 = ax1.plot(time, S*1e3, 'r', label='Substrate')
lns2 = ax1.plot(time, x1*1e3, 'k', label='Growing biomass')
# ax2 = ax1.twinx()
lns3 = ax1.plot(time, x2*1e3, 'b', label='Non-growing biomass')
lns4 = ax1.plot(time, (x1+x2)*1e3, 'g', label='Total biomass')
# yscale("log")
ax1.set_ylabel(r'Concentration [$\mu$gL$^{-1}$]')
ax1.set_xlabel(r'Time [h]')
# added these three lines
lns = lns1+lns2+lns3+lns4
labs = [l.get_label() for l in lns]
ax1.legend(lns, labs, loc=0, ncol=4)
ax1.grid(True)
ax1.set_ylim(-5,750)
fig.savefig("inter_feed_het.svg")
# y_true = np.array([88, 628, 42])
# ax1.scatter([120, 120], y_true[:2], c=['r', 'k'])
# ax2.scatter([120], y_true[2], c='b')
In [11]:
def dy_dt(t, y):
"""
y1 = dS/dt
y2 = dx/dt
INPUTS:
t: time
y: the time varying concentrations: substrate (S) and biomass (x)
"""
V = 2000 # ml 1
Y_g = 0.028 # g/g
mu_max = 0.11 # /h
Ks = 0.237 #mg/L
m=0.22 # g/g-h
kd=m*Y_g
F = 0.023*V# ml/hr
P =40 # period
D =5 # width of pulse
sig = t % P < D
S_in = sig*30
F = 0.023*V*sig # ml/hr
S = y[0]
x = y[1]
y = np.zeros((2,1))
qs = m * x
p = (mu_max*S/(Ks+S)) * x
if F <=0:
y[0]=S - (1/Y_g)*p- qs
y[1]= p
if y[0]<=0:
y[0]=S - (1/Y_g)*p
y[1]= p-kd*x
else:
y[0] = F/V*(S_in - S) - (1/Y_g)*p- qs
y[1] = -x*F/V +p
return y
r = integrate.ode(dy_dt).set_integrator('vode', method='bdf')
ICs = [0.080, 0.40]
t_0 = 0.0
r.set_initial_value(ICs, t_0)
t_final = 800.0
dt = 2
# Create vectors to store the solutions in;
# add extra space for intial condition
n_steps = np.int((t_final - t_0)/dt) + 1
time = np.zeros(n_steps)
S = np.zeros(n_steps)
x = np.zeros(n_steps)
S[0], x[0] = ICs
k = 1
while r.successful() and r.t < t_final:
r.integrate(r.t + dt)
time[k] = r.t
S[k] = r.y[0]
x[k] = r.y[1]
k += 1
# Clear figure window from previous simulation
sns.set(style="white")
sns.set(font_scale=1.2)
fig, ax1 = plt.subplots(1,1, figsize=(11,4))
ax1.plot(time, S*1e3, 'r', label='Substrate')
ax1.plot(time, x*1e3, 'k', label='Biomass')
ax1.set_ylabel(r'Concentration [$\mu$gL$^{-1}$]')
ax1.set_xlabel(r'Time [h]')
ax1.legend(loc='best', ncol=2)
ax1.grid(True)
ax1.set_ylim(-5,500)
fig.savefig("inter_feed_homo.svg")
Sinusoidal variation in inflow rate¶
In [12]:
def dy_dt(t, y):
"""
y1 = dS/dt
y2 = dx/dt
INPUTS:
t: time
y: the time varying concentrations: substrate (S) and biomass (x)
"""
V = 2000 # ml
Y_g = 0.028 # g/g
mu_max = 0.11 # /h
Ks = 0.237 #mg/L
m=0.22 # g/g-h
m2=0.10 # g/g-h
D = 0.006
F = D*V+0.6*D*V*sin(2*np.pi/24*t -np.pi/5)
S_in=30
S = y[0]
x1 = y[1]
x2 = y[2]
y = np.zeros((3,1))
gamma, delta,a = 0.006,0.002,0.028 #####################
qs = m*x1+m2*x2
p = (mu_max*S/(Ks+S)) * x1
mu= mu_max*S/(Ks+S)
K = (1-mu/mu_max)*gamma*x1-mu/mu_max*delta*x2
y[0] = F/V*(S_in - S) - (1/Y_g)*p -qs-((F/V)/mu_max)*delta*x2*a
y[1] = -x1*F/V +p -K
y[2] = -x2*F/V +K
return y
r = integrate.ode(dy_dt).set_integrator('vode', method='bdf')
# Part 1:
ICs = [0.076, 0.43, 0.1]
# Part 2: use the steady-state solution
# ICs = [33.33333333333, 53.3333333333]
t_0 = 0.0
r.set_initial_value(ICs, t_0)
t_final = 400.0
dt = 1
# Create vectors to store the solutions in;
# add extra space for intial condition
n_steps = np.int((t_final - t_0)/dt) + 1
time = np.zeros(n_steps)
S = np.zeros(n_steps)
x1 = np.zeros(n_steps)
x2 = np.zeros(n_steps)
S[0], x1[0], x2[0] = ICs
k = 1
while r.successful() and r.t < t_final:
r.integrate(r.t + dt)
time[k] = r.t
S[k] = r.y[0]
x1[k] = r.y[1]
x2[k] = r.y[2]
k += 1
sns.set(style="white")
sns.set(font_scale=1.2)
fig, ax1 = plt.subplots(1,1, figsize=(11,4))
lns1 = ax1.plot(time, S*1e3, 'r', label='Substrate')
lns2 = ax1.plot(time, x1*1e3, 'k', label='Growing biomass')
# ax2 = ax1.twinx()
lns3 = ax1.plot(time, x2*1e3, 'b', label='Non-growning biomass')
lns4 = ax1.plot(time, (x1+x2)*1e3, 'g', label='Total biomass')
# yscale("log")
ax1.set_ylabel(r'Concentration [$\mu$gL$^{-1}$]')
ax1.set_xlabel(r'Time [h]')
ax1.set_ylim(0,750)
# added these three lines
lns = lns1+lns2+lns3+lns4
labs = [l.get_label() for l in lns]
ax1.legend(lns, labs, loc=0)
ax1.grid(True)
ax1.set_ylim(-5,750)
fig.savefig("fluc_rapid.svg")
In [13]:
def dy_dt(t, y):
"""
y1 = dS/dt
y2 = dx/dt
INPUTS:
t: time
y: the time varying concentrations: substrate (S) and biomass (x)
"""
V = 2000 # ml
Y_g = 0.028 # g/g
mu_max = 0.11 # /h
Ks = 0.237 #mg/L
m=0.22 # g/g-h
kd=m*Y_g
D = 0.006
F = D*V+0.6*D*V*sin(2*np.pi/24*t -np.pi/5)
S_in = 30 ## mg/L
S = y[0]
x = y[1]
y = np.zeros((2,1))
qs = m * x
p = (mu_max*S/(Ks+S)) * x
if F <=0:
y[0]=S - (1/Y_g)*p- qs
y[1]= p
if y[0]<=0:
y[0]=S - (1/Y_g)*p
y[1]= p-kd*x
else:
y[0] = F/V*(S_in - S) - (1/Y_g)*p- qs
y[1] = -x*F/V +p
return y
r = integrate.ode(dy_dt).set_integrator('vode', method='bdf')
# Part 1:
ICs = [0.080, 0.40]
# Part 2: use the steady-state solution
# ICs = [33.33333333333, 53.3333333333]
t_0 = 0.0
r.set_initial_value(ICs, t_0)
t_final = 400.0
dt = 2
# Create vectors to store the solutions in;
# add extra space for intial condition
n_steps = np.int((t_final - t_0)/dt) + 1
time = np.zeros(n_steps)
S = np.zeros(n_steps)
x = np.zeros(n_steps)
S[0], x[0] = ICs
k = 1
while r.successful() and r.t < t_final:
r.integrate(r.t + dt)
time[k] = r.t
S[k] = r.y[0]
x[k] = r.y[1]
k += 1
# Clear figure window from previous simulation
sns.set(style="white")
sns.set(font_scale=1.2)
fig, ax1 = plt.subplots(1,1, figsize=(11,4))
ax1.plot(time, S*1e3, 'r', label='Substrate')
ax1.plot(time, x*1e3, 'k', label='Biomass')
ax1.set_ylabel(r'Concentration [$\mu$gL$^{-1}$]')
ax1.set_xlabel(r'Time [h]')
ax1.legend(loc='best')
ax1.grid(True)
ax1.set_ylim(-5,750)
fig.savefig("fluc_rapid_homo.svg")
