面向智能血糖管理的餐前胰島素劑量貝葉斯學習優化方法

史大威 蔡德恒 劉 蔚 王軍政 紀立農

糖尿病是一種影響人類健康的嚴重慢性疾病,主要包括1 型(胰島素完全不分泌)和2 型(胰島素分泌相對不足).長期良好的血糖控制是避免出現糖尿病慢性并發癥(如糖尿病視網膜病變、糖尿病、腎病)的關鍵,同時也是降低心血管死亡風險的重要因素.中國擁有全世界最多的糖尿病患者,在預防和控制糖尿病方面面臨巨大挑戰[1-3].因而,研究開發適合中國糖尿病患者的先進治療技術具有重要意義.

胰島素治療是控制血糖達標的重要治療手段,其中每日多次胰島素皮下注射和胰島素泵治療是實施胰島素強化治療的兩種主要策略[4].隨著血糖監測等硬件設備的發展,人工胰腺系統成為一種新型的治療方式,它是一類閉環給藥智能系統[5-7],可全天候控制患者血糖.然而,在上述治療方法中,餐后血糖管理是關鍵步驟[8].現階段主要采用前饋控制方法來控制餐后血糖水平,即患者在就餐前提前輸注大劑量胰島素,以控制餐后血糖的波動.

目前,餐前胰島素劑量主要根據標準計算方法確定,其常用于患者的臨床治療中[9-10].具體地,標準計算方法形式如下:

式中,CHO 為攝入膳食的碳水化合物含量;Gc為餐前血糖水平;Gsp為血糖控制目標;uIOB為體內殘余活性胰島素(Insulin on board,IOB)量;CR和CF 為個體化生理參數,分別為碳水化合物系數和胰島素敏感系數.現已有不少研究聚焦于提升餐前胰島素劑量決策性能,例如文獻[11]利用血糖測量信息,通過計算餐后血糖最小值與目標值之間的差值來調整CR 和CF 系數,改善餐后血糖管理.文獻[12]提出一種基于胰島素靈敏度指標優化CR系數的方法,該指標利用患者的血糖監測和胰島素輸注數據估計得到.文獻[13]基于患者歷史血糖管理數據,利用案例推理(Case based reasoning,CBR)方法調整CR 和CF 系數.另外,也有一些研究致力于改善人工胰腺中的餐后血糖控制,例如文獻[14]結合伯格曼最小模型和無跡卡爾曼濾波自動確定餐前胰島素劑量,無需患者參與.文獻[15]結合CBR和R2R (Run-to-run)方法自適應確定餐前胰島素劑量.文獻[16]提出一種貝葉斯優化輔助學習方法,可利用患者歷史血糖監測和胰島素輸注信息調整CR系數,提升餐后血糖控制水平.

基于以上分析發現,現有方法大多通過調整CR和CF 系數來改善餐后血糖管理.這種方式不能有效地將餐后血糖動力學考慮進餐前劑量決策.隨著血糖監測設備的普及和提升,數據驅動控制方法有望為餐前劑量決策提供一種嶄新的方式.數據驅動控制無需系統模型知識,僅利用系統的輸入/輸出數據,即可完成建模和控制[17],近年來被廣泛應用于復雜的工程領域.文獻[18]采用神經網絡,提出數據驅動的赤鐵礦磨礦過程運行優化控制方法.文獻[19]針對機理建模預測和控制性能差的問題,建立基于多元時間序列的高爐鐵水硅含量數據驅動預測模型.文獻[20]提出了一種基于迭代學習模型預測控制(Iterative learning model predictive control,ILMPC)的閉環血糖控制方法,該方法利用患者歷史數據迭代更新控制目標,改善血糖控制,并通過臨床試驗[21-22]驗證了算法的有效性與安全性.與該方法所研究的實時閉環血糖控制問題不同,本文重點關注有限血糖和胰島素數據訓練樣本情形下單次餐前胰島素劑量決策問題,旨在通過單次劑量決策改善餐后血糖控制;在研究方法方面,與文獻[20-22]采用的ILMPC 方法不同,本文基于患者有限歷史數據,建立非參數化的餐后血糖高斯過程預測模型,利用風險敏感(Risk sensitive,RS)優化框架,完成考慮血糖預測不確定性的胰島素劑量魯棒優化決策,并通過引入專家經驗輔助決策,確保學習優化過程的安全性,以達到改善餐后血糖管理水平的目的.文獻[23]結合高斯過程和風險敏感控制,提出一種數據驅動餐前胰島素劑量決策方法.然而該方法采用離線建模的方式,所需樣本數量大,不適用于患者短期住院或小樣本情況下餐前劑量的決策.在小樣本情況下,數據驅動建模無法有效地學習患者血糖代謝規律,難以確保餐前劑量的安全和有效決策.

