Reputation: 395
I have a model with 10 equations that describe a fed-batch bioreactor. Basically, every 24h "food" (Glucose and other components) is added to the system. My problem is that this feeding procedure is currently being modeled as the flow-rate (L/H) over two time steps (delta_T), instead of a single discrete food addition (delta_T = 0).
This is what the glucose equation looks like:
e4 = m.Intermediate(**(Fi/V)*SG** - (Fo/V)*G + (-mu/YXG - mG)*XV)
m.Equation(G.dt() == e4)
Where G
is the Glucose Concentration in the bioreactor (mM), Fi
is input feed-rate (L/h), V
is the bioreactor volume (L), and SG
is the glucose concentration in the feed (mM).
I managed to make the system feasible by calling this delta_T = 0.2 hours
, in other words, the Glucose concentration continuously (instead of instantly) rises from G1
at time t1
to G2
at time t1 + 0.2h
. If I try to lower this delta_T
, the system shows a very bizarre behavior.
Time array looks like this: [..., 19.5, 20.0, 20.5, 21.0, 21.5, 22.0, 22.5, 23.0, 23.5, 24.0, 24.2, 24.5, 25.0, ...].
It changes in steps of 0.5h, and every 24h, when I add glucose to the bioreactor, I forced the next time step to be only 0.2 long, instead of 0.5. I want this delta to be 0.
My feed-rate looks like this:
[..., 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0,...]
Is there anyway that I can make this feeding process discrete? The full code is shown below. Thanks!!
import numpy as np
from gekko import GEKKO
import matplotlib.pyplot as plt
import pandas as pd
m = GEKKO(remote=False) # create GEKKO model
#constants 3L fed-batch
KdQ = 0.001 #degree of degradation of glutamine (1/h)
mG = 1.1e-10 #glucose maintenance coefficient (mmol/cell/hour)
YAQ = 0.90 #yield of ammonia from glutamine
YLG = 2 #yield of lactate from glucose
YXG = 2.2e8 #yield of cells from glucose (cells/mmol)
YXQ = 1.5e9 #yield of cells from glutamine (cells/mmol)
KL = 150 #lactate saturation constant (mM)
KA = 40 #ammonia saturation constant (mM)
Kdmax = 0.01 #maximum death rate (1/h)
mumax = 0.044 #maximum growth rate (1/h)
KG = 1 #glucose saturation constant (mM)
KQ = 0.22 #glutamine saturation constant (mM)
mQ = 0 #glutamine maintenance coefficient (mmol/cell/hour)
kmu = 0.01 #intrinsic death rate (1/h)
Klysis = 2e-2 #rate of cell lysis (1/h)
Ci_star = 100 #inhibitor saturation concentration (mM)
qi = 2.5e-10 #specific inhibitor production rate (1/h)
N_HOURS = 150 #Number of hours of the experiment
TIME_STEP = 0.5
feed_small_delta_t = 0.2
#create time array. It will be from 0 to N_HOURS, with time step TIME_STEP,
#and every 24h, there will be a feed_small_delta_t
t = []
for i in range(int(N_HOURS/TIME_STEP +1)):
t_value = i*TIME_STEP
t.append(t_value)
if t_value%24 == 0:
t.append(t_value + feed_small_delta_t)
m.time = t
#Create input feed-rate array
Fi = np.zeros(len(t))
for i in range(1,len(t)):
if t[i]%(24) == 0:
Fi[i] = 0.1/feed_small_delta_t #it is divided by feed_small_delta_t so volume added is constant.
#Flow, volume and concentration
Fi = m.Param(Fi) #input feed-rate (L/h)
Fo = 0 #output feed-rate (L/h)
V = 3 #volume (L)
SG = 653 #glucose concentration in the feed (mM)
SQ = 58.8 #glutamine concentration in the feed (mM)
XTMM = m.Var(value=2,lb=-0.0000,name='XT') #total cell density (MMcells/L)
XVMM = m.Var(value=2,lb=-0.0000, name='XV') #viable cell density (MMcells/L)
XDMM = m.Var(value=0,lb=-0.0000,name='XD') #dead cell density (MMcells/L)
G = m.Var(value = 20,lb=-0.0000,name='G') #glucose concentration (mM)
Q = m.Var(value = 3.8,lb=-0.0000, name='Q') #glutamine concentration (mM)
L = m.Var(value=0.1,lb=-0.0000,name='L') #lactate concentration (mM)
A = m.Var(value=1.8,lb=-0.0000,name='A') #ammonia concentration (mM)
Ci = m.Var(lb=-0.0000, name='Ci') #inhibitor concentration (mM)
mu = m.Var(lb=-0.0000, name='mu') #growth rate (1/h)
Kd = m.Var(lb=-0.0000, name='Kd') #death rate(1/h)
# scale back to cells/L from million cells/L
XT = m.Intermediate(XTMM*1e8)
XV = m.Intermediate(XVMM*1e8)
XD = m.Intermediate(XDMM*1e8)
e1 = m.Intermediate((mu*XV - Klysis*XD - XT*Fo/V)/1e8)
e2 = m.Intermediate(((mu - Kd)*XV - XV*Fo/V)/1e8)
e3 = m.Intermediate((Kd*XV - Klysis*XD - XV*Fo/V)/1e8)
e4 = m.Intermediate((Fi/V)*SG - (Fo/V)*G + (-mu/YXG - mG)*XV)
e5 = m.Intermediate((Fi/V)*SQ - (Fo/V)*Q + (-mu/YXQ - mQ)*XV - KdQ*Q)
e6 = m.Intermediate(-YLG*(-mu/YXG -mG)*XV-(Fo/V)*L)
e7 = m.Intermediate(-YAQ*(-mu/YXQ - mQ)*XV +KdQ*Q-(Fo/V)*A)
e8 = m.Intermediate(qi*XV - (Fo/V)*Ci)
e9a = m.Intermediate((Ci_star*(KG+G)*(KQ+Q)*(L/KL + 1)*(A/KA + 1)))
e9b = m.Intermediate((mumax*G*Q*(Ci_star-Ci)))
e10a = m.Intermediate((mu+kmu))
e10b = m.Intermediate(Kdmax*kmu)
# create GEEKO equations
m.Equation(XTMM.dt() == e1)
m.Equation(XVMM.dt() == e2)
m.Equation(XDMM.dt() == e3)
m.Equation(G.dt() == e4)
m.Equation(Q.dt() == e5)
m.Equation(L.dt() == e6)
m.Equation(A.dt() == e7)
m.Equation(Ci.dt() == e8)
m.Equation(e9a * mu == e9b)
m.Equation(e10a * Kd == e10b)
# solve ODE
m.options.IMODE = 4
m.options.SOLVER = 1
m.options.NODES = 3
m.options.COLDSTART = 2
#m.open_folder()
m.solve(display=False)
plt.figure()
plt.subplot(3,1,1)
plt.plot(m.time, XV.value,label='XV')
plt.plot(m.time, XT.value,label='XT')
plt.plot(m.time, XD.value,label='XD')
plt.legend()
plt.subplot(3,1,2)
plt.plot(m.time, G.value,label='G')
plt.plot(m.time, Q.value,label='Q')
plt.plot(m.time, L.value,label='L')
plt.plot(m.time, A.value,label='A')
plt.plot(m.time, Ci.value,label='Ci')
plt.legend()
plt.subplot(3,1,3)
plt.plot(m.time, mu.value,label='mu')
plt.plot(m.time, Kd.value,label='Kd')
plt.legend()
plt.xlabel('Time (hr)')
plt.figure()
plt.plot(m.time, e1.value,'r-.',label='eqn1')
plt.plot(m.time, e2.value,'g:',label='eqn2')
plt.plot(m.time, e3.value,'b:',label='eqn3')
plt.plot(m.time, e4.value,'b--',label='eqn4')
plt.plot(m.time, e5.value,'y:',label='eqn5')
plt.plot(m.time, e6.value,'m--',label='eqn6')
plt.plot(m.time, e7.value,'b-.',label='eqn7')
plt.plot(m.time, e8.value,'g--',label='eqn8')
plt.plot(m.time, e9a.value,'r:',label='eqn9a')
plt.plot(m.time, e9b.value,'r--',label='eqn9b')
plt.plot(m.time, e10a.value,'k:',label='eqn10a')
plt.plot(m.time, e10b.value,'k--',label='eqn10b')
plt.legend()
plt.show()
Upvotes: 2
Views: 237
Reputation: 14346
Your strategy to feed glucose as a pulse is a good method to have a discontinuous input. The problem with a discrete jump in glucose concentration is that there is a glucose derivative term as equation 4: m.Equation(G.dt() == e4)
. If the dG/dt
term changes over a very short amount of time then the derivative term gets very large.
One strategy to deal with the large derivatives at discrete points is to use m.options.NODES=2
to avoid problems with the additional internal nodes with orthogonal collocation on finite elements. With no internal nodes, you may need to increase the number of time points to maintain accuracy for the integration. This allows a very short-duration impulse input of glucose to the batch reactor such as 3.6 seconds
for the addition.
feed_small_delta_t = 0.001 # 3.6 seconds
The index for the feed input is off by one so Fi[i+1]
should be where the impulse is applied, not Fi[i]
.
Fi = np.zeros(len(t))
for i in range(1,len(t)):
if t[i]%(24) == 0:
Fi[i+1] = 0.1/feed_small_delta_t
This gives similar results to what you had before but with a shorter input impulse for the discrete event of adding additional sugar to the batch.
Full Script with Modifications
import numpy as np
from gekko import GEKKO
import matplotlib.pyplot as plt
import pandas as pd
m = GEKKO(remote=False) # create GEKKO model
#constants 3L fed-batch
KdQ = 0.001 #degree of degradation of glutamine (1/h)
mG = 1.1e-10 #glucose maintenance coefficient (mmol/cell/hour)
YAQ = 0.90 #yield of ammonia from glutamine
YLG = 2 #yield of lactate from glucose
YXG = 2.2e8 #yield of cells from glucose (cells/mmol)
YXQ = 1.5e9 #yield of cells from glutamine (cells/mmol)
KL = 150 #lactate saturation constant (mM)
KA = 40 #ammonia saturation constant (mM)
Kdmax = 0.01 #maximum death rate (1/h)
mumax = 0.044 #maximum growth rate (1/h)
KG = 1 #glucose saturation constant (mM)
KQ = 0.22 #glutamine saturation constant (mM)
mQ = 0 #glutamine maintenance coefficient (mmol/cell/hour)
kmu = 0.01 #intrinsic death rate (1/h)
Klysis = 2e-2 #rate of cell lysis (1/h)
Ci_star = 100 #inhibitor saturation concentration (mM)
qi = 2.5e-10 #specific inhibitor production rate (1/h)
N_HOURS = 150 #Number of hours of the experiment
TIME_STEP = 0.5
feed_small_delta_t = 0.001
#create time array. It will be from 0 to N_HOURS, with time step TIME_STEP,
#and every 24h, there will be a feed_small_delta_t
t = []
for i in range(int(N_HOURS/TIME_STEP +1)):
t_value = i*TIME_STEP
t.append(t_value)
if t_value%24 == 0:
t.append(t_value + feed_small_delta_t)
m.time = t
#Create input feed-rate array
Fi = np.zeros(len(t))
for i in range(1,len(t)):
if t[i]%(24) == 0:
Fi[i+1] = 0.1/feed_small_delta_t #it is divided by feed_small_delta_t so volume added is constant.
#Flow, volume and concentration
Fi = m.Param(Fi) #input feed-rate (L/h)
Fo = 0 #output feed-rate (L/h)
V = 3 #volume (L)
SG = 653 #glucose concentration in the feed (mM)
SQ = 58.8 #glutamine concentration in the feed (mM)
XTMM = m.Var(value=2,lb=-0.0000,name='XT') #total cell density (MMcells/L)
XVMM = m.Var(value=2,lb=-0.0000, name='XV') #viable cell density (MMcells/L)
XDMM = m.Var(value=0,lb=-0.0000,name='XD') #dead cell density (MMcells/L)
G = m.Var(value = 20,lb=-0.0000,name='G') #glucose concentration (mM)
Q = m.Var(value = 3.8,lb=-0.0000, name='Q') #glutamine concentration (mM)
L = m.Var(value=0.1,lb=-0.0000,name='L') #lactate concentration (mM)
A = m.Var(value=1.8,lb=-0.0000,name='A') #ammonia concentration (mM)
Ci = m.Var(lb=-0.0000, name='Ci') #inhibitor concentration (mM)
mu = m.Var(lb=-0.0000, name='mu') #growth rate (1/h)
Kd = m.Var(lb=-0.0000, name='Kd') #death rate(1/h)
# scale back to cells/L from million cells/L
XT = m.Intermediate(XTMM*1e8)
XV = m.Intermediate(XVMM*1e8)
XD = m.Intermediate(XDMM*1e8)
e1 = m.Intermediate((mu*XV - Klysis*XD - XT*Fo/V)/1e8)
e2 = m.Intermediate(((mu - Kd)*XV - XV*Fo/V)/1e8)
e3 = m.Intermediate((Kd*XV - Klysis*XD - XV*Fo/V)/1e8)
e4 = m.Intermediate((Fi/V)*SG - (Fo/V)*G + (-mu/YXG - mG)*XV)
e5 = m.Intermediate((Fi/V)*SQ - (Fo/V)*Q + (-mu/YXQ - mQ)*XV - KdQ*Q)
e6 = m.Intermediate(-YLG*(-mu/YXG -mG)*XV-(Fo/V)*L)
e7 = m.Intermediate(-YAQ*(-mu/YXQ - mQ)*XV +KdQ*Q-(Fo/V)*A)
e8 = m.Intermediate(qi*XV - (Fo/V)*Ci)
e9a = m.Intermediate((Ci_star*(KG+G)*(KQ+Q)*(L/KL + 1)*(A/KA + 1)))
e9b = m.Intermediate((mumax*G*Q*(Ci_star-Ci)))
e10a = m.Intermediate((mu+kmu))
e10b = m.Intermediate(Kdmax*kmu)
# create GEEKO equations
m.Equation(XTMM.dt() == e1)
m.Equation(XVMM.dt() == e2)
m.Equation(XDMM.dt() == e3)
m.Equation(G.dt() == e4)
m.Equation(Q.dt() == e5)
m.Equation(L.dt() == e6)
m.Equation(A.dt() == e7)
m.Equation(Ci.dt() == e8)
m.Equation(e9a * mu == e9b)
m.Equation(e10a * Kd == e10b)
# solve ODE
m.options.IMODE = 4
m.options.SOLVER = 1
m.options.NODES = 2
m.options.COLDSTART = 2
#m.open_folder()
m.solve(display=False)
plt.figure()
plt.subplot(3,1,1)
plt.plot(m.time, XV.value,label='XV')
plt.plot(m.time, XT.value,label='XT')
plt.plot(m.time, XD.value,label='XD')
plt.legend()
plt.subplot(3,1,2)
plt.plot(m.time, G.value,label='G')
plt.plot(m.time, Q.value,label='Q')
plt.plot(m.time, L.value,label='L')
plt.plot(m.time, A.value,label='A')
plt.plot(m.time, Ci.value,label='Ci')
plt.legend()
plt.subplot(3,1,3)
plt.plot(m.time, mu.value,label='mu')
plt.plot(m.time, Kd.value,label='Kd')
plt.legend()
plt.xlabel('Time (hr)')
plt.figure()
plt.plot(m.time, e1.value,'r-.',label='eqn1')
plt.plot(m.time, e2.value,'g:',label='eqn2')
plt.plot(m.time, e3.value,'b:',label='eqn3')
plt.plot(m.time, e4.value,'b--',label='eqn4')
plt.plot(m.time, e5.value,'y:',label='eqn5')
plt.plot(m.time, e6.value,'m--',label='eqn6')
plt.plot(m.time, e7.value,'b-.',label='eqn7')
plt.plot(m.time, e8.value,'g--',label='eqn8')
plt.plot(m.time, e9a.value,'r:',label='eqn9a')
plt.plot(m.time, e9b.value,'r--',label='eqn9b')
plt.plot(m.time, e10a.value,'k:',label='eqn10a')
plt.plot(m.time, e10b.value,'k--',label='eqn10b')
plt.legend()
plt.show()
Upvotes: 1