GEKKO分批补料生物反应器的ODE方程的不可行系统

Question

我对GEKKO还是生物反应器建模都是新手，所以我可能会遗漏一些明显的东西。

我有一个由10个ODE组成的系统，用于描述补料分批生物反应器。给出了所有常数。下图显示了该模型的预期行为（摘自论文）。但是，我发现的唯一可行解决方案是活细胞密度（XV）= 0，并且在整个时间t内都保持为0，或者如果时间T确实很小（<20）。如果下边界> = 0或初始值设置为XV且t> 20，则系统将变得不可行。

多次检查方程式和常数。我尝试给变量赋初始值，但是也没有用。我只能想到两个问题：我没有正确地初始化变量，或者我没有正确地使用GEKKO。有任何想法吗？谢谢！！

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)
使用完全相同的模型的文章：

[1）使用GEKKO的“哺乳动物细胞培养模型”的硕士学位论文

链接：

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

2）描述该方程的原始论文：“过程模型的比较和跨生物反应器规模的转移性和哺乳动物细胞生物过程的操作模式“

链接：https://sci-hub.tw/10.1002/btpr.1664

3）具有使用该模型的控制系统的论文：“使用非线性模型预测控制器对分批补料的哺乳动物细胞生物过程的葡萄糖浓度控制”

链接：https://sci-hub.tw/https://doi.org/10.1016/j.jprocont.2014.02.007

我对GEKKO以及对生物反应器建模都是新手，所以我可能缺少明显的东西。我有一个由10个ODE组成的系统，用于描述补料分批生物反应器。给出了所有常数。图片...

Answer 1

有几个问题：

Answer 2

最后，这不是编程问题，而是阅读方程式并正确翻译它们的问题。