Я пытаюсь согласовать свои экспериментальные данные с очень известной моделью в фармакокинетике. Система уравнений довольно сложна:
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])
Подгонка прошла гладко, но подогнанная кривая неудовлетворена.
Ребята, у вас есть другие комментарии, предложения?