針對以上問題,本文提出臨床經驗輔助的餐前劑量貝葉斯學習優化方法,通過評估數據驅動建立預測模型的有效性,決定是否實施臨床經驗決策輔助,避免小樣本情況下貝葉斯優化決策失誤.同時,實時收集樣本,更新預測模型,形成餐前劑量自適應優化決策框架.其中考慮到中國現有糖尿病患者大多采用固定的餐前胰島素劑量[24],缺乏系統的決策規則,本文結合醫生臨床決策經驗,建立基于患者歷史血糖管理數據的系統性臨床經驗決策規則.另外,為了在小樣本情況下,快速確定合適的餐前劑量,利用迭代更新的思想,設計包含上次劑量懲罰項的代價函數,約束貝葉斯優化基于上次劑量進行求解.本文使用美國食品藥品監督管理局接受的UVA/Padova T1DM 血糖代謝仿真器[25]對所設計的方法進行仿真驗證,同時利用1 型糖尿病患者臨床治療數據,基于顧問模式分析方法[26]評估所設計方法的決策性能.實驗結果表明,本文設計的方法可為餐前胰島素劑量智能決策及臨床試驗提供方法基礎和技術支持,也為中國糖尿病患者血糖管理水平的有效改善,提供了精準醫學治療手段.

1 基于高斯過程的餐前劑量優化設計

首先,利用高斯過程建立餐后血糖預測模型;然后,利用預測的餐后血糖設計懲罰不對稱風險敏感代價函數,同時設計包含上次劑量的懲罰項;最后,形成決策餐前劑量的約束優化問題.設計方法的總體框架如圖1 所示.

圖1 餐前胰島素劑量貝葉斯學習優化方法框架圖Fig.1 Block diagram of the Bayesian learning based optimization method for meal bolus

1.1 高斯過程

高斯過程是在機器學習領域發展形成的概率化建模方法[27],與傳統建立確定的參數化模型的方法不同,其假設系統在任意的輸入x∈Rn下,輸出f(x)∈R服從聯合高斯分布,記為:

式中,m(x)和k(x,x′)分別表示均值函數和協方差函數,可根據對系統的先驗知識進行設計,其中待定的參數稱為超參數.在實際觀測中,系統的輸出可能會受到噪聲污染,假設噪聲變量服從ω～N(0,)且觀測方程為y=f(x)+ω,則觀測輸出y仍是一個高斯過程:

式中,δ(x,x′)為克羅內克函數,當且僅當x=x′時,為1;其他,為0.假定給定N個訓練輸入樣本X=[x1,x2,···,xN]T對應的輸出觀測樣本Y=[y1,y2,···,yN]T,y的后驗超參數可通過最大似然法確定[27].

確定超參數后的高斯過程可預測新的輸入x?對應的觀測值y?.由聯合高斯分布的性質可知,y?服從高斯分布y?～N(m?,),滿足:

式中,K?=[k(x?,x1),···,k(x?,xN)]T和K是協方差矩陣,其元素Kij=k(xi,xj),I為單位矩陣.可見,對于一個新的輸入,高斯過程可給出一個確定的預測值(式(4))和對該預測值不確定性的評估,以方差進行衡量(式(5)).

1.2 餐后血糖預測

現有研究大多采用具有生理學意義的房室模型描述人體血糖代謝過程[25,28].這些模型由多個耦合的微分方程組成,變量繁多且包含大量個體化參數,有時僅利用血糖測量數據難以完成參數辨識.另外,考慮到人體血糖變化的連續性,在已知過去的血糖值、胰島素輸注劑量和碳水化合物攝入量的條件下,當前血糖值必定在一定合理范圍內變化,因此在上述信息作為輸入向量條件下,可利用高斯分布去近似刻畫當前血糖的變化.基于以上分析,本文利用高斯過程對餐后血糖動力學進行建模.特別地,通過自回歸輸入/輸出信號作為模型的回歸變量,高斯過程被廣泛應用于非線性控制系統的建模[29-30].

