Python GEKKO求解Sedaghat胰岛素通路ODE结果异常排查

API 探路者

2026-8-29

问题描述

我尝试使用Python的GEKKO库实现文献提出的Sedaghat胰岛素信号通路仿真模型，本次采用无反馈结构的系统版本，模型方程与常数取值均参考该文献附录A。当前仿真结果中前6组状态变量的输出曲线表现正常，但x16~x21状态变量的计算结果存在明显异常，与同文献中图6、图7的参考结果不匹配。

我已逐一核对所有建模方程与常数取值、尝试为状态变量添加取值上下界、调整m.options.IMODE与m.options.NODES求解参数，以上操作均未解决异常问题，恳请提供相关排查思路与解决建议。

复现代码

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

def insulin_pathway_CM(insulin, time_interval): 

    m = GEKKO(remote=False)
    m.time = np.linspace(0, time_interval-1, time_interval)

    # 变量初始化
    x1 = m.Param(value=insulin)
    x2 = m.Var(9*1e-13)
    x3 = m.Var()
    x4 = m.Var()
    x5 = m.Var()
    x6 = m.Var(1e-13)
    x7 = m.Var()
    x8 = m.Var()
    x9 = m.Var(1e-12)
    x10 = m.Var()
    x11 = m.Var(1e-13)
    x12 = m.Var()
    x13 = m.Var(0.0031)
    x14 = m.Var(0.994)
    x15 = m.Var(0.0029)
    x16 = m.Var(1)
    x17 = m.Var()
    x18 = m.Var(1)
    x19 = m.Var()
    x20 = m.Var(0.96)
    x21 = m.Var(0.04)
  
    # 常数初始化
    km1 = 0.2
    k1 = 6*1e7
    km2 = 100*km1
    k2 = k1
    km3 =  km1
    k3 = 2500
    km4 = 0.003
    k4 = km4/9
    km4a = 2.1*1e-4
    k4a = 2.1*1e-3
    km5 = 1.67*1e-18
    k5 = 10*km5    # 测试后做的简化
                   # 所有关注取值范围内x6均高于1e-13
    k6 = 0.461
    k7 = 4.16
    km7 = (2.5/7.45)*k7
    km8 = 10
    k8 = km8 * (5/70.775) * 1e12
    km9 = (94/3.1)*1.39
    PI3K = 5*1e-15
    k9basal = (0.31/99.4)*km9
    k9 = m.Var() 
    m.Equation(k9 == (1.39 - k9basal)*(x12/PI3K) + k9basal)
    km10 = 2.77
    k10 = (3.1/2.9)*km10
    km11 = 6.9314718
    k11 = m.Var()
    m.Equation(k11 == (0.1*km11)*(x13-0.31)/2.79)
    km12 = 6.93147
    k12 = m.Var()
    m.Equation(k12 == (0.1*km12)*(x13-0.31)/2.79)
    km13 = 0.167
    k13 = (4/96)*km13
    km14 = 0.001155
    k14 = 96*km14
    
    Effect = m.Var()
    APequil = 100/11
    m.Equation(Effect == (0.2*x17 + 0.8*x19)/APequil)
    k13a = ((4/6) - (4/96))*km13*Effect
    
    SHIP = 1
    PTEN = 1
    PTP = 1
    IRp = 8.97*1e-13

    # 模型方程
    m.Equation(x2.dt() == km1*x3 + km3*PTP*x5 - k1*x1*x2 + km4*x6 - k4*x2)
    m.Equation(x3.dt() == k1*x1*x2 - km1*x3 - k3*x3)
    m.Equation(x4.dt() == k2*x1*x5 -km2*x4 + km4a*x7 -k4a*x4)
    m.Equation(x5.dt() == k3*x3 + km2*x4 - k2*x1*x5 - km3*PTP*x5 + km4a*x8 - k4a*x5)
    m.Equation(x6.dt() == k5 - km5*x6 + k6*PTP*(x7 + x8) + k4*x2 - km4*x6)
    m.Equation(x7.dt() == k4a*x4 - km4a*x7 - k6*PTP*x7)
    m.Equation(x8.dt() == k4a*x5 - km4a*x8 - k6*PTP*x8)
    m.Equation(x9.dt() == km7*PTP*x10 - k7*x9*(x4+x5)/IRp)
    m.Equation(x10.dt() == k7*x9*(x4+x5)/IRp + km8*x12 - (km7*PTP + k8*x11)*x10)
    m.Equation(x11.dt() == km8*x12 - k8*x10*x11)
    m.Equation(x12.dt() == k8*x10*x11 - km8*x12)
    m.Equation(x13.dt() == k9*x14 + k10*x15 - (km9*PTEN + km10*SHIP)*x13)
    m.Equation(x14.dt() == km9*PTEN*x13 - k9*x14) 
    m.Equation(x15.dt() == km10*SHIP*x13 - k10*x15)
    m.Equation(x16.dt() == km11*x17 - k11*x16) 
    m.Equation(x17.dt() == k11*x16 - km11*x17)  
    m.Equation(x18.dt() == km12*x19 - k12*x18)  
    m.Equation(x19.dt() == k12*x18 - km12*x19)
    m.Equation(x20.dt() == km13*x21 - (k13 + k13a)*x20 + k14 - km14*x20)  
    m.Equation(x21.dt() == (k13 + k13a)*x20 - km13*x21) 
  
    m.options.IMODE = 7
    m.options.OTOL  = 1e-8
    m.options.RTOL  = 1e-8
    m.options.NODES = 3
    m.solve(disp=False)
  
    # 绘图
    x45 = []
    for i in range(len(x4.value)):
        x45.append(x4[i]+x5[i])
  
    fig, axs = plt.subplots(4, 2, figsize=(20, 20))
    axs[0, 0].set_title('Free surface receptors (x2)')
    axs[0, 0].plot(m.time, x2)
    axs[0, 1].set_title('Once and twice bound phosphorylated receptors (x4:orange, x5:green)')
    axs[0, 1].plot(m.time, x4, 'tab:orange')
    axs[0, 1].plot(m.time, x5, 'tab:green')
    axs[1, 0].set_title('Surface phosphorylated receptors (x4 + x5)')
    axs[1, 0].plot(m.time, x45, 'tab:red')
    axs[1, 1].set_title('Unphosphorylated and tyrosine phosphorylated IRS-1 (x9:orange, x10:green)')
    axs[1, 1].plot(m.time, x9, 'tab:orange')
    axs[1, 1].plot(m.time, x10, 'tab:green')
    axs[2, 0].set_title('Activated PI 3-kinase (x12)')
    axs[2, 0].plot(m.time, x12, 'tab:orange')
    axs[2, 1].set_title('PI(3,4,5)P3 and PI(3,4,5)P2 (x13:orange, x15:green)')
    axs[2, 1].plot(m.time, x13, 'tab:orange')
    axs[2, 1].plot(m.time, x15, 'tab:green')
    axs[3, 0].set_title('Activated PKC-ζ (x19)')
    axs[3, 0].plot(m.time, x19, 'tab:red')
    axs[3, 1].set_title('Percentage of cell surface GLUT4 (x21)')
    axs[3, 1].plot(m.time, x21, 'tab:cyan')
    return

