利用易于获取的患者数据对1型糖尿病的未来血糖进行个性化模型识别和预测的贝叶斯框架的制作方法

利用易于获取的患者数据对1型糖尿病的未来血糖进行个性化模型识别和预测的贝叶斯框架
1.相关申请的交叉引用
2.本技术要求于2021年1月25日提交的美国临时专利申请第63/141345号的权益。本技术涉及均于2020年12月4日提交的美国申请第17/112870号、第17/112828号和第17/112814号。上述申请的内容全文以引用方式并入本文，并且因此明确成为本说明书的一部分。


背景技术：

3.1型糖尿病(t1d)是一种由于胰腺中β细胞的进行性破坏引起的慢性自身免疫疾病，其使机体不能产生内源性胰岛素。因此，血糖(bg)浓度往往会超过高血糖阈值(bg》180mg/dl)，如果这种情况频繁出现且长期存在，可能会导致几种严重的心血管长期并发症，以及肾病和神经病。为了降低血糖(bg)水平，必须每天数次施用外源性胰岛素。不幸的是，过量的外源性胰岛素给药可能使患者出现低血糖症，即bg《70mg/dl，由于低血糖可能产生昏厥、头晕、昏迷甚至死亡，因此即使在短期内这也是危险的。
4.有效的t1d治疗依靠于通过传统的手指针刺设备或更现代的微创型动态血糖监测(cgm)传感器进行的频繁的血糖(bg)监测，而且并不容易实现。事实上，从患者的角度来看，t1d管理是一个终身学习过程，从而了解数个日常因素(例如，生病、饮食和体能活动)如何影响bg水平以及如何使用干预措施(例如，减少碳水化合物摄入量当然以及施用胰岛素)将bg保持在安全范围内。
5.在这方面，研究界已经做出了许多努力来提供能够帮助t1d患者的新工具。其中，通过使主动治疗决策基于预期的未来葡萄糖水平而不是当前的葡萄糖水平，能够实时预测未来bg浓度的基于cgm的算法具有显著改善t1d治疗效果的潜力。


技术实现要素：