為此,本文利用自回歸模型描述餐后血糖動力學,每一時刻t+i,i={1,2,···,n} 都由帶有加性噪聲的非線性函數表示:

式中,t表示攝入膳食的時刻,ω表示高斯白噪聲,P為血糖測量值,l表示自回歸模型的階數,u為餐前劑量,d為攝入的碳水化合物CHO 含量.為了盡可能地包含血糖態勢信息,本文考慮l=7 和采樣時間T=15 min,時間長度為2 h.另外,考慮n=8,時間長度為2 h 的餐后血糖軌跡.

高斯過程的應用分為建模和預測兩個階段.在建模階段,基于模型(6),為每種餐型構建預測模型,餐型包括早餐、午餐和晚餐.具體地,針對某一餐型,分別利用高斯過程對餐后每一時刻的血糖動力學進行建模.各時刻的建模方式類似,以t+1 時刻為例,收集zt=[Pt-7,···,Pt,u,d]T作為訓練輸入,其中血糖測量值利用最小/最大歸一化法,進行歸一化,同時為減小預測不確定性,收集差值?Pt=Pt+1-Pt作為對應的訓練輸出[29].均值函數考慮為線性形式,協方差函數采用常用的平方指數形式:

在預測階段,給定攝入的CHO 含量和注射的胰島素劑量,根據餐前血糖水平,利用上階段訓練得到的高斯過程,逐時刻預測餐后血糖水平,得到餐后2 h (8 步)的血糖變化軌跡概率分布.其中預測輸入中的餐后血糖值采用血糖預測值.由于建模階段,訓練輸出采用差值的形式,因此在得到原始的預測值后,需要將此值與上一時刻的血糖測量(預測)值相加,確定為當前時刻的血糖預測值,預測不確定性(方差)保持不變.考慮到糖尿病患者一般具有規律的飲食習慣(如相似的用餐時間和膳食攝入量),因此,可將膳食攝入對餐后血糖變化的影響,近似成為一個常值擾動,利用高斯過程魯棒性建模的特點,建立無需膳食信息的餐后血糖預測模型,進一步減輕患者估算CHO 攝入量的負擔,實現無需膳食信息的劑量優化決策.

1.3 不對稱風險敏感代價函數

血糖控制中存在高低血糖不對稱風險;與長時間高血糖所引發的并發癥不同,即時的嚴重低血糖會直接造成患者死亡.因此在決策劑量糾正患者餐后高血糖的同時,需要格外小心低血糖事件的發生.現有研究中一個可行的解決方案是構建不對稱代價函數[31-32].受文獻[32]啟發,本文對高于和低于血糖目標值的偏離量實施不對稱懲罰;此外,考慮到預測信息存在的不確定性,本文以風險敏感的形式來設計代價函數.記高斯過程提供的8 步血糖預測量分別為Gi～N(mi,),i∈ {1,2,···, 8},并組成一個向量狀態G=[G1,G2,···,G8]T,該狀態描述了餐后2 h 血糖變化軌跡的概率分布.進一步地,根據風險敏感代價函數的定義[33],設計不對稱風險敏感代價函數如下:

式中,E 為數學期望;Gr為餐后血糖控制目標;Q為懲罰矩陣,其為對角形式,用于對預測的高低血糖值施加不對稱懲罰權值;γ <0 表示風險敏感系數,決定著優化者對預測不確定性的態度[34-35],γ <0使得優化者傾向于選擇減少預測不確定性的控制量.基于上述不對稱風險敏感代價函數,優化者能夠充分利用數據樣本中的決策經驗,同時形成自我決策能力.

對于懲罰矩陣Q,記Q的第i個對角元素為Qi,設計不對稱懲罰如下:

式中,1(·)為示性函數,Gri為Gr的第i個元素,0為基準懲罰系數,Γ:={α,β,c1,c2} 為決定懲罰強度的元素集合.記h(Gi,Γ)=c1/(1+exp{α(β-(Gri-Gi))}+c2),可得在Gi-Gri <0 時,h(Gi,Γ)的變化曲線如圖2 所示.由圖2 可知,懲罰強度的上/下界分別由c1+c2和c2決定,中間值由β決定,增長速度由α控制.總之,在高于血糖目標值時,設計常值的懲罰系數可使得優化者對高血糖保持相對合理的響應;在低于血糖目標值時,設計隨偏離量增大而指數式增大的懲罰系數可使得優化者在偏離血糖目標值不多時保持保守的響應,但同時能及時快速地對大偏離量做出響應,保證安全性.

