Python GEKKO ODE 意外结果

Question

我一直在尝试实现 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)

Answer 1

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 年提出了一个描述血糖和胰岛素动力学的最小模型。