GEKKO Infeasible system of ODE equations of a fed-batch Bioreactor

Question

I am new to GEKKO and also to modeling bioreactors, so I might be missing something obvious.

I have a system of 10 ODEs that describe a fed-batch bioreactor. All constants are given. The picture below shows the expected behavior of this model (extracted from a paper). However, the only feasible solution I found is when Viable Cells Density (XV) = 0, and stays 0 for all time t, or if time T is really small (<20). If a lower boundary >= 0 or initial value is set to XV and t > 20, the system becomes infeasible.

Equations and constants were checked multiple times. I tried giving initial values to my variables, but it didn't work either. I can only think of two problems: I am not initiating variables properly, or I am not using GEKKO properly. Any ideas? Thanks!!

import numpy as np
from gekko import GEKKO
import matplotlib.pyplot as plt

m = GEKKO(remote=False)    # create GEKKO model

#constants 3L continuous fed-batch
KdQ = 0.001        #degree of degradation of glutamine (1/h)
mG = 1.1*10**-10   #glucose maintenance coefficient (mmol/cell/hour)
YAQ = 0.90         #yield of ammonia from glutamine
YLG = 2            #yield of lactate from glucose
YXG = 2.2*10**8    #yield of cells from glucose (cells/mmol)
YXQ = 1.5*10**9    #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 = 2*10**-2  #rate of cell lysis (1/h)
Ci_star = 100      #inhibitor saturation concentration (mM)
qi = 2.5*10**-10   #specific inhibitor production rate (1/h)

#Flow, volume and concentration
Fo = 0.001         #feed-rate (L/h)
Fi = 0.001         #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)

# create GEKKO parameter
t = np.linspace(0,120,121)
m.time = t

XT = m.Var(name='XT')            #total cell density (cells/L)
XV = m.Var(lb=0, name='XV')      #viable cell density (cells/L)
XD = m.Var(name='XD')            #dead cell density (cells/L)
G = m.Var(value = 30, name='G')  #glucose concentration (mM)
Q = m.Var(value = 4.5, name='Q') #glutamine concentration (mM)
L = m.Var(name='L')              #lactate concentration (mM)
A = m.Var(name='A')              #ammonia concentration (mM)
Ci = m.Var(name='Ci')            #inhibitor concentration (mM)
mu = m.Var(name='mu')            #growth rate (1/h)
Kd = m.Var(name='Kd')            #death rate(1/h)

# create GEEKO equations
m.Equation(XT.dt() == mu*XV - Klysis*XD - XT*Fo/V)
m.Equation(XV.dt() == (mu - Kd)*XV - XV*Fo/V)
m.Equation(XD.dt() == Kd*XV - Klysis*XD - XV*Fo/V)
m.Equation(G.dt() == (Fi/V)*SG - (Fo/V)*G + (-mu/YXG - mG)*XV)
m.Equation(Q.dt() == (Fi/V)*SQ - (Fo/V)*Q + (-mu/YXQ - mQ)*XV - KdQ*Q)
m.Equation(L.dt() == -YLG*(-mu/YXG -mG)*XV-(Fo/V)*L)
m.Equation(A.dt() == -YAQ*(-mu/YXQ - mQ)*XV +KdQ*Q-(Fo/V)*A)
m.Equation(Ci.dt() == qi*XV - (Fo/V)*Ci)
m.Equation(mu.dt() == (mumax*G*Q*(Ci_star-Ci)) / (Ci_star*(KG+G)*(KQ+Q)*(L/KL + 1)*(A/KA + 1)))
m.Equation(Kd.dt() == Kdmax*(kmu/(mu+kmu)))

# solve ODE
m.options.IMODE = 4
m.open_folder()
m.solve(display = False)

plt.plot(m.time, XV.value)

Articles that used the exact same model:

1) Master thesis using GEKKO "MODELING OF MAMMALIAN CELL CULTURE" link:

https://search.proquest.com/openview/e4df2d115cbc48ec63235a64b352249c/1.pdf?pq-origsite=gscholar&cbl=18750&diss=y

2) Original paper describing the equations: "Process Model Comparison and Transferability Across Bioreactor Scales and Modes of Operation for a Mammalian Cell Bioprocess"

link: https://sci-hub.tw/10.1002/btpr.1664

3) Paper with a control sytem using that model: "Glucose concentration control of a fed-batch mammalian cell bioprocess using a nonlinear model predictive controller"

link: https://sci-hub.tw/https://doi.org/10.1016/j.jprocont.2014.02.007

Answer 1

There are a couple issues:

• the last two equations are algebraic, not differential. It should be mu==... not mu.dt()==...
• a few of the equations have the potential for divide by zero. Equations such as x.dt() = z/y can be replaced with y * x.dt()==z so that the equation becomes 0 * x.dt() == z if y approaches zero.
• some of the initial conditions aren't set so they use a default value of 0. This is likely creating a zero solution.

I put in some different values and used m.options.COLDSTART=2 to help it find an initial solution. I also used Intermediates to help visualize any terms that are getting big. I put the cell concentrations in units of Millions of cells per liter to help with scaling.

import numpy as np
from gekko import GEKKO
import matplotlib.pyplot as plt

m = GEKKO(remote=False)    # create GEKKO model

#constants 3L continuous 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)

#Flow, volume and concentration
Fo = 0.001         #feed-rate (L/h)
Fi = 0.001         #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)

# create GEKKO parameter
t = np.linspace(0,50,121)
m.time = t

XTMM = m.Var(value=1,name='XT')            #total cell density (MMcells/L)
XVMM = m.Var(value=1,lb=0, name='XV')      #viable cell density (MMcells/L)
XDMM = m.Var(value=1.0,name='XD')          #dead cell density (MMcells/L)
G = m.Var(value = 20, name='G')            #glucose concentration (mM)
Q = m.Var(value = 4.5, name='Q')           #glutamine concentration (mM)
L = m.Var(value=1,name='L')                #lactate concentration (mM)
A = m.Var(value=1.6,name='A')              #ammonia concentration (mM)
Ci = m.Var(value=0.1,name='Ci')            #inhibitor concentration (mM)
mu = m.Var(value=0.1,name='mu')            #growth rate (1/h)
Kd = m.Var(value=0.5,name='Kd')            #death rate(1/h)

# scale back to cells/L from million cells/L
XT = m.Intermediate(XTMM*1e7)
XV = m.Intermediate(XVMM*1e7)
XD = m.Intermediate(XDMM*1e7)

e1 = m.Intermediate((mu*XV - Klysis*XD - XT*Fo/V)/1e7)
e2 = m.Intermediate(((mu - Kd)*XV - XV*Fo/V)/1e7)
e3 = m.Intermediate((Kd*XV - Klysis*XD - XV*Fo/V)/1e7)
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.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()

Answer 2

In the end it is not a programming problem, but a problem reading the equations and correctly translating them.

mu and Kd are not dynamical variables, they are ordinary functions of the state vector (which then only has dimension 8). For such intermediate variables Gekko has the construction function m.Intermediate

mu = m.Intermediate((mumax*G*Q*(Ci_star-Ci)) / (Ci_star*(KG+G)*(KQ+Q)*(L/KL + 1)*(A/KA + 1)), name='mu') #growth rate (1/h)
Kd = m.Intermediate(Kdmax*(kmu/(mu+kmu)))    #death rate(1/h)

With this change and initial value 5e8 for XT and XV, the script gives the plots below that look about what can be found in your cited papers.