圖2 懲罰強度變化示意圖Fig.2 Schematic diagram of the penalty changes

根據上述設計原則,本文利用美國食品藥品監督管理局接受的血糖代謝仿真器確定參數集合 Γ 的數值.本文方法的參數確定選用仿真器中的“平均成人患者”完成,性能驗證利用10 個虛擬成人患者完成.在參數確定過程中,首先選取平均成人患者進行8 h 的仿真,模擬仿真從9:00 開始,假定患者在12:00 進食70 g 碳水化合物,觀察不同餐前胰島素劑量下患者餐后血糖的變化,得到發生低血糖事件與餐后血糖偏離控制目標Gr程度之間的大致關系.基于此關系,首先根據圖2 選取確定合適的β和c2,使得預測血糖在控制目標附近時,懲罰程度大致不變;然后,設計α和c1確定懲罰增長速度和最大懲罰強度,使得預測血糖偏離控制目標的程度超過β時,懲罰急速增大,對預測低血糖高敏感性;最后,利用所得的參數組進行不同含量碳水化合物進食測試,確保參數組安全和有效.

1.4 優化問題

基于上述不對稱風險敏感代價函數LA,本文將餐前劑量決策設計成帶約束的隨機優化問題,具體形式如下:

式中,u為待優化的餐前劑量,up為上一次相同餐型餐前輸注的劑量,R為權重系數.設計懲罰項R(u-up)2引入迭代更新的方式,使得在小樣本情況下,優化者能基于上次劑量,快速確定合適的劑量.其中R可用于調整對上次劑量的信任程度,可使得優化者偏向于選擇上次劑量附近的劑量.

2 臨床經驗輔助的貝葉斯優化決策

由式(11)可知,隨機變量的性質和不對稱懲罰矩陣的設計使得難以計算該代價函數對待優化劑量u的梯度.因而,本文結合貝葉斯優化[36]和蒙特卡洛法求解上述隨機優化問題.貝葉斯優化首先利用高斯過程擬合代價函數(11),然后利用采集函數確定下次評估的劑量,同時利用蒙特卡洛法計算對應的代價函數觀測值,收集觀測樣本更新高斯過程,重復上述操作,最后確定優化問題的解.另外,考慮到在小樣本情況下,數據驅動建模利用少量的輸入/輸出數據,無法準確學習患者餐后血糖代謝規律,本文設計模型有效性在線評估機制,評估血糖預測模型的有效性,進而實施臨床經驗決策輔助,避免小樣本情況下貝葉斯優化決策失誤;同時實時收集樣本,更新預測模型,形成自適應優化決策.

2.1 貝葉斯優化

2.1.1 代價函數建模

首先,利用高斯過程擬合式(11)中決策變量u與代價函數L(u)之間的映射關系.均值函數考慮為常數形式,協方差函數采用平方指數形式.給定N個代價函數的觀測樣本,D1:N={u1:N,L(u1:N)},通過最大化后驗概率確定高斯過程的超參數.接著,利用預測式(4)、式(5),預測新劑量u?對應的代價函數值,該值將被用于構建采集函數(見第2.1.2 節)來確定下一次評估的劑量.在確定評估的劑量后,利用蒙特卡洛法估算對應的代價函數觀測值.具體地,首先,給定一劑量,根據第1.2 節方法計算餐后血糖變化軌跡的概率分布;然后,獨立生成1000 條餐后血糖軌跡樣本;最后,計算這些樣本對應的代價函數的均值,該均值即為該劑量下代價函數的觀測值.

2.1.2 采集函數

采集函數是貝葉斯優化中最重要的一個環節,其直接決定優化方向,影響求解速度.具體地,基于建模階段提供的預測信息,采集函數通過權衡探索新的搜索空間和利用現有的可行區域確定下一次評估點.目前,相關研究已提出了多種采集函數[36],本文采用改善期望(Expected improvement,EI)采集函數[37]來確定每次的評估劑量.

為最小化代價函數(11),改善期望采集函數αEI(u?),設計如下:

