微循環血管不斷地調節自身結構，以適應組織的功能性需求。微循環結構自調節模型能夠仿真這一過程，為生理研究提供輔助，但目前缺少合適的模型參數設置方法，限制了模型的進一步應用。本文提出一種改進的量子粒子群優化算法用于設置模型參數，并在真實的大鼠腸系膜微循環血管網絡上進行仿真實驗。仿真結果表明，該方法的參數優化能力優于標準粒子群算法、標準量子粒子群算法和相關文獻報道的 Downhill 算法，可使微循環結構自調節模型的仿真更接近動物實驗數據，并顯著提高模型的有效性。
引用本文： 潘清, 姚家良, 王若帆, 曹平, 寧鋼民, 方路平. 基于改進量子粒子群算法的微循環結構自調節模型參數優化研究. 生物醫學工程學雜志, 2017, 34(5): 784-789. doi: 10.7507/1001-5515.201607030 復制
引言
在人體循環系統中，微循環為組織提供氧氣和營養物質，并帶走新陳代謝產物，對于人體正常的生理功能起到重要的作用。微循環的結構性自調節是指血管系統響應組織和器官的功能需求調節血管的舒張與收縮。它確保血液充分灌注器官，使其正常運作，是微循環完成其生理功能的主要機制之一[1]。微循環結構自調節功能異常與高血壓、冠心病等疾病的發生與發展密切相關[2-4]。
生理實驗和臨床試驗是研究循環系統生理病理特性的傳統方法。然而，由于微循環的空間尺度較小，傳統方法難以定量、持續地觀察其生理狀態，造成了微循環系統研究的困難。數學模型與仿真技術能夠有效輔助微循環系統的生理功能研究。一些模型重點關注真實微血管網絡中的血流動力學特性[5-6]，另一些模型則研究簡單血管網絡響應各類外部刺激而發生結構改變的過程[7-8]，但鮮有考慮真實微血管網絡中血管結構自調節過程的數學模型。Pries 等[9]提出了一個描述真實微循環網絡中血管結構在血流動力學與新陳代謝因素刺激下進行自我調節的數學模型。該模型研究了血管剪切力、血管透壁壓、新陳代謝信號、傳導信號四個刺激因素對血管管徑調節的共同作用。基于此模型，Pries 等[10]發現了腫瘤組織中動靜脈短路病理性擴張的一種信號傳遞機制，Vorderwülbecke 等[11]研究了血流剪切力介導的內皮間聯蛋白生成對于微循環灌注的作用，這些研究表明該模型能夠比較全面地反映微循環網絡中的結構自調節過程，從而輔助生理、病理研究。
該模型中，四類刺激信號的比例關系等參數難以基于動物實驗數據確定。Pries 團隊[9-11]采用了 Do-wnhill Simplex 最優化方法確定模型參數。該方法盡管合理地設置了參數，但存在著初值依靠經驗、優化限于局部等限制，因而難以推廣至其它微循環網絡的研究，并且無法確定參數設定的最優性。為克服上述局限性，同時考慮所研究的參數優化問題具有高維度、非線性、不連續等特點，本研究采用具有量子行為的粒子群優化（particle swarm opti-mization，PSO）算法[12]。在其基礎上，本研究改進了其壓縮-擴展系數設置，加入了粒子選擇算子，使算法具有更強的全局優化能力，以期獲得全局最優的模型參數，提高模型的可靠性、有效性和適用性。
1 微循環結構自調節模型
微循環結構自調節模型主要包含兩部分：① 血流動力學與氧分壓模型；② 結構自調節模型。模型結構框圖如圖 1 所示。首先，基于動物實驗獲得的大鼠腸系膜微循環網絡的形態、拓撲和邊界條件數據（德國柏林 Charité 醫學中心提供）建立模型，計算血流動力學參數（包括血壓、血流速度等）和新陳代謝參數（氧分壓等）。隨后，結構自調節模型根據血流動力學和新陳代謝水平計算相應的刺激信號，調節血管管徑。最后，基于調節后的血管結構，再次仿真血流動力學和新陳代謝參數。上述過程迭代進行，直至血管結構調節收斂。以下分別介紹兩部分模型的建模方法。
圖1
微循環結構自調節模型框圖
Figure1.
Block diagram of the microcirculatory structural adaptation model
1.1 血流動力學與氧分壓模型
根據泊肅葉公式，每段血管的阻力 R、壓降 ?P 和血流量 Q 的關系可表示為
 |
