梁露花,黃葉真,普思靜,鄭浩東,劉欣然,項(xiàng)榮武

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

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

1 材料與方法

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

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

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

1.4WBPBPK和集總的生理藥動(dòng)學(xué)模型的評(píng)估 對(duì)集總前茶堿的WBPBPK和利用SA集總得到的適當(dāng)集總簡(jiǎn)化的茶堿生理藥動(dòng)模型進(jìn)行評(píng)估。簡(jiǎn)化模型和原始模型之間茶堿濃度-時(shí)間曲線下面積的ARD%<0.002%和模擬動(dòng)脈濃度-時(shí)間的分布圖中原始模型和集總模型的曲線擬合程度,若二者比較良好則認(rèn)為簡(jiǎn)化模型可以較好地替代原始PBPK模型,并且可以應(yīng)用到臨床[8]。

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

2 結(jié)果

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

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

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

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

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

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

3 討論

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