# 文献对应输入参数
time_interval = 60
insulin = np.zeros(60)
insulin[:15] = 1e-7

insulin_pathway_CM(insulin, time_interval)

排查思路与解决建议

补全所有状态变量的初始值：代码中x3、x4、x5、x7、x8、x10、x12、x17、x19未手动赋值，GEKKO默认取0作为初值，但对照文献的基础稳态计算规则，胰岛素输入为0时上述变量的稳态值不为0，尤其是x17、x19作为下游通路的核心变量，初值为0会直接导致k11、k12、Effect、k13a的计算从错误起点出发，触发数值刚性问题。需要先求解零输入下的系统稳态，将得到的所有状态变量值作为初始值传入对应m.Var()。
修正速率常数的非负约束：k11、k12的计算式包含(x13-0.31)/2.79项，当仿真过程中x13计算值低于0.31时，k11、k12会得到负值，而反应速率常数不能为负，会直接导致x16~x19的微分方程通量方向错误。需要给k11、k12添加下限为0的约束，或按照文献的激活阈值逻辑修正表达式，保证速率常数始终非负。
处理方程组刚性问题：代码中k8的取值量级达到1e12，和其余反应速率常数差10个数量级以上，属于典型的刚性生化系统，默认求解参数容易出现数值振荡。可将m.options.IMODE改为9，先自动求解稳态再启动动态仿真，同时将NODES调整为4提升配置点精度，初始调试阶段可将OTOL、RTOL放宽到1e-6，趋势正确后再逐步收紧容差。
验证守恒循环约束：x16+x17、x18+x19、x20+x21属于共价修饰守恒循环，正常仿真中三者总和应始终等于初始值之和，求解后可打印这三个总和的时序变化，快速定位是否存在方程项漏写、符号写反导致的通量不守恒问题。

内容的提问来源于stack exchange，提问作者nefeli