在血管分叉或匯聚節點處，根據流量守恒可得
 |
其中 P0 為該節點的血壓值，下標 1–3 表示連接該節點的三段血管。血管阻力可根據
 |
計算得到。其中，
為血液粘滯度，L 為血管長度，r 為血管半徑。對血管網絡中的每個節點列出式(2)，并設置合理的邊界條件，即可求解血管網絡中每段血管的血壓和流量水平。
所研究的大鼠腸系膜血管網絡具有 36 段邊界血管（5 段輸出血管，31 段輸入血管）。其中 1 段為主輸入血管，1 段為主輸出血管，其它為次級邊界。對于主輸入血管，根據 Charité 醫學中心提供的動物實驗數據將其邊界流速設置為 10.064 mm/s。對于主輸出血管，根據文獻數據將其邊界血壓設為 13.95 mm Hg[13]。對于其它次級邊界，基于此前建立的 PSO 算法設置邊界流速值。設置流程如圖 2 所示，具體算法說明詳見文獻[14]。
圖2
血流動力學模型邊界條件設置流程
Figure2.
Procedure of setting the boundary condition for the hemodynamic model
微循環血流的非牛頓現象顯著。模型考慮了血液粘滯度關于血細胞比容和血管管徑的定量關系。首先，模型在血管分叉處應用血漿與血細胞不均勻分配的經驗公式，在血管匯合處應用血細胞數量守恒原則計算整個血管網絡中血細胞比容的分布。隨后，應用 Fahraeus-Lindqvist 效應的經驗公式，依據每段血管的管徑和血細胞比容計算血液粘滯度，用于計算血管阻力[15]。上述過程迭代進行，直至血液粘滯度變化收斂。具體計算流程及公式詳見文獻[6]。
得到每段血管的血流量后，根據血流攜帶氧氣的能力和平均組織氧氣消耗速率計算氧分壓在血管網絡中的分布，具體公式詳見文獻[9]。
1.2 結構自調節機制
微循環血管結構在多種因素的刺激下會發生自我調節，以適應周圍環境的變化。在模型中，考慮血流剪切信號、血管透壁壓信號、組織新陳代謝信號和傳導信號對血管管徑的調節作用[9]。以下分別介紹各信號的計算方法。
（1）血流剪切信號
血流剪切力
的計算公式為
 |
其中，Q 為血流量，
為血液粘滯度，D 為血管管徑。每段血管中血流剪切信號 S
的計算公式為
 |
其中，
為血流剪切力參考值。由于剪切力數值變化范圍較大，故對其進行對數處理，使其與其它信號的數值變化范圍相當。
（2）血管透壁壓信號
血流剪切信號假設微動脈與微靜脈對血流剪切的響應程度是一致的，因而單純基于該信號仿真會導致動靜脈管徑較為接近，與實際情況不符。由于在動物實驗中發現血管網絡中剪切力隨透壁壓的上升而增大[16]，因而考慮通過剪切力相關的血管透壁壓信號對血流剪切信號進行修正。每段血管中透壁壓信號 Sp 的計算公式為
 |
其中，
表示受血管透壁壓影響的血流剪切力水平。基于大鼠腸系膜實驗數據[16]，剪切力與透壁壓的關系式為
 |
該修正可使微動脈處（高透壁壓）的調節信號較小，微靜脈處（低透壁壓）的調節信號較大，使調節后的網絡具有較細的微動脈和較粗的微靜脈，與實際生理情況一致。
（3）新陳代謝信號
新陳代謝信號的傳遞介質為血流，因此它由上游向下游傳遞。新陳代謝信號傳遞量
的計算公式為
 |
式中，上標 up 和 down 分別表示上游與下游血管，
為血管長度，
為氧分壓，
為氧分壓參考值，根據文獻設為 93.2 mm Hg[9]。新陳代謝信號
與傳遞量的關系為
 |
其中 Qref 是血管段的參考血流量。
（4）傳導信號
組織產生的氧氣需求信號需以傳導信號（conduc-ted signal）的形式向上游傳遞。該傳遞的介質主要是血管內皮細胞間的間連蛋白[17]。傳導信號傳遞量 Jc 的傳導公式為
 |
式中，Lref 為血管長度參考值。傳導信號 Sc 的計算公式為
 |
其中，Jref 為傳遞量參考值。
上述四類刺激信號作用的總效果為[9]
 |