式中,Φ (·)為標準正態分布的累積分布函數,?(·)為標準正態分布的概率密度函數.最大化αEI,即可確定下一次評估劑量u?的值.

式中,約束條件由設計的不對稱風險敏感代價函數引入,使得采集函數能實時利用更新的血糖管理樣本,基于上次劑量up快速確定下一次評估的劑量.

貝葉斯優化的求解步驟見算法1,在進行M次順序優化求解后,得到優化問題(11)～(14)的解,記為ub.為保證安全性,利用體內殘余活性胰島素約束修正該劑量,IOB 是指上次給藥一段時間后在體內仍然殘留的有降糖作用的胰島素,用來約束當前胰島素輸注,以防止過量注射胰島素導致低血糖,其計算方法可參見文獻[31].記IOB 約束為uIOB,則最終確定的餐前胰島素劑量為:

算法1.餐前劑量貝葉斯優化決策算法

2.2 預測模型有效性評估

餐后血糖預測模型影響著貝葉斯優化的決策方向.考慮到只有符合生理代謝規律的預測模型才能保證貝葉斯優化能根據預測的高低血糖正確地決策劑量,因而本文依據胰島素具有降糖作用這一生理代謝規律設計模型有效性在線評估機制,防止貝葉斯優化決策失誤.具體地,該機制利用訓練得到的高斯過程,根據此次餐前血糖監測水平,基于上次相同餐型的餐前劑量,分別預測增加 ?u、保持不變和減少 ?u劑量情況下的餐后血糖軌跡,并依次記為,其都為由8步血糖預測值組成的向量,進而設計如下的判斷條件,比較確定模型的有效性:

式中,J為判斷結果,若為1,則認為預測模型有效,采用貝葉斯優化決策推薦的劑量;若為0,則認為預測模型無效,采用臨床經驗輔助決策推薦的劑量(見第2.3 節),避免貝葉斯優化決策失誤.“<” 表示一一比較兩向量相同位置上元素的大小,若都滿足小于條件,則返回邏輯1;否則,返回邏輯0.

2.3 臨床經驗輔助決策

臨床經驗輔助決策可約束在預測模型不可信時的貝葉斯優化,推薦合適的餐前劑量,保證安全性.根據發布的中國1 型糖尿病胰島素治療指南[24],國內患者普遍規律飲食,缺乏估算膳食CHO 含量的經驗,通常采用固定的餐前劑量,之后根據血糖監測結果進行調整.一般地,餐前劑量基于全天胰島素總量的40%～60%,根據早、中、晚三餐按1/3、1/3、1/3 或1/5、2/5、2/5 分配.這種方式不能有效地利用餐后血糖管理歷史數據,同時沒有形成系統的劑量實時調整方法,容易受到醫生或患者主觀意向的影響.鑒于此,本文將醫生臨床決策經驗量化,形成系統的臨床經驗決策規則.具體地,利用歷史數據,分3 檔評估上次相同餐型餐后血糖的控制效果,再結合考慮本次餐前血糖態勢(包括血糖水平和變化趨勢),基于上次劑量調整確定本次劑量,所設計的臨床經驗決策規則見圖3.

圖3 臨床經驗決策規則Fig.3 Clinical experience based decision rules

圖3 中Gh和Gl分別表示評估餐后血糖控制水平的高低閾值,表示上一次相同餐型餐后到下一餐餐前血糖監測的最小值,表示此次餐前半小時的血糖監測平均值,表示上一次相同餐型餐前半小時的血糖監測平均值,?G表示血糖均值的增量,上升、平穩和下降表示此次餐前血糖水平的變化趨勢.臨床經驗輔助決策,首先通過查詢上次相同餐型餐后血糖控制效果和此次餐前血糖態勢符合的決策規則,推薦圖中對應的劑量ub;然后,同樣利用IOB 約束,修正此次劑量.

3 結果

本文方法參數利用美國食品藥品監督管理局接受的UVA/Padova T1DM 血糖代謝仿真器調試得到,同時根據臨床經驗,考慮其是否符合人體血糖代謝實際情況,最后加以確認,具體數值見表1.進一步地,設計不同仿真測試方案,利用仿真器中的成年患者組(內含10 個模擬成年人)驗證本文方法的性能.另外,收集1 型糖尿病患者的臨床治療數據,利用顧問模式分析方法[26],評估本文方法的決策性能.

表1 方法參數Table 1 The parameters of the proposed method

基于以下場景對本文方法進行仿真.模擬仿真假定患者采用每日多次胰島素皮下注射治療方案,即注射基礎胰島素和餐前大劑量胰島素.仿真從第1 天的5:00 開始,持續5 天,共120 h.考慮到糖尿病患者通常具有規律的日常活動(如類似的用餐時間和CHO 攝入量),但為了模擬生活方式存在的隨機性,假設患者早餐、午餐、晚餐的CHO 攝入量服從正態分布(均值和標準差分別為[50,75,75] g 和[3,4,4] g);用餐時間分別在[7:00,9:00]、[11:00,13:00]、[18:00,20:00]之間,服從均勻分布.假定5天的胰島素基礎率保持正常值不變,前3 天的餐前胰島素劑量由標準計算方法調整得到:

式中,Gsp通常選為140 mg/dL;r為劑量調整系數,可用于仿真前期不同的劑量決策情況;后2 天的餐前胰島素劑量利用本文方法確定.

對于仿真器而言,虛擬患者最優的CR 和CF參數是已知的,標準計算方法可在仿真器中決策近似最優的餐前劑量,實現良好的餐后血糖控制.因此,為仿真現實生活中存在的劑量決策情況(偏多、偏少和正常),本文將標準計算方法決策的劑量乘以不同調整系數 (r=1.5、r=0.5 和r=1)確定為前期決策的劑量,生成本文方法前期所需的訓練樣本,后2 天的餐前劑量利用本文方法確定,并與保持調整系數不變的前期決策方法對比,以此說明本文方法的學習與改善血糖控制的能力.特別地,當本文方法能實現與標準計算方法 (r=1)類似的控制性能時,可在一定程度上說明,利用本文方法可實現無需CR 和CF 參數的餐前胰島素劑量優化決策.在下文圖表中,本文提出的方法用“本文方法”表示,不同調整系數的標準計算方法用“標準方法”表示,同時在表頭和圖注標明調整系數.

由于仿真實驗中,每天的用餐時間和CHO 攝入量都類似,本文方法分別收集前3 天的早餐、午餐和晚餐樣本,各3 個,為每個餐型構建無需膳食信息的餐后血糖預測模型,并在后2 天依次收集樣本進行更新.仿真結果主要使用正常血糖70～180 mg/dL 時間比率、低血糖低于70 mg/dL 時間比率、嚴重低血糖低于54 mg/dL 時間比率、高血糖高于250 mg/dL 時間比率和血糖濃度平均值大小等,如表2、表3 所示的7 種評價指標,對方法性能進行評估.表中數據以中位數(四分位距)的形式展現,p值由仿真的患者數量計算得到,p<0.05 被認為具有統計學差異性,并在表中以粗體標出.同時,餐后血糖調節及劑量決策比較結果見圖4～圖6,圖中曲線包括5%、25%、75%和95%的四分位數曲線和中位數曲線,相關討論見第3.1.1 節.

表2 前期餐前劑量正常情形下 r=1 血糖控制評估結果 (%)Table 2 Evaluation results of the glucose control in the case of normal boluses in the early stage when r=1 (%)

表3 前期餐前劑量偏少/偏多情形下血糖控制評估結果(%)Table 3 Evaluation results of the glucose control in the case of under-estimated/over-estimated boluses in the early stage (%)

圖4 前期劑量正常 r=1 情形下,餐后血糖調節和劑量決策比較結果Fig.4 Comparison of postprandial glucose regulation and dosage decision in the case of normal boluses in the early stage when r=1

為進一步說明本文方法的學習能力和魯棒決策能力,本文利用仿真測試實際生活場景下本文方法的血糖控制性能,并與R2R 方法[38]進行比較.具體地,模擬仿真從第1 天的5:00 開始,持續5 天,共120 h.考慮到患者在實際住院調整餐前胰島素劑量過程中,每天要求攝入固定碳水化合物含量的膳食,餐前劑量由固定的碳水化合物含量確定,但患者執行始終存在偏差(少食或過食).因此為模擬實際生活場景,假設患者早餐、午餐和晚餐的CHO 攝入量分別在標準量[50,75,75] g 附近隨機變動,變動范圍為[-15,15] g;用餐時間分別在[7:00,9:00]、[11:00,13:00]、[18:00,20:00],服從均勻分布.餐前劑量由標準計算方法按CHO 攝入標準量確定,但三餐分別有著不適當的CR 參數,CF 參數始終保持正常不變.在后續4 天,利用R2R 方法優化三餐的CR 參數,其中R2R 方法形式如下:

式中,k表示第k次R2R 優化,下標i∈{1,2,3} 分別對應早餐、午餐和晚餐,表示CR 初始值的倒數,k1、k2和k3表示R2R 增益,Tb(k)、Ta(k)和Gm為餐后血糖控制評價指標,分別表示低血糖低于70 mg/dL 時間比率,高血糖高于180 mg/dL時間比率和血糖濃度平均值,其統計時長最長為5 h,或從本次餐后到下次餐前.式中參數數值與文獻[38]一致.為確保比較的公平性,本文方法收集R2R 方法前3 天的數據進行建模,并在后續兩天進行決策與R2R 方法比較,比較結果如表4 和圖7 所示,相關討論見第3.1.2 節.

表4 膳食擾動下血糖控制評估結果(%)Table 4 Evaluation results of the glucose control in the case of meal disturbance (%)

顧問模式分析方法可將臨床試驗中獲得的血糖監測數據重新提供給本文方法,使得本文方法在對應的用餐時刻,根據餐前血糖監測數據,確定相應的胰島素劑量,并可與同時刻醫生建議的胰島素劑量進行比較,評估方法的決策性能.本文方法決策的胰島素劑量只是用于比較,并不會對歷史血糖監測數據造成因果性影響.評估結果見圖8、圖9,相關討論見第3.2 節.

3.1 仿真驗證

3.1.1 前期不同劑量決策情況

1)前期餐前劑量正常情形

取劑量調整系數r=1,仿真前期餐前劑量正常的情形.將后2 天的血糖控制效果進行比較,統計結果如表2 所示.后3 天的餐后血糖調節及劑量決策比較結果如圖4 所示.由表2 和圖4 可看出,在前期訓練樣本良好時,本文方法可在無需膳食信息下實現與標準方法類似的血糖控制效果,表現在正常血糖70～180 mg/dL 時間比率(92.3% 對93.6%,p=0.188)和血糖平均值(132.5 對133.2,p=1.000)等評價指標,兩方法均未出現低血糖情況(0 對0,p=1.000).

2)前期餐前劑量偏低情形

取劑量調整系數r=0.5,仿真前期餐前劑量偏低的情形.該情形用于測試本文方法在前期樣本餐前劑量偏低時,是否可快速地確定合適的餐前劑量,安全有效地改善患者出現的餐后高血糖情況.將仿真的后2 天血糖控制效果進行比較,比較結果如表3 和圖5 所示.可以看出,在前期劑量決策偏少時,與繼續采用劑量偏少的決策方法相比,在后2天的決策中,本文方法可基于前1 天的劑量,快速確定合適的劑量,有效緩解繼續高血糖情況,這可以從正常血糖70～180 mg/dL 時間比率(77.1%對64.8%,p=0.004)、高血糖高于250 mg/dL 時間比率(0% 對1%,p=0.031)和血糖平均值(148.7對162.1,p=0.004)等評價指標看出.

圖5 前期餐前劑量偏少 r=0.5 情形下,餐后血糖調節和劑量決策比較結果Fig.5 Comparison of postprandial glucose regulation and dosage decision in the case of under-estimated boluses in the early stage whenr=0.5

3)前期餐前劑量偏高情形

取劑量調整系數r=1.5,仿真前期餐前胰島素劑量偏高的情形.該情形用于測試本文方法在前期樣本餐前劑量偏高時,是否可以快速地確定合適的餐前劑量,安全有效地改善患者出現的餐后低血糖情況.將仿真的后2 天血糖控制效果進行比較,比較結果如表3 和圖6 所示.可以看出,在前期劑量決策偏多時,與上一種情況相似,在后2 天的決策中,本文方法可基于前一天的劑量,快速確定合適的劑量,有效緩解繼續采用劑量偏多的決策方式所帶來的低血糖情況,體現在低血糖低于54 mg/dL時間比率(0%對1.6%,p=0.031)、低于70 mg/dL時間比率(0% 對6.5%,p=0.008)和血糖平均值(125.6 對114.9,p=0.002)等評價指標.

