CN115472308A - 一种全血药物谷浓度预测系统 - Google Patents

一种全血药物谷浓度预测系统 Download PDF

Info

Publication number
CN115472308A
CN115472308A CN202211139186.XA CN202211139186A CN115472308A CN 115472308 A CN115472308 A CN 115472308A CN 202211139186 A CN202211139186 A CN 202211139186A CN 115472308 A CN115472308 A CN 115472308A
Authority
CN
China
Prior art keywords
prediction
model
pharmacokinetic
module
theta
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
CN202211139186.XA
Other languages
English (en)
Inventor
罗子峻
吴瑶洁
钱蔚
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Shanghai Weixin Technology Co ltd
Original Assignee
Shanghai Weixin Technology Co ltd
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Shanghai Weixin Technology Co ltd filed Critical Shanghai Weixin Technology Co ltd
Priority to CN202211139186.XA priority Critical patent/CN115472308A/zh
Publication of CN115472308A publication Critical patent/CN115472308A/zh
Pending legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16H—HEALTHCARE INFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR THE HANDLING OR PROCESSING OF MEDICAL OR HEALTHCARE DATA
    • G16H70/00—ICT specially adapted for the handling or processing of medical references
    • G16H70/40—ICT specially adapted for the handling or processing of medical references relating to drugs, e.g. their side effects or intended usage
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16C—COMPUTATIONAL CHEMISTRY; CHEMOINFORMATICS; COMPUTATIONAL MATERIALS SCIENCE
    • G16C10/00—Computational theoretical chemistry, i.e. ICT specially adapted for theoretical aspects of quantum chemistry, molecular mechanics, molecular dynamics or the like
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16C—COMPUTATIONAL CHEMISTRY; CHEMOINFORMATICS; COMPUTATIONAL MATERIALS SCIENCE
    • G16C20/00—Chemoinformatics, i.e. ICT specially adapted for the handling of physicochemical or structural data of chemical particles, elements, compounds or mixtures
    • G16C20/50—Molecular design, e.g. of drugs
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16C—COMPUTATIONAL CHEMISTRY; CHEMOINFORMATICS; COMPUTATIONAL MATERIALS SCIENCE
    • G16C20/00—Chemoinformatics, i.e. ICT specially adapted for the handling of physicochemical or structural data of chemical particles, elements, compounds or mixtures
    • G16C20/70—Machine learning, data mining or chemometrics
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16H—HEALTHCARE INFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR THE HANDLING OR PROCESSING OF MEDICAL OR HEALTHCARE DATA
    • G16H50/00—ICT specially adapted for medical diagnosis, medical simulation or medical data mining; ICT specially adapted for detecting, monitoring or modelling epidemics or pandemics
    • G16H50/50—ICT specially adapted for medical diagnosis, medical simulation or medical data mining; ICT specially adapted for detecting, monitoring or modelling epidemics or pandemics for simulation or modelling of medical disorders

Landscapes

  • Engineering & Computer Science (AREA)
  • Health & Medical Sciences (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • General Health & Medical Sciences (AREA)
  • Computing Systems (AREA)
  • Chemical & Material Sciences (AREA)
  • Medical Informatics (AREA)
  • Theoretical Computer Science (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Public Health (AREA)
  • Crystallography & Structural Chemistry (AREA)
  • Epidemiology (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • Pharmacology & Pharmacy (AREA)
  • Physics & Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Databases & Information Systems (AREA)
  • Primary Health Care (AREA)
  • Medicinal Chemistry (AREA)
  • Artificial Intelligence (AREA)
  • Toxicology (AREA)
  • Software Systems (AREA)
  • Evolutionary Computation (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Biomedical Technology (AREA)
  • Pathology (AREA)
  • Investigating Or Analysing Biological Materials (AREA)

Abstract

本发明属于药物检测技术领域，具体涉及一种全血药物谷浓度预测系统，包括原始数据采集模块、药动学计算模块、机器学习回归预测模块、预测结果输出模块以及模型决策树解释模块。克服了现有技术的不足，通过患者信息预测一定时间之后的血药浓度谷浓度，改变临床对剂量决策困难并且需要频繁采血的现状，降低患者感染风险以及经济负担，提高患者治疗依从性以及满意度，最终实现个体化精准用药。

Description

一种全血药物谷浓度预测系统
技术领域
本发明属于药物检测技术领域，具体涉及一种全血药物谷浓度预测系统。
背景技术
他克莫司(Tacrolimus)，又名FK506，是一种大环内酯类抗生素，为一种强力的新型免疫抑制剂，主要通过抑制白介素-2(IL-2)的释放，全面抑制T淋巴细胞的作用，较环孢素(CsA)强100倍。由于其良好的免疫抑制效果，成为了肝脏、心脏、肾脏及骨髓移植患者的首选免疫抑制药物，移植后排斥反应对传统免疫抑制方案耐药者，也可选用该药物。
他克莫司的其中一种制剂他克莫司胶囊在全球范围内作为移植患者一线用药广泛应用，随着临床实践的增加，医疗人员以及药物研究者逐渐对他克莫司有了更多的认识。由于免疫抑制药物本身对剂量的要求较其他的药物更加精准，并且在临床上他克莫司使用过量展现出的肾毒性，让其成为治疗药物监测(TDM)名单之一。
群体药代动力学是基础药代动力学与统计学结合的科学，主要研究药物在体内代谢的群体规律、药代动力学参数的统计分布以及影响因素。群体药动学当中提出变异的概念，用于描述实际观测值与基础药动学模型估计值的差异。群体药代动力学的变异包括确定性的和随机性的两类，确定性变异通常指的是由已知因素造成的药代动力学差异，随机性变异指的是遵循某种分布的均值为0的随机差异。在建模方法上，药代动力学金标准为非线性混合效应模型法，其混合效应模型指的是固定效应与随机效应。在建模过程中，由影响因素表达的确定型变异与药动学参数共同建立固定效应模型，在固定效应模型基础之上，添加随机的个体间差异以及残差变异，建立随机效应模型。模型在计算结束之后能够得到估计的药动学参数以及变异情况，将估计的药动学参数代入到药动学公式当中，结合随机变异的分布，使用最大后验估计贝叶斯法能够计算出观测目标(大部分为血药浓度)的点估计以及全分布。
现有技术结合药代动力学以及统计学，使用较少的样本量即可建立较为稳健的模型，且对影响因素的分析较为充分。
使用群体药动学在小样本当中可建立稳健的模型，造成偏差较大，方差较小的情况，在个体预测上误差较大；并且群体药动学建模需要很强的专业领域知识，大部分的从业人员无法自行开展研究，造成应用模型指导用药困难
发明内容
本发明的目的在于提供一种全血药物谷浓度预测系统。
为解决上述问题，本发明所采取的技术方案如下：
一种全血药物谷浓度预测系统，包括原始数据采集模块、药动学计算模块、机器学习回归预测模块、预测结果输出模块以及模型决策树解释模块；
所述原始数据采集模块用于自医院系统中获取预测所需的相关数据；
所述药动学计算模块通过非线性混合效应模型法计算样本群体药动学参数以及群体预测值，同时通过最大后验贝叶斯计算点估计得到个体预测值；
所述机器学习回归预测模块通过手动输入需要预测的血药浓度时间，将其与药动学计算模块输出的结果一并输入到训练好的机器学习回归预测模型当中，计算输出群体药动学血药浓度点估计残差；
所述预测结果输出模块用于根据机器学习回归预测模块输出的结果对多个药物谷浓度预测数值进行理论药时曲线函数图像绘制，叠加群体参数曲线以及个体参数曲线，单个谷浓度预测直接显示数值；
所述模型决策树解释模块用于根据模型预测时计算出的SHAP值绘制出力图。
进一步，所述预测所需的相关数据包括患者的体征数据、相关检测结果以及基因监测结果，其中对于缺失的数据采用预设的典型值填充，对于采集的数据按照特征工程处理并暂存为Python字典格式数据一份，以及后续预测使用的事件格式数据一份后进入系统流程下一步，并同时存入研究数据库。
进一步，所述非线性混合效应模型法在拟合数据建模过程使用的是扩展的最小二乘法，其目标函数公式如下：
Figure BDA0003852730960000031
其中O(θ,ξ,σ2)为非线性混合效应模型法中关于θ,ξ,σ2的目标函数；σ2为残差变异的平方；θ为药动学参数；ξ为指数矢量参数，通常设置为0到2之间；n为样本量；yi为个体血药浓度实测值；xi为个体自变量(例如剂量、时间等)；f(xi,θ)为使用θ,xi计算血药浓度的药动学公式；u(xi,θ,ξ)为关于xi,θ,ξ的权重函数，计算为f(xi,θ)ξ；
通过定义药动学参数以及公式，代入到上述式中，进行ELS运算使目标函数O(θ,ξ,σ2)尽可能小(获得极小值)，同时获得对应的θ,ξ,σ2的集合。
进一步，所述药动学公式如下：
Figure BDA0003852730960000032
其中
Figure BDA0003852730960000033
为血药浓度估计值；ke为消除速率常数；X0为单次剂量；Clast doseadj为上次调整剂量时的血药浓度实测值；τ为给药间隔；n为自上次调整剂量开始计算的给药次数。
进一步，所述药动学计算模块使用非线性混合效应模型法建立药动学模型，所述药动学模型的具体建模方法包括：
(1)仅采用最基础的个体自变量建立无其他信息的基础药动学模型，即仅采用时间、剂量对血药浓度进行建模，得到基础模型如下表示：
CL/F＝TV(θ1)×EXP(η1)
V/F＝TV(θ2)×EXP(η2)
其中，CL/F和V/F分别为表观清除率生物利用度比及表观分布容积生物利用度比；TV(θ1)和TV(θ2)分别为CL/F与V/F的群体典型值，η1与η2代表参数的个体间变异；
(2)随后将固定效应模型参数加入服从正态或对数正态分布的随机效应，对残差变异使用加法常系数混合随机效应，并认为随机效应同样存在于个体间变异中，得到如下公式：
Figure BDA0003852730960000041
Figure BDA0003852730960000042
其中ηi服从均值为0，方差为
Figure BDA0003852730960000043
的对数正态分布；εij服从均值为0，方差为
Figure BDA0003852730960000044
的正态分布。将上述模型结构代入到非线性混合效应模型法的对应参数当中，进行运算后，得到ω,σ2的集合。
进一步，所述机器学习回归预测模块的具体算法包括：
(1)预测函数
设存在数据集D＝(xi,yi)(|D|＝n,xi∈Rm,yi∈R)，m为特征数量，n为样本数量；另存在一个由K棵树组成的模型
Figure BDA0003852730960000051
其中
Figure BDA0003852730960000052
表示第k棵树对样本xi的预测结果；T表示该树叶子节点的数量；ωt为该树第t个叶子节点的权重，可理解为预测值；q(x)函数用于搜索样本xi所属叶子节点的下标；
(2)损失函数与树的分裂
定义损失函数
Figure BDA0003852730960000053
其中惩罚项
Figure BDA0003852730960000054
用于防止树的结构过于复杂造成过拟合；
Figure BDA0003852730960000055
为Pseudo-Huber损失函数；
设
Figure BDA0003852730960000056
是xi在第t次迭代的预测值，有
Figure BDA0003852730960000057
将其进行泰勒二次展开得到
Figure BDA0003852730960000058
由于
Figure BDA0003852730960000059
为常数，对函数无影响，因此可省略；
定义被分至第j个叶子节点的样本对应下标集合Ij＝{iq(xi)＝j}，则
Figure BDA00038527309600000510
对ωj求导并计算最小值，得到
Figure BDA0003852730960000061
令I＝IL∪IR，得到衡量节点划分与否的公式，
Figure BDA0003852730960000062
公式每次都用不同的特征计算，得分最大的特征作为当前树节点的划分点；
(3)分割算法
定义Dk＝{(x1k,h1),(x2k,h2)…(xnk,hn)}表示训练样本第k个特征的值以及对应的二阶梯度；以此定义排名函数
Figure BDA0003852730960000063
并结合如下公式用来查找候选分割点，
Figure BDA0003852730960000064
以下是近似分割点查找算法伪代码：
For k＝1 to m do
查找第k个特征的候选分割点记为Sk＝{sk1,sk2…skl}；
可按树进行计算(Global)也可按每次分割进行计算；
End for
For k＝1 to m do
Figure BDA0003852730960000065
Figure BDA0003852730960000066
End for
算法输入为原始数据的K列特征、群体药动学预测模块估计的残差以及标签列，即共有K+1列特征以及1列标签，进行交叉验证以及网格搜索调参训练后得到XGBoost的残差预测模型，该模型预测结果加上步骤二群体药动学预测模块点估计值得到总预测结果。
本发明与现有技术相比较，具有以下有益效果：
本发明所述一种全血药物谷浓度预测系统，通过患者信息预测一定时间之后的血药浓度谷浓度，改变临床对剂量决策困难并且需要频繁采血的现状，降低患者感染风险以及经济负担，提高患者治疗依从性以及满意度，最终实现个体化精准用药。
附图说明
图1为一种全血药物谷浓度预测系统的应用场景流程图。
图2为具体实施过程中计算XGBoost回归器对应的特征SHAP均值图。
图3为具体实施过程中机器学习模型真实值与预测值的散点图。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
实施例
本发明所述一种全血药物谷浓度预测系统，包括原始数据采集模块、药动学计算模块、机器学习回归预测模块、预测结果输出模块以及模型决策树解释模块。
原始数据采集模块用于自医院系统中获取预测所需的相关数据，包括患者的体征数据、相关检测结果以及基因监测结果，其中对于缺失的数据采用预设的典型值填充，对于采集的数据按照特征工程处理并暂存为Python字典格式数据一份，以及后续预测使用的事件格式数据一份后进入系统流程下一步，并同时存入研究数据库。
药动学计算模块通过非线性混合效应模型法计算样本群体药动学参数以及群体预测值，同时通过最大后验贝叶斯计算点估计得到个体预测值。
1.回归方法
非线性混合效应模型法在拟合数据建模过程使用的是扩展的最小二乘法(ELS)，目标函数公式如下：
Figure BDA0003852730960000081
其中O(θ,ξ,σ2)为NONMEM定义的关于θ,ξ,σ2的目标函数；σ2为残差变异的平方；θ为药动学参数(固定参数)；ξ为指数矢量参数，通常设置为0到2之间；n为样本量；yi为个体血药浓度实测值；xi为个体自变量(例如剂量、时间等)；f(xi,θ)为使用θ,xi计算血药浓度的药动学公式；u(xi,θ,ξ)为关于xi,θ,ξ的权重函数，计算为f(xi,θ)ξ。
通过定义药动学参数以及公式，代入到(1)式中，进行ELS运算使目标函数O(θ,ξ,σ2)尽可能小(获得极小值)。同时获得对应的θ,ξ,σ2的集合。
2.药动学模型
为方便应用，本方法调整后的一房室一级消除药代动力学公式如下：
Figure BDA0003852730960000082
其中
Figure BDA0003852730960000083
为血药浓度估计值；ke为消除速率常数；X0为单次剂量；Clast doseadj为上次调整剂量时的血药浓度实测值；τ为给药间隔；n为自上次调整剂量开始计算的给药次数。
3.药动学基础模型
在本方法中，仅采用最基础的个体自变量建立无其他信息的基础药动学模型，即仅采用时间、剂量对血药浓度进行建模，得到基础模型如下表示：
CL/F＝TV(θ1)×EXP(η1) (3)
V/F＝TV(θ2)×EXP(η2) (4)
其中CL/F和V/F分别为表观清除率生物利用度比及表观分布容积生物利用度比；TV(θ1)和TV(θ2)分别为CL/F与V/F的群体典型值，η1与η2代表参数的个体间变异。
4.药动学统计学模型
随后将固定效应模型参数加入服从正态或对数正态分布的随机效应，本方法对残差变异使用加法常系数(乘法)混合随机效应，并认为随机效应同样存在于个体间变异中，得到如下公式：
Figure BDA0003852730960000091
Figure BDA0003852730960000092
其中ηi服从均值为0，方差为
Figure BDA0003852730960000093
的对数正态分布；εij服从均值为0，方差为
Figure BDA0003852730960000094
的正态分布。将上述模型结构代入到非线性混合效应模型法的对应参数当中，进行运算后，得到ω,σ2的集合。
至此完成药动学模型设置的工作，除了(5)(6)外，吸收速率常数ka，记为θ3。以下是非线性混合效应模型法估计的固定参数。
Figure BDA0003852730960000101
5.基于NUTS采样器的最大后验估计
在此使用Hoffman和Gelman在2011年发表的NUTS算法进行重采样以进行最大后验估计[1]。算法如下：
给定θ0，δ，L，M，Madapt：
设定ε0＝FindReasonableEpsilon(θ),μ＝log(10ε0),
Figure BDA0003852730960000102
γ＝0.05，
t0＝10，κ＝0.75.
for m＝1 to M do
在N(0，I)采样r0.
在
Figure BDA0003852730960000103
重采样u.
初始化θ-＝θ+＝θm＝θm-1,r-＝r+＝r0,j＝0,n＝1,s＝1.
While s＝1 do
在B(-1，1)中选择vj.
if vj＝-1 then
θ-，r-，-，-，θ′，n′，s′α，nα←BulidTree(θ-,r-,u,vj,j,εm-1θm-1,r0)
else
-，-θ+，r+，θ′，n′，s′α，nα←BulidTree(θ+,r+,u,vj,j,εm-1θm-1,r0)
end if
随机在U(0,1)中取1个值赋值给临时变量π
if s′＝1 and
Figure BDA0003852730960000111
then
θm←θ′.
end if
n←n+n′.
s←s′and[(θ+-θ-)·r-≥0]and[(θ+-θ-)·r+≥0].
j←j+1
end while
if m≤Madapt then
设定
Figure BDA0003852730960000112
设定
Figure BDA0003852730960000113
else
设定
Figure BDA0003852730960000114
end if
end for
function BulidTree(θ,r,u,vj,j,ε,θ0,r0)
if j＝0 then
θ′,r′←Leapfrog(θ,r,vε).
Figure BDA0003852730960000121
Figure BDA0003852730960000122
return
Figure BDA0003852730960000123
else
θ-,r-,θ+,r+,θ′,n′,s′,α′,n′α←BuildTree(θ,r,u,v,j-1,ε,θ0,r0).
if s′＝1 then
if v＝-1 then
θ-,r-,-,-,θ″,n″,s″,α″,n″α←BuildTree(θ-,r-,u,v,j-1,ε,θ0,r0).
else
-,-,θ+,r+,θ″,n″,s″,α″,n″α←BuildTree(θ+,r+,u,v,j-1,ε,θ0,r0).
end if
随机在U(0,1)中取1个值赋值给临时变量λ.
if
Figure BDA0003852730960000124
then
θ′←θ″.
end if
α′←α′+α″,n′α←n′α+n″α.
s←s′and[(θ+-θ-)·r-≥0]and[(θ+-θ-)·r+≥0].
n′←n′+n″.
end if
returnθ-,r-,θ+,r+,θ′,n′,s′,α′,n′α.
end if
function Leapfrog(θ,r,ε)
Figure BDA0003852730960000131
Figure BDA0003852730960000132
Figure BDA0003852730960000133
return
Figure BDA0003852730960000134
function FindReasonableEpsilon(θ)
初始化ε＝1,r～(0,I).
θ′,r′←Leapfrog(θ,r,ε).
Figure BDA0003852730960000135
while
Figure BDA0003852730960000136
do
ε＝2aε.
θ′,r′←Leapfrog(θ,r,ε).
end while
returnε.
其中θ为我们输入的D维参数；L(θ)为θ的对数联合概率密度；r为θ状态转移的动量；
Figure BDA0003852730960000141
为对θ的梯度；I表示单位矩阵；N(μ，∑)表示多元正态分布，均值为μ和协方差矩阵为∑。
另外，定义平均接收概率以及期望计算公式如下：
Figure BDA0003852730960000142
其中
Figure BDA0003852730960000143
为在马尔可夫链中第t次迭代最后加倍时所搜索的所有状态的集合；θt-1和rt,0是马尔可夫链第t次迭代的初始位置和(重采样)动量。HNUTS可以理解为HMC在最后的加倍迭代中所搜索的位置-动量状态的平均接受概率。
此时能够使用Ht≡δ-HNUTS,x≡logε应用到状态更新公式迫使hNUTS＝δ，假定我们需要找到一个参数组合x∈R，状态更新公式如下：
Figure BDA0003852730960000144
其中μ是一个自由选择的点，迭代中的xt向它收缩；γ＞0是一个自由参数控制向μ收缩的量；t0≥0是一个自由参数用于稳定算法的初始迭代；ηt≡t-κ是遵循
Figure BDA0003852730960000145
的步长计划，且定义
Figure BDA0003852730960000146
在上述算法内，定义模型、θ0、先验概率分布p(θ)，输入超参数目标平均接受概率δ以及迭代次数Madapt，运行算法可得到满足先验分布p(θ)的抽样样本，结合最大后验点估计的公式
Figure BDA0003852730960000151
可求出θ的点估计值。将θ的点估计值代入到药动学模型的公式中，输入该点对应的x，即可算出对血药浓度的点估计并计算估计的残差，残差以列合并到原数据集当中。
机器学习回归预测模块通过手动输入需要预测的血药浓度时间，将其与药动学计算模块输出的结果一并输入到训练好的机器学习回归预测模型当中，计算输出群体药动学血药浓度点估计残差。
机器学习回归预测模块的具体算法包括：
1.预测函数
设存在数据集D＝(xi,yi)(|D|＝n,xi∈Rm,yi∈R)，m为特征数量，n为样本数量。另存在一个由K棵树组成的模型
Figure BDA0003852730960000152
其中
Figure BDA0003852730960000153
表示第k棵树对样本xi的预测结果；T表示该树叶子节点的数量；ωt为该树第t个叶子节点的权重，可理解为预测值；q(x)函数用于搜索样本xi所属叶子节点的下标。
2.损失函数与树的分裂
定义损失函数
Figure BDA0003852730960000154
其中惩罚项
Figure BDA0003852730960000155
用于防止树的结构过于复杂造成过拟合；
Figure BDA0003852730960000161
为Pseudo-Huber损失函数。设
Figure BDA0003852730960000162
是xi在第t次迭代的预测值，有
Figure BDA0003852730960000163
将其进行泰勒二次展开得到
Figure BDA0003852730960000164
由于
Figure BDA0003852730960000165
为常数，对函数无影响，因此可省略。
定义被分至第j个叶子节点的样本对应下标集合Ij＝{i|q(xi)＝j}，则
Figure BDA0003852730960000166
对ωj求导并计算最小值，得到
Figure BDA0003852730960000167
令I＝IL∪IR，得到衡量节点划分与否的公式，
Figure BDA0003852730960000168
该公式每次都用不同的特征计算，得分最大的特征作为当前树节点的划分点。
3.分割算法
定义Dk＝{(x1k,h1),(x2k,h2)…(xnk,hn)}表示训练样本第k个特征的值以及对应的二阶梯度。以此定义排名函数
Figure BDA0003852730960000171
并结合如下公式用来查找候选分割点，
Figure BDA0003852730960000172
以下是近似分割点查找算法伪代码：
For k＝1 to m do
查找第k个特征的候选分割点记为Sk＝{sk1,sk2…skl}.
可按树进行计算(Global)也可按每次分割进行计算。
End for
For k＝1 to m do
Figure BDA0003852730960000173
Figure BDA0003852730960000174
End for
以上为XGBoost机器学习算法内容，算法输入为原始数据的K列特征、群体药动学预测模块估计的残差以及标签列，即共有K+1列特征以及1列标签。进行交叉验证以及网格搜索调参训练后得到XGBoost的残差预测模型，该模型预测结果加上步骤二群体药动学预测模块点估计值得到总预测结果。
预测结果输出模块用于根据机器学习回归预测模块输出的结果对多个药物谷浓度预测数值进行理论药时曲线函数图像绘制，叠加群体参数曲线以及个体参数曲线，单个谷浓度预测直接显示数值。
模型决策树解释模块用于根据模型预测时计算出的SHAP值绘制出力图。
具体实施
在医院采集到肾移植术后应用他克莫司患者的原始数据，并处理为事件格式如下：
ID DATE TIME AMT II ADDL CONC LNDV CMT MDV
1 1 6:00 4 24 9 . . 1 1
1 1 18:00 4 24 9 . . 1 1
1 3 6:00 . . . 7.4 2.00148 2 0
1 4 6:00 . . . 7.9 2.066862759 2 0
1 6 6:00 . . . 7.7 2.041220329 2 0
1 10 6:00 . . . 5.9 1.774952351 2 0
1 11 6:00 . . . 5.9 1.774952351 2 0
将事件格式数据输入到步骤二当中，计算得到群体药动学模型估计值。
得到样本的群体药动学估计值
Figure BDA0003852730960000181
后，扩展每个样本Xi＝{x1,x2…xk}使变为如下形式：
X′i＝{x1,x2…xk+1}作为步骤三机器学习的特征矩阵，且群体药动学估计值与实测值的残差
Figure BDA0003852730960000182
作为步骤三机器学习的标签。
使用Python中的xgboost框架将X,yresidual输入到定义好的XGBoost回归器当中，并通过10*10嵌套交叉验证以及网格搜索优化超参数，调整的超参数如下表：
Figure BDA0003852730960000183
Figure BDA0003852730960000191
其中boosting_type选择gbdt指的是使用梯度提升树方法；objective的值为reg:pseudohubererror即
Figure BDA0003852730960000192
n_estimators为426则表示树模型估计器有426个即有426棵树；eta为学习步长，在每次boosting之后收缩比例通过它增加权重，用于降低每棵树对整体模型的影响，将更多的空间留给后面的；max_depth为3，可以看到由于特征不多(本例子仅仅包含12个特征)，树的深度不会太大，防止过拟合；min_child_weight是一个阈值，在树分裂过程中，一个叶子节点的实例权重之和如果小于min_child_weight设置的值，则放弃继续分裂，这里的值为2.5，较为保守；lambda为权重的L2正则化项，这里的值为2.11，较为保守；subsample为对样本随机采样比率，实际训练样本N*＝N×subsample，每次树的迭代都会进行一次随机采样；colsample_bytree为构造每棵树时对列进行抽样的比率，构建每一棵树时对列进行一次抽样。在实验过程中的10*10嵌套交叉验证过程是将数据集分成10等份作为数据子集，并分为内外两层嵌套，内层嵌套使用其中9个子集进行训练，1个子集进行验证，同一超参数对应的模型但不同的数据子集进行验证需要进行10次，使用10次验证的模型评估指标均值(RMSE)进行网格搜素，选择到最优超参数后将该超参数进行外层嵌套验证，即为正常的10折交叉验证，完成模型评估。
完成训练调参工作后，得到局部最优模型。使用上述超参数的XGBoost回归器预测结果
Figure BDA0003852730960000201
加上群体药动学点估计值
Figure BDA0003852730960000202
得到最终预测值；在模型决策树解释模块计算XGBoost回归器对应的特征SHAP均值如图2所示。
2.3、本发明带来的有益效果：
通过结合基础群体药动学模型以及机器学习模型，利用药动学原理提供的先验知识以及梯度提升树的高预测效能，有效提高了模型的预测精度，同时保留了更强的可解释性。另外在建模时使用嵌套交叉验证能够在小数据集当中获得更加鲁棒的模型，更适合临床上得到的数据。
与现有技术使用集成软件进行建模不同，采用Python进行模型开发，能够广泛接入到其他各类系统当中，具有很强的工程泛用性。
利用可视化模块能够更直观地展示模型结果，便于决策者快速捕捉到数据提供的信息。
采用本系统进行肾移植患者口服他克莫司胶囊后的全血药物谷浓度真实世界数据(740个有效样本)试验，取得效果明显优于现有技术下的研究成果，其真实值与预测值散点图如图3所示：
Figure BDA0003852730960000203
对于本领域技术人员而言，显然本发明不限于上述示范性实施例的细节，而且在不背离本发明的精神或基本特征的情况下，能够以其他的具体形式实现本发明。因此，无论从哪一点来看，均应将实施例看作是示范性的，而且是非限制性的，本发明的范围由所附权利要求而不是上述说明限定，因此旨在将落在权利要求的等同要件的含义和范围内的所有变化囊括在本发明内。不应将权利要求中的任何附图标记视为限制所涉及的权利要求。

Claims (6)

1.一种全血药物谷浓度预测系统，其特征在于：包括原始数据采集模块、药动学计算模块、机器学习回归预测模块、预测结果输出模块以及模型决策树解释模块；
所述原始数据采集模块用于自医院系统中获取预测所需的相关数据；
所述药动学计算模块通过非线性混合效应模型法计算样本群体药动学参数以及群体预测值，同时通过最大后验贝叶斯计算点估计得到个体预测值；
所述机器学习回归预测模块通过手动输入需要预测的血药浓度时间，将其与药动学计算模块输出的结果一并输入到训练好的机器学习回归预测模型当中，计算输出群体药动学血药浓度点估计残差；
所述预测结果输出模块用于根据机器学习回归预测模块输出的结果对多个药物谷浓度预测数值进行理论药时曲线函数图像绘制，叠加群体参数曲线以及个体参数曲线，单个谷浓度预测直接显示数值；
所述模型决策树解释模块用于根据模型预测时计算出的SHAP值绘制出力图。
2.根据权利要求1所述的一种全血药物谷浓度预测系统，其特征在于：所述预测所需的相关数据包括患者的体征数据、相关检测结果以及基因监测结果，其中对于缺失的数据采用预设的典型值填充，对于采集的数据按照特征工程处理并暂存为Python字典格式数据一份，以及后续预测使用的事件格式数据一份后进入系统流程下一步，并同时存入研究数据库。
3.根据权利要求1所述的一种全血药物谷浓度预测系统，其特征在于：所述非线性混合效应模型法在拟合数据建模过程使用的是扩展的最小二乘法，其目标函数公式如下：
Figure FDA0003852730950000011
其中O(θ,ξ,σ2)为非线性混合效应模型法中关于θ,ξ,σ2的目标函数；σ2为残差变异的平方；θ为药动学参数；ξ为指数矢量参数，通常设置为0到2之间；n为样本量；yi为个体血药浓度实测值；xi为个体自变量(例如剂量、时间等)；f(xi,θ)为使用θ,xi计算血药浓度的药动学公式；u(xi,θ,ξ)为关于xi,θ,ξ的权重函数，计算为f(xi,θ)ξ；
通过定义药动学参数以及公式，代入到上述式中，进行ELS运算使目标函数O(θ,ξ,σ2)尽可能小(获得极小值)，同时获得对应的θ,ξ,σ2的集合。
4.根据权利要求3所述的一种全血药物谷浓度预测系统，其特征在于：所述药动学公式如下：
Figure FDA0003852730950000021
其中
Figure FDA0003852730950000022
为血药浓度估计值；ke为消除速率常数；X0为单次剂量；Clastdoseadj为上次调整剂量时的血药浓度实测值；τ为给药间隔；n为自上次调整剂量开始计算的给药次数。
5.根据权利要求4所述的一种全血药物谷浓度预测系统，其特征在于：所述药动学计算模块使用非线性混合效应模型法建立药动学模型，所述药动学模型的具体建模方法包括：
(1)仅采用最基础的个体自变量建立无其他信息的基础药动学模型，即仅采用时间、剂量对血药浓度进行建模，得到基础模型如下表示：
CL/F＝TV(θ1)×EXP(η1)
V/F＝TV(θ2)×EXP(η2)
其中，CL/F和V/F分别为表观清除率生物利用度比及表观分布容积生物利用度比；TV(θ1)和TV(θ2)分别为CL/F与V/F的群体典型值，η1与η2代表参数的个体间变异；
(2)随后将固定效应模型参数加入服从正态或对数正态分布的随机效应，对残差变异使用加法常系数混合随机效应，并认为随机效应同样存在于个体间变异中，得到如下公式：
Figure FDA0003852730950000031
Figure FDA0003852730950000032
其中ηi服从均值为0，方差为
Figure FDA0003852730950000033
的对数正态分布；εij服从均值为0，方差为
Figure FDA0003852730950000034
的正态分布。将上述模型结构代入到非线性混合效应模型法的对应参数当中，进行运算后，得到ω,σ2的集合。
6.根据权利要求1所述的一种全血药物谷浓度预测系统，其特征在于：所述机器学习回归预测模块的具体算法包括：
(1)预测函数
设存在数据集D＝(xi,yi)(|D|＝n,xi∈Rm,yi∈R)，m为特征数量，n为样本数量；另存在一个由K棵树组成的模型
Figure FDA0003852730950000035
其中
Figure FDA0003852730950000036
表示第k棵树对样本xi的预测结果；T表示该树叶子节点的数量；ωt为该树第t个叶子节点的权重，可理解为预测值；q(x)函数用于搜索样本xi所属叶子节点的下标；
(2)损失函数与树的分裂
定义损失函数
Figure FDA0003852730950000041
其中惩罚项
Figure FDA0003852730950000042
用于防止树的结构过于复杂造成过拟合；
Figure FDA0003852730950000043
为Pseudo-Huber损失函数；
设
Figure FDA0003852730950000044
是xi在第t次迭代的预测值，有
Figure FDA0003852730950000045
将其进行泰勒二次展开得到
Figure FDA0003852730950000046
由于
Figure FDA0003852730950000047
为常数，对函数无影响，因此可省略；
定义被分至第j个叶子节点的样本对应下标集合Ij＝{i|q(xi)＝j}，则
Figure FDA0003852730950000048
对ωj求导并计算最小值，得到
Figure FDA0003852730950000049
令I＝IL∪IR，得到衡量节点划分与否的公式，
Figure FDA0003852730950000051
公式每次都用不同的特征计算，得分最大的特征作为当前树节点的划分点；
(3)分割算法
定义Dk＝{(x1k,h1),(x2k,h2)...(xnk,hn)}表示训练样本第k个特征的值以及对应的二阶梯度；以此定义排名函数
Figure FDA0003852730950000052
并结合如下公式用来查找候选分割点，
Figure FDA0003852730950000053
以下是近似分割点查找算法伪代码：
For k＝1to m do
查找第k个特征的候选分割点记为Sk＝{sk1,sk2...skl}；
可按树进行计算(Global)也可按每次分割进行计算；
End for
For k＝1to m do
Figure FDA0003852730950000054
Figure FDA0003852730950000055
End for
算法输入为原始数据的K列特征、群体药动学预测模块估计的残差以及标签列，即共有K+1列特征以及1列标签，进行交叉验证以及网格搜索调参训练后得到XGBoost的残差预测模型，该模型预测结果加上步骤二群体药动学预测模块点估计值得到总预测结果。
CN202211139186.XA 2022-09-19 2022-09-19 一种全血药物谷浓度预测系统 Pending CN115472308A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202211139186.XA CN115472308A (zh) 2022-09-19 2022-09-19 一种全血药物谷浓度预测系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202211139186.XA CN115472308A (zh) 2022-09-19 2022-09-19 一种全血药物谷浓度预测系统

Publications (1)

Publication Number Publication Date
CN115472308A true CN115472308A (zh) 2022-12-13

Family

ID=84332557

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202211139186.XA Pending CN115472308A (zh) 2022-09-19 2022-09-19 一种全血药物谷浓度预测系统

Country Status (1)

Country Link
CN (1) CN115472308A (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116504413A (zh) * 2023-04-11 2023-07-28 长沙市妇幼保健院(长沙市妇幼保健计划生育服务中心) 人工智能麻醉管理系统

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116504413A (zh) * 2023-04-11 2023-07-28 长沙市妇幼保健院(长沙市妇幼保健计划生育服务中心) 人工智能麻醉管理系统
CN116504413B (zh) * 2023-04-11 2024-04-26 长沙市妇幼保健院(长沙市妇幼保健计划生育服务中心) 人工智能麻醉管理系统

Similar Documents

Publication Publication Date Title
CN109033738B (zh) 一种基于深度学习的药物活性预测方法
CN109994200B (zh) 一种基于相似度融合的多组学癌症数据整合分析方法
CN110957002B (zh) 一种基于协同矩阵分解的药物靶点相互作用关系预测方法
CN111090764B (zh) 基于多任务学习和图卷积神经网络的影像分类方法及装置
Boyko et al. Use of machine learning in the forecast of clinical consequences of cancer diseases
CN112001218A (zh) 一种基于卷积神经网络的三维颗粒类别检测方法及系统
CN106055922A (zh) 一种基于基因表达数据的混合网络基因筛选方法
CN105139430A (zh) 一种基于图熵的医学图像聚类方法
CN115472308A (zh) 一种全血药物谷浓度预测系统
CN116580848A (zh) 一种基于多头注意力机制的分析癌症多组学数据方法
Dutta Detecting Lung Cancer Using Machine Learning Techniques.
CN116959585B (zh) 基于深度学习的全基因组预测方法
CN113283465B (zh) 一种弥散张量成像数据分析方法及装置
Nurdin et al. Performance comparison of hybrid CNN-XGBoost and CNN-LightGBM methods in pneumonia detection
Hamid New location model based on automatic trimming and smoothing approaches
CN114818969A (zh) 变点位置及类别的检测方法及装置
CN108304546B (zh) 一种基于内容相似度和Softmax分类器的医学图像检索方法
CN113296947A (zh) 基于改进XGBoost模型的资源需求预测方法
CN112784886A (zh) 一种基于多层最大生成树图核的脑图像分类方法
CN113971984A (zh) 分类模型构建方法及装置、电子设备、存储介质
Wahyu et al. The application of particle swarm optimization using Naive Bayes method for predicting heart disease
CN116226629B (zh) 一种基于特征贡献的多模型特征选择方法及系统
Zhao et al. Class-specific Prompts in Vision Transformer for Continual Learning of New Diseases
Zhang et al. An automatic detection model of pulmonary nodules based on deep belief network
CN116597902B (zh) 基于药物敏感性数据的多组学生物标志物筛选方法和装置

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination