6

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!!

Expected model behavior

Equations

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

eyllanesc
  • 235,170
  • 19
  • 170
  • 241
Felipe Mello
  • 395
  • 2
  • 8
  • 2
    You have 10 equations but only set 2 values, are the other initial values automatically set to 0? Did you try to reduce the time interval? You would see that, trying with initial XV=5, the solution diverges around t=20 – Lutz Lehmann Dec 10 '19 at 10:03
  • 2
    mu grows like `0.04*t` and XV like `5*exp(0.02*t^2)`, which is about 3e5 for t=24. Is that the expected dynamic? The reported problem leading to non-convergence is probably that the occurring error scale for these large values deviates too much from the default tolerances. – Lutz Lehmann Dec 10 '19 at 10:47
  • 1
    Thanks for anwering @Dr.LutzLehmann ! I added a picture with the model expected behavior. As you can see, the simulation is longer than 100h, and XV is on a 10^9 scale, so having 3e5 for t=24 is totally acceptable. All variables have initial value close to 0 (except for XV). I tried giving it initial values, but I am not sure if it is needed. I have 10 equations and 10 variables. Shouldn't it be enough? Do you think that loosening tolerances would help? – Felipe Mello Dec 10 '19 at 14:56
  • 1
    Here are some additional tips on model scaling and initialization: https://apmonitor.com/do/index.php/Main/ModelInitialization Systems biology models often have very nonlinear dynamics that need some extra care. I also like to make smaller time steps and solve with IMODE=7 for these types of problems. I'll take a more in-depth look soon. – John Hedengren Dec 10 '19 at 15:29
  • 2
    It is unclear if GEKKO is the right tool for this, I assume this is only the first step leading to parameter estimation or similar. Note that the collocation solver transforms every time step into a system of 10 equations which are then all solved simultaneously using the external optimization tool. While that adds flexibility for adding optimization tasks on top of the ODE, it is doing much more work than any conventional ODE integrator does. Using `odeint` still gives errors over the full interval, but it is much faster over small intervals. – Lutz Lehmann Dec 10 '19 at 15:48
  • 2
    A universal recommendation is to rescale the variables so that the solver ideally mainly sees values in the range 0.01 to 1000. This means removing 1e9 from XV,XD,XT and where necessary changing the constants appropriately. Then the control heuristic based on error estimators of the solver is more within the range of the test examples used in its development. – Lutz Lehmann Dec 10 '19 at 15:49
  • 2
    Could you cite in some way the source for the ODE system? Either as image or as link? I find similar articles, but they seem to have different topics and different conventions in naming and scaling the quantities. – Lutz Lehmann Dec 10 '19 at 15:50
  • 1
    @Dr.LutzLehmann, thanks for your replies. I edited my question and added 3 articles that used this exact same model. One of them is a master thesis that used GEKKO to model it, but the code is proprietary, so I don't know how he did it. – Felipe Mello Dec 10 '19 at 16:25
  • 1
    Thank you very much @JohnHedengren. I will check that link. Let me know if I can help in any way. – Felipe Mello Dec 10 '19 at 16:26
  • 1
    Robert Jackson (MS Thesis author) https://www.linkedin.com/in/robert-jackson-0ab075187/ may also be willing to share the model and other specifics of his implementation. – John Hedengren Dec 10 '19 at 20:53

2 Answers2

1

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.

enter image description here

John Hedengren
  • 12,068
  • 1
  • 21
  • 25
Lutz Lehmann
  • 25,219
  • 2
  • 22
  • 51
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.

results

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()
John Hedengren
  • 12,068
  • 1
  • 21
  • 25