我一直在尝试实现 ODE 模型来模拟胰岛素信号通路，正如 supplementary material 中所介绍的那样的 this paper , 使用 python's GEKKO .
实现的模型变体是“Md3”，具有以下方程、常量和初始值:
即使这篇论文没有提供补充代码的结果，人们也希望结果是曲线，而不是我的代码产生的结果:
我已经检查了方程式和常量值，尝试为变量添加下限和上限，并尝试使用 m.options.IMODE 和 m.options.NODES参数，但这些似乎没有帮助。
如有任何建议，我们将不胜感激。
def insulin_pathway_CM(time_interval, insulin=0):
m = GEKKO()
t = np.linspace(0, time_interval-1, time_interval)
m.time = np.insert(t,1,[1e-5,1e-4,1e-3,1e-2,0.1])
# initialization of variables
ins = m.Param(value=insulin)
IR = m.Var(10)
IRp = m.Var()
IRS = m.Var(10)
IRSp = m.Var()
PKB = m.Var(10)
PKBp = m.Var()
GLUT4_C = m.Var(10)
GLUT4_M = m.Var()
# initialization of constants
k1aBasic = 1863.78
k1f = 38668.300000000003
k1b = 40229000
k2f = 401394
k2b = 36704300
k4f = 336782
k4b = 187399
k5Basic = 85530.899999999994
k5f = 11264.799999999999
k5b = 26389900
# equations
m.Equation(IR.dt() == k1b*IRp - (k1f*ins*1000*IR + k1aBasic*IR))
m.Equation(IRp.dt() == -k1b*IRp + k1f*ins*1000*IR + k1aBasic*IR)
m.Equation(IRS.dt() == k2b*IRSp - (k2f*IRp*IRS))
m.Equation(IRSp.dt() == -k2b*IRSp + k2f*IRp*IRS)
m.Equation(PKB.dt() == k4b*PKBp - k4f*IRSp*PKB)
m.Equation(PKBp.dt() == -k4b*PKBp + k4f*IRSp*PKB)
m.Equation(GLUT4_C.dt() == k5b*GLUT4_M - (k5f*PKBp*GLUT4_C + k5Basic*GLUT4_C))
m.Equation(GLUT4_M.dt() == -k5b*GLUT4_M + (k5f*PKBp*GLUT4_C + k5Basic*GLUT4_C))
m.options.IMODE = 7
m.options.OTOL = 1e-8
m.options.RTOL = 1e-8
m.options.NODES = 3
m.solve(disp=False)
# plotting
fig, axs = plt.subplots(4, 2, figsize=(20, 20))
axs[0, 0].plot(m.time, IR)
axs[0, 0].set_title('IR(t)')
axs[0, 1].plot(m.time, IRp)
axs[0, 1].set_title('IRp(t)')
axs[1, 0].plot(m.time, IRS, 'tab:orange')
axs[1, 0].set_title('IRS(t)')
axs[1, 1].plot(m.time, IRSp, 'tab:green')
axs[1, 1].set_title('IRSp(t)')
axs[2, 0].plot(m.time, PKB, 'tab:red')
axs[2, 0].set_title('PKB(t)')
axs[2, 1].plot(m.time, PKBp, 'tab:purple')
axs[2, 1].set_title('PKBp(t)')
axs[3, 0].plot(m.time, GLUT4_C, 'tab:olive')
axs[3, 0].set_title('GLUT4_C(t)')
axs[3, 1].plot(m.time, GLUT4_M, 'tab:cyan')
axs[3, 1].set_title('GLUT4_M(t)')
plt.figure()
plt.show()
return
time_interval = 60
insulin = 100
insulin_pathway_CM(time_interval, insulin)
最佳答案
Lutz Lehmann 是正确的。结束时间为 1e-5 的导数图表明大部分 Action 发生在短时间内。
尝试缩短时间跨度。
m.time = np.linspace(0,1e-5,100)
这显示了快速动态。
方程式可能存在错误，例如单位问题(以天或年为单位的时间？)或者作者忘记为患者包含某种类型的volume 因素，例如血容量或体积。
# equations
V1 = 1e9 # hypothetical volume 1
V2 = 10 # hypothetical volume 2
m.Equation(V1*IR.dt() == k1b*IRp - (k1f*ins*1000*IR + k1aBasic*IR))
m.Equation(V1*IRp.dt() == -k1b*IRp + k1f*ins*1000*IR + k1aBasic*IR)
m.Equation(V1*IRS.dt() == k2b*IRSp - (k2f*IRp*IRS))
m.Equation(V1*IRSp.dt() == -k2b*IRSp + k2f*IRp*IRS)
m.Equation(V2*PKB.dt() == k4b*PKBp - k4f*IRSp*PKB)
m.Equation(V2*PKBp.dt() == -k4b*PKBp + k4f*IRSp*PKB)
m.Equation(V2*GLUT4_C.dt() == k5b*GLUT4_M - (k5f*PKBp*GLUT4_C + k5Basic*GLUT4_C))
m.Equation(V2*GLUT4_M.dt() == -k5b*GLUT4_M + (k5f*PKBp*GLUT4_C + k5Basic*GLUT4_C))
这里是原始时间尺度(0-60 小时)下更合理的动态。
from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt
def insulin_pathway_CM(time_interval, insulin=0):
m = GEKKO()
t = np.linspace(0, time_interval-1, time_interval)
m.time = np.insert(t,1,[1e-5,1e-4,1e-3,1e-2,0.1])
# initialization of variables
ins = m.Param(value=insulin)
IR = m.Var(10)
IRp = m.Var()
IRS = m.Var(10)
IRSp = m.Var()
PKB = m.Var(10)
PKBp = m.Var()
GLUT4_C = m.Var(10)
GLUT4_M = m.Var()
# initialization of constants
k1aBasic = 1863.78
k1f = 38668.300000000003
k1b = 40229000
k2f = 401394
k2b = 36704300
k4f = 336782
k4b = 187399
k5Basic = 85530.899999999994
k5f = 11264.799999999999
k5b = 26389900
# calculate derivatives
d = m.Array(m.Var,8)
m.Equation(d[0] == k1b*IRp - (k1f*ins*1000*IR + k1aBasic*IR))
m.Equation(d[1] == -k1b*IRp + k1f*ins*1000*IR + k1aBasic*IR)
m.Equation(d[2] == k2b*IRSp - (k2f*IRp*IRS))
m.Equation(d[3] == -k2b*IRSp + k2f*IRp*IRS)
m.Equation(d[4] == k4b*PKBp - k4f*IRSp*PKB)
m.Equation(d[5] == -k4b*PKBp + k4f*IRSp*PKB)
m.Equation(d[6] == k5b*GLUT4_M - (k5f*PKBp*GLUT4_C + k5Basic*GLUT4_C))
m.Equation(d[7] == -k5b*GLUT4_M + (k5f*PKBp*GLUT4_C + k5Basic*GLUT4_C))
# equations
V1 = 1e9 # hypothetical volume 1
V2 = 10 # hypothetical volume 2
m.Equation(V1*IR.dt() == k1b*IRp - (k1f*ins*1000*IR + k1aBasic*IR))
m.Equation(V1*IRp.dt() == -k1b*IRp + k1f*ins*1000*IR + k1aBasic*IR)
m.Equation(V1*IRS.dt() == k2b*IRSp - (k2f*IRp*IRS))
m.Equation(V1*IRSp.dt() == -k2b*IRSp + k2f*IRp*IRS)
m.Equation(V2*PKB.dt() == k4b*PKBp - k4f*IRSp*PKB)
m.Equation(V2*PKBp.dt() == -k4b*PKBp + k4f*IRSp*PKB)
m.Equation(V2*GLUT4_C.dt() == k5b*GLUT4_M - (k5f*PKBp*GLUT4_C + k5Basic*GLUT4_C))
m.Equation(V2*GLUT4_M.dt() == -k5b*GLUT4_M + (k5f*PKBp*GLUT4_C + k5Basic*GLUT4_C))
m.options.IMODE = 4
m.options.OTOL = 1e-8
m.options.RTOL = 1e-8
m.options.NODES = 3
m.solve(disp=True)
# plotting
fig, axs = plt.subplots(4, 2, figsize=(20, 20))
axs[0, 0].plot(m.time, IR)
axs[0, 0].set_title('IR(t)')
axs[0, 1].plot(m.time, IRp)
axs[0, 1].set_title('IRp(t)')
axs[1, 0].plot(m.time, IRS, 'tab:orange')
axs[1, 0].set_title('IRS(t)')
axs[1, 1].plot(m.time, IRSp, 'tab:green')
axs[1, 1].set_title('IRSp(t)')
axs[2, 0].plot(m.time, PKB, 'tab:red')
axs[2, 0].set_title('PKB(t)')
axs[2, 1].plot(m.time, PKBp, 'tab:purple')
axs[2, 1].set_title('PKBp(t)')
axs[3, 0].plot(m.time, GLUT4_C, 'tab:olive')
axs[3, 0].set_title('GLUT4_C(t)')
axs[3, 1].plot(m.time, GLUT4_M, 'tab:cyan')
axs[3, 1].set_title('GLUT4_M(t)')
plt.figure()
for i in range(8):
plt.plot(m.time,d[i].value)
plt.show()
return
time_interval = 60
insulin = 100
insulin_pathway_CM(time_interval, insulin)
作为引用，这里有一个 simplified blood glucose response model (类似于伯格曼模型)供大家引用。 Richard Bergman 和 Claudio Cobelli 于 1979 年提出了一个描述血糖和胰岛素动力学的最小模型。
关于Python GEKKO ODE 意外结果，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/73435108/