2017, Vol. 52
2. 中南大学湘雅医院 医学工程中心, 湖南 长沙 410008;
3. 中南大学数学与统计学院, 湖南 长沙 410083
2. Center of Medical Engineering, Xiangya Hospital, Central South University, Changsha 410008, China;
3. School of Mathematics and Statistics, Central South University, Changsha 410083, China
针对短半衰期静脉麻醉药物的多协变量群体药代建模, 普遍采用的是非线性混合效应建模法(nonlinear mixed effects modeling, NONMEM)[1, 2]。其原理是在药动学建模过程中引入协变量和误差项, 通过最大似然法求解群体药代动力学参数, 达到对个体间差异的解析、影响药物代谢的分析, 以及血药浓度预测等目的[3, 4]。然而, NONMEM方法需要采用半对数图对模型进行判断, 再根据结果校正模型, 工作繁杂, 且在建模过程中存在较多主观判断。其次, NONMEM方法采用的扩展最小二乘和最大似然法对初值的要求都较高, 初值选取不慎极易导致拟合失败。
灰色系统理论由邓聚龙教授于1982年首创[5], 在经济、工程、管理等领域都已得到广泛应用[6, 7]。该理论是控制理论在社会、经济系统的应用, 是自动控制科学与数学方法结合的初步尝试。它能对不确定系统特征值的变化进行预测, 其建模优势是样本点少, 不必有较好的分布规律。但早期灰色模型预测存在精度较低、背景值模糊等缺陷, 因此许多学者对其优化方式进行了广泛的研究[8, 9]。
智能优化算法是一种近年新兴的优化方法, 是目前受到关注最多的优化研究领域之一。它利用群体中的个体之间的信息交互和合作来实现寻优的目的。其不但具有较强的全局寻优能力, 且易于实现。近年来在生物信息[10, 11]、医学[12]等领域也得到广泛应用。差分进化(differential evolution, DE)算法是由Storn等[13]于1995年提出的一种基于实数编码、具有保优思维的全局优化算法, 其原理简单, 控制参数少, 无需设置初值和计算导数, 且易于理解和实现。因此, DE作为一种高效的并行搜索算法, 对其进行理论和应用研究具有重要的学术意义和工程价值[14, 15]。
异丙酚(propofol), 化学名为2、6-双异丙基酚, 是一种具有高度脂溶性的快速强效全身静脉麻醉药物[16]。异丙酚的药物代谢过程颇具特性, 包括一个快速分布阶段(α消除的半衰期约为1.8~8.3 min)和一个快速代谢消除阶段(β消除的半衰期约为34~64 min)。异丙酚在肝脏经酶的降解生成无活性的代谢产物, 这些代谢产物很快经肾脏排出体外, 且对肝肾功能无不良影响。以上这些特性, 都说明异丙酚明显优于目前常用的静脉麻醉药[17]。近年来, 异丙酚已在临床麻醉中得到广泛应用, 尤其是在麻醉药物计算机辅助设计-靶控输注技术(target controlled infusion, TCI)广泛应用于临床, 以及麻醉自动反馈控制系统的不断深入研究之后[18], 医生和学者们愈加重视对异丙酚等静脉麻醉药物和其群体药代学的研究。
本文采用差分进化优化算法, 结合多因素灰色模型时序性、小样本的优势, 以实测数据(含患者的年龄、身高、瘦体重、体表面积、剂量、输注率、采样时间等信息)为基础, 对异丙酚群体药代模型进行拟合分析, 构建血药浓度预测模型, 指导麻醉临床用药和研究血药浓度与患者体征的非线性关系。
材料与方法灰色模型 在灰色模型(grey model, GM)中, GM (1, 1) 是目前应用最为普遍的模型, 其通过单变量的时间序列单次累加处理, 对生成的序列建立一阶微分方程来揭示其内在的发展规律。随着研究的发展, 有学者提出了多变量灰色模型(multi-variable grey model, MGM), 它能够较好地反映系统中各变量之间相互关联、共同作用的关系, 能从系统的角度对各变量进行统一描述。MGM (1, m)模型是GM (1, 1) 模型在m个变量情况下的扩展, 通过建立m个一阶微分方程组成微分方程组进行模拟预测, 不同于GM (1, 1) 模型的单个一阶微分方程。自灰色多变量MGM (1, m)模型提出以来, 许多学者利用该模型应用于实际的社会系统中进行预测[19, 20]。
多因素灰色模型MGM (1, m) 设原始非负数据序列为
| $\begin{array}{l} x_i^{(1)}(k) = \sum\nolimits_{j = 1}^k {x_i^{(0)}(k)}, \\ i = 1,2, \ldots \ldots ,n;k = 1,2, \ldots \ldots ,N \end{array}$ | (1) |
对此生成序列建立m元一阶常微分方程, 多因素灰色模型的矩阵形式为:
| $\frac{{d{x^{(1)}}}}{{dt}} = A{x^{(1)}} + B$ | (2) |
其中
对MGM (1, m)模型求解, 将微分方程离散化, 并由梯形积分公式及矩形积分公式得:
| $x_i^{(0)}(k) = \sum\nolimits_{j = 1}^n {\frac{{{a_{ij}}}}{2}\left( {x_j^{(1)}(k) + x_j^{(1)}(k - 1)} \right)} + {b_i}$ | (3) |
其中, i=1, 2, ……, n; k=1, 2, ……, m。
通过求解方程(3) 得模型预测解为:
| $\hat x_i^{(1)}(k) = {e^{\hat A(k - 1)}}\left( {x_i^{(1)}(1) + {{\hat A}^{ - 1}}\hat B} \right) - {{\hat A}^{ - 1}}\hat B$ | (4) |
还原预测模型为:
| $\left\{ {\begin{array}{*{20}{c}} {\hat x_i^{(0)}(k) = \hat x_i^{(1)}(k)}\\ {\hat x_i^{(0)}(k) = \hat x_i^{(1)}(k) - \hat x_i^{(1)}(k - 1)} \end{array}} \right.$ | (5) |
其中, i=1, 2, ……, n; k=2, ……, m。
差分进化优化MGM模型 DE算法通过进化的方式来改善种群中候选解的质量, 利用当前种群个体间的差异信息来产生下一代种群。
(1) 种群初始化
生成初始群体, 设D为个体维数, 随机产生满足约束条件的NP个种群规模的向量
| ${x_{ji,0}} = ran{d_j}({b_{j,U}} - {b_{j,L}})$ | (6) |
其中, r为区间(0, 1) 上的随机数, 1≤i≤NP, 1≤j≤D, bj, U和bj, L分别为第j个变量xj的上界和下界。
(2) 变异
对每一个目标向量xi(G), 在{1, 2, ……, NP}范围内随机选择3个互异整数r1, r2, r3, 且使得i≠rj (1≤j≤3), 可用公式表达:
| $\begin{array}{l} {y_i}(G) = {x_{{r_1}}}(G) + F \cdot ({x_{{r_2}}}(G) + {x_{{r_3}}}(G)),\\ i = 1,2, \ldots \ldots ,NP \end{array}$ | (7) |
其中, 缩放因子F∈(0, 2) 为一常数, 控制偏差向量的放大程度(xr2(G)-xr3(G))。
(3) 交叉
为提高种群多样性, 引入离散杂交算子, 并通过如下方式得到试验向量:
| $\begin{array}{l} {u_{ij}}(G + 1) = \left\{ {\begin{array}{*{20}{c}} {{y_{ij}}(G + 1),{\rm{ if\;}}ran{d_j} \le CR\;{\rm{ or }}\;j = ran{d_i}}\\ {{x_{ij}}(G),{\rm{else}}} \end{array}} \right.,\\ i = 1,2, \ldots \ldots ,NP;j = 1,2, \ldots \ldots ,D \end{array}$ | (8) |
式中, randi∈{1, 2, ……, D}为随机序列, randj∈(0, 1) 用来确保yij(G+1) 从uij(G+1) 至少获取一个分量, 交叉因子R∈[0, 1]。
(4) 选择
差分进化的选择机制是基于局部竞争基于最小优化问题, 选择算子可由下式描述:
| $\begin{array}{l} {x_i}(G + 1) = \left\{ {\begin{array}{*{20}{c}} {{u_i}(G + 1),{\rm{if \;}}{F_x}({u_i}(G + 1) ＜ {F_x}({x_i}(G))}\\ {{x_i}(G),{\rm{else}}} \end{array},} \right.\\ i = 1,2, \ldots \ldots ,NP \end{array}$ | (9) |
针对MGM(1, m)模型, 提出模型λ-MGM(1, m), 其表达式为:
| $\begin{array}{l} {z^{\left( 1 \right)}}\left( k \right) = \lambda {x^{\left( 1 \right)}}\left( k \right) + \left( {1 - \lambda } \right){x^{(1)}}\left( {k - 1} \right),\\ k = 2, \ldots \ldots ,n;\lambda \in \left[ {0,1} \right] \end{array}$ | (10) |
其中, λ直接影响计算A和B的值。
DE优化MGM流程图如图 1所示。
|
Figure 1 Flow chart of DE-MGM (1, m). DE-MGM: Differential evolution-multi-variable grey model |
数据库建立 本研究所选用的血药浓度实测数据自来于World SIVA Open TCI Initiative数据库中Schnider的实验[21]。Schnider实验24例成人(13男, 11女), 每人两组实验数据, 共48组, 年龄、身高、体重和瘦体重分布如表 1, 患者接受异丙酚团注1 h后, 分别进行25、50、100和200 μg·kg-1·min-1持续输注。选用其实验总数据中的1 055对数据作为构建模型的总数据集, 其中761对(约占总样本数的70%)作为训练数据, 294对数据(约占总样本数的30%)作为测试数据。为减少在样本验证过程中出现外推现象, 将样本集中的最大、最小的样本作为训练样本, 保证训练样本覆盖检验样本的极值。
|
Table 1 Evaluation of data base structure. LBM: Lean body mass |
此外, 多数研究认为异丙酚符合药代三室开放模型, 其中, 体重、性别、年龄和瘦体重(lean body mass, LBM)均为异丙酚药代学的重要影响因素[22, 23]。因此, 基于之前对异丙酚的研究[16, 17], 选取每个患者的年龄、瘦体重、身高、体表面积、血药浓度采样时间、总剂量、注射率作为自变量, 该患者在采样时刻的异丙酚血药浓度为因变量, 构建时序DE-MGM异丙酚血药浓度预测模型。同时, 因数据数量级差别较大, 为避免造成误差较大, 并加快收敛, 对输入与输出数据进行归一化处理。建模环境为Matlab软件, 含grey model_tool box工具箱等。最后, 将DE-MGM预测结果与Schnider实验中采用NONMEM方法所得的预测效果进行对比。
结果 1 预测结果的相关性NONMEM、DE-MGM的预测趋势相似, NONMEM、DE-MGM预测值与实测值均呈强正相关性。DE-MGM预测结果与实测浓度的相关性分析见图 2A, NONMEM网络预测结果与实测浓度的相关性分析见图 2B。
|
Figure 2 Correlation analysis of prediction in DE-MGM (A) and NONMEM (B). NONMEM: Nonlinear mixed effects modeling |
在进行验证的14组数据中, DE-MGM算法最优个体的相对误差为7.27%, 最差个体的相对误差为20.51%, NONMEM算法最优个体的相对误差为11.13%, 最差个体的相对误差为39.75%。DE-MGM、NONMEM在麻醉期最优预测和最差预测个体的比较见图 3。由于血药浓度数据数量级差别较大, 因此, 纵轴采用指数形式表达, 横坐标为时间秒。
|
Figure 3 Comparison of best and worst prediction during the anesthesia |
DE-MGM, NONMEM两种模型的预测误差较为集中, 绝大部分均在[-1 1]区间。DE-MGM预测浓度、NONMEM预测浓度与其对应的预测误差的相关性分别见图 4。
|
Figure 4 Correlation analysis of the error in DE-MGM (A) and NONMEM (B) |
本研究选用较为常用的分析和评价模型预测误差指标[21]。CM为模型预测浓度, CP为实测血药浓度。PE (performance error, PE)表示执行误差, 又称相对误差, 常用于检验药代模型的指标。
| ${\rm{PE}} = \frac{{{C_{\rm{M}}} - {C_{\rm{P}}}}}{{{C_{\rm{P}}}}} \times 100$ | (11) |
DE-MGM算法预测误差明显优于NONME算法。MDE-MGM与NONMEM预测误差统计结果对比见表 2。其中, 偏离性(median performance error, MDPE):执行误差的中位数, MDPEi = median{PEij, j=1, 2, ……, Ni}, 代表系统偏离预期浓度的误差。精确度(median absolute performance error, MDAPE):执行误差绝对值的中位数, MDAPEi = median{|PE|ij, j=1, 2, ……, Ni}, 是所有实测浓度与预计浓度的误差。摆动性(wobble):用中位绝对偏差值(MDADPE)表示, 该参数衡量误差的变异程度, Wobblei = median {|PEij-MDPEi|, j=1, 2, ……, Ni}, 即代表执行误差的易变性。
|
|
Table 2 Comparison of prediction error of DE-MGM and NONMEM |
NONMEM、DE-MGM预测值与实测值的相关系数均在0.95以上, 预测质量较好。其中, DE-MGM预测值的相关系数为0.99, 优于NONMEM (0.97)。
然而, NONMEM预测值与Y=X线之间偏角较大, 随着浓度升高, 偏离趋势愈加明显。而DE-MGM预测浓度一直与实测浓度明显相关性更高。DE-MGM, NONMEM两种模型的预测误差均在可接受范围。其中, NONMEM预测误差随实测浓度的增大而增大的趋势较为明显, 预测误差相对分散, 而DE-MGM误差相对集中。
从表 2中不难看出, 无论是偏离性、精确度、摆动性, DE-MGM建模法预测效果均优于NONMEM, 其中摆动性两者较为接近。从标准差和置信区间不难看出, NONMEM预测离散度较大, 而经DE算法优化后的MGM其预测误差各指标的离散程度均优于NONMEM方法, 模型预测误差较稳定, 波动性小。
同时, 不难看出DE-MGM的预测存在着一定的波动性, 预测值与实际值不能完全一致。此外, 基于DE-MGM算法的建模仍无法解决智能算法在药物血药浓度预测中缺乏生理理论模型的问题。
结论本文利用差分进化算法优化灰色模型, 将两者的优势结合, 实现了灰色模型参数的优化。用灰色预测模型的平均误差作为适应度值, 弥补了传统算法当中的部分缺陷。
以复杂的非线性、时变性较强的短效静脉麻醉药物异丙酚的实测血药浓度临床数据为基础进行模型仿真, 结果表明, 相比传统的NONMEM方法, DE-MGM预测结果各误差更小, 精度更高, 反映出较高的非线性拟合能力。
DE-MGM方法能够迅速地捕捉血药浓度和患者体征、给药方案、时间等关系, 为非线性群体药代药效学研究分析提供了新方法, 可以进一步有效地指导临床麻醉用药, 特别针对短半衰期的速效阿片类药物的群体药代药效学研究分析, 为预测术中麻醉深度提供技术支持。
| [1] | Bonate PL. Pharmacokinetic-pharmacodynamic Modeling and Simulation[M]. 2 ed. New York: Springer Science + Business Media, Inc, 2001. |
| [2] | Lin R, Lin W, Wang C, et al. Population pharmacokinetic/pharmacodynamic modeling of warfarin by nonlinear mixed effects model[J]. Acta Pharm Sin (药学学报), 2015, 50: 1280–1284. |
| [3] | Bergstrand M, Karlsson MO. Handling data below the limit of quantification in mixed effect models[J]. AAPS J, 2009, 11: 371–380. DOI:10.1208/s12248-009-9112-5 |
| [4] | Wang L, Cao J, Ramsay JO, et al. Estimating mixed-effects differential equation models[J]. Stat Comput, 2014, 24: 111–121. DOI:10.1007/s11222-012-9357-1 |
| [5] | Deng JL. Introduction to grey system theory[J]. J Grey Syst-UK, 1989, 1: 1–24. |
| [6] | Kayacan E, Ulutas B, Kaynak O. Grey system the-ory-based models in time series prediction[J]. Expert Syst Appl, 2010, 37: 1784–1789. DOI:10.1016/j.eswa.2009.07.064 |
| [7] | Bezuglov A, Comert G. Short-term freeway traffic parameter prediction:application of grey system theory models[J]. Expert Syst Appl, 2016, 62: 284–292. DOI:10.1016/j.eswa.2016.06.032 |
| [8] | Bahrami S, Hooshmand RA, Parastegari M. Short term electric load forecasting by wavelet transform and grey model improved by PSO (particle swarm optimization) algorithm[J]. Energy, 2014, 72: 434–442. DOI:10.1016/j.energy.2014.05.065 |
| [9] | Zhang L, Zheng Y, Wang K, et al. An optimized Nash nonlinear grey Bernoulli model based on particle swarm optimization and its application in prediction for the incidence of hepatitis B in Xinjiang, China[J]. Comput Biol Med, 2014, 49: 67–73. DOI:10.1016/j.compbiomed.2014.02.008 |
| [10] | Baytas IM, Lin K, Wang F, et al. Phenotree:interactive visual analytics for hierarchical phenotyping from large-scale electronic health records[J]. IEEE Trans Multimedia, 2016, 18: 2257–2270. DOI:10.1109/TMM.2016.2614225 |
| [11] | Liu X, Sun F, Jin Y, et al. Application of near infrared spectroscopy combined with particles warm optimization based least square support vactor machine to rapid quantitative analysis of Corni Fructus[J]. Acta Pharm Sin (药学学报), 2015, 50: 1645–1651. |
| [12] | Zhang MQ, Wang H, Xiong K. Is the neocortex a novel reservoir for adult mammalian neurogenesis?[J]. Neural Regen Res, 2011, 6: 1334–1341. |
| [13] | Storn R, Price K. Differential evolution-a simple and efficient heuristic for global optimization over continuous spaces[J]. J Global Optim, 1997, 11: 341–359. DOI:10.1023/A:1008202821328 |
| [14] | Wang WJ, Xie B, Pan KJ, et al. Electromagnetic imaging of heterogeneous media based on the accelerated differential evolution algorithm[J]. Prog Geophys (地球物理学进展), 2010, 25: 2002–2008. |
| [15] | De Falco I. Differential evolution for automatic rule extraction from medical databases[J]. Appl Soft Comput, 2013, 13: 1265–1283. DOI:10.1016/j.asoc.2012.10.022 |
| [16] | Sebel PS, Lowdon JD. Propofol:a new intravenous anesthetic[J]. Anesthesiology, 1989, 71: 260–277. DOI:10.1097/00000542-198908000-00015 |
| [17] | Bandschapp O, Filitz J, Ihmsen H, et al. Analgesic and antihyperalgesic properties of propofol in a human pain model[J]. Anesthesiology, 2010, 113: 421–428. DOI:10.1097/ALN.0b013e3181e33ac8 |
| [18] | Ting CH, Arnott RH, Linkens DA, et al. Migrating from target-controlled infusion to closed-loop control in general anaesthesia[J]. Comput Methods Programs Biomed, 2004, 75: 127–139. DOI:10.1016/j.cmpb.2003.11.005 |
| [19] | Xiong P, Dang Y, Wu X, et al. Combined model based on optimized multi-variable grey model and multiple linear regression[J]. J Syst Eng Electron, 2011, 22: 615–620. DOI:10.3969/j.issn.1004-4132.2011.04.010 |
| [20] | Guo X, Liu S, Wu L, et al. A multi-variable grey model with a self-memory component and its application on engineering prediction[J]. Eng Appl Artif Intel, 2015, 42: 82–93. DOI:10.1016/j.engappai.2015.03.014 |
| [21] | Schnider TW, Minto CF, Gambus PL, et al. The influence of method of administration and covariates on the pharma-cokinetics of propofol in adult volunteers[J]. Anesthesiology, 1998, 88: 1170–1182. DOI:10.1097/00000542-199805000-00006 |
| [22] | Choong E, Loryan I, Lindqvist M, et al. Sex difference in formation of propofol metabolites:a replication study[J]. Basic Clin Pharmacol, 2013, 113: 126–131. DOI:10.1111/bcpt.2013.113.issue-2 |
| [23] | Vuyk J, Oostwouder CJ, Vletter AA, et al. Gender differ-ences in the pharmacokinetics of propofol in elderly patients during and after continuous infusion[J]. Brit J Anaesth, 2001, 86: 183–188. DOI:10.1093/bja/86.2.183 |