Created on Thu Sep 6 13:34:32 2018 By Simoneta Cano de las Heras
Adapted by Mariana Albino, October 2024
Aerobic fermentation model (S. cerevisiae)#
This model describes an aerobic batch S. cerevisiae fermentation. It uses Monod-Herbert expressions to describe the kinetics of the process. It can take into account the metabolic heat produced by the cells, in which case a PID controller is used to maintain a stable temperature
Package import#
This portion of the code handles the import of all the relevant python packages.
#import the necesary packages
from scipy.integrate import odeint
#Package for plotting
import math
#Package for the use of vectors and matrix
import numpy as np
import array as arr
from matplotlib.figure import Figure
import sys
import os
import matplotlib.pyplot as plt
from matplotlib.ticker import FormatStrFormatter
import glob
from random import sample
import random
import time
import plotly.graph_objects as go
import plotly
import json
import pandas as pd
---------------------------------------------------------------------------
ModuleNotFoundError Traceback (most recent call last)
Cell In[1], line 3
1 #import the necesary packages
----> 3 from scipy.integrate import odeint
4 #Package for plotting
5 import math
ModuleNotFoundError: No module named 'scipy'
Model definition#
The model is defined by using a class. This has several adavantages, one being that different classes can be defined independtly, e.g., one class represents the model definition while another will handle parameter optimization with experimental data. The same “optimization” class can be used with different model classes, so the code becomes more easily reusable.
The class Monod_Herbert includes several functions, each with a specific function.
init This function initialises the model class by defining all the relavant parameters and initial conditions.
rxn This function includes all the model equations. The model uses matrix notation to define the ODEs
solve This function generates the timesteps for solving the ODEs. It includes an approach for when the “CONTROL” option is turned OFF and an approach that solves a PDI controller for temperature control, for when the “CONTROL” feature is ON.
create_plot This function stores the results of the simulation on a dataframe which is used for plotting all relevant variables.
class Monod_Herbert:
#initialise the model
def __init__(self, Control=False): #change to Control=True if you want to see the effect of Temperature control
# define value of model parameters
self.Y_XS = 0.8
self.Y_OX = 1.05
self.y_x = 0.5
self.mu_max = 1.1
self.Ks = 0.17
self.kd = 0.08
self.kla = 1004
self.O_sat = 0.0755
#define initial conditions
self.G0 = 18
self.O0 = 0.0755
self.X0 = 0.01
self.V0 = 40
#define parameters for control, default every 1/24 hours:
self.t_end = 30
self.t_start = 0
self.Control = False
self.coolingOn = True
self.steps = (self.t_end - self.t_start)*24
self.T0 = 30
self.K_p = 2.31e+01
self.K_i = 3.03e-01
self.K_d = -3.58e-03
self.Tset = 30
self.u_max = 150
self.u_min = 1
#define the stoichiometric matrix
def rxn(self, C,t, u):
#when there is no control, k has no effect
k=1
#when cooling is off than u = 0
if self.coolingOn == False:
u = 0 #flow of cooling water
#define temperature controller
if self.Control == True :
#Cardinal temperature model with inflection: Salvado et al 2011 "Temperature Adaptation Markedly Determines Evolution within the Genus Saccharomyces"
#Strain S. cerevisiae PE35 M
#How is the growth rate of the organism affected by the changes in temperature
Topt = 30
Tmax = 45.48
Tmin = 5.04
T = C[4]
if T < Tmin or T > Tmax:
k = 0 #growth rate
else:
D = (T-Tmax)*(T-Tmin)**2
E = (Topt-Tmin)*((Topt-Tmin)*(T-Topt)-(Topt-Tmax)*(Topt+Tmin-2*T))
k = D/E
#number of components
n = 3
m = 3
#initialize the stoichiometric matrix, s
s = np.zeros((m,n))
s[0,0] = -1/self.Y_XS
s[0,1] = -1/self.Y_OX
s[0,2] = 1
s[1,0] = 0
s[1,1] = 1/self.y_x
s[1,2] = -1
s[2,0] = 0
s[2,1] = self.kla
s[2,2] = 0
#initialize the rate vector
rho = np.zeros((m,1))
##initialize the overall conversion vector
r=np.zeros((n,1))
rho[0,0] = self.mu_max*(C[0]/(C[0]+self.Ks))*C[2]
rho[1,0] = self.kd*C[2]
rho[2,0] = self.kla*(self.O_sat - C[1])
#Developing the matrix, the overall conversion rate is stoichiometric *rates
r[0,0] = (s[0,0]*rho[0,0])+(s[1,0]*rho[1,0])+(s[2,0]*rho[2,0])
r[1,0] = (s[0,1]*rho[0,0])+(s[1,1]*rho[1,0])+(s[2,1]*rho[2,0])
r[2,0] = (s[0,2]*rho[0,0])+(s[1,2]*rho[1,0])+(s[2,2]*rho[2,0])
#Solving the mass balances
dSdt = r[0,0]
dOdt = r[1,0]
dXdt = r[2,0]
dVdt = 0
if self.Control == True :
'''
dHrxn heat produced by cells estimated by yeast heat combustion coeficcient dhc0 = -21.2 kJ/g
dHrxn = dGdt*V*dhc0(G)-dEdt*V*dhc0(E)-dXdt*V*dhc0(X)
(when cooling is working) Q = - dHrxn -W ,
dT = V[L] * 1000 g/L / 4.1868 [J/gK]*dE [kJ]*1000 J/KJ
dhc0(EtOH) = -1366.8 kJ/gmol/46 g/gmol [KJ/g]
dhc0(Glc) = -2805 kJ/gmol/180g/gmol [KJ/g]
'''
#Metabolic heat: [W]=[J/s], dhc0 from book "Bioprocess Engineering Principles" (Pauline M. Doran) : Appendix Table C.8
dHrxndt = dXdt*C[4]*(-21200) #[J/s]
#Shaft work 1 W/L1
W = -1*C[4] #[J/S] negative because exothermic
#Mass flow cooling water
M=u/3600*1000 #[kg/s]
#Define Tin = 5 C, Tout=TReactor
#heat capacity water = 4190 J/kgK
Tin = 5
#Estimate water at outlet same as Temp in reactor
Tout = C[4]
cpc = 4190
#Calculate Q from Eq 9.47
Q=-M*cpc*(Tout-Tin) # J/s
#Calculate Temperature change
dTdt = -1*(dHrxndt - Q + W)/(C[4]*1000*4.1868) #[K/s]
else:
dTdt = 0
return [dSdt, dOdt, dXdt, dVdt, dTdt]
def solve(self):
t = np.linspace(self.t_start, self.t_end, self.steps) #generation of the time-points
#solve if Control is OFF:
if self.Control == False :
u = 0
C0 = [self.G0, self.O0, self.X0,self.V0, self.T0] #initial conditions vector
C = odeint(self.rxn, C0, t, rtol = 1e-7, mxstep= 500000, args=(u,)) #solve ODEs
#solve for Control ON
else:
"""
PID Temperature Control:
"""
# storage for recording values
C = np.ones([len(t), 5])
C0 = [self.G0, self.O0, self.X0,self.V0,self.T0]
self.ctrl_output = np.zeros(len(t)) # controller output
e = np.zeros(len(t)) # error
ie = np.zeros(len(t)) # integral of the error
dpv = np.zeros(len(t)) # derivative of the pv
P = np.zeros(len(t)) # proportional
I = np.zeros(len(t)) # integral
D = np.zeros(len(t)) # derivative
for i in range(len(t)-1):
#PID control of cooling water
dt = t[i+1]-t[i]
#Error
e[i] = C[i,4] - self.Tset
#print(e[i])
if i >= 1:
dpv[i] = (C[i,4]-C[i-1,4])/dt
ie[i] = ie[i-1] + e[i] * dt
#print(ie)
P[i]=self.K_p*e[i]
I[i]=self.K_i*ie[i]
D[i]=self.K_d*dpv[i]
self.ctrl_output[i]=P[i]+I[i]+D[i]
u=self.ctrl_output[i]
#print(u)
if u>self.u_max:
u=self.u_max
ie[i] = ie[i] - e[i]*dt # anti-reset windup
if u < self.u_min:
u =self.u_min
ie[i] = ie[i] - e[i]*dt # anti-reset windup
#print(u)
if i > 0 and e[i] * e[i-1] < 0:
ie[i] = 0
#time for solving ODE
ts = [t[i],t[i+1]]
#solve ODE from last timepoint to new timepoint with old values
y = odeint(self.rxn, C0, ts, rtol = 1e-7, mxstep= 500000, args=(u,))
#update C0
C0 = y[-1]
#merge y to C
C[i+1]=y[-1]
return t, C
#generate the plot of model variables
def create_plot(self, t, C):
G = C[:, 0]
O = C[:, 1]
B = C[:, 2]
T = C[:, 4]
df = pd.DataFrame({'t': t, 'Substrate': G, 'Biomass': B, 'Oxygen': O, 'Temperature' :T})
fig = go.Figure()
fig.add_trace(go.Scatter(x=df['t'], y=df['Substrate'], name='Substrate'))
fig.add_trace(go.Scatter(x=df['t'], y=df['Biomass'], name='Biomass'))
fig.add_trace(go.Scatter(x=df['t'], y=df['Oxygen'], name='Oxygen'))
fig.update_layout( xaxis_title='time (h)',
yaxis_title='Concentration (g/L)')
figT =go.Figure()
figT.add_trace(go.Scatter(x=df['t'], y=df['Temperature'], name='Temperature'))
figT.update_layout(xaxis_title='time (h)',
yaxis_title='Temperature (C)')
return fig, figT
Extract results and make plots#
This portion of the code extracts the results from the Monod_Herbert class and creates the figures to display the simulation results.
model = Monod_Herbert() # Instantiate the Monod_Herbert class
t, C = model.solve() # Call the solve method to get the simulation results
fig = model.create_plot(t, C)[0] # Call create_plot with the simulation results
figT = model.create_plot(t, C)[1] # Call create_plot with the simulation results
fig.show()
figT.show()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[3], line 2
1 model = Monod_Herbert() # Instantiate the Monod_Herbert class
----> 2 t, C = model.solve() # Call the solve method to get the simulation results
3 fig = model.create_plot(t, C)[0] # Call create_plot with the simulation results
4 figT = model.create_plot(t, C)[1] # Call create_plot with the simulation results
Cell In[2], line 140, in Monod_Herbert.solve(self)
138 def solve(self):
--> 140 t = np.linspace(self.t_start, self.t_end, self.steps) #generation of the time-points
142 #solve if Control is OFF:
143 if self.Control == False :
NameError: name 'np' is not defined
Using the model#
With this code you are able to test the behaviour of the process in a more simple case when we assume temperature is constant since the generation of metabolic heat is ignored. You can also turn ON the control option and evaluate the dymamics in Temperature. Furthermore, you can try to test different initila conditions, e.g., initial substrate concentrations (remember, these should be realistic values), and study the impact on the dynamics of Biomass.