ChE-311 "Biochemical Engineering"

Ipython demonstration for the exercises in analogy to the matlab scripts.

Batch Growth Model

The biological model assumes growth on xylene only. So the second Monod equation can be used for the specific growth $\mu$ and death rates $\alpha$ to predict substrate concentration $s$ and the biomass production for the given saturation (Monod) constant $K_s$: $$\mu=\mu_{max}\frac{s}{s+K_s}.$$ The rate of substrate usage is given by the third Monod equation: $$Y_{x/s}=\frac{r_x}{-r_s}$$ From the rates the concentration change can be obtained directly (first Monod equation): $$\frac{dx}{dt}=\left(\mu - \alpha\right) x.$$ The first order linear ODE is solved below numerically but it is also easily solved analytically by separation of variables: $$x(t)=x_0\cdot\exp^{(\mu-\alpha)t}$$ and indeed by simulation we find this exponential behaviour limited by the substrate concentration. However, for detailed and complex biological systems often analytical solutions are not available.

What is the effect of the concentration causing death due to xylene toxicity?

In [184]:
##Biochemical Model Function
import unyt as u
def BioChem_Bmodel(x,t):
    # Model Parameters Values
    MWt_Xyl = 106.2 *u.g/u.mol    # Molecular Weight of Xylene
    Ks_l    = 0.01                # Ks for our cells grown on Xylene (mM)
    max_g   = 0.008 /u.hr         # Maximum Theoretical Growthrate (h-1)
    Yxx     = 2.5                 # Yield of Biomass on Xylene (--)
    Ks_d    = 58.4                # Ks for death due to Xylene Toxicity (mM)
    max_d   = 0.25  /u.hr         # Maximum Theoretical Deathrate (h-1) 


    # Growth Rate (2. Monod eq.) 
    growth = max_g * (x[1]) / (x[1] + Ks_l)
    # Death Rate 
    death  = max_d * (x[1]) / (x[1] + Ks_d)

    # Change of Substrate Concentrateion Based on Rate (3. Monod eq.)
    # Use if you are having trouble with substrate < 0 values
    # if x[1] < 0:
    #     rxylene = 0
    # else:
    rxylene = (1/MWt_Xyl) * growth / Yxx

    # Differential for Biomass and Xylene Concentration
    biomass  = (growth - death) * x[0]   # biomass [=] Biomass (mg/L)
    xxylene  = - rxylene * x[0]          # xxylene [=] s, Xylene Concentration (mM)
    return [biomass, xxylene]




##load some modules for plotting and math for integration
%matplotlib inline
import matplotlib.pyplot as plt
import numpy
from scipy.integrate import odeint

##Integration Time Scale
tspan = numpy.linspace(0,340,num=500)

##Initial Conditions
Init_Cond = numpy.array([82.096,1.0]) #Xylene Concentration (mM), Biomass (mg/L)

##Call ODE Solver
solu = odeint(BioChem_Bmodel,Init_Cond,tspan)

##Plotting With Two Axes
x_bio, x_xyl = solu.transpose() #Xylene Concentration (mM), Biomass (mg/L)

fig, ax1 = plt.subplots()
color = 'tab:red'
ax1.set_xlabel('Time [hr]')
ax1.set_ylabel('Cell Concentration [mg/liter]', color=color)
ax1.plot(tspan, x_bio, color=color)
#Analytical solution for unlimited exponential behavior
initial_growth   =  (x_bio[1]-x_bio[0])/Init_Cond[0]/(tspan[1]-tspan[0])
plt.ylim(top=285)
ax1.plot(tspan, Init_Cond[0]*numpy.exp(initial_growth*tspan), color='red', linestyle=':', label='Exponential Growth at Initial Rate')
handles, labels = ax1.get_legend_handles_labels()
ax1.legend(handles, labels)
ax1.tick_params(axis='y', labelcolor=color)

ax2 = ax1.twinx()  # instantiate a second axes that shares the same x-axis
color = 'tab:blue'
ax2.set_ylabel('Substrate Concentration [mM]', color=color)
ax2.plot(tspan, x_xyl, color=color)
ax2.tick_params(axis='y', labelcolor=color)
plt.ylim(top=1.1) #to avoid cutting the label

fig.tight_layout()  # otherwise the right y-label is slightly clipped
plt.show()

Fed-Batch Cultivation Model