圖6 前期餐前劑量偏多 r=1.5 情形下,餐后血糖調節和劑量決策比較結果Fig.6 Comparison of postprandial glucose regulation and dosage decision in the case of over-estimated boluses in the early stage when r=1.5

3.1.2 膳食擾動場景

由表4 和圖7 可看出,R2R 方法容易受到膳食實際攝入量與標準量不匹配的影響,決策過多或過少的餐前劑量,導致餐后低血糖事件或持續高血糖事件的發生.由于患者每天攝入類似CHO 含量的膳食,本文方法利用歷史數據建立包含預測不確定性的高斯過程預測模型,有效降低了膳食信息對決策性能的影響,從而在改善血糖控制性能的同時,對膳食擾動具有更好的魯棒性.這可以從低血糖情況的改善(0%對3.6%,p=0.016)和血糖平均值的提升(131.3 對123.8,p=0.006)等評價指標的改善看出.

圖7 膳食擾動下,本文方法與R2R 方法餐后血糖調節和劑量決策比較結果Fig.7 Comparison of postprandial glucose regulation and dosage decision in the case of meal disturbance between the proposed method and the R2R method

3.2 臨床數據驗證

本文收集在北京大學人民醫院采用每日多次胰島素皮下注射治療方案的某位1 型糖尿病患者臨床數據,并利用顧問模式評估方法性能.在患者治療過程中,餐前胰島素劑量仍由醫生決策確定,本文方法決策的劑量只用于后續的比較分析,并未實際應用于患者治療中.患者住院期間采用瞬感血糖監測技術(Flash glucose monitoring,FGM),每15 min監測一次皮下血糖濃度.該數據收集通過北京大學人民醫院倫理審查委員會的批準,并獲得所有參與者的書面知情同意[39].由于患者住院時間普遍在一周左右且在住院期間飲食規律,因此利用前3 天的血糖監測數據和醫生餐前決策劑量訓練得到無需膳食信息的餐后血糖預測模型,并利用后4 天臨床數據驗證系統性能.患者前3 天的歷史血糖管理數據如圖8 所示,后4 天的顧問模式分析結果如圖9 所示,圖中醫生決策的劑量以“醫生建議”表示,本文方法決策的劑量中由臨床經驗決策規則決策的以“決策規則”表示,由貝葉斯優化決策的以“貝葉斯優化”表示.

圖8 前期血糖管理臨床數據Fig.8 The clinical data of the glucose management in the early stage

圖9 基于臨床數據的餐前劑量決策評估結果Fig.9 The evaluation results for the meal bolus decision based on the clinical data

由圖8 和圖9 可以看出,本文方法可基于上次相同餐型的劑量,快速確定合適的餐前劑量.其中臨床經驗輔助決策考慮到上次相同餐型的餐后血糖控制水平高于高閾值Gh且本次餐前血糖水平較高,依據邏輯規則推薦高于醫生建議的劑量,有效地利用患者血糖歷史管理數據進行劑量決策.由圖9可以看出,在臨床經驗輔助決策下,本文方法不斷收集樣本實現貝葉斯學習優化,持續學習患者血糖代謝規律,使得貝葉斯優化決策模塊可基于預測的餐后高血糖,優化確定高于醫生建議的劑量.這些高于醫生建議的劑量,有望緩解患者餐后的高血糖情況.

4 結束語

本文設計一種臨床經驗輔助的餐前胰島素劑量貝葉斯學習優化方法.該方法通過評估高斯過程預測模型是否符合血糖代謝規律,決定是否實施臨床經驗輔助決策,避免小樣本情況下貝葉斯優化決策的失誤.同時,不斷收集樣本,持續更新預測模型,實現自適應優化決策.仿真和臨床數據驗證結果表明,本文所設計的臨床經驗決策規則,可有效地利用患者血糖歷史管理數據輔助決策,使得貝葉斯優化決策可持續有效地學習患者血糖代謝規律,根據預測的高/低血糖水平,基于上次相同餐型的劑量,快速確定合適的劑量.本文方法輔助醫生決策的臨床試驗申請已通過北京大學人民醫院倫理審查(批件號:2020PH338-01),未來將通過隨機對照臨床試驗,進一步評估方法的安全性和有效性.