使用 lmfit 的 2 组织 3 隔室模型拟合数据

Question

我正在尝试将我的实验数据与一个非常知名的药代动力学模型进行拟合。方程组相当复杂：

dC1/dt = k1*Cp - (k2+k3)*C1 + k4*C2
dC2/dt = k3*C1 - k4*C2
Ctissue = (1-vB)*(C1 + C2) + vB*Cp

vB 是一个常数，Cp 是一个数组（通过测量已知的因变量），k1、k2、k3、k4 是不同隔间之间的动力学常数，并且是我希望从拟合中获得的参数。 Ctissue 是我想用真实数据拟合的。 C1 和 C2 是我应该能够在执行拟合后计算的两个数组。有一个商业软件（PKIN）可以计算这个方程组，所以我确信这是可能的，但我不知道如何用 python 实现它。

这是我的代码

tini = np.array([  15.,   45.,   75.,  120.,  180.,  240.,  300.,  360.,  450.,
        570.,  690.,  810.,  930., 1080., 1260., 1440., 1650., 1890.,
       2130., 2400., 2700., 3000., 3300., 3525.])

Ctissue = np.array([  1.00229754,  25.06505484,  60.0265695 ,  82.87576127,
        68.07901198,  67.10795788,  81.42071546,  81.05644343,
       100.6740041 ,  90.43091176, 111.7861611 , 110.3851624 ,
       116.4682562 , 126.7390119 , 133.8460856 , 132.8657165 ,
       145.3951029 , 141.4012821 , 156.7317122 , 159.8293774 ,
       163.609847  , 175.7823822 , 168.5340708 , 171.5013387 ])

Cp = np.array([ 13.99461153, 559.5563251 , 914.7457277 , 782.4498718 ,
       574.7527458 , 521.4668956 , 412.9772775 , 421.5475443 ,
       403.2700613 , 368.6237412 , 355.8405377 , 340.0395723 ,
       306.9848032 , 295.0192494 , 295.0294368 , 240.9861338 ,
       245.9420067 , 217.3042524 , 229.6231028 , 196.4563327 ,
       190.8358096 , 190.161142  , 182.2021123 , 169.1384708 ])

vB = 0.05

# initial conditions
x10 = 0.1
x20 = 0.1
y0 = [x10, x20]
guess = [0.1,0.1,0.1,0.1]

import scipy as sp
import numpy as np
import pandas as pd
import matplotlib as mpl
from matplotlib import pyplot as plt
import math as m
from scipy.integrate import odeint
from lmfit import minimize, Parameters, Parameter, report_fit
from scipy.interpolate import interp1d

def myCp( t ):
    cp = interp1d( tini, Cp )
    if np.all(t < tini[0]):  
        out = Cp[0] 
    elif np.all(t > tini[-1]):
        out = 0 
    else:
        out = cp( t )
    return out

def f(y, t, paras):
#define differential equations
    x1 = y[0]
    x2 = y[1]

    try:
        k1 = paras['k1'].value
        k2 = paras['k2'].value
        k3 = paras['k3'].value 
        k4 = paras['k4'].value

    except KeyError:
        k1, k2, k3, k4 = paras
    f1 = k1*myCp( t ) - (k2+k3)*x1 + k4*x2
    f2 = k3*x1 - k4*x2
    return [f1, f2]

def g(t, x0, paras):
    x = odeint(f, x0, t, args=(paras,))
    return x

def tis2comp3(t, paras):
    x0 = params['x10'].value, params['x20'].value
    model = g(t, x0, paras)
    x1_model = model[:, 0]
    x2_model = model[:, 1]
    Ct = (1-vB)*(x1_model + x2_model) + vB*myCp( t )
    return Ct

def residual(paras, t, data):
    Ct = tis2comp3(t, params)
    return (Ct - data).ravel()

# set parameters
params = Parameters()
params.add('x10', value=x10, vary=False)
params.add('x20', value=x20, vary=False)
params.add('k1', value=guess[0], min=0.0001, max=2.)
params.add('k2', value=guess[1], min=0.0001, max=2.)
params.add('k3', value=guess[2], min=0.0001, max=2.)
params.add('k4', value=guess[3], min=0.0001, max=2.)

# fit model
result = minimize(residual, params, args=(tini, Ctissue), method='leastsq')  # leastsq nelder
# check results of the fit
xfit = np.linspace(15., 3525., 100)
yfit = tis2comp3(xfit, result.params)

#plot the final optimization results
figopt = plt.figure(figsize=(10,6))
lineini, = plt.plot(tini,Ctissue, 'b', linestyle='none', marker='o', markersize=7, label='data')
lineopt, = plt.plot(xfit,yfit, 'r-', label='optimized curve')
plt.legend(handles=[lineini,lineopt]) 

拟合运行顺利，但拟合曲线不满意。 各位有没有其他的cmet，建议？

Answer 1

您的问题似乎是您将使用Cp 的离散图片与准连续图片混合使用odeint 和f 如果您在拟合迭代中查看f 的第一个输出，您将看到第一个，即f1，是一个数组，而第二个是一个数字。所以存在概念上的错误。

将您的 f 更改为如下内容：

import numpy as np
from matplotlib import pyplot as plt
from scipy.integrate import odeint
from scipy.interpolate import interp1d


tini = np.array([  
        15.,   45.,   75.,  120.,  180.,  240.,  300.,  360.,  450.,
        570.,  690.,  810.,  930., 1080., 1260., 1440., 1650., 1890.,
       2130., 2400., 2700., 3000., 3300., 3525.])

Ctissue = np.array([  
        1.00229754,   25.06505484,  60.0265695 ,  82.87576127,
        68.07901198,  67.10795788,  81.42071546,  81.05644343,
       100.6740041,   90.43091176, 111.7861611 , 110.3851624 ,
       116.4682562,  126.7390119,  133.8460856 , 132.8657165 ,
       145.3951029,  141.4012821,  156.7317122 , 159.8293774 ,
       163.609847,   175.7823822,  168.5340708 , 171.5013387 ])

Cp = np.array([ 
        13.99461153, 559.5563251 , 914.7457277 , 782.4498718 ,
       574.7527458 , 521.4668956 , 412.9772775 , 421.5475443 ,
       403.2700613 , 368.6237412 , 355.8405377 , 340.0395723 ,
       306.9848032 , 295.0192494 , 295.0294368 , 240.9861338 ,
       245.9420067 , 217.3042524 , 229.6231028 , 196.4563327 ,
       190.8358096 , 190.161142  , 182.2021123 , 169.1384708 ])

vB = 0.05

def myCp( t ):
    cp = interp1d( tini, Cp )
    if t < tini[0]: # does this makes sense 
        out = Cp[0] # may require to be refined
    elif t > tini[-1]:
        out = Cp[-1] # same here.
    else:
        out = cp( t )
    return out

def f( y, t, paras ):
#define differential equations
    x1 = y[0]
    x2 = y[1]
    k1, k2, k3, k4 = paras
    f1 = k1 * myCp( t ) - ( k2 + k3 ) * x1 + k4 * x2
    f2 = k3 * x1 - k4 * x2
    return [ f1, f2 ]
paras=[.1, .12, .14, .15 ]
sol = odeint( f, [ .1, .2], tini, args=( paras, ) )

print sol[ :, 0 ]
print sol[ :, 1 ]

应该可以。

当odeint 检查超出t 限制的值时，您必须弄清楚合理的推断是什么样的。