python - 使用 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]) 

拟合运行顺利，但拟合曲线不满意。大家还有其他意见、建议吗？

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 限制的值时，您必须弄清楚合理的外推是什么样的。