Python GEKKO ODE résultats erronés

Question

J'ai essayé de mettre en œuvre le modèle Sedaghat pour la simulation de la voie de signalisation de l'insuline, telle qu'elle est présentée. aquí en utilisant la fonction GEKKO . Le système modélisé est celui sans rétroaction et les équations du modèle et les valeurs des constantes se trouvent dans l'annexe A. Alors que les 6 premiers diagrammes de mes résultats semblent bien ressortir, les autres états du modèle (x16-x21) semblent un peu problématiques (fig. 6,7 du même document utilisé pour la comparaison).

J'ai vérifié les équations et les valeurs constantes, j'ai essayé d'ajouter des limites inférieures et supérieures aux variables et j'ai expérimenté avec la fonction m.options.IMODE y m.options.NODES mais cela n'a pas semblé aider.

Tout conseil serait très apprécié.

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)

    # initialization of variables 
    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)

    # initialization of constants
    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    # simplification made after testing  
                   # x6 seems to be above 1e-13 for all values of interest
    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

    # equations
    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)

    # plotting
    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

# input as in the paper
time_interval = 60
insulin = np.zeros(60)
insulin[:15] = 1e-7

insulin_pathway_CM(insulin, time_interval)

Answer 1

Voici quelques suggestions :

• Remplacer certains des états par Intermediate type. Cela donne la même solution mais le calcul est plus rapide.

  #k9 = m.Var() 
  #m.Equation(k9 == (1.39 - k9basal)*(x12/PI3K) + k9basal)
  k9 = m.Intermediate(1.39 - k9basal)*(x12/PI3K) + k9basal)

• Ajouter quelques petits pas de temps au début en cas de dynamique rapide précoce qui ne suit pas bien avec un pas de temps de 1. Cela n'a pas eu d'effet, mais c'est une bonne chose à essayer.

  m = GEKKO(remote=False)
  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])
  insulin = np.insert(insulin,1,[insulin[0]]*5)

• Examinez toutes les solutions non intuitives telles que x17 y x19 comme des nombres négatifs.

Équations pour x16 à x21 ne semblent pas influencer les autres parties du modèle et la seule entrée est k11 qui dépend également de x13 , x14 et x15 .

    km11 = 6.9314718
    k11 = m.Intermediate((0.1*km11)*(x13-0.31)/2.79)
    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) 

Il peut y avoir une erreur de modèle dans ces équations, voire dans le document original (cela arrive parfois). Se concentrer sur les paires x16 y x17 il apparaît que x17 est juste le changement dans x16 . Il en va de même pour x18 y x19 . À première vue, je ne voyais pas de différences évidentes dans les équations entre le code source Python et l'article.

Voici le script complet avec quelques changements avec les Intermédiaires, 2 parcelles ajoutées, et les points de temps supplémentaires au début.

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)
    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])
    print(m.time)
    insulin = np.insert(insulin,1,[insulin[0]]*5)

    # initialization of variables 
    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)

    # initialization of constants
    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    # simplification made after testing  
                   # x6 seems to be above 1e-13 for all values of interest
    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.Intermediate((1.39 - k9basal)*(x12/PI3K) + k9basal)
    km10 = 2.77
    k10 = (3.1/2.9)*km10
    km11 = 6.9314718
    k11 = m.Intermediate((0.1*km11)*(x13-0.31)/2.79)
    km12 = 6.93147
    k12 = m.Intermediate((0.1*km12)*(x13-0.31)/2.79)
    km13 = 0.167
    k13 = (4/96)*km13
    km14 = 0.001155
    k14 = 96*km14

    APequil = 100/11
    Effect = m.Intermediate((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

    # equations
    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-10
    m.options.RTOL  = 1e-10
    m.options.NODES = 3
    m.solve(disp=False)

    # plotting
    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')

    plt.figure()
    for ni in range(16,22):
        nx = 'x'+str(ni); x = eval(nx)
        plt.subplot(3,2,ni-15)
        plt.plot(m.time,x.value,label=nx); plt.legend()

    plt.figure()
    plt.subplot(2,1,1)
    plt.plot(m.time,k11.value,label='k11')
    plt.subplot(2,1,2)
    plt.plot(m.time,x13.value,label='x13')

    return

# input as in the paper
time_interval = 60
insulin = np.zeros(60)
insulin[:15] = 1e-7

insulin_pathway_CM(insulin, time_interval)

plt.show()