于广华1,裔照国2,陈国忠1,王宁1 (1.盐城卫生职业技术学院,江苏盐城
224006;2.盐城市第三人民医院,江苏盐城
224001)
[摘要] 目的:应用Matlab程序模拟药动学模型。方法:分析药动学模型的构建,用Matlab的命令(函数)建立M函数与文 件,做药动学模型的模拟,还可使用曲线拟合工具箱以及Simulink仿真等工具箱。结果:Matlab可做房室模型、生理动力学模型等的模拟及仿真,模拟简捷、快速、稳定并可视化。结论:Matlab功能强大,程序设计自由,源代码开放,可构建药动学专用软件,是药物动力学研究的有力工具。[关键词] 药动学;模型;Matlab;模拟 [中图分类号]R969.1 [文献标识码]B [文章编号]100125213(2008)2221956204
药动学系采用数学模型,定量描述药物体内过程的量时变化。模型建立、实验数据拟合、参 数估算以及药物体内过程的模拟,会涉及复杂的数学运算。如何选择正确的处理方
式,做数据的正确处理,从而得到准确有效的结论,是药动学研究中的重要问题,计算机虚拟技术尤为关注。这些研究与计算技术的应用密切相关,同时也是向纵深发展的需要,因而离不开计算软件的辅助。国内外研究者开发了许多通用及专用计算软件[122],用于药动学模型的模拟。Matlab是MathWorks公司推出的一套高性能的科学计算软件,其强大的数学运算能力、绘图功能、数据可视化及高度集成的语言,在科学与工程领域的应用越来越广,是科学计算通行的权威工具[324]。本文介绍Matlab7.1在药动学模型模拟中的应用。1 Matlab的模型拟合
Matlab以数组和矩阵为基础,具有丰富的库函数,介绍了许多线性(非线性)拟合方法,如定义非线性方程命令in2line和拟合曲线方程命令nlinfit、curvefit、lsqcurvefit、lsqnon2lin、多项式拟合函数polyfit等,可用于药动学线性(非线性)曲线的拟合。插值函数interpl、3次样条插值函数spline,求常(偏)微分方程的初值函数ode23、ode45,符号极限函数limit、导数函数diff、拉普拉斯变换函数laplace等等,这些函数在药动学模型的微积分中运算非常有用,如用极值函数fmin和fmins求血药浓度的谷浓度或峰浓度,函数quad、quad8计算药时曲线下的面积AUC。回归分析是常用的统计方法,函数polyfit可用于线性回归,ieastsq函数可做非线性回归,
回归分析各参数估计量的分布、回归方程的显著性检验、回归系数的假设检验、置信区间估计以及回归模型预测,Matlab都可以得到完美的实现。
危机管理理论Matlab提供了一系列绘图及图形控制函数(命令),例如线性坐标函数plot、对数坐标函数semilogy、半对数坐标函数loglog以及各种二维、三维统计分析图函数、图函数等。Matlab在图形窗口中提供了可视化的图形编辑工具,以完成各种图形对象的编辑处理,可供制作药时曲线、剂量—效应曲线、药物效应—时间曲线以及数据统计图等图表,实现数据的可视化,使得药动学数据便于分析与讨论。
面向具体应用的强大工具箱和众多模块集是Matlab的重大特之一,如曲线拟合工具箱(CurveFittingTool,cftool)、数理统计工具箱(Toolbox)、分析与综合工具箱(A2nalysisandSynthesisToolbox)、优化工具箱(OptimizationToolbox)、仿真工具箱(Simulink)、神经网络工具箱(NeuralNetworkToolbox,ANN)、小波分析(waveletToolbox)、符号计算工具箱(SymbolicMathToolbox)等等,是药动学模型的模拟及运算的有力工具。2 药动学模型的拟合2.1 M程序法2.1.1 模型的建立 药动学是以药物体内过程的流向速率构建相应的模型,描述药物的体内经时变化,常用的模型有:房室模型、非线性动力学模型、统计矩模型、生理动力学模型等[5]。如单剂量静注(iv)给药的二室模型为
图中,D为剂量,C1、C2分别为中央室及周边室的血药浓度,
Vc、Vp分别为中央室及周边室的表观分布容积,K12为中央室到周边室的转运速率常数,K21为周边室到中央室的转运速率,Kel为药物从中央室的消除速率常数。ied
故建立微分方程组
dC1
dt
=-(K12+Kel)C1+K21C2dC2
dt
=K12C1-K21C2血药浓度C与时间t符合模型:C=Ae-αt+Beβt
并且有
K21=Aβ+B
αA+B,Kel=αβK21
,K12=α+β-K21-Kel
2.1.2 模型的拟合 建立M函数[iv2cfun.m],并保存。
functiondydt=iv2cfun(t,y)
obs=[27.10;155.80;305.40;604.00;903.40;1202.95;1802.75;2402.20;3601.90;7201.56];%输入时间t(min)-血药浓度C(mg・L-1)数据。静注给药200mg
cent=inline(’b(1)3exp(-b(2)3x)+b(3)3exp(-b(4)3x)’,’b’,’x’);%定义模型函数
b0=[100.0110.001];%给出初值
beta=nlinfit(obs(:,1),obs(:,2),cent,b0);%非线性
拟合,求各系数值
x=0:5:720;
k21=(beta(1)3beta(4)+beta(3)3beta(2))/(beta(1)+beta(3));
kel=beta(2)3beta(4)/k21;
k12=beta(2)+beta(4)-k21-kel;
dydt=[-(k12+kel)3y(1)+k213y(2);k123y(1)-k21
3y(2)];
建立M文件[iv2cfit.m],拟合模型,求有关参数。
clear,clc
obs=[27.10;155.80;305.40;604.00;903.40;1202.95;1802.75;2402.20;3601.90;7201.56];
cent=inline(’b(1)3exp(-b(2)3x)+b(3)3exp(-b(4)3x)’,’b’,’x’);
b0=[100.0110.001];
beta=nlinfit(obs(:,1),obs(:,2),cent,b0)x=0:5:720;
k21=[beta(1)3beta(4)+beta(3)3beta(2)]/[beta(1)+beta(3)]%显示药动学参数计算结果
kel=beta(2)3beta(4)/k21k12=beta(2)+beta(4)-k21-kel
Vc=200/[beta(1)+beta(3)]%剂量200mg
2012豪门盛宴
figure(1)%绘制药时数据的拟合曲线图
plot(x,cent(beta,x));axis([0720010]),holdon;scatter(obs(:,1),obs(:,2),5,’r’,’filled’);title(’静注给药二室模型药时数据的拟合曲线’),xlabel[’时间(min)’],ylabel[’血药浓度(mg・L-1)’];
figure(2)%模拟中央室及周边室血药浓度经时动态变
化
ode45(’iv2cfun’,[0720],[100]);%微分方程组数值的
动态解
title(’静注给药二室模型中央室及周边室血药浓度经时
中国教育学刊
变化’),xlabel(’时间(min)’),ylabel(’血药浓度(μg・mL-1)’);运行M文件,显示计算结果、药时曲线图和模拟效果,
计算结果为:
beta=4.39260.01852.76130.0009
© 1994-2009 China Academic Journal Electronic Publishing House. All rights reserved. wwwki
k21=0.0077kel=0.0021k12=0.0096Vc=27.9567即C=4.3926e-0.0185t
+2.7613e
-0.0009t
结果未按有效数字处理,以便比较(下同)。由中央室及
周边室浓度经时变化的模拟可见,单剂量静注给药后,中央
室浓度呈持续下降趋势,而周边室浓度先迅速上升然后逐渐
下降,最终两室的浓度都趋近于零。同理,如分别定义为一
室、三室模型函数,则可作一室、三室模型的拟合。所以,各
类药动学模型方程均可用inline命令和拟合曲线方程命令,
建立M函数和M文件来模拟。
2.1.3 模型拟合的评价 在拟合药动学模型特征时,由于
实验数据的误差,有时难以用一般方法推断模型的种类,因而常用残差平方和SSE、拟合度R2、AIC值等来推断。本数据用M程序法与药动学软件3P87/97(中国药理学会数学药理专业委员会编制)、Kinetica4.4.1(ThermoElectronCorporation编制)和PKSolver1.0(中国药科大学药剂教研室编制,最新开发的基于Excel的药动学软件)拟合结果的比较,见表1。可见,4种方法的结果一致,拟合效果良好,而M程序法与PKSolver法最为吻合。表1 M程序法与3P87/97、Kinetica和PKSolver软件二室模型拟合的比较Tab1 ComparisonoftheMprogramfittingoftwocompartmentmodelwith3P87/97,KineticaandPKSolver
拟合方法
模型参数值
判据值
A(mg・L-1)
B(mg・L
-1)
α(1/min)β(1/min)
SSE
R2AIC3M程序法4.39262.76130.01850.000
90.20120.9988-8.035263P87/974.35642.73710.018010.000840.20520.9939-8.8517Kinetica4.37822.781120.018610.0008740.201370.9988-8.02608PKSolver
4.39264
2.761270.01847
0.000
85
0.201
19
0.998
8-
8.03526
注:3权重为1
2.2 曲线拟合工具箱(cftool)法 cftool具有强大的曲线拟合功能,其核心为Powell’sDoglegMethod,是一种结合了
Leven2berg2Marquardt法和牛顿最速下降法的混合算法,在
线性和非线性拟合中是相当出的,拟合结果有很高的可信度。cftool提供了众多的数学模型,其指数(Exponential)模型中的单指数方程(y=a3exp(b3x))、双指数方程(y=a3exp(b3x)+c3exp(d3x)),可分别用于静注一室模型和静注二室模型的拟合。其他药动学方程可用“CustomEqua2
公差配合与测量技术tions”,做自定义拟合。从“Fitting”对话框的“Results”中可
以看到拟合参数值,“Analysis”工具可以得到SSE、R2
、RMSE(标准误差)等数据统计信息,供分析拟合结果的准确性。在用工具箱拟合时,需注意拟合参数和权值的选择,权值是与数据相联系的一个权向量,可通过“fitoptions”进行设置。本数据cftool拟合结果为:
腓特烈大帝C=4.393e
-0.01847t
+2.761e
-0.0008529t
,与M程序法完全一致。
2.3 药动学模型Simulink仿真 Matlab组件中的Simulink
豫公网安备41160202000603
站长QQ:729038198
关于我们
投诉建议