As for the batch process Monod behavior governs the growth ($\mu=\mu_{max}\frac{S}{S+K_S}$) and yield ($Y_{X/S}=\frac{r_X}{-r_S}$) but by a feed schedule can add substrate by a flow rate $F$. The new aspect for the fed-batch is that the volume changes ($F=\frac{dV}{dt}$ or integrated $V=V_0+Ft$) and with it concentrations of substrate, products and biomass are volume dependent, in other words: a dilution effect. Mathematically the time dependence of volume shows in the chain rule with intensive (independen of system size) concentration variables ($S$, $P$ & $X$) and the extensive volume $V$ (scaling with system size): $$\begin{align*} \frac{d(SV)}{dt}&=V\frac{dS}{dt}+S\frac{dV}{dt}=V\frac{dS}{dt}+SF(t)\\ \frac{d(XV)}{dt}&=V\frac{dX}{dt}+X\frac{dV}{dt}=V\frac{dX}{dt}+XF(t). \end{align*}$$ The first term has the volume dependence of the intensive variables, i.e. it explicitly contains the time dependence of the batch process (the rates from the mass balance) scaled by the volume. The second term the has the flow rate dependence. Also, the mass balance that holds at every instance: $$\begin{align*} \frac{d(SV)}{dt}&=FS_\text{feed}-V\frac{r_x(X,S)}{Y_{x/s}}\\ \frac{d(XV)}{dt}&=Vr_X. \end{align*}$$ By equating these and rearranging for the time dependence one obtains: $$\begin{align*} \frac{dS}{dt}&=\frac{F(t)(S_\text{feed}-S)}{V}-\frac{r_x(X,S)}{Y_{x/s}}\\ \frac{dX}{dt}&=-\frac{F(t)X}{V}+r_X. \end{align*}$$

These concentrations and the volume will be integrated numerically. Concentrations of products $P$ behave in analogue (proportional) to the biomass.

The feeding schedule in this model is exponential up to the limit of the machine. Hence, as discussed for exercise 6 3) the growth initially supports exponential growth at a defined growth rate. With the change to a linear feed the growth rate reduces and the biomass ultimately reaches an upper limited as the available substrate ($s\ll s_0$) is used as maintainance energy. In this quasi-steady-state of the total biomass the specific growth rate $\mu$ equals the dilution rate $D$ as becomes clear form the time drivative of the total cell concentration: $$\begin{align*} X&=x\cdot V [g]\\ \frac{dX}{dt}&=\frac{V\frac{dx}{dt}-x\frac{dV}{dt}}{V^2}\text{ with }F=\frac{dV}{dt}\text{ and }\mu=\frac{dx}{dt}/x\\ \frac{dX}{dt}&=\left(\mu-D\right)\cdot X. \end{align*}$$

In [26]:
#load some modules for plotting and math
%matplotlib inline
import matplotlib.pyplot as plt
import numpy
from scipy.integrate import odeint
import unyt as u

##Biochemical Model Function
def BioChem_FBmodel(y,t):
    x=y[0]
    s=y[1]
    V=y[2]*u.dm**3
    ## Model parameters
    qSmax   = 1.4   /u.hr #maximum growth rate
    Ks      = 0.1         #saturation or Monode constant
    qm      = 0.04  /u.hr #death rate
    Yem     = 0.5   *u.g/u.g  #yield
    Si      = 500         #Substrate initial
    F0      = 0.03  *u.dm**3/u.hr #initial feed rate
    SFR     = 0.2   *u.hr #rate of the exponential feed
    Fmax    = 1.5   *u.dm**3/u.hr #maximum possible feed
    
    ## Exponential Feed Schedule
    F = F0*numpy.exp(SFR*t)     #dm**3/hr
    if F > Fmax:
        F = Fmax
    
    qS   = qSmax*s/(s+Ks) #growth rate of your cells at given substrate concentration (2. Monod eq.) /hr
    My   = (qS-qm)*Yem          #biomass change /hr

    dy = numpy.zeros(3)
    dy[0] = -F/V*x+My*x       
    dy[1] = F/V *(Si-s)-qS*x   # /hr *  g/dm**3" (dimensions "(mass)/(length)**3") and "dimensionless"
    dy[2] = F
    return dy


## Integration Time Scale
tspan = numpy.linspace(0, 50, 300) *u.hr
## Initial Values
x_i   = 0.7           #mM
s_i   = 0.1           #mM
V_i   = 40  *u.dm**3  #L
y   = [x_i,s_i,V_i]

## Call ODE Solver
sol = odeint(BioChem_FBmodel,y,tspan)
x, s, v = sol.transpose()
V=v*u.dm**3  #only to get the units right, does not work with integration of function

## Growth and Dillution Rate 
mu=numpy.gradient(x)/(tspan[1]-tspan[0])/x                          #/hr
dil=numpy.gradient(V)/(tspan[1]-tspan[0])/V                         #/hr
X_total=numpy.array([x[i]/V[i] for i in range(x.shape[0])])         #/hr
deltaDilX_total=numpy.gradient(X_total)/(tspan[1]-tspan[0])/X_total
print(mu[1])