6.在第一方面，提供了一种预测个体患者的未来血糖浓度的方法，该方法包括：选择个体化非线性葡萄糖-胰岛素动力学生理模型，该选择的模型具有多个模型参数，该模型参数的值待确定；估计所述多个模型参数中的每个模型参数的值，通过将参数估计技术应用于所述个体患者的先验信息和数据以获得后验信息，所述模型参数的第一子集具有按照先验群体数据估计的值，并且所述模型参数的第二子集具有针对所述个体患者个性化的值；以及使用模型参数中的每个模型参数的估计值将非线性预测技术应用于选择的模型，以获得个体患者在未来某一时间的预测血糖浓度。
7.在第一方面的实施方案中，参数估计技术是马尔可夫链蒙特卡尔(mcmc)贝叶斯估计(bayesian estimator)。
8.在第一方面的实施方案中，参数估计技术选自由最大后验法、最大似然法和预测误差最小化法组成的组。
9.在第一方面的实施方案中，非线性预测技术是粒子滤波。
10.在第一方面的实施方案中，非线性预测技术选自由扩展卡尔曼滤波技术和无迹卡尔曼滤波技术组成的组。
11.在第一方面的实施方案中，选择的个体化非线性生理模型包括皮下胰岛素吸收子模型、口服葡萄糖吸收子模型和葡萄糖-胰岛素动力学模型。
12.在第一方面的实施方案中，选择的个体化非线性生理模型使用个体患者的过去血糖浓度水平、碳水化合物摄入量信息和外源性胰岛素数据作为输入。
13.在第一方面的实施方案中，该皮下胰岛素吸收子模型包括第一模型参数和第二模型参数，第一模型参数指定胰岛素分布的值，第二模型参数指定胰岛素在第一房室中出现的延迟，第一和第二模型参数位于具有按照先验群体数据估计的值的模型参数的第一子集中。
14.在第一方面的实施方案中，该皮下胰岛素吸收子模型包括第三模型参数和第四模型参数，第三模型参数指定从第一房室到第二房室的扩散速率常数，第四模型参数指定从第二房室到血浆的皮下胰岛素吸收速率常数，第三和第四模型参数位于具有针对个体患者个性化的值的模型参数的第二子集中。
15.在第一方面的一个实施方案中，选择的个体化非线性生理模型是最大生理模型，其中忽略了模型参数的选择的参数。
16.在第一方面的实施方案中，使用实时动态血糖监测仪(cgm)数据、碳水化合物摄入量信息和外源性胰岛素数据实时更新预测血糖浓度。
17.在第一方面的实施方案中，该方法还包括用残差模型增强选择的个体化非线性生理模型，从而描述残差建模误差。
18.在第一方面的实施方案中，残差建模误差被建模成阶数为(5,5,1)的自回归综合移动平均(arima)。
19.在第一方面的实施方案中，残差模型具有使用预测误差法(pem)估计的模型参数。
20.在第一方面的实施方案中，在多个不同时间对预测血糖浓度进行预测。
附图说明
21.图1是高层次示出本文所述的基于个体化模型的预测方案中涉及的步骤的流程图。
22.图2示出了皮下胰岛素吸收子系统方案的一个示例。
23.图3示出了口服葡萄糖吸收子系统方案的一个示例。
24.图4示出了葡萄糖-胰岛素动力学子系统方案的一个示例。
25.图5示出了残差模型方案的一个示例。
26.图6a示出了对于15min多至120min的不同ph使用arx、arimax和pf获得的测定系数(cod)；并且图6b示出了对于15min多至120min的不同ph使用arx、arimax和pf获得的均方根误差(rmse)。
具体实施方式
27.i.引言
28.尽管已经提出了用于血糖(bg)实时预测的几种不同的方法论，但是由于改变bg水
平的几种未知的干扰因素以及葡萄糖生理学中的高患者间/患者内可变性，该问题仍然未解决。现有预测技术依靠着描述某组输入特征与bg之间的关系的一种数学模型。选择葡萄糖-胰岛素调节的合适模型是一个重要步骤，并且选项范围包括从完全数据驱动的“黑盒”模型到基于代谢生理学的机械/半机械非线性描述的模型。与黑盒模型相比，生理模型具有保证正确捕捉以下的葡萄糖-胰岛素代谢机制的明显优点。在t1d文献中可获得的非线性生理模型中，所谓的最小模型被证明过于僵化和简单化以至于无法给予有效的葡萄糖预测。另一方面，大型模型有数个方程和多个参数，这使得它们难以识别并且更适合计算机模拟而不是基于模型的预测。
29.图1是高层次示出本文所述的基于个体化模型的预测方案中涉及的步骤的流程图。如图所示，方案中的第一步骤10为从上述各种可用选项中选择合适的模型结构。选择的模型是葡萄糖-胰岛素调节的个体化非线性生理模型。如下所述，该选择的模型具有模型参数，该模型参数可以根据易获得的患者记录诸如碳水化合物摄入量和胰岛素输注数据进行个体化，以模拟个体的特定生理机能。
30.在图1所示的预测方案的第二步骤20中，使用参数估计技术估计从针对个体患者待个性化的模型参数中选择的参数。在一个实施方案中，该使用的参数估计技术采用了由马尔可夫链蒙特卡尔(mcmc)实施的贝叶斯方法，以在存在复杂非线性动态的情况下估计大量的参数。接下来，在第三步骤30中，在基于例如粒子滤波方法论的非线性预测方案中使用所获得的个性化模型。
31.如将在下文更详细进一步描述的，使用包括由以下文献中描述的模拟模型生成的100个虚拟t1d患者的数据的计算机模拟临床试验测试本文所述的预测方案：r.visentin,e.campos-m.schiavon,d.lv,m.vettoretti,m.breton,b.p.kovatchev,c.dalla man和c.cobelli,“the uva/padova type 1diabetes simulator goes from single meal to single day,”j.diabetes sci.technol.，第12卷，编号2，第273
–
281页，2018，该文献全文以引用并入本文。为了进行比较，考虑了两个基线黑盒模型。所选择的用于评估的度量是均方根误差、测定系数和时间增益。结果表明，所提出的方法在预测未来葡萄糖浓度方面优于两种基线方法论，特别是在时间增益方面。
32.ii.葡萄糖-胰岛素调节的生理模型
33.该部分描述了用于个性化葡萄糖预测的非线性葡萄糖-胰岛素动态生理模型的一个示例。该示例从如上述引用的visentin等人的参考文献中描述的最近版本的uva/padova t1d模拟器(t1ds)中实施的最大生理模型开始，并且使该最大生理模型尽可能地简化，以减少用于个体化的待估计的参数的数量，同时保持其实现良好的葡萄糖预测的能力。
34.所提出的模型具有两个输入，即胰岛素输注i(t)和碳水化合物摄入量cho(t)，以及一个输出，即间质葡萄糖浓度ig(t)。该模型由描述皮下胰岛素吸收、口服葡萄糖吸收和葡萄糖-胰岛素动力学的三个子系统组成。
35.a.皮下胰岛素吸收子系统
36.皮下胰岛素吸收模型是以下文献中描述的t1ds模型中所用吸收模型的略微简化版本：m.schiavon,c.dalla man和c.cobelli,“modeling subcutaneous absorption of fast-acting insulin in type 1diabetes,”ieee trans.biomed.eng.，第65卷，编号9，第2079
–
2086页，2018，该模型在图2中示出。该模型由三个房室组成。将外源性胰岛素i输注到
第一房室，其在第一房室中延迟β后出现。在具有非单体状态的胰岛素的第一房室中，胰岛素转化为单体状态，然后扩散到血浆中。模型方程是：
[0037][0038]
其中，i
sc1
(mu/kg)和i
sc2
(mu/kg)分别表示非单体状态和单体状态的胰岛素；i
p
(mu/l)是血浆胰岛素浓度；kd(min-1
)是从第一房室向第二房室的扩散速率常数；k
a2
(min-1
)是从第二房室到血浆的皮下胰岛素吸收速率常数；ke(min-1
)是部分清除率；vi(l/kg)是胰岛素分布体积；β(min)是第一房室中胰岛素出现的延迟。相对于t1ds中使用的模型，这种简化在于忽略了一部分的非单体胰岛素可以直接到达血浆的事实。该吸收路径在图1中以灰色示出并且与速率常数k
a1
相关联。简化的原因是胰岛素分数和发生这种情况的患者的数量都很小。在表1中报告的关于模型参数的先验信息已从文献中获得。具体地，将vi和β设定为群体值，即分别为0.126l/kg和8min。此外，kd被限制为kd≥k
a2
，因为这两种组合是可互换的。未知模型参数为θ
ins
＝[k
a2
，kd]。
[0039]
b.口服葡萄糖吸收子系统
[0040]
来自以下文献的口服葡萄糖吸收子系统模型是t1ds中所用模型的简化版本：c.dalla man,m.camilleri和c.cobelli,“a system model of oral glucose absorption:validation on gold standard data,”ieee trans.biomed.eng.，第53卷，编号12，第2472
–
2478页，2006，并且该模型在图3中示出。该模型将胃肠道描述为三房室系统：前两个房室负责胃中的食物(固体和研磨状态)，而第三个房室将吸收了cho的小肠上端建模。模型方程是：
[0041][0042]
其中，q
sto1
(mg/kg)和q
sto2
(mg/kg)分别是胃中固态和液态的葡萄糖的量；q
gut
(mg/kg)是肠道中的葡萄糖浓度；k
gri
(min-1
)是研磨速率常数；k
empt
(min-1
)是胃排空速率常数；k
abs
(min-1
)是肠道吸收速率常数；cho(mg/kg/min)是碳水化合物摄入率。根据模型(2)估计血浆中葡萄糖出现的速率ra(mg/kg/min)，其为：
[0043]
ra(t)＝f.k
abs
·qgut
(t)
ꢀꢀꢀ
（3）
[0044]
其中，f(无量纲)是血浆中吸收的肠道内容物的分数。
[0045]
相对于t1ds中使用的模型，简化之处在于假定胃排空率不变，从而忽略了其对胃内容物的依赖性(在图3中以灰色示出)。
[0046]
从文献中获得了模型(2)的先验信息，并且已在表1中详述。特别地，设定f等于0.9并且限制k
gri
＝k
empt
。此外，将k
abs
限制为k
abs
≤k
empt
，因为这两种组合是可互换的。因此，未知模型参数是θo
ral
＝[k
abs
,k
empt
]。
[0047]
c.葡萄糖-胰岛素动力学子系统
[0048]
关于在t1ds模拟器中使用的模型，我们使用了一种更简洁的描述，即葡萄糖-胰岛
素动力学子系统。该子系统基于一个已知的二房室模型，该二房室模型描述了血浆胰岛素作用和葡萄糖出现率对血浆葡萄糖浓度的影响，该影响在以下文献中介绍：r.n.bergman,c.r.bowden和c.cobelli,“quantitative estimation of insulin sensitivity,”am.j.physiol.，第236卷，编号6，第e667
–
e677页，1979。该模型进一步配有第三房室以描述葡萄糖从血浆到间质的运输，在间质中由传感器测量葡萄糖。该模型在图4中示出。模型方程是：
[0049][0050]
其中，g(mg/dl)是血浆葡萄糖浓度，x(min-1
)是胰岛素对葡萄糖处置和生产的作用；sg(min-1
)是描述葡萄糖本身促进葡萄糖处置和抑制葡萄糖产生的能力的葡萄糖有效性；gb(mg/dl)是血浆中的基础葡萄糖浓度；vg(dl/kg)是葡萄糖分布体积；p2(min-1
)是胰岛素作用动力学的速率常数；si(ml/μu
·
min)是胰岛素敏感性；i
pb
(mu/l)是血浆中的基础胰岛素浓度；ig(mg/dl)是间质葡萄糖浓度；α(min)是血浆与间质葡萄糖浓度房室之间的延迟。
[0051]
已知具有恒定单位ρ，的上述模型很难捕捉低血糖，这可能是由于对胰岛素作用的描述不充分，该胰岛素作用被证明在葡萄糖降低到某个阈值以下时会增加。为此，按照以下文献中提出的原理在上述模型中引入术语ρ(g)：c.dalla man,f.micheletto,d.lv,m.breton,b.kovatchev和c.cobelli,“the uva/padova type 1diabetes simulator:new features,”j.diabetes sci.technol.，第8卷，编号1，第26
–
34页，2017：
[0052][0053]
其中，g
th
是低血糖阈值(设定为60mg/dl)，并且r1(无量纲)和r2(无量纲)是没有直接进行生理解释的两个模型参数。
[0054]
为了解释患者特定的日内胰岛素敏感性变异性，并对所谓的黎明现象进行建模，参数si和gb被视为在一天内随时间变化：
[0055][0056][0057]
在表1中详细报告的关于参数分布的先验信息系从文献中获得。详细地，r1、r2和vg已固定为群体值，即分别为1.44dl/kg、0.81dl/kg和1.45dl/kg。葡萄糖-胰岛素子系统的未
知模型参数是知模型参数是
[0058]
d.整体生理模型
[0059]
综上所述，通过组合目前介绍的三个子模型得到整体生理模型。
[0060][0061]
其中，x
phy
(t)是状态向量，其定义如下：
[0062]
x
phy
(t)：＝[x
ins
，|x
oral
，|x
glu
]
t
[0063]
＝[i
sc1
，i
sc2
，i
p
，|q
sto1
，q
sto2
，q
gut
，|g，x，ig]
t
；
[0064]uphy
(t)：＝[i(t)，cho(t)]是输入向量；
[0065]fphy
(
·
)是组合(1)-(5)所获得的状态更新函数。
[0066]fphy
(
·
)取决于一组未知参数
[0067]
θ
phy
：＝[θ
ins
，θ
oral
，θ
glu
]
[0068]
该组参数的估计将在下一部分中讨论。
[0069]
表1：模型参数的先验信息
[0070][0071]
n(μ,σ)代表均值μ和标准偏差σ的正态分布。logn(μ,σ)代表均值μ和标准偏差σ的对数正态分布。gamma(a,b)代表具有形状参数a和尺度参数b的gamma分布。
[0072]
iii.通过马尔可夫链蒙特卡尔理论对个性化模型参数进行离线贝叶斯估计
[0073]
通过使用训练数据y:＝{cgm(tk),tk＝k
·
ts,k＝1,
…
,d}和u:＝{u
phy
(tk),tk＝k
·
ts,k＝1,
…
,d}识别每个患者的未知模型参数θ
phy
，从而进行模型个性化，其中d是可用数据点的数量。
[0074]
通过采用贝叶斯方法进行识别，并且具体地，在这项工作中，θ
phy
通过其后验均值进行估计，该后验均值定义为
[0075][0076]
已知后验均值是θ
phy
的最小变异数不偏估计量。
[0077]
根据贝叶斯定理，获得后验密度函数p
θ|y,u
(θ|y,u)，其为：
[0078]
[0079]
其中，p
y|θ,i
(y|θ,u)是似然函数，即在给定参数向量θ和输入u的情况下观察到某个y的概率。
[0080]
由于(9)中的积分难以分析，因此必须借助mcmc对它近似求解。特别地，通过创建马尔可夫链，根据后验分布p
θ|y,u
(θ|y,u)生成n个样本θi,i＝1,
…
,n，该马尔可夫链的静态分布正好是该后验分布(目标分布)。然后，这些样本θi用于执行蒙特卡洛积分以获得θ
phy
的点估计值：
[0081][0082]
为了构建这样一个链，使用了单分量梅特罗波利斯－黑斯廷斯(scmh)算法。
[0083]
a.实施细节
[0084]
将θ
phy
划分成五组θ
phy
:＝[θ1,θ2,θ3,θ4,θ5]，即]，即θ3:＝[sid],θ4:＝[p2,k
a2
,kd],θ5:＝[k
empt
,k
abs
]。
[0085]
选择这种划分方案是因为它改善了mc混合并且允许打破si与p2之间的相关性，从文献中可知它很重要。
[0086]
该算法的迭代i由五个步骤p＝1,
…
,5组成，并且每个步骤通过批准/拒绝从提议密度函数q
p
(
·
|
·
)中提取的样本φ
p
来更新θ
phy
的第p个分区θ
p
。具体地，如shcm程序所规定的，批准发生的概率为α
[0087][0088]
鉴于其它分量θ-p
取值θ-p
＝θ
i,-p
，π(θ
p
|θ
i,-p
)与θ
p
的后验成正比：
[0089]
π(θ
p
|θ
i,-p
)＝p
y|θ,u
(y|θ
p
,θ
i,-p
,u)p
θ
(θ
p
|θ
i,-p
,u)
[0090]
θ
i,-p
包括除了θ
p
之外的θ
phy
的所有其它分量。准确的说，θ
i,-p
包含在算法的当前阶段可用的每个分量的最新版本：θ
i,-p
＝[θ
i,1
,
…
,θ
i,p-1
,θ
i-1,p+1
,θ
i-1,5
]。当在迭代i处理第p个分量时，到p-1的分量已经更新，而从p+1到5的其它分量还没有更新，所以使用在先前迭代i-1中计算的它们的值。
[0091]
对于它所涉及的建议分布，使用高斯分布，它以在前一个的链迭代中由θ
p
假设的值为中心
[0092][0093]
其中，σ
p
是调节链的接受率的调整参数。将σ
p
设定为对角矩阵，其分量是分区p,的每一元素的条件标准偏差的估计值乘以如h.haario、e.saksman和j.taminen，“adaptive proposal distribution for random walk metropolis algorithm”，computational statistics，第14卷，第375-395页，1999年中所提的比例因子通过运行两个探索性mcmc(niter＝600迭代)计算该估计值，并且每1500次迭代更新一次该估计值，从而实现自适应shcm。
[0094]
最后，通过已知的raftery-lewis标准验证mcmc的收敛性，提供确保马尔可夫链表
示目标后验分布所必需的迭代次数。
[0095]
在算法(1)中总结了自适应单分量梅特罗波利斯－黑斯廷斯算法。
[0096][0097]
iv.残差模型
[0098]
作为最后一步，通过使用图5的残差模型方案以随机误差e(t)破坏生理模型ig(t)的输出，获得测量的输出cgm(t)。
[0099]
cgm(t)＝ig(t)+e(t)。
ꢀꢀꢀꢀꢀ
(12)
[0100]
e(t)(mg/dl)描述了残留建模误差和cgm传感器引入的测量误差。建模成阶数为(5,5,1)的自回归综合移动平均(arima)是一个有色随机过程：
[0101][0102]
其中，{ai}
i＝1,
…
,6
和{bi}
i＝1,
…
,5
是模型的系数，并且∈(t)(mg/dl)是具有零均值和恒定标准偏差sd
∈
的白噪声。由于考虑到arima，分母的零点之一在1中。
[0103]
选择模型阶数以保持模型尽可能简单但仍然能够很好地描述残差e(t)。使用标准预测误差法(pem)估计模型参数。
[0104]
a.增强模型
[0105]
用残差模型的状态空间表示增强(8)，然后得到最终模型。
[0106][0107]
其中，通过使用具有tk＝k
·
ts的后向欧拉法对f
phy
(x
phy
,u
phy
,t,θ
phy
)进行离散化获
得fd，并且进一步假定其被模型误差v_phy破坏，该误差是一个具有协方差的i.i.d高斯过程。用e即与描述残差(13)的arima过程的状态空间表示相关联的状态向量增强该模型，其中矩阵a、b和c如g.e.p.box、g.m.jenkins、g.c.reinsel，和g.m.ljung，“timeseries analysis:forecasting and control”，第5版，hoboken，nj:wileyand sons，2015年中所定义。
[0108][0109]
c＝[1 0 0 0 0]
[0110]
最后，通过定义增强状态x(tk):＝[x
phy
(tk)|e(tk)]
t
和增强错误v(tk):＝[v
phy
(tk)|b
·
∈(tk)]
t
，将方程(14a)和(14b)进行组合。
[0111]
增强模型更新方程读数
[0112]
x(t
k+1
)＝f(x(t
k+1
),u(t
k+1
),t,θ)+v(tk)
ꢀꢀꢀ
(15)
[0113]
v.通过粒子滤波进行实时葡萄糖预测
[0114]
至此，我们着重估计t1d受试者的生理模型的参数，以捕捉患者特异性动力学。可以使用从回顾性患者数据中获得的可用信息“离线”完成该过程。
[0115]
在该部分中，将呈现如何使用这种个性化模型实时预测未来葡萄糖浓度。该任务必须“在线”执行，因此借助了顺序算法，即在每个时间步长tk，当新的测量值y(tk)＝cgm(tk)变得可用时，更新模型状态x(tk)的当前估计值并且使用它来推断未来葡萄糖浓度。特别地，采用粒子滤波(pf)，其是能够处理模型的非线性结构的一种顺序贝叶斯预测技术。
[0116]
pf基于后验概率函数p(x(tk)|y(t
1:k
),u(t
1:k
))的递归更新，其中y(t
1:k
)是变量y(t1),
…
,y(tk)的简写，并且u(t
1:k
)表示u(t1),
…
,u(tk)。
[0117]
通过两个基本步骤，即单步超前预测和测量值更新，执行p(x(t
k-1
)|y(t
1:k-1
),u(t
1:k-1
))的递归更新。
[0118]
单步超前预测步骤假设在时间t
k-1
可获得后验概率p(x(t
k-1
)|y(t
1:k-1
),u(t
1:k-1
))，并且使用这种后验概率推断
[0119]
p(x(tk)|y(t
1:k-1
),u(t
1:k
))
[0120]
然后，当新的测量值y(tk)在时间tk变得可用时，在测量值更新步骤中使用这种测量值计算后验概率
[0121]
p(x(tk)|y(t
1:k
),u(t
1:k
))。
[0122]
然后，针对数据集中的每个可用的测量值，重复这两个步骤。pf通过使用起作用的概率函数的采样近似值执行这些步骤：
[0123][0124]
其中，是一组p点，称为“粒子”，以p(x(t
k-1
)|y(t
1:k-1
),u(t
1:k-1
))为支撑。每个粒子都与权重∑
pwp
(t
k-1
)＝1相关联，并且
[0125][0126]
其中，是表示p(x(tk)|y(t
1:k-1
),u(t
1:k
))的p粒子，每一个p粒子与权重∑
pw*p
(tk)＝1相关联。
[0127]
当处理诸如(15)的状态空间模型时，pf的单步超前预测步骤特别方便。
[0128]
可以看出
[0129]
p(x(tk)|y(t
1:k-1
),u(t
1:k
))＝∫p(x(tk)|x(t
k-1
))p(x(t
k-1
)|y(t
1:k-1
),u(t
1:k-1
))dx(t
k-1
)
[0130]
其中，通过(15)的状态更新方程充分描述p(x(tk)|x(t
k-1
))。
[0131]
关于测量值更新步骤，可以证明它是成立的p(x(tk)|y(t
1:k
),u(t
1:k
))
∝
p(y(tk)|x(tk),u(t
1:k
))p(x(tk)|y(t
1:k-1
),u(t
1:k
))
[0132]
其中，p(y(tk)|x(tk),u(t
1:k
))是由(14c)完全指定的似然函数。
[0133]
在实施过程中这些量的详细计算方式参见“实施细节”部分。
[0134]
附加的结果为，后验概率p(x(tk)|y(t
1:k
),u(t
1:k
))进一步被pf用来计算后验概率
[0135][0136]
其描述了i步超前到ph步超前及时预测的状态分布(其中ph为预测时域)。
[0137]
特别地，对于如(15)的状态空间模型，可以示出
[0138][0139]
其由方程(14c)完全指定。
[0140]
最后，可以使用后验的期望值导出在时间t_(k+i),i＝1,
…
,ph的未来cgm值的点估计值：
[0141][0142]
a.实施细节
[0143]
在下文中提出了通过pf实施的数值方案，以执行单步超前预测、测量值更新和多步超前预测。
[0144]
单步超前预测步骤。回顾一下，在时间t
k-1
，p(x(t
k-1
)|y(t
1:k-1
),u(t
1:k-1
))可以以
具有相关权重∑
p
w(t
k-1
)p＝1的一组p粒子定义的采样形式获得，使得
[0145][0146]
pf通过从p(x(tk)|x(t
k-1
))抽取一组新的粒子来执行单步超前预测步骤：
[0147]
x
p
(tk)～p(x
p
(tk)|x
p
(t
k-1
))。
ꢀꢀꢀꢀ
(16)
[0148]
由(15)指定该概率：
[0149]
p(x
p
(tk)|x
p
(t
k-1
))～n(f(x
p
(t
k-1
),u,t
k-1
,θ),σv)。
[0150]
鉴于此，为了抽取一组新的粒子，只需让每个粒子x
p
(t
k-1
)根据模型(15)进行演变并且以噪声v的实现来破坏它。
[0151]
测量值更新步骤。pf算法将每个第p个粒子x
p
(tk)的权重w
*p
(tk)设定为在x
p
(tk)
[0152]w*p
(tk)＝p(y(tk)|x
p
(tk),u(t
1:k
))。上估计的似然函数。
ꢀꢀꢀꢀꢀ
(17)
[0153]
特别地，根据方程(14c)和随机建模误差e的统计，p(y(tk)|x
p
(tk),u(t
1:k
))定义为：
[0154]
p(y(tk)|x(tk),u(t
1:k
))＝n(y(tk)-y
p
(tk),sd
∈
)。
ꢀꢀꢀꢀꢀꢀ
(18)
[0155]
其中使用(14c)获得y
p
(tk)。
[0156]
然后对权重进行归一化，使得∑
pw*p
(tk)＝1。
[0157]
这提供了后验密度的采样形式表示
[0158][0159]
再采样步骤。为了提高pf的准确性，通过更新一组粒子来完成测量值更新步骤。具体地，用一组新的p粒子，即由p(x(tk)|y(t
1:k-1
),u(t
1:k
))的采样表示生成的取代使得pr(x
*p
(tk)＝x
p
(tk))＝w
*p
(tk)。通过以下文献中所述的公认的重采样算法执行该步骤：g.kitagawa,“monte carlo filter and smoother for non-gaussian non-linear state space models,”j.comput.graph.statist.，第5卷，编号1，第1
–
25页，1996。
[0160]
因此，所有的新粒子与同一权重w
*p
(tk)＝1/p相关联，从而使p(x(tk)|y(t
1:k-1
),u(t
1:k
))的近似简化成
[0161]
[0162]
多步超前预测。可按如下方式获得多步超前预测值。首先，获得采样形式的概率该采样形式从p(x(tk)|y(t
1:k
),u(t
1:k
))开始，并且如在单步超前预测步骤中那样，i步超前传播p粒子
[0163]
然后，使用方程(14c)计算一组粒子并且将其用于获得采样形式的
[0164]
最后，获得葡萄糖i步超前的点估计y(t
k+i
)作为对p(y(t
k+i
)|y(t
1:k
),u(t
1:k+i
))的采样形式计算的平均值，即：
[0165][0166]
所实施的pf概括在算法(2)中。
[0167][0168]
vi.方法论评估
[0169]
a.模拟数据集
[0170]
t1ds是已被fda批准的一种软件工具，其替代闭环胰岛素输送策略和开环疗法的临床前试验。它提供了100个虚拟成年受试者的群体，每个受试者的特征在于生理参数的向量是不相同，从而捕获t1d群体的受试者间/受试者内可变性。在本文所述的该特定示例中，
使用了模拟器的最新版本，其能够更好地解释患者间/患者内的可变性(即，患者生理和患者行为随时间的变化)。使用t1ds为100个虚拟成人产生两周的合成数据，模拟所谓的传感器增强型胰岛素泵疗法(sap)。在该治疗方案中，通过胰岛素泵向患者施用胰岛素，该胰岛素的施用由两部分组成：第一部分是基础胰岛素输注速率，其是由临床医生针对每个患者调整的每日固定模式；第二部分由在每餐时施用的多个胰岛素剂量组成以抵消相应的葡萄糖增加。注射的胰岛素量与估计的碳水化合物摄入量(cho)和餐前血糖测量值成正比。
[0171]
在模拟场景中，每个受试者每天吃3餐，在时间隔为[06:30、08:30](早餐)、[12:00、14:00](午餐)和[19:30、21:30](晚餐)的情况下以均匀的概率发生，cho消耗量分别均匀分布在[20,60]g、[50,90]g和[50,90]g中。实施一个简单的安全方案以治疗低血糖症，即在患者测得的葡萄糖低于70mg/dl的情况下，每15min消耗15g速效cho。最后，根据以下文献中的建议对cgm传感器的测量误差进行建模：a.facchinetti,s.del favero,g.sparacino和c.cobelli,“model of glucose sensor error components:identification and assessment for new dexcom g4 generation devices,”med.biol.eng.comput.，第53卷，编号12，第1259
–
1269页，2015。将获得的计算机数据集划分成训练集和测试集，其中将第一周的数据分配给前者并且将第二周的数据分配给后者。
[0172]
b.用于比较的预测方法
[0173]
考虑了两个最先进的黑盒模型作为预测准确性的基线，在相同的训练数据上训练阶数为5的具有外源输入的自回归(arx)模型和阶数为(5,5,1)的具有外源输入的自回归综合移动平均(arimax)模型，两者都具有两个输入，胰岛素输注i(t)和碳水化合物摄入量cho(t)，以及一个输出，葡萄糖浓度cgm(t)。
[0174]
c.评估度量
[0175]
通过使用不同的ph将所获得的葡萄糖预测值与实际cgm值进行比较，在测试集中进行预测准确性的评估。在该项工作中，考虑了三个度量：均方根误差(rmse)、测定系数(cod)和时间增益(tg)。在文献中经常使用的这些度量定义如下：
[0176][0177]
其中，||
·
||2表示l2范数。rmse越大，预测准确性越差。
[0178][0179]
其中，是血糖测量时间序列的样本均值。cod取决于针对信号样本方差归一化的rmse。在完美预测的情况下(rmse＝0mg/dl)，cod最多为cod＝100％，并且随着rmse的增加而降低。
[0180][0181]
其中，延迟(s1(tk),s2(tk))[min]对两个信号s1(tk)和s2(tk)之间的延迟进行量化，并且基于互相关(xcorr)，用于测量两个信号的相似性。具体地，该延迟定义为时域偏移τ，其使s1与s2的τ偏移版本之间的互相关(xcorr)最大化：
[0182]
延迟(s1(tk),s2(tk))＝argmax
τ
[xcorr(s1(tk),s2(t
k-τ)]
1.47mg/dl(-6％)、-3.85mg/dl(-15％)、-7.83mg/dl(-23％)和-11.78mg/dl(-28％)。当从平均δcod方面评估预测准确性时，同样的考虑也是成立的，在ph＝30、60、90和120min时，平均cod分别提高了2.19％(2％)、9.86％(14％)、20.30％(40％)和32.99％(118％)。针对δtg，在ph＝30、60、90和120min时pf将结果分别显著提高15.0min(100％)、40.0min(200％)、67.5min(300％)、95.0min(380％)。
[0190]
将pf与arimax进行比较，前者在ph》30min时实现了明显更好的rmse和cod，ph越高，提高程度越大。详细地，在ph＝60、90和120min时，通过pf将(p《0.05)δrmse＝rmse
pf-rmse
arimax
分别显著地提高了-0.54mg/dl(-2％)、-2.78mg/dl(-10％)和-5.81mg/dl(-16％)，并且在ph＝60、90和120min时，通过pf将(p《0.05)δcod分别显著地提高了3.17％(4％)、10.62％(18％)和20.70％(51％)。在ph＝30min时，性能没有显著统计学差异：p值(rmse，ph＝30min)＝0.22，p值(cod，ph＝30min)＝0.28。
[0191]
最后，在ph＝30、60、90和120min时，通过pf将(p《0.05)δtg分别显著提高了5.0min(20％)、10.0min(20％)、10.0min(13％)和10.0min(9％)。值得注意的是，平均而言，pf允许实现tg等于ph。
[0192]
vii.另选地实施方案
[0193]
使用特定个体化非线性生理模型(诸如uva/padova t1d模拟器之类的最大模型)、特定参数估计技术(mcmc)和特定非线性预测方案(粒子滤波方法论)的示例性示例来说明结合图1的流程图在此呈现的基于个体化模型的预测方案。更一般地，基于个体化模型的预测方案可以采用多种另选的参数估计技术，并且多种另选的非线性预测方案可以与不同的个体化非线性生理模型组合使用。例如，在一些实施方案中，所采用的参数估计技术可以是最大似然法或预测误差最小化法。同样地，在一些实施方案中，所采用的非线性预测方案可以是扩展卡尔曼滤波技术或无迹卡尔曼滤波技术。
[0194]
viii.结论
[0195]
t1d中葡萄糖浓度的个性化实时预测对于实现有效的主动治疗非常重要。已提出了一种使用个性化生理模型来实现该目标的方法论。在一些实施方案中，在包括两个主要阶段的贝叶斯框架中处理该灰盒非线性模型。在第一阶段中，我们使用mcmc例如识别t1d受试者的生理模型的患者特异性参数。然后，在第二阶段中，通过粒子滤波或其他非线性预测方案使用个性化模型推断未来葡萄糖数据。第一阶段可以离线进行，回顾性地处理患者历史数据。当新的患者数据变得可用时，第二阶段可以在线运行。
[0196]
上述的方法的各种操作可通过能够执行这些操作的任何合适装置来执行，例如各种硬件和/或软件组件、电路和/或模块。通常，图中所示的任何操作都可以由能够执行这些操作的相应功能装置来执行。
[0197]
可以使用目的在于执行本文所述的功能的通用处理器、数字信号处理器(dsp)、专用集成电路(asic)、现场可编程门阵列信号(fpga)或其他可编程逻辑器件(pld)、离散门或晶体管逻辑、离散硬件组件或它们的任何组合来实施或执行结合本公开描述的各种说明性逻辑块、模块和电路。通用处理器可为微处理器，但在另选地实施方案中，该处理器可为任何市售处理器、控制器、微控制器或状态机。处理器还可实施为计算装置的组合，例如，dsp与微处理器的组合、多个微处理器的组合、与dsp核心结合的一个或多个微处理器的组合，或任何其它此类配置。
[0198]
在一个或多个方面，可以使用硬件、软件、固件或它们的任何组合来实施所述的功能。如果用软件实施，这些功能可以作为一个或多个指令或代码存储在或传输到非临时性计算机可读介质上。作为示例而非限制地，这种非暂时性计算机可读介质可包括ram、rom、eeprom、cd-rom或其它光盘存储装置、磁盘存储装置或其它磁性存储装置。
[0199]
本文公开的方法包括用于实现所述方法的一个或多个步骤或动作。在不脱离权利要求书的范围的情况下，方法步骤和/或动作可彼此互换。换句话讲，除非指定步骤或动作的特定顺序，否则可在不脱离权利要求书的范围的情况下修改特定步骤和/或动作的顺序和/或使用。
[0200]
某些方面可包括用于执行本文所呈现的操作的计算机程序产品。例如，这种计算机程序产品可包括其上存储(和/或编码)的有指令的计算机可读介质，这些指令可由一个或多个处理器执行以进行本文所述的操作。对于某些方面，计算机程序产品可包括封装材料。
[0201]
此外，应当理解的是，用于执行本文所述的方法和技术的模块和/或其他合适装置可以由用户终端和/或基站下载和/或以其他方式获得，如果适用的话。例如，这种设备可以耦合到服务器，以便于用于执行本文所述的方法的装置的转移。或者，可经由存储装置(例如，ram、rom、诸如压缩光盘(cd)或软盘等物理储存介质)来提供本文所述的各种方法，使得用户终端和/或基站可在将存储装置耦合到或提供给设备时获得各种方法。此外，可以利用用于向设备提供本文所述的方法和技术的任何其他合适的技术。
[0202]
应当理解的是，权利要求不限于上述的精准配件和组件。在不脱离权利要求书的范围的情况下，可对上述的方法及设备的布置、操作及细节作出各种修改、改变及变化。
[0203]
除非另有定义，否则所有术语(包括技术和科学术语)对于本领域普通技术人员来说具有普通和惯用的含义，并且不限于特殊或自定义的含义，除非在此明确定义。应当注意的是，当描述本公开的某些特征或方面时，特定术语的使用不应被视为暗示该术语在本文中被重新定义以被限制为包括与该术语相关联的本公开的特征或方面的任何特定特性。除非另有明确说明，本技术中使用的术语和短语及它们的变体，特别是所附权利要求中使用的术语和短语及它们的变体，应当被理解为开放性的，而不是限制性的。作为前述的示例，术语“包括”应理解为是指“包括但不限于(including,without limitation和including but not limited to)”等；本文使用的术语“包含”与“包括”、“含有”或“特征在于”同义，并且具有包容性或开放性，并且不排除额外的未列举的要素或方法步骤；术语“具有”应解释为“至少具有”；术语“包括”应解释为“包括但不限于”；术语“示例”用于提供所讨论事项的示例性实例，而不是详尽地或限制性其列表；形容词如“已知的”、“正常的”、“标准的”以及具有类似含义的术语不应被理解为将所述的事项限制在某一特定时间段或截止某一特定时间可获得的事项，而是应当被理解为涵盖现在或未来任何时间可知的或可获得的已知、正常或标准技术；并且术语如“优选地”、“优选”、“期望的”、或“可取的”以及类似含义的字词的使用不应理解为暗示某些特征对于发明的结构或功能来说是关键的、必要的或甚至是重要的，而是仅旨在强调在发明的特定实施方案中可以使用或可以不使用的另选地或附加特征。同样，除非另有明确说明，否则与连接词“和”连接的一组事项不应解读为要求这些事项中的每一个事项都存在于该组中，而是应解读为“和/或”。类似地，除非另有明确说明，否则与连接词“或”连接的一组事项不应解读为要求该组事项之间的相互排斥，而是应解读为“和/或”。
[0204]
在本说明书中使用的表示成分的量、反应条件等的所有数字应理解为在所有情况下都由术语“约”修饰。因此，除非有相反的说明，否则本文中列出的数值参数是近似值，其可以根据所求获得的期望性质而变化。最低程度上说，应当根据有效数字的数目和普通舍入法解读每一个数值参数，而不是试图将均等原则的应用限制于在要求本技术优先权的任何申请中任何权利要求的范围。
[0205]
本文所引用的所有参考文献全文以引用方式并入本文。通过引用并入的出版物和专利或专利申请与本说明书中包含的公开内容在一定程度上相矛盾，本说明书旨在取代任何此类矛盾的材料和/或获得其优先权。
[0206]
本文包含的小标题仅供参考，并且以帮助查找各个部分。这些小标题并不是为了限制针对其描述的概念的范围。此类概念可在整个说明书中具有适用性。
[0207]
此外，尽管为了清晰易懂，已经通过图示和示例的方式相当详细地描述了前述内容，但对于本领域的技术人员来说显而易见的是可以进行某些改变和修改。因此，本说明书和示例不应被理解为将发明的范围限制于本文所述的具体实施方案和示例，而是还涵盖所有落入本发明的真实范围和精神的修改和替换。

技术特征：
1.一种预测个体患者的未来血糖浓度的方法，所述方法包括：选择个体化非线性葡萄糖-胰岛素动力学生理模型，所述选择的模型具有多个模型参数，所述模型参数的值待确定；估计所述多个模型参数中的每个模型参数的值，通过将参数估计技术应用于所述个体患者的先验信息和数据以获得后验信息，所述模型参数的第一子集具有按照先验群体数据估计的值，并且所述模型参数的第二子集具有针对所述个体患者个性化的值；以及使用所述模型参数中的每个模型参数的估计值将非线性预测技术应用于所述选择的模型，以获得所述个体患者在未来某一时间的预测血糖浓度。2.根据权利要求1所述的方法，其中所述参数估计技术是马尔可夫链蒙特卡尔(mcmc)贝叶斯估计。3.根据权利要求1所述的方法，其中所述参数估计技术选自由最大后验法、最大似然法和预测误差最小化法组成的组。4.根据权利要求1所述的方法，其中所述非线性预测技术是粒子滤波。5.根据权利要求2所述的方法，其中所述非线性预测技术是粒子滤波。6.根据权利要求1所述的方法，其中所述非线性预测技术选自由扩展卡尔曼滤波技术和无迹卡尔曼滤波技术组成的组。7.根据权利要求2所述的方法，其中所述非线性预测技术选自由扩展卡尔曼滤波技术和无迹卡尔曼滤波技术组成的组。8.根据权利要求1所述的方法，其中选择的个体化非线性生理模型包括皮下胰岛素吸收子模型、口服葡萄糖吸收子模型和葡萄糖-胰岛素动力学模型。9.根据权利要求8所述的方法，其中所述选择的个体化非线性生理模型使用所述个体患者的过去血糖浓度水平、碳水化合物摄入量信息和外源性胰岛素数据作为输入。10.根据权利要求8所述的方法，其中所述皮下胰岛素吸收子模型包括第一模型参数和第二模型参数，所述第一模型参数指定胰岛素分布的wo 2022/159163a1值，所述第二模型参数指定胰岛素在第一房室中出现的延迟，所述第一和第二模型参数位于具有按照所述先验群体数据估计的值的所述模型参数的第一子集中。11.根据权利要求10所述的方法，其中所述皮下胰岛素吸收子模型包括第三模型参数和第四模型参数，所述第三模型参数指定从第一房室到第二房室的扩散速率常数，所述第四模型参数指定从第二房室到血浆的皮下胰岛素吸收速率常数，所述第三和第四模型参数位于具有针对所述个体患者个性化的值的所述模型参数的第二子集中。12.根据权利要求1所述的方法，其中所述选择的个体化非线性生理模型是最大生理模型，其中忽略了所述模型参数的选择的参数。13.根据权利要求1所述的方法，其中使用实时动态血糖监测仪(cgm)数据、碳水化合物摄入量信息和外源性胰岛素数据实时更新所述预测血糖浓度。14.根据权利要求1所述的方法，所述方法还包括用残差模型增强所述选择的个体化非线性生理模型，从而描述残差建模误差。15.根据权利要求15所述的方法，其中所述残差建模误差被建模成阶数为(5,5,1)的自回归综合移动平均(arima)。16.根据权利要求14所述的方法，其中所述残差模型具有使用预测误差法(pem)估计的
模型参数。17.根据权利要求1所述的方法，其中在多个不同时间对所述预测血糖浓度预测。

技术总结
一种预测个体患者的未来血糖浓度的方法，该方法包括：选择个体化非线性葡萄糖-胰岛素动力学生理模型，该选择的模型具有多个模型参数，该模型参数的值待确定；估计所述多个模型参数中的每个模型参数的值，通过将参数估计技术应用于所述个体患者的先验信息和数据以获得后验信息，所述模型参数的第一子集具有按照先验群体数据估计的值，并且所述模型参数的第二子集具有针对所述个体患者个性化的值；以及使用模型参数中的每个模型参数的估计值将非线性预测技术应用于选择的模型，以获得个体患者在未来某一时间的预测血糖浓度。者在未来某一时间的预测血糖浓度。者在未来某一时间的预测血糖浓度。


技术研发人员：G
受保护的技术使用者：德克斯康公司
技术研发日：2021.11.12
技术公布日：2023/8/24