基于模擬退火算法的適當集總法簡化茶堿生理藥動學模型*

梁露花,黃葉真,普思靜,鄭浩東,劉欣然,項榮武

(1.沈陽藥科大學醫(yī)療器械學院生物醫(yī)藥信息教研室,沈陽 110016;2.遼寧省醫(yī)藥大數(shù)據(jù)與人工智能工程技術(shù)研究中心,沈陽 110016;3.沈陽藥科大學制藥工程學院,沈陽 110016;4.沈陽藥科大學藥學院,沈陽 110016)

全身生理藥動學模型 (whole-body pharmacokinetics model,WBPBPK)是高維復雜的,通常是非線性的,表示為一個常微分方程(ordinary differential equation,ODE)系統(tǒng),其中每個常微分方程表示一個組織或器官腔室的狀態(tài)[1]。在臨床研究中,通常希望建立的WBPBPK能吻合動物或臨床數(shù)據(jù),進而改善藥物相關的特異性參數(shù),但是,由于倫理限制,通常不能直接測定組織濃度,許多參數(shù)(例如組織-血液分配系數(shù))不能僅僅通過血漿數(shù)據(jù)得到[2]。因此,根據(jù)模型參數(shù)先驗的信息,在高維系統(tǒng)中,統(tǒng)計模型可能無法被識別,將WBPBPK擬合到臨床數(shù)據(jù)通常會產(chǎn)生數(shù)值上不穩(wěn)定的分析。故一個常見的策略是,將全身生理藥動學系統(tǒng)進行降階,把高維復雜的模型簡化為更少的狀態(tài),從而減少生理藥動學模型中ODE和參數(shù)的數(shù)量,實現(xiàn)對WBPBPK的簡化,以穩(wěn)定臨床數(shù)據(jù)[3]。目前常用簡化生理藥動學模型的4種方法來搜索最優(yōu)集總模型,包括完全枚舉法(參考方法)、非自適應隨機搜索(non-adaptive random search,NARS)、碎石圖加非自適應隨機搜索(scree plot plus NARS)、模擬退火算法(simulated annealing algorithm,SA)[4]。筆者通過已有的茶堿成人WBPBPK,運用SA對全身模型進行適當集總簡化。對簡化模型和原始模型預測的血藥濃度-時間分布進行比較,以保留對系統(tǒng)的機制解釋。

1 材料與方法

1.1數(shù)據(jù)來源 模型數(shù)據(jù)來源于KRAUSS等[5]提供的非線性混合效應模型用戶指南的第五部分12例30歲的健康志愿者的身高、體質(zhì)量等基本參數(shù),將這些基本參數(shù)輸入到PK-sim數(shù)據(jù)庫,得到志愿者的系統(tǒng)特定的解剖學參數(shù)。在開放系統(tǒng)藥理學(www.open-systems-pharmacology.org)中找到茶堿的吸收、分布、代謝和排泄(absorption distribution, metabolism, and excretion,ADME)參數(shù),放入PK-sim數(shù)據(jù)庫(版本10.0,是Open Systems藥理學套件的一部分)中查看,得到藥物特定的參數(shù)。按照12例志愿者不同的體質(zhì)量對其分別給予不同劑量,為了計算機模型的簡單性,文獻[5]給出計算其給藥劑量的平均值(單位:mg),這并不會影響最終集總模型的通用性。

1.2茶堿的原始模型 從BJ?RKMAN 等[6]已發(fā)表的文獻中確定茶堿的WBPBPK。茶堿的PBPK模型在MATLAB?(版本 R2016)中被編碼為常微分方程,使用MATLAB矩陣指數(shù)求解其常微分方程。茶堿的WBPBPK是由14個隔室組成的PBPK模型[7]。為了方便將茶堿的原始PBPK模型轉(zhuǎn)化為矩陣的形式,首先就是列出原始微速率常數(shù)矩陣K。把每個隔室進行單獨的編號,其中靜脈血為輸出狀態(tài),做為非集總狀態(tài),不對其進行編號,剩余的13個隔室對其進行編號,肺為1號、動脈血為2號、大腦為3號、心臟為4號、腎臟為5號、肝臟為6號、肌肉為7號、皮膚為8號、脂肪為9號、酮體為10號、腸道(小腸)為11號、脾臟和胰臟合并為12號、胃為13號。

1.3模擬退火算法適當集總以簡化模型 茶堿的WBPBPK為14個隔室的高維復雜模型,為了保證原始系統(tǒng)動力學不會改變,選用一種對于線性或線性化模型來說是一種潛在的更靈活的技術(shù)——集總,其中狀態(tài)的合并是根據(jù)模型的性能特征來確定的,集總技術(shù)參考文獻[7]。在本文中SA被用來在尋找集總矩陣時最小化相對差異的絕對值(absolute value of the relative difference,ARD%),使用SA對M矩陣進行最優(yōu)搜索[8]。SA的變量設置為:初始溫度為104℃,負責在搜索過程中監(jiān)管接受新的能量狀態(tài)的概率;溫度下降速率為每周期0.999,即第一次循環(huán)之后溫度從104℃下降到104×0.999 ℃,線性冷卻方案繼續(xù)每個循環(huán);最大迭代次數(shù)為2×104。直到尋找到最優(yōu)M矩陣停止。

1.4WBPBPK和集總的生理藥動學模型的評估 對集總前茶堿的WBPBPK和利用SA集總得到的適當集總簡化的茶堿生理藥動模型進行評估。簡化模型和原始模型之間茶堿濃度-時間曲線下面積的ARD%<0.002%和模擬動脈濃度-時間的分布圖中原始模型和集總模型的曲線擬合程度,若二者比較良好則認為簡化模型可以較好地替代原始PBPK模型,并且可以應用到臨床[8]。

1.5軟件 PK-sim(版本10.0,是Open Systems藥理學套件的一部分,www.open-systems-pharmacology.org)用于查找系統(tǒng)特定解剖學參數(shù)和藥物的特異性參數(shù)。MATLAB?(版本 R2016)用于SA進行茶堿PBPK模型的簡化,尋找最優(yōu)集總矩陣M和進行生理藥動學建模。

2 結(jié)果

2.1SA適當集總算法簡化的生理藥動學模型 將SA適當集總算法代碼,分別保存為m文件,輸入12例志愿者生理解剖學參數(shù)和茶堿的相關參數(shù)進行運行,最終得到最優(yōu)集總M矩陣為集總成7態(tài)的全身藥動學模型,最終簡化模型的集總矩陣具體如下:

將最優(yōu)集總M矩陣轉(zhuǎn)化為具體的隔室,得到茶堿集總后的簡化PBPK模型集總成7態(tài)時,大多數(shù)集總狀態(tài)下的預測與原始系統(tǒng)預測相當一致,動脈血液單獨為一態(tài);第2個偽狀態(tài)是肌肉、酮體、腸道;大腦、合并的脾臟和胰臟、胃合并為一個偽狀態(tài);腎臟程單獨的一態(tài);肝臟為單獨的一態(tài);肺為單獨的一態(tài);第7個偽狀態(tài)是心臟、皮膚、脂肪合并為一態(tài)。第二偽狀態(tài)是肌肉、酮體、腸道的集總,第三偽狀態(tài)是大腦和內(nèi)臟器官的集總。然而,預測這2個集總隔室中的藥物量對茶堿沒有特別的臨床意義。

2.2適當集總簡化PBPK模型的評估 運用MATLAB進行SA搜索集總得到一個7態(tài)的M最優(yōu)集總矩陣,此矩陣結(jié)果輸出的血藥濃度-時間曲線下面積為332.47 ng·h·mL-1,ARD%為0.000 190 82%。在對茶堿PBPK模型簡化之前人為規(guī)定集總模型和原始模型之間茶堿濃度-時間曲線下面積的最大差異為0.002%,即ARD%<0.002%,表示原始模型和集總模型血藥濃度-時間偏差不大,此時認為簡化后的茶堿模型有意義。額外運行不同初始溫度(即分別為103和105℃)在集總數(shù)量方面給出了相同的結(jié)果狀態(tài)。圖1溫度變化趨勢和經(jīng)過迭代后對最優(yōu)ARD%解的收斂顯示了模擬退火算法對數(shù)尺度的溫度變化的趨勢,和經(jīng)過迭代后對最優(yōu)ARD%解的收斂。隨著溫度指數(shù)的下降即退火,在連續(xù)迭代中形成ARD%的卷積曲面,在經(jīng)過20 000次數(shù)的迭代后,ARD%達到了穩(wěn)定,表明迭代收斂到最佳狀態(tài)。上坡跳躍隨著溫度的降低而減小,最終使搜索收斂到最優(yōu)解。

圖1 溫度變化趨勢(左)和經(jīng)過迭代后對最優(yōu)ARD%解的收斂(右)

2.3WBPBPK和集總的生理藥動學模型的評估 在MATLAB中使用模擬退火算法得到最優(yōu)M矩陣之后,繼續(xù)在MATLAB中調(diào)用茶堿的集總模型參數(shù),模擬動脈血藥濃度-時間分布圖(圖2),并與原始的茶堿PBPK模型的模擬動脈血藥濃度-時間分布圖進行比較,可以看出集總后隔室與原始系統(tǒng)的偏差幾乎不會影響血漿中藥物濃度的時間過程,認為集總模型與原始模型相關生理意義繼續(xù)保持,此次集總的茶堿簡化成人模型有生理學意義。

圖2 原始模型和7態(tài)集總模型模擬動脈濃度-時間分布

3 討論

本研究通過查閱茶堿的相關文獻得到關于該藥物的相關參數(shù)、臨床實驗數(shù)據(jù)以及茶堿的WBPBPK模型,運用SA適當集總算法對PBPK模型進行簡化,對簡化后的模型進行評估。盡管簡化后模型與原始模型比較模擬濃度-時間曲線偏差較小。但研究中仍然存在局限性,本研究所使用的參數(shù)為12例30歲男性歐洲人生理參數(shù),只考慮了肝、腎的清除率忽略了肺、皮膚等其他組織器官的清除率。總體看來,本文利用SA適當集總算法對茶堿的成人模型的簡化相對有意義,經(jīng)評估與原始模型預測血藥濃度的差異性較小,可以用于表征潛在生理和ADME相關參數(shù)以及藥物特定參數(shù)的個體間變異性。后續(xù)還需要更多的藥物例子來證明這種簡化方法在藥物開發(fā)中的優(yōu)勢。