## plotting
##simple plot
#plt.plot(tspan, sol)
#plt.title('Volume, Substrate and Biomass')
#plt.xlabel('Time [hr]')
#plt.ylabel('Volume [liter], Concentration [mM]')
#plt.show()
##plot of Volume, Dilution and Growth Rate
fig, ax1 = plt.subplots()
color = 'tab:red'
ax1.set_xlabel('Time [hr]')
ax1.set_ylabel('Volume, [liter]', color=color)
ax1.plot(tspan, V, color=color)
ax1.tick_params(axis='y', labelcolor=color)
# instantiate a second axes that shares the same x-axis
ax2 = ax1.twinx()  
color = 'tab:blue'
ax2.set_ylabel('Dilution and Growth Rate $\mu$ [/hr]', color=color) 
ax2.plot(tspan, mu, color=color, label='growth')
ax2.plot(tspan, dil, color=color, linestyle='--', label='dilution')
ax2.plot(tspan, deltaDilX_total, linestyle=':', label='difference')
ax2.tick_params(axis='y', labelcolor=color)
handles, labels = ax2.get_legend_handles_labels()
ax2.legend(handles, labels)
fig.tight_layout()  # otherwise the right y-label is slightly clipped
plt.show()
##second plot for Biomass and Substrate
fig, ax1 = plt.subplots()
color = 'tab:red'
ax1.set_xlabel('Time [hr]')
ax1.set_ylabel('Biomass Concentration, [mM/liter]', color=color)
ax1.plot(tspan, x, color=color)
ax1.tick_params(axis='y', labelcolor=color)
# instantiate a second axes that shares the same x-axis
ax2 = ax1.twinx()  
color = 'tab:blue'
ax2.set_ylabel('Substrate Concentration [mM]', color=color) 
ax2.plot(tspan, s, color=color)
ax2.tick_params(axis='y', labelcolor=color)
fig.tight_layout()  # otherwise the right y-label is slightly clipped
plt.show()

0.2982891329592172 1/hr
In [ ]:

...for teaching:

In [30]:
#export notebook as html
import os
os.system('jupyter nbconvert --to html Batch_and_Fed_Batch.ipynb')
Out[30]:
0

Below the same as shown above but more complicated in terms of software because the interactive widget for demonstration, requires running kernel as in JupyterLab, does not work in html. You can also just edit the code above. To start run the cell and click the "Run Interact" button!

In [28]:
##load some modules for plotting, math and for the interactive part (slider)
%matplotlib inline
import matplotlib.pyplot as plt
import numpy
from scipy.integrate import odeint
from ipywidgets import widgets
from IPython.display import display
import unyt as u

#Biochemical Model Function
def BioChem_Bmodel(x,t,Ks_d):
    # Model Parameters Values
    MWt_Xyl = 106.2 *u.g/u.mol    # Molecular Weight of Xylene
    Ks_l    = 0.01                # Ks for our cells grown on Xylene (mM)
    max_d   = 0.25  /u.hr         # Maximum Theoretical Deathrate (h-1) 
    max_g   = 0.008 /u.hr         # Maximum Theoretical Growthrate (h-1)
    Yxx     = 2.5                 # Yield of Biomass on Xylene (--)
    #Ks_d    = 58.4                # Ks for death due to Xylene Toxicity (mM)


    # Growth Rate (2. Monod eq.)
    growth = max_g * (x[1]) / (x[1] + Ks_l)
    # Death Rate
    death  = max_d * (x[1]) / (x[1] + Ks_d)

    # Change of Substrate Concentrateion Based on Rate (3. Monod eq.)
    # Use if you are having trouble with substrate < 0 values
    # if x[1] < 0:
    #     rxylene = 0
    # else:
    rxylene = (1/MWt_Xyl) * growth / Yxx

    # Differential for Biomass and Xylene Concentration
    biomass  = (growth - death) * x[0]   # biomass [=] Biomass (mg/L)
    xxylene  = - rxylene * x[0]          # xxylene [=] s, Xylene Concentration (mM)
    return [biomass, xxylene]

##Integration Time Scale
tspan = numpy.linspace(0,340,num=300)

##Initial Conditions
Init_Cond = numpy.array([82.096,1.0]) #Biomass (mg/L), Xylene Concentration (mM)

#function to update with call of widget defined below
def f(Ks_d_move):
    ##Call ODE Solver
    solu = odeint(BioChem_Bmodel,Init_Cond,tspan,args=(Ks_d_move,))
    x_bio, x_xyl = solu.transpose() #Xylene Concentration (mM), Biomass (mg/L)
    ##plotting
    fig, ax1 = plt.subplots()
    color = 'tab:red'
    ax1.set_xlabel('Time [hr]')
    ax1.set_ylabel('Cell Concentration [mg/liter]', color=color)
    ax1.plot(tspan, x_bio, color=color)
    #Analytical solution for unlimited exponential behavior
    initial_growth   =  (x_bio[1]-x_bio[0])/Init_Cond[0]/(tspan[1]-tspan[0])
    plt.ylim(top=x_bio[-1]+10)
    ax1.plot(tspan, Init_Cond[0]*numpy.exp(initial_growth*tspan), color='red', linestyle=':', label='Exponential Growth at Initial Rate')
    handles, labels = ax1.get_legend_handles_labels()
    ax1.legend(handles, labels)
    ax1.tick_params(axis='y', labelcolor=color)
    #instantiate a second axes that shares the same x-axis
    ax2 = ax1.twinx()  
    color = 'tab:blue'
    ax2.set_ylabel('Substrate Concentration [mM]', color=color)
    ax2.plot(tspan, x_xyl, color=color)
    ax2.tick_params(axis='y', labelcolor=color)

interactive_plot = widgets.interact_manual(f, Ks_d_move=widgets.FloatSlider(min=20, max=150, value=58.4, description='Xyl. Tox. Lim.'))
In [ ]: