Here, we are going to learn about the mathematical model of the aerobic growth of Saccharomyces cerevisiae

First, we need to create a class, as we will be using object-oriented modelling (see introductory video in biovl.com) and we initialize all the parameters in the function def init()

In [2]:
class Scerevisiae_Aero:
def __init__(self):
self.Yox_XG = 0.8 #
self.Yred_XG = 0.05
self.Yox_XE =  0.72
self.Y_OG = 1.067
self.Y_EG = 0.5
self.Y_OE = 1.5
self.q_g = 3.5
self.q_o = 0.37
self.q_e = 0.32
self.t_lag = 4.66
self.Kg = 0.17
self.Ke = 0.56
self.Ko = 0.0001
self.Ki = 0.31
self.O_sat = 0.00755
self.kla = 1004

self.G0 = 50
self.E0 = 0.0
self.O0 = 0.00755
self.X0 = 0.4
self.V0 = 20
self.T0 = 30

self.t_end = 50
self.t_start = 2
self.steps = (self.t_end - self.t_start)*24
self.Cin = 50
self.Fconst = 0.5
self.t_feed_start = 5

self.u= 75 #mass flow of cooling water

def rxn(self, C,t):
#number of components
n = 4
m = 4

#initialize the stoichiometric matrix, s
s = np.zeros((m,n))
s[0,0] = -1
s[0,1] = 0
s[0,2] = -self.Y_OG
s[0,3] = self.Yox_XG

s[1,0] = -1
s[1,1] = self.Y_EG
s[1,2] = 0
s[1,3] = self.Yred_XG

s[2,0] = 0
s[2,1] = -1
s[2,2] = -self.Y_OE
s[2,3] = self.Yox_XE

s[3,0] = 0
s[3,1] = 0
s[3,2] = 1
s[3,3] = 0

#initialize the rate vector
rho = np.zeros((4,1))

##initialize the overall conversion vector
r=np.zeros((4,1))

#Developing the matrix, the overall conversion rate is stoichiometric *rates
rho[0,0] = ((1/self.Y_OG)*min(self.q_o*(C[2]/(C[2]+self.Ko)),self.Y_OG*(self.q_g*(C[0]/(C[0]+self.Kg)))))*C[3]
rho[1,0] = ((1-math.exp(-t/self.t_lag))*((self.q_g*(C[0]/(C[0]+self.Kg)))-(1/self.Y_OG)*min(self.q_o*(C[2]/(C[2]+self.Ko)),self.Y_OG*(self.q_g*(C[0]/(C[0]+self.Kg))))))*C[3]
rho[2,0] = ((1/self.Y_OE)*min(self.q_o*(C[2]/(C[2]+self.Ko))-(1/self.Y_OG)*min(self.q_o*(C[2]/(C[2]+self.Ko)),self.Y_OG*(self.q_g*(C[0]/(C[0]+self.Kg)))),self.Y_OE*(self.q_e*(C[1]/(C[1]+self.Ke))*(self.Ki/(C[0]+self.Ki)))))*C[3]
rho[3,0] = self.kla*(self.O_sat - C[2])

#Volume balance
if (t>=self.t_feed_start):
Fin = self.Fconst
Fout = self.Fconst
else:
Fin = 0
Fout = 0
F = Fin - Fout

#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])+(s[3,0]*rho[3,0])
r[1,0] = (s[0,1]*rho[0,0])+(s[1,1]*rho[1,0])+(s[2,1]*rho[2,0])+(s[3,1]*rho[3,0])
r[2,0] = (s[0,2]*rho[0,0])+(s[1,2]*rho[1,0])+(s[2,2]*rho[2,0])+(s[3,2]*rho[3,0])
r[3,0] = (s[0,3]*rho[0,0])+(s[1,3]*rho[1,0])+(s[2,3]*rho[2,0])+(s[3,3]*rho[3,0])

#Solving the mass balances
dGdt = r[0,0] -F/C[4]*C[0] + Fin/C[4]*self.Cin - Fout/C[4]*C[0]
dEdt = r[1,0] -F/C[4]*C[1] - Fout/C[4]*C[1]
dOdt = r[2,0]
dXdt = r[3,0] -F/C[4]*C[3] - Fout/C[4]*C[3]
dVdt = F

#Solving the heat balance
#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]  + dGdt*C[4]*(15580)- dEdt*C[4]*(29710)
#Shaft work 1 W/L1
W = -1*C[4] #[J/S] negative because exothermic
#Cooling just an initial value (constant cooling to see what happens)
#dQdt = -0.03*C[4]*(-21200) #[J/S]
#velocity of cooling water: u [m3/h] -->controlled by PID
#Mass flow cooling water
M= self.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[5]
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]

return [dGdt,dEdt,dOdt,dXdt,dVdt, dTdt]

def solve(self):
#time
t = np.linspace(self.t_start, self.t_end, self.steps)
C0 = [self.G0, self.E0, self.O0, self.X0, self.V0, self.T0]
C = odeint(self.rxn, C0, t, rtol = 1e-7, mxstep= 500000)
return C,t



Now, we do write the decomposed reactions, using a matrix notation

Once we have implemented the model, we need to solve the differential equations of [dGdt,dEdt,dOdt,dXdt,dVdt] through a solver. In this case the solver will be scipy.odeint

Function def solve(self) starts with initial concentrations (C0) and through the self.rxn (= def rxn(self, C,t)), calculate the concentration for the time step 1, then it goes back to def solve and using the concentrations for the time step 1, calculates the concentrations for the time step 2, etc.

However, for this to work, we need to import some packages

In [3]:
#Package for solving the ODE
from scipy.integrate import odeint
#Package for plotting
import math
#Package for the use of vectors and matrix
import numpy as np
import sys
import os
import matplotlib.pyplot as plt

In [4]:
model = Scerevisiae_Aero()
solution_model = model.solve()

In [5]:
C = solution_model[0]
G = C[:,0]
E = C[:,1]
X = C[:,3]

t = solution_model[1]

In [6]:
plt.plot(t, G, label = 'Glucose')
plt.plot(t, E,label = 'Ethanol')
plt.plot(t, X,label = 'Biomass')

plt.xlabel('t (h)')
plt.ylabel('Concentration (g/L)')
plt.legend()

plt.show

Out[6]:
<function matplotlib.pyplot.show(*args, **kw)>

# How to... modify the model¶

First, we learn how we change variables and recalculate the model

In [17]:
model = Scerevisiae_Aero()
print('Old value')
print(model.Yox_XG)

model_new = Scerevisiae_Aero()

model_new.Yox_XG = 0.02
print('New value')
print(model_new.Yox_XG)

solution_model_new = model_new.solve()

0.8
0.02

In [18]:
C = solution_model_new[0]
G = C[:,0]
E = C[:,1]
X = C[:,3]

t = solution_model_new[1]

plt.plot(t, G, label = 'Glucose')
plt.plot(t, E,label = 'Ethanol')
plt.plot(t, X,label = 'Biomass')

plt.xlabel('t (h)')
plt.ylabel('Concentration (g/L)')
plt.legend()

plt.show

Out[18]:
<function matplotlib.pyplot.show(*args, **kw)>

What differences do you see?

In [ ]:



## Question for you¶

Now that you know how to change a variable and recalculate the model; propose a new value for a variable with the aim to increase the consumption of ethanol.

Plot it so you can see the effect

In [ ]: