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

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

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

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

1 材料與方法

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

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

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

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

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

2 結果

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

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

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

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

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

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

3 討論

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