其中，參數 kp、km 和 kc 分別是透壁壓信號、新陳代謝信號和傳導信號相對于剪切力信號的作用程度。參數 ks 表達了沒有增加信號時血管的基本收縮趨勢。每段血管的管徑調節量 ?D 是上述四個刺激信號作用的總和 Stot 與仿真步長 ?t 和初始管徑 D 的乘積，即
 |
當血管網絡自調節過程到達平衡時，管徑不再隨著時間變化（?D = 0），對應于 Stot = 0，即這些影響管徑結構狀態的刺激信號相互協同達到平衡。在仿真中，當 Stot < 1.0 × 10 –5 時認為自調節過程收斂。
2 模型參數優化
在模型中，共有 8 個參數需要通過優化確定，包括 4 個自調節因素的比例系數（kp，km，kc 和 ks）和 4 個計算中涉及的參考值（
，Qref，Lref 和 Jref）。本節將介紹參數的優化方法。
2.1 最優化適應度函數
如 1.2 節所述，當模型所考慮的四個自調節因素達到平衡時，微循環系統的結構與血流動力學水平達到了一個新的狀態。本研究對參數優化的目標是使模型仿真計算得到的血流流速與動物實驗觀察值的誤差最小。因此，最優化的適應度函數即為血流流速誤差
 |
式中
為仿真計算得到的第 i 段血管的血流速度，
為動物實驗測量的流速，n 為血管網絡中血管的段數。
2.2 標準 PSO 算法
PSO 算法是一種受生物行為啟發的智能優化算法。它初始化一群隨機粒子（即隨機解），通過迭代尋找最優解。在本文中，一個粒子即為上述 8 個需優化的參數組成的向量。每個粒子的特征包括它所在的位置 X 以及它飛行的速度 V。應用每個粒子可計算得到一個適應度函數值。粒子的位置和飛行速度根據適應度函數值更新，更新公式為
 |
式中下標 i 和 d 分別表示第 i 個粒子和粒子中的第 d 個參數，t 表示當前的迭代輪次，
、
分別表示單個解在迭代過程中最好的位置和整個基本解群組中全局最佳的位置。算法中下次迭代計算的粒子位置
由當前粒子位置
和粒子速度
求和得到，而粒子速度
與前一次粒子速度
、當前位置與本粒子最佳位置之差
、當前位置和全局最佳位置之差
三者有關。ω 為慣性系數，c1、c2 為學習因子，r1、r2 是 0 到 1 之間的隨機數。
2.3 量子 PSO 算法及其改進
2.3.1 標準量子 PSO 算法 標準 PSO 在優化多峰函數時，可能會發生早熟現象，無法跳出局部最優解。量子粒子群（quantum-behaved PSO，QPSO）是一種深度改進的 PSO 算法。它借鑒了量子力學中波函數的描述，認為粒子出現在解空間的某一個位置滿足一定的 Delta 勢阱函數，用概率描述的方式表達粒子在解集合的分布。QPSO 的粒子更新公式為
 |
其中 M 為粒子個數，α 和 u 為 0 到 1 之間的隨機數。β 為壓縮擴展系數，其計算公式為
 |
其中 t 表示當前的迭代輪次，Niter 為仿真所設的最大迭代次數。
2.3.2 改進的量子粒子群算法 鑒于本文所研究的優化問題具有較強的非線性和不連續性的特點，對 QPSO 進行了改進，提出了一種帶選擇算子的隨機壓縮擴展 QPSO 算法。首先，為增強算法的全局搜索能力，提出一種隨機壓縮擴展系數的算子，將 β 設置為
 |
為保證 β 小于 1.78（保證算法收斂[12]），將基礎壓縮擴展系數 βbase 設為 0.78，設 γ 為 0 到 1 之間的隨機數。壓縮擴展系數的隨機化有助于增強算法在全局范圍內尋優的能力。
其次，考慮到參數的生理意義，需要對參數施加一定的約束條件。常用的處理方法是當粒子越過邊界時，直接將粒子賦值為邊界上的值。但該方法容易使粒子陷入邊界周邊難以跳出。本文采用邊界變異的方法，使粒子會以一定的概率重新回到整個搜索域中。具體變異方法為
 |
(19)
式中的
是粒子當前位置，
、
是粒子中第 d 個參數的上下界，r 是 0 到 1 之間的隨機數，c 是變異常量，一般取 0 到 0.5 之間。
最后，引入選擇算子。在計算適應度函數 EV 后，將粒子位置按照 EV 的值進行排序，淘汰優化結果較差（EV 值較大）的一部分粒子（隨機選取 20%～50%），以相同數量的較優粒子填補被淘汰的粒子的位置，以期在下一輪迭代時變異出更好的粒子，增強算法的局部搜索性能。
2.4 仿真方法
仿真環境為 Intel i5-4460 處理器，8 GB 內存，1 TB 硬盤，Windows 10（64 位）操作系統。血流動力學模型的線性方程組求解采用 Intel Math Kernel Library （MKL）11.3。
首先，對比標準 PSO、標準 QPSO 和改進的 QPSO 算法的優化性能。具體做法為：采用每種算法優化 10 次。各算法均采用 20 個隨機粒子作為初始種群。標準 PSO 算法權重 ω 參數值設為 0.7，c1、c2 分別設為 1.8 和 2.3。其次，使用改進的 QPSO 優化模型參數進行 70 次參數優化，并分析仿真結果。
在仿真中，若系數組合不合適，血管網絡自調節過程可能無法收斂，因而無法計算適應度函數值。此時我們將適應度函數值設為一個極大值（100），以確保優化能夠繼續進行。
3 仿真結果及分析
3.1 算法對比
不同優化算法的性能對比如表 1 所示。改進的 QPSO 具有最強的收斂能力。在所進行的 10 組測試中，標準 PSO 算法有 5 組無法收斂，標準 QPSO 算法有 1 組無法收斂，而改進的 QPSO 算法全部得到了收斂的結果。從 EV 和最佳 EV 結果來看，改進的 QPSO 得到的結果也小于標準 PSO 和標準 QPSO 得到的結果，表明改進的 QPSO 能夠使模型的仿真結果與動物實驗結果更為接近。
表1
不同優化算法性能對比
Table1.
Comparison among different optimization algorithms
3.2 改進 QPSO 算法優化結果
為獲得最優的模型參數，本文基于改進的 QPSO 算法進行了 70 次優化，優化的平均 EV 值為 0.67，共有 8 組優化的結果優于文獻報道的使用 Downhill 方法得到的最優 EV 值 0.57[11]，其中最優結果為 0.51，表明經過改進的 QPSO 算法優化后，模型仿真的血流速度更接近動物實驗觀測值。每次仿真平均的迭代次數為 421 次。其中，小于 1 000 次迭代的仿真有 63 次，表明優化具有很好的收斂性。
基于優化的參數進行仿真使調節后的血管結構相比于文獻的報道而言更接近實際測量值。仿真得到的平均管徑誤差為 12.56%，顯著低于基于文獻參數得到的 24.75%（t 檢驗，P<0.001）。圖 3 以偽彩圖的形式呈現了血管網絡中血壓的分布。血管形態及拓撲根據動物實驗測量的數據繪制。血管的顏色與 10～90 mm Hg 范圍內相應區段的血壓值對應。仿真得到的血壓在血管網絡入口處為 88.83 mm Hg，在微動脈和毛細血管中迅速下降，在微靜脈中維持在較低的水平，符合一般的生理認識。整個網絡中血壓下降了 74.88 mm Hg，與文獻報道的大鼠腸系膜微循環中的血壓下降程度相符[13]。
圖3
大鼠腸系膜血管網絡中血壓分布圖
Figure3.
The distribution of blood pressure in the rat mesenteric network
3.3 參數優化結果的分析
根據公式（12），四個刺激信號的比例關系可表示為
 |
圖 4 呈現了文獻報道以及兩組最優仿真中各部分刺激信號的比例關系。從圖中可以發現兩點：① 本文優化的結果與文獻報道的結果有一定差異，體現在 Sp 占比減小，Sm、
和 Sc 占比增大。② 本文的兩組最優仿真之間也存在一定差異，最優仿真對應于更小的 Sp 比重和更大的 Sm、Sτ 比重。這一結果說明在模型中減小透壁壓刺激的比重，增加其它刺激因素的作用，有助于使仿真更接近于真實情況。這表明在該模型的原參數設置中透壁壓信號的作用可能被高估了，這需要通過后續動物實驗進一步驗證。
圖4
仿真和文獻中各刺激信號的比例關系
Figure4.
Comparison of the proportion of each stimulus between the simulation and the literature
4 結論
本文介紹了在建立微循環網絡結構自適應調節模型中，采用改進的 QPSO 算法解決模型參數的優化問題。結果顯示該優化方法優于傳統算法，能使模型的仿真結果與動物實驗數據更為接近。但模型仿真與動物實驗的結果仍存在一定差異，這可能是由于當前模型中每段血管對刺激信號的響應程度是一致的，而實際生理情況下，每段血管的反應可能有所不同。這還需要今后結合大量動物實驗加以改進。
引言
在人體循環系統中，微循環為組織提供氧氣和營養物質，并帶走新陳代謝產物，對于人體正常的生理功能起到重要的作用。微循環的結構性自調節是指血管系統響應組織和器官的功能需求調節血管的舒張與收縮。它確保血液充分灌注器官，使其正常運作，是微循環完成其生理功能的主要機制之一[1]。微循環結構自調節功能異常與高血壓、冠心病等疾病的發生與發展密切相關[2-4]。
生理實驗和臨床試驗是研究循環系統生理病理特性的傳統方法。然而，由于微循環的空間尺度較小，傳統方法難以定量、持續地觀察其生理狀態，造成了微循環系統研究的困難。數學模型與仿真技術能夠有效輔助微循環系統的生理功能研究。一些模型重點關注真實微血管網絡中的血流動力學特性[5-6]，另一些模型則研究簡單血管網絡響應各類外部刺激而發生結構改變的過程[7-8]，但鮮有考慮真實微血管網絡中血管結構自調節過程的數學模型。Pries 等[9]提出了一個描述真實微循環網絡中血管結構在血流動力學與新陳代謝因素刺激下進行自我調節的數學模型。該模型研究了血管剪切力、血管透壁壓、新陳代謝信號、傳導信號四個刺激因素對血管管徑調節的共同作用。基于此模型，Pries 等[10]發現了腫瘤組織中動靜脈短路病理性擴張的一種信號傳遞機制，Vorderwülbecke 等[11]研究了血流剪切力介導的內皮間聯蛋白生成對于微循環灌注的作用，這些研究表明該模型能夠比較全面地反映微循環網絡中的結構自調節過程，從而輔助生理、病理研究。
該模型中，四類刺激信號的比例關系等參數難以基于動物實驗數據確定。Pries 團隊[9-11]采用了 Do-wnhill Simplex 最優化方法確定模型參數。該方法盡管合理地設置了參數，但存在著初值依靠經驗、優化限于局部等限制，因而難以推廣至其它微循環網絡的研究，并且無法確定參數設定的最優性。為克服上述局限性，同時考慮所研究的參數優化問題具有高維度、非線性、不連續等特點，本研究采用具有量子行為的粒子群優化（particle swarm opti-mization，PSO）算法[12]。在其基礎上，本研究改進了其壓縮-擴展系數設置，加入了粒子選擇算子，使算法具有更強的全局優化能力，以期獲得全局最優的模型參數，提高模型的可靠性、有效性和適用性。
1 微循環結構自調節模型
微循環結構自調節模型主要包含兩部分：① 血流動力學與氧分壓模型；② 結構自調節模型。模型結構框圖如圖 1 所示。首先，基于動物實驗獲得的大鼠腸系膜微循環網絡的形態、拓撲和邊界條件數據（德國柏林 Charité 醫學中心提供）建立模型，計算血流動力學參數（包括血壓、血流速度等）和新陳代謝參數（氧分壓等）。隨后，結構自調節模型根據血流動力學和新陳代謝水平計算相應的刺激信號，調節血管管徑。最后，基于調節后的血管結構，再次仿真血流動力學和新陳代謝參數。上述過程迭代進行，直至血管結構調節收斂。以下分別介紹兩部分模型的建模方法。
圖1
微循環結構自調節模型框圖
Figure1.
Block diagram of the microcirculatory structural adaptation model
1.1 血流動力學與氧分壓模型
根據泊肅葉公式，每段血管的阻力 R、壓降 ?P 和血流量 Q 的關系可表示為
|
在血管分叉或匯聚節點處，根據流量守恒可得
|
其中 P0 為該節點的血壓值，下標 1–3 表示連接該節點的三段血管。血管阻力可根據
|
計算得到。其中，
為血液粘滯度，L 為血管長度，r 為血管半徑。對血管網絡中的每個節點列出式(2)，并設置合理的邊界條件，即可求解血管網絡中每段血管的血壓和流量水平。
所研究的大鼠腸系膜血管網絡具有 36 段邊界血管（5 段輸出血管，31 段輸入血管）。其中 1 段為主輸入血管，1 段為主輸出血管，其它為次級邊界。對于主輸入血管，根據 Charité 醫學中心提供的動物實驗數據將其邊界流速設置為 10.064 mm/s。對于主輸出血管，根據文獻數據將其邊界血壓設為 13.95 mm Hg[13]。對于其它次級邊界，基于此前建立的 PSO 算法設置邊界流速值。設置流程如圖 2 所示，具體算法說明詳見文獻[14]。
圖2
血流動力學模型邊界條件設置流程
Figure2.
Procedure of setting the boundary condition for the hemodynamic model
微循環血流的非牛頓現象顯著。模型考慮了血液粘滯度關于血細胞比容和血管管徑的定量關系。首先，模型在血管分叉處應用血漿與血細胞不均勻分配的經驗公式，在血管匯合處應用血細胞數量守恒原則計算整個血管網絡中血細胞比容的分布。隨后，應用 Fahraeus-Lindqvist 效應的經驗公式，依據每段血管的管徑和血細胞比容計算血液粘滯度，用于計算血管阻力[15]。上述過程迭代進行，直至血液粘滯度變化收斂。具體計算流程及公式詳見文獻[6]。
得到每段血管的血流量后，根據血流攜帶氧氣的能力和平均組織氧氣消耗速率計算氧分壓在血管網絡中的分布，具體公式詳見文獻[9]。
1.2 結構自調節機制
微循環血管結構在多種因素的刺激下會發生自我調節，以適應周圍環境的變化。在模型中，考慮血流剪切信號、血管透壁壓信號、組織新陳代謝信號和傳導信號對血管管徑的調節作用[9]。以下分別介紹各信號的計算方法。
（1）血流剪切信號
血流剪切力
的計算公式為
|
其中，Q 為血流量，
為血液粘滯度，D 為血管管徑。每段血管中血流剪切信號 S
的計算公式為
|
其中，
為血流剪切力參考值。由于剪切力數值變化范圍較大，故對其進行對數處理，使其與其它信號的數值變化范圍相當。
（2）血管透壁壓信號
血流剪切信號假設微動脈與微靜脈對血流剪切的響應程度是一致的，因而單純基于該信號仿真會導致動靜脈管徑較為接近，與實際情況不符。由于在動物實驗中發現血管網絡中剪切力隨透壁壓的上升而增大[16]，因而考慮通過剪切力相關的血管透壁壓信號對血流剪切信號進行修正。每段血管中透壁壓信號 Sp 的計算公式為
|
其中，
表示受血管透壁壓影響的血流剪切力水平。基于大鼠腸系膜實驗數據[16]，剪切力與透壁壓的關系式為
|
該修正可使微動脈處（高透壁壓）的調節信號較小，微靜脈處（低透壁壓）的調節信號較大，使調節后的網絡具有較細的微動脈和較粗的微靜脈，與實際生理情況一致。
（3）新陳代謝信號
新陳代謝信號的傳遞介質為血流，因此它由上游向下游傳遞。新陳代謝信號傳遞量
的計算公式為
|
式中，上標 up 和 down 分別表示上游與下游血管，
為血管長度，
為氧分壓，
為氧分壓參考值，根據文獻設為 93.2 mm Hg[9]。新陳代謝信號
與傳遞量的關系為
|
其中 Qref 是血管段的參考血流量。
（4）傳導信號
組織產生的氧氣需求信號需以傳導信號（conduc-ted signal）的形式向上游傳遞。該傳遞的介質主要是血管內皮細胞間的間連蛋白[17]。傳導信號傳遞量 Jc 的傳導公式為
|
式中，Lref 為血管長度參考值。傳導信號 Sc 的計算公式為
|
其中，Jref 為傳遞量參考值。
上述四類刺激信號作用的總效果為[9]
|
其中，參數 kp、km 和 kc 分別是透壁壓信號、新陳代謝信號和傳導信號相對于剪切力信號的作用程度。參數 ks 表達了沒有增加信號時血管的基本收縮趨勢。每段血管的管徑調節量 ?D 是上述四個刺激信號作用的總和 Stot 與仿真步長 ?t 和初始管徑 D 的乘積，即
|
當血管網絡自調節過程到達平衡時，管徑不再隨著時間變化（?D = 0），對應于 Stot = 0，即這些影響管徑結構狀態的刺激信號相互協同達到平衡。在仿真中，當 Stot < 1.0 × 10 –5 時認為自調節過程收斂。
2 模型參數優化
在模型中，共有 8 個參數需要通過優化確定，包括 4 個自調節因素的比例系數（kp，km，kc 和 ks）和 4 個計算中涉及的參考值（
，Qref，Lref 和 Jref）。本節將介紹參數的優化方法。
2.1 最優化適應度函數
如 1.2 節所述，當模型所考慮的四個自調節因素達到平衡時，微循環系統的結構與血流動力學水平達到了一個新的狀態。本研究對參數優化的目標是使模型仿真計算得到的血流流速與動物實驗觀察值的誤差最小。因此，最優化的適應度函數即為血流流速誤差
|
式中
為仿真計算得到的第 i 段血管的血流速度，
為動物實驗測量的流速，n 為血管網絡中血管的段數。
2.2 標準 PSO 算法
PSO 算法是一種受生物行為啟發的智能優化算法。它初始化一群隨機粒子（即隨機解），通過迭代尋找最優解。在本文中，一個粒子即為上述 8 個需優化的參數組成的向量。每個粒子的特征包括它所在的位置 X 以及它飛行的速度 V。應用每個粒子可計算得到一個適應度函數值。粒子的位置和飛行速度根據適應度函數值更新，更新公式為
|
式中下標 i 和 d 分別表示第 i 個粒子和粒子中的第 d 個參數，t 表示當前的迭代輪次，
、
分別表示單個解在迭代過程中最好的位置和整個基本解群組中全局最佳的位置。算法中下次迭代計算的粒子位置
由當前粒子位置
和粒子速度
求和得到，而粒子速度
與前一次粒子速度
、當前位置與本粒子最佳位置之差
、當前位置和全局最佳位置之差
三者有關。ω 為慣性系數，c1、c2 為學習因子，r1、r2 是 0 到 1 之間的隨機數。
2.3 量子 PSO 算法及其改進
2.3.1 標準量子 PSO 算法 標準 PSO 在優化多峰函數時，可能會發生早熟現象，無法跳出局部最優解。量子粒子群（quantum-behaved PSO，QPSO）是一種深度改進的 PSO 算法。它借鑒了量子力學中波函數的描述，認為粒子出現在解空間的某一個位置滿足一定的 Delta 勢阱函數，用概率描述的方式表達粒子在解集合的分布。QPSO 的粒子更新公式為
|
其中 M 為粒子個數，α 和 u 為 0 到 1 之間的隨機數。β 為壓縮擴展系數，其計算公式為
|
其中 t 表示當前的迭代輪次，Niter 為仿真所設的最大迭代次數。
2.3.2 改進的量子粒子群算法 鑒于本文所研究的優化問題具有較強的非線性和不連續性的特點，對 QPSO 進行了改進，提出了一種帶選擇算子的隨機壓縮擴展 QPSO 算法。首先，為增強算法的全局搜索能力，提出一種隨機壓縮擴展系數的算子，將 β 設置為
|
為保證 β 小于 1.78（保證算法收斂[12]），將基礎壓縮擴展系數 βbase 設為 0.78，設 γ 為 0 到 1 之間的隨機數。壓縮擴展系數的隨機化有助于增強算法在全局范圍內尋優的能力。
其次，考慮到參數的生理意義，需要對參數施加一定的約束條件。常用的處理方法是當粒子越過邊界時，直接將粒子賦值為邊界上的值。但該方法容易使粒子陷入邊界周邊難以跳出。本文采用邊界變異的方法，使粒子會以一定的概率重新回到整個搜索域中。具體變異方法為
|
(19)
式中的
是粒子當前位置，
、
是粒子中第 d 個參數的上下界，r 是 0 到 1 之間的隨機數，c 是變異常量，一般取 0 到 0.5 之間。
最后，引入選擇算子。在計算適應度函數 EV 后，將粒子位置按照 EV 的值進行排序，淘汰優化結果較差（EV 值較大）的一部分粒子（隨機選取 20%～50%），以相同數量的較優粒子填補被淘汰的粒子的位置，以期在下一輪迭代時變異出更好的粒子，增強算法的局部搜索性能。
2.4 仿真方法
仿真環境為 Intel i5-4460 處理器，8 GB 內存，1 TB 硬盤，Windows 10（64 位）操作系統。血流動力學模型的線性方程組求解采用 Intel Math Kernel Library （MKL）11.3。
首先，對比標準 PSO、標準 QPSO 和改進的 QPSO 算法的優化性能。具體做法為：采用每種算法優化 10 次。各算法均采用 20 個隨機粒子作為初始種群。標準 PSO 算法權重 ω 參數值設為 0.7，c1、c2 分別設為 1.8 和 2.3。其次，使用改進的 QPSO 優化模型參數進行 70 次參數優化，并分析仿真結果。
在仿真中，若系數組合不合適，血管網絡自調節過程可能無法收斂，因而無法計算適應度函數值。此時我們將適應度函數值設為一個極大值（100），以確保優化能夠繼續進行。
3 仿真結果及分析
3.1 算法對比
不同優化算法的性能對比如表 1 所示。改進的 QPSO 具有最強的收斂能力。在所進行的 10 組測試中，標準 PSO 算法有 5 組無法收斂，標準 QPSO 算法有 1 組無法收斂，而改進的 QPSO 算法全部得到了收斂的結果。從 EV 和最佳 EV 結果來看，改進的 QPSO 得到的結果也小于標準 PSO 和標準 QPSO 得到的結果，表明改進的 QPSO 能夠使模型的仿真結果與動物實驗結果更為接近。
表1
不同優化算法性能對比
Table1.
Comparison among different optimization algorithms
3.2 改進 QPSO 算法優化結果
為獲得最優的模型參數，本文基于改進的 QPSO 算法進行了 70 次優化，優化的平均 EV 值為 0.67，共有 8 組優化的結果優于文獻報道的使用 Downhill 方法得到的最優 EV 值 0.57[11]，其中最優結果為 0.51，表明經過改進的 QPSO 算法優化后，模型仿真的血流速度更接近動物實驗觀測值。每次仿真平均的迭代次數為 421 次。其中，小于 1 000 次迭代的仿真有 63 次，表明優化具有很好的收斂性。
基于優化的參數進行仿真使調節后的血管結構相比于文獻的報道而言更接近實際測量值。仿真得到的平均管徑誤差為 12.56%，顯著低于基于文獻參數得到的 24.75%（t 檢驗，P<0.001）。圖 3 以偽彩圖的形式呈現了血管網絡中血壓的分布。血管形態及拓撲根據動物實驗測量的數據繪制。血管的顏色與 10～90 mm Hg 范圍內相應區段的血壓值對應。仿真得到的血壓在血管網絡入口處為 88.83 mm Hg，在微動脈和毛細血管中迅速下降，在微靜脈中維持在較低的水平，符合一般的生理認識。整個網絡中血壓下降了 74.88 mm Hg，與文獻報道的大鼠腸系膜微循環中的血壓下降程度相符[13]。
圖3
大鼠腸系膜血管網絡中血壓分布圖
Figure3.
The distribution of blood pressure in the rat mesenteric network
3.3 參數優化結果的分析
根據公式（12），四個刺激信號的比例關系可表示為
|
圖 4 呈現了文獻報道以及兩組最優仿真中各部分刺激信號的比例關系。從圖中可以發現兩點：① 本文優化的結果與文獻報道的結果有一定差異，體現在 Sp 占比減小，Sm、
和 Sc 占比增大。② 本文的兩組最優仿真之間也存在一定差異，最優仿真對應于更小的 Sp 比重和更大的 Sm、Sτ 比重。這一結果說明在模型中減小透壁壓刺激的比重，增加其它刺激因素的作用，有助于使仿真更接近于真實情況。這表明在該模型的原參數設置中透壁壓信號的作用可能被高估了，這需要通過后續動物實驗進一步驗證。
圖4
仿真和文獻中各刺激信號的比例關系
Figure4.
Comparison of the proportion of each stimulus between the simulation and the literature
4 結論
本文介紹了在建立微循環網絡結構自適應調節模型中，采用改進的 QPSO 算法解決模型參數的優化問題。結果顯示該優化方法優于傳統算法，能使模型的仿真結果與動物實驗數據更為接近。但模型仿真與動物實驗的結果仍存在一定差異，這可能是由于當前模型中每段血管對刺激信號的響應程度是一致的，而實際生理情況下，每段血管的反應可能有所不同。這還需要今后結合大量動物實驗加以改進。
