本文提出一種基于左心輔助的血液循環系統控制模型，該模型由左心血液循環系統和旋轉式心臟泵模型耦合而成。考慮心臟泵的水力特性，根據泵相似性原理提出了新型旋轉式心臟泵的數學模型。泵模型耦合心血管系統后采用七階集總參數非線性空間狀態方程表示，控制變量為電機轉速。本文采用泵心臟舒張期最小流量變化斜率作為反饋機制，確保最大灌注量的同時避免抽吸。通過Matlab對模型進行仿真，結果顯示在開環控制下，心衰患者的血液循環系統各項血液動力特性均有改善。當轉速為9 000 r/min時，經左心輔助的患者每搏心輸出量由最初34 mL/s提高至正常水平82 mL/s，左心室壓力-容積環向左偏移并且面積減小，說明輔助裝置能夠明顯改善泵血不足并幫助心室卸載。若持續提高泵轉速，當轉速達到12 800 r/min時，會發生抽吸異常狀況，這是因為轉速過高導致靜脈回流量不足而形成的。經采用反饋控制后，系統可有效避免抽吸現象。本文研究結果表明，基于左心輔助的血液循環系統的控制模型能在確保足量灌注量的同時避免抽吸異常狀態。該模型反映了左心輔助裝置與心血管系統的交互作用，為左心輔助心衰治療和生理控制策略設計提供理論依據。
引用本文： 王芳群, 徐慶, 吳振海, 溫太陽, 吉敬華, 賀照明. 基于左心輔助的血液循環系統的控制研究. 生物醫學工程學雜志, 2016, 33(6): 1075-1083. doi: 10.7507/1001-5515.20160172 復制
引言
隨著生活水平提高，各種心臟病尤其是心衰發病率逐年上升，給患者和社會帶來巨大的經濟負擔。心衰，即心力衰竭，心臟泵血功能衰退是其晚期常見的一組臨床綜合癥狀。據統計，全球超過2 000萬患者受心衰威脅，我國成年人慢性心衰患病率高達0.9%，每年的心衰患者超過450萬[1]。目前，心臟移植是公認的治療心衰晚期的最佳方案，但心臟供體不足限制了其臨床應用。2014年，美國僅完成2 000多例的心臟移植手術，在中國僅為200多例[2]。左心輔助裝置(left ventricular assist device,LVAD)通過輔助左心室泵血功能幫助人體完成血液循環，已成為另一種高效的治療方案。研究左心輔助裝置對血液循環系統的血液動力學影響，實現對血流的生理控制是國際心血管領域關注和研究的熱點。
人體血液循環系統模型是一個非線性時變分布參數模型，建模和仿真是研究這一棘手問題的有效手段。血液循環系統的建模研究起源于19世紀后期Qtto Frank提出著名的彈性腔模型，在此基礎上形成了經典的Windkessel四元件模型[3]。此后，計算機技術促進了血液循環系統建模研究的迅速發展。1996年Zacek等[4]利用15個電子管嚴格按照人體結構搭建出血液循環系統，該模型能提供合理的生理數據。20世紀末清華大學Lu等[5]建立了一系列血液循環系統數字模型研究血液動力學。Simaan等[6]建立五階集總參數心血管網絡模型，能很好地模擬人體正常生理情況下的血流動力學特性。馮敏[7]建立周身血液循環系統集總參數模型，仿真二尖瓣關閉不全時血液動力特性曲線。陳泓[8]建立一個集總參數數學模型，并仿真研究心衰及高血壓病理情況。溫太陽等[9]用時變電阻和電感模擬二尖瓣建立改進型左心血液循環系統，文中進行的仿真實驗能很好地模擬二尖瓣返流現象。
心臟泵數學模型與血液循環系統模型的交互作用也是心衰疾病的研究重點。由于心臟輔助裝置是血液循環系統的一部分，通過對耦合系統的研究可以了解左心輔助裝置的運行狀態及對系統的輔助性能，對心臟泵的設計尤其是控制策略的制定具有重要意義。Choi等[10]通過最小二乘法擬合得到最早的軸流式心臟泵數學模型，該模型為泵壓差關于泵流量和轉速的多元微分方程，被學者廣泛用于心臟輔助裝置控制建模中。在此基礎上Faragallah等[11]通過能量守恒定律推導出改進型泵模型，該模型將機械功率和電功率相互轉化，用電機電流代替轉速作為控制變量，能大幅度簡化實際控制策略的設計難度。
基于以上研究，本文將在課題組已建立的改進型血液循環系統基礎上，通過水力特性和泵的相似性原理提出新型的心臟泵模型，建立優化的基于左心輔助的血液循環系統模型。通過控制泵的轉速，研究不同工況下系統的血液動力學特性，并且采用泵最小流量(心臟舒張期泵流量)變化斜率作為反饋機制，能在提供最大泵輸出量的同時避免抽吸，為心室輔助裝置的設計優化及應用提供理論依據。
1 耦合心臟泵的血液循環系統
左心室是心臟泵血的主要動力源，當左心室功能嚴重衰退、泵血不足時，需左心輔助裝置輔助。隨著旋轉電機制造工藝的進步，適于長期植入人體的左心輔助裝置設計和控制策略研究也不斷完善。本文左心輔助裝置類型為軸流式心臟泵，以并聯方式將軸流泵橋接于左心室和主動脈之間，輔助左心室共同泵血，如圖 1所示。
圖1
左心輔助裝置連接心臟循環示意圖
Figure1.
Cross section of the heart with the cardiovascular-pump model superimposed
1.1 心臟泵模型
用于耦合心血管模型的旋轉泵數學模型最早由Chio等[10]提出，該數學模型是由泵壓差與泵流量以及轉速的關系擬合得到，并未考慮泵的壓差和流量之間的關系，即水力特性關系。為了準確描述心臟泵的水力特性及轉速與壓差、流量之間的關系，本文從泵理想工作點出發，根據泵的水力特性和旋轉泵的相似性原理進行推導。
泵的理想特性曲線為直線(理想狀況下，不考慮任何損失)。當泵的轉速為ω0時，泵的實際揚程H0和理想揚程Ht之間存在以下關系，其中ΔH為水力損失
| $~{{H}_{0}}={{H}_{t}}-\Delta H$ |
由于理想的水力特性關系中，理想揚程Ht和流量Q0之間為線性關系，因此有
| ${{H}_{t}}={{\beta }_{0}}{{Q}_{0}}+b$ |
對于旋轉式心臟泵，在設計工況點附近存在沖擊水力損失，由于該水力損失與流量Q0和設計流量Qopt之間的偏離值ΔQ的平方成正比[12]，同時考慮水力慣性，則綜合得到轉速為ω0時的揚程公式：
| ${{H}_{0}}={{\beta }_{0}}{{Q}_{0}}-k{{({{Q}_{0}}-{{Q}_{opt}})}^{2}}+{{\beta }_{1}}\frac{d{{Q}_{0}}}{dt}+b$ |
根據泵相似理論，某一轉速ω與已知轉速ω0之間的流量和揚程關系如下：
| $\frac{{{Q}_{0}}}{Q}=\frac{{{\omega }_{0}}}{\omega },\frac{{{H}_{0}}}{H}={{(\frac{{{\omega }_{0}}}{\omega })}^{2}}$ |
將式(4)代入式(3)，推廣可得新型旋轉泵數學模型：
| $H=\frac{{{\beta }_{0}}\omega }{{{\omega }_{0}}}\cdot Q-k{{(Q-{{Q}_{opt}}\cdot \frac{\omega }{{{\omega }_{0}}})}^{2}}+b\cdot {{(\frac{\omega }{{{\omega }_{0}}})}^{2}}+\frac{{{\beta }_{1}}\omega }{{{\omega }_{0}}}\cdot \frac{dQ}{dt}$ |
可改寫為：
| $H={{R}_{P}}Q+{{L}_{P}}\frac{dQ}{dt}+{{\beta }_{2}}{{\omega }^{2}}$ |
其中，。
該模型從泵的水力特性曲線出發，同時考慮水力損失，通過控制轉速ω，能夠實現不同工況下對血液循環系統的交互作用。
1.2 耦合模型
運用基爾霍夫定律模擬血液動力學，搭建心血管集總參數電路模型是研究左心輔助裝置與心血管交互作用的有效方法。根據臨床常見情況，假設右心室和肺循環健康正常，則左心室和輔助裝置才是影響血液動力學的主要因素。模型前負荷和肺循環可簡化為左心室的順應性，后負荷可由Windkessel四元件表征[13-14]。在此基礎上，本課題組基于時變電阻二尖瓣模型建立改進型左心血液循環系統，用于模擬健康狀況下血液的動力特性，并仿真模擬二尖瓣關閉期間的反流狀況。旋轉泵和左心血液循環系統耦合后等效電路如圖 2所示，由電阻、電感及電容搭建而成，分別代表血液的黏滯性、慣性以及器官的順應性，用節點電流和回路電壓等效心臟各腔室和血管的血壓和血流，其中H為泵壓差。表 1為等效電路圖系統狀態變量，分別為：左心室壓(left ventricular pressure,LVP)，左心房壓(left atrial pressure,LAP)，動脈壓(arterial pressure,AP)，主動脈壓(aorta pressure,AOP )；模型參數如表 2所示。
圖2
基于左心輔助裝置的左心血液循環系統等效電路圖
Figure2.
Combined cardiovascular and LVAD model
表1
系統狀態變量
Table1.
State variables of the system
表2
模型參數
Table2.
Model parameters
左心室是負責心臟泵血工作的最主要器官，同時也是該電路模型主要動力源，C(t)用以表征左心室壓力-容積變化。本文采用Suga等[15]提出的時變彈性函數的方法對左心室進行建模：
| $E\left( t \right)=\frac{LVP\left( t \right)}{LVV\left( t \right)-{{V}_{0}}}$ |
其中，LVP(t)為左心室壓力，LVV(t)為左心室容積(left ventricular volume,LVV)，V0為左心室在收縮末期血液排空時的初始體積。左心室時變彈性函數又稱時變倒電容，其獨立于心臟的前、后負荷，工程上一般如式(6)和式(7)所示：
| $E\left( t \right)=({{E}_{MAX}}-{{E}_{MIN}}){{E}_{n}}({{t}_{n}})+{{E}_{MIN}}$ |
| ${{E}_{n}}({{t}_{n}})=1.55\left( \frac{{{(\frac{{{t}_{n}}}{0.7})}^{1.9}}}{1+{{(\frac{{{t}_{n}}}{0.7})}^{1.9}}} \right)\left( \frac{1}{1+{{({{t}_{n}}1.17)}^{21.9}}} \right)$ |
為了能夠模擬二尖瓣工作中的返流現象，課題組提出了基于實驗結果的改進型二尖瓣模型[9]。該模型采用時變電阻RM和電感LM代替傳統二極管瓣膜模型，能夠模擬二尖瓣返流現象。如式(8)所示，為RM的表達式，其中Td為心室彈性曲線E(t)達到最大值時間點時的取值，為0.325 s，T1為二尖瓣開始閉合時間點取0.28×Td，T2為二尖瓣完全閉合時間點取0.53×Td，T3為二尖瓣張開時間點取1.3×Td：
| ${{R}_{M}}\left( t \right)=\left\{ \begin{align} &\begin{matrix} 0.004&if~0\le t<T \\ \end{matrix}{{~}_{1}} \\ &\begin{matrix} 1E+11\times {{t}^{12.82}}&if~{{T}_{1}}\le t<{{T}_{2}} \\ \end{matrix}~ \\ &{{\begin{matrix} 15&if~{{T}_{2}}\le t<T \\ \end{matrix}}_{3}} \\ &\begin{matrix} 0.004&f~{{T}_{3}}\le t\le {{t}_{c}} \\ \end{matrix}i \\ \end{align} \right.$ |
為模擬左心室抽 吸，引入時變非線性電阻Rk[16]，表達式如式(9)所示，其中α和X1抽吸函數的權重和發生閾值分別為3.5和1:
| ${{R}_{k}}=\left\{ \begin{matrix} 0&if~{{x}_{1}}\left( t \right)>\overline{{{x}_{1}}} \\ \alpha ({{x}_{1}}\left( t \right)-\overline{{{x}_{1}}})~&if~{{x}_{1}}\left( t \right)\le \overline{{{x}_{1}}} \\ \end{matrix} \right.$ |
根據基爾霍夫定律，耦合左心輔助裝置的心血管循環系統可列寫為七階狀態方程，如式(10)所示，控制變量為u(t)=ω(t)2：
| $\dot{x}={{A}_{p}}\left( t \right)\times {{x}_{i}}+{{B}_{p\left( t \right)}}\times {{C}_{p}}\left( x \right)+b\cdot u(t)$ |
其中
| $b={{[\begin{matrix} 0&0&0&0&0&0&-{{\beta }_{2}}/{{L}^{*}} \\ \end{matrix}]}^{T}}$ |
| ${B_{p\left( t \right)}} = \left( {\begin{array}{*{20}{c}} 0&{ - \frac{1}{{C\left( t \right)}}}\\ 0&0\\ 0&0\\ 0&{\frac{1}{{{C_A}}}}\\ 0&0\\ 0&0\\ 0&0 \end{array}} \right){C_p} = \left( {\begin{array}{*{20}{c}} {\frac{1}{{{R_M}\left( t \right)}}\left( {{x_2} - {x_1}} \right)}\\ {\frac{1}{{{R_A}D}}({x_1} - {x_4})} \end{array}} \right)$ |
| ${A_{p\left( t \right) = }}\left( {\begin{array}{*{20}{c}} { - \frac{{C\prime \left( t \right)}}{{C\left( t \right)}}}&0&0&0&0&{\frac{1}{{C\left( t \right)}}}&{ - \frac{1}{{C\left( t \right)}}}\\ 0&{ - \frac{1}{{{R_S}{C_R}}}}&{\frac{1}{{{R_S}{C_R}}}}&0&0&{ - \frac{1}{{{C_R}}}}&{\frac{1}{{{C_R}}}}\\ 0&{\frac{1}{{{R_S}{C_S}}}}&{ - \frac{1}{{{R_S}{C_S}}}}&0&{\frac{1}{{{C_S}}}}&0&{ - \frac{1}{{{C_S}}}}\\ 0&0&0&0&{ - \frac{1}{{{C_A}}}}&0&{\frac{1}{{{C_A}}}}\\ 0&0&{ - \frac{1}{{{L_S}}}}&{\frac{1}{{{L_S}}}}&{ - \frac{{{R_C}}}{{{L_S}}}}&0&0\\ { - \frac{1}{{{L_M}}}}&{\frac{1}{{{L_M}}}}&0&0&0&{ - \frac{{{R_M}\left( t \right)}}{{{L_M}}}}&0\\ {\frac{1}{{{L^*}}}}&0&0&{ - \frac{1}{{{L^*}}}}&0&0&{ - \frac{{{R^*}}}{{{L^*}}}} \end{array}} \right)$ |
其中，L*=LI+LO+LP，R*=RI+RO+Rk+Rp，D(ξ)=
2 反饋控制
左心輔助裝置控制器的主要目標是根據患者的需求，輔助輸送血液以滿足患者生理要求的心搏出量，同時還需避免發生抽吸異常現象。抽吸是由于靜脈回流量不足或泵導管與心室壁接觸所產生的異常狀況，會導致嚴重的溶血和組織破壞等嚴重醫療事故。對于病變的心臟，目前還沒有能夠完善解決抽吸的辦法。
為了保持最大的心輸出量同時及時檢測和避免抽吸異常，本文采用根據泵最小流量(心臟舒張期)變化斜率作為反饋的控制策略[17]。如圖 3所示為本文采用的泵轉速反饋控制器。該控制器由三部分組成：第一部分提取泵流量每一個波動周期內的最小值；第二部分計算由最小量組成包絡線的斜率dx7/dt，此處x7為泵搏出量x7最小值所組成的包絡線；第三部分為更新轉速模塊，如式(11)所示：
| $\omega \left( k+1 \right)=\omega \left( k \right)+c\frac{d{{{\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{x}}}_{7}}}{dt}$ |
式中c為增益常量(c值應該小于泵轉速轉換速率)，常量大小可由醫生根據患者狀況進行確定，過小會導致較長的調節過程，過高可能導致泵轉速調節過頭甚至引發抽吸現象。當控制器開始工作時，左心輔助裝置啟動轉速低于抽吸點，此時dx7/dt＞0，根據式(12)可知，轉速保持上升，當接近抽吸點時dx7/dt近似為零，此時泵轉速可維持低于抽吸轉速。
圖3
泵轉速反饋控制器
Figure3.
Pump speed feedback controller
3 數值仿真
3.1 數值方法
為了評估基于左心輔助的血液循環系統，本文采用Matlab軟件對模型進行了仿真研究，模擬心臟泵在不同工況下對系統的血液動力學特性的影響和控制器性能。通過改變心室收縮期壓容比(EMAX)，對血液循環系統進行健康狀態和擴張型心肌病(dilate cardial myopathy,DCM)狀態(因心肌擴張而導致心室收縮無力、射血不足，是心衰后期主要癥狀)模擬仿真[18]，為帶有左心輔助裝置治療的模型提供對照。并對心肌擴張狀態心血管系統附加心臟泵，觀察其在開環加速控制下的血液動力特性以及抽吸現象。最后，課題組通過進行反饋控制，評估該反饋控制器對抽吸現象的預防。
3.2 數值結果
3.2.1 開環左心輔助血液循環動力學特性
如圖 4所示，為健康狀態和心肌擴張狀態以及左心輔助裝置輔助下的左心室壓力-容積環。左心室健康狀態下的壓力-容積環飽滿；而當心肌擴張時，為了維持必要心搏出量，心室容積增大而壓力下降，導致壓力-容積環整體左移。可見，心衰導致左心室壓也下降。臨床上環路面積代表心臟每搏功，自然心臟泵血功能受損，心臟每搏功下降，無法提供足夠的灌注量，因而射血出現問題需要裝置輔助泵血。對心肌擴張患者采用左心輔助裝置輔助治療后，其血液動力學行為發生改變且受轉速影響。對比8 000、9 000、10 000 r/min的壓力容積曲線，可以得出，隨著轉速升高，環面積減小并左移，即左心室壓和心室容積均減小。原因在于隨著心臟泵轉速的增加，心臟泵的輔助程度加大，從而對自然心臟的左心室起到卸載作用，減輕患者左心室負荷，降低左心室的壓力和每搏輸出量。
圖4
左心室壓容關系曲線
Figure4.
Left ventricular pressure-volume loops for health and DCM and LVAD supported models
表 3為左心輔助裝置不同轉速下血液循環系統的主要血液動力學參數值。當轉速為8 000 r/min時系統流量(即總搏出量)為4.26 L/min，主動脈收縮壓為80.5 mm Hg，當轉速為10 000 r/min時搏出量達到7.23 L/min，主動脈收縮壓達到82.6 mm Hg，滿足人體正常生理灌注，可見適當程度輔助有助于心衰患者治療。
表3
基于左心輔助裝置的血液循環系統動力學參數
Table3.
Hemodynamic parameter for DCM with LVAD
圖 5為健康狀態、心肌擴張狀態以及采用左心輔助裝置下的血液動力學參數曲線。維持足量的泵血是心臟最重要的功能，心肌擴張狀態下每搏輸出量不足34 mL，在進行左心輔助治療過程中(設置轉速為9 000 r/min)，平均每搏輸出量上升至82 mL，達到健康狀態的心輸出量要求。健康狀態左心室壓峰值為118 mm Hg，主動脈壓在78~118 mm Hg之間變化，而心肌擴張狀態下左心室壓峰值降至83 mm Hg，主動脈壓也隨之下降，在左心輔助后主動脈壓回升到77~100 mm Hg。
圖5
血液動力特性參數
Figure5.
Hemodynamic parameters
可見，在心肌擴張狀態下，左心室功能下降導致泵血不足，心室壓和每搏輸出量均降低。進行左心輔助后心搏出量和平均主動脈壓均有明顯回升。由于心臟泵幫助左心卸載，使得主動脈量上升但左心室壓降低，可在確保灌注量的同時分擔左心室工作壓力，此過程中，主動脈壓高于左心室壓，主動脈瓣保持關閉狀態。
3.2.2 抽吸及反饋控制仿真
抽吸現象是進行左心輔助過程中嚴重的醫療事故，可能給患者帶來生命危險，其主要原因在于泵轉速過高。將心臟泵初步轉速設為8 000 r/min并逐步提高泵轉速，驗證心肌擴張狀態中開環響應下血液動力特性，并模擬抽吸異常。如圖 6所示，隨著泵轉速逐漸升高，泵凈流量隨著轉速增加而上升，當t＞5 s時，左心室壓持續低于主動脈壓，主動脈保持關閉狀態。在48 s時，轉速為12 800 r/min，泵流量波形開始劇烈振蕩，此后泵平均流量不再隨著轉速的增加而上升。這是由于心臟泵轉速過高，使得一個心動周期內泵和左心室的總搏出量大于靜脈回流量，導致發生了抽吸異常情況。
圖6
左心輔助裝置線性加速下血液循環動力學變化及抽空異常現象
Figure6.
Pump flow and pressure of linearly increasing pump speed,and suction
圖 7為采用反饋控制的泵搏出量和心臟腔室壓力，由于起始收縮期泵最小流量斜率dx7/dt較大，經過反饋控制泵流量迅速上升。隨著轉速上升接近12 800 r/min，泵流量下包絡線接近水平，即其斜率dx7/dt近似為零，轉速穩定下來。結果顯示該反饋控制能夠增加并維持轉速在抽吸點以下，在保持足量灌注的同時避免了抽吸現象。
圖7
反饋控制下血液循環系統流量及壓力變化
Figure7.
Pump flow and pressure with feedback control
4 結論
本文根據旋轉泵的相似性原理和水力動力學建立了一種新型旋轉心臟泵模型，并耦合課題組提出的左心血液循環系統模型，采用基于泵最小流量變化斜率的反饋策略，提出了基于旋轉式左心輔助裝置的血液循環系統的改進型控制模型，對左心輔助裝置與血液循環系統的交互作用進行仿真研究，并對過度輔助造成的抽吸異常現象及可采取的避免措施進行了模擬。本文分別針對健康狀態、心肌擴張狀態及采用心臟泵左心輔助后的血液循環系統的動力學特性進行了仿真研究。結果表明，進行左心輔助后，心肌擴張模型患者心搏出量回升至健康水平(人體靜息時每搏心輸出量約為70 mL)，并且主動脈壓也有一定提升，同時對左心室壓起卸載作用，使得左心室壓下降，證明左心輔助裝置對心衰患者能起到良好的輔助卸載作用，且對左心室輔助程度與心臟泵轉速成正相關。但臨床上衰弱的心臟對左心輔助裝置的高轉速特別敏感，在輔助灌注的過程中常伴隨著抽吸異常的危險。本文模擬實驗中，隨著轉速的不斷增加，若不經過反饋控制，到達抽吸點后將會發生嚴重的抽吸異常；而當我們采用預防抽吸的反饋控制后，左心輔助裝置能在提供足量搏出量的同時很好地避免抽吸情況。
該改進型模型能夠較好地模擬心臟泵對血液循環系統動力學特性的影響，可為復雜的人體生理時變系統控制策略研究提供理論參考。但本文基于假設左心為心臟功能主要載體，未考慮右心室和肺循環系統等，為了系統的簡易性也忽略了其對血液循環系統的影響。另外，控制模型的數值結果的有效性還需要借助體外物理實驗加以驗證。 左心輔助過程中，隨著轉速的上升，可能會使左心室壓持續低于主動脈壓，這會導致主動脈瓣一直保持關閉狀態，不利于心臟后期康復。因此，更加符合生理和臨床需求的多目標控制策略仍需深入研究和改進。
引言
隨著生活水平提高，各種心臟病尤其是心衰發病率逐年上升，給患者和社會帶來巨大的經濟負擔。心衰，即心力衰竭，心臟泵血功能衰退是其晚期常見的一組臨床綜合癥狀。據統計，全球超過2 000萬患者受心衰威脅，我國成年人慢性心衰患病率高達0.9%，每年的心衰患者超過450萬[1]。目前，心臟移植是公認的治療心衰晚期的最佳方案，但心臟供體不足限制了其臨床應用。2014年，美國僅完成2 000多例的心臟移植手術，在中國僅為200多例[2]。左心輔助裝置(left ventricular assist device,LVAD)通過輔助左心室泵血功能幫助人體完成血液循環，已成為另一種高效的治療方案。研究左心輔助裝置對血液循環系統的血液動力學影響，實現對血流的生理控制是國際心血管領域關注和研究的熱點。
人體血液循環系統模型是一個非線性時變分布參數模型，建模和仿真是研究這一棘手問題的有效手段。血液循環系統的建模研究起源于19世紀后期Qtto Frank提出著名的彈性腔模型，在此基礎上形成了經典的Windkessel四元件模型[3]。此后，計算機技術促進了血液循環系統建模研究的迅速發展。1996年Zacek等[4]利用15個電子管嚴格按照人體結構搭建出血液循環系統，該模型能提供合理的生理數據。20世紀末清華大學Lu等[5]建立了一系列血液循環系統數字模型研究血液動力學。Simaan等[6]建立五階集總參數心血管網絡模型，能很好地模擬人體正常生理情況下的血流動力學特性。馮敏[7]建立周身血液循環系統集總參數模型，仿真二尖瓣關閉不全時血液動力特性曲線。陳泓[8]建立一個集總參數數學模型，并仿真研究心衰及高血壓病理情況。溫太陽等[9]用時變電阻和電感模擬二尖瓣建立改進型左心血液循環系統，文中進行的仿真實驗能很好地模擬二尖瓣返流現象。
心臟泵數學模型與血液循環系統模型的交互作用也是心衰疾病的研究重點。由于心臟輔助裝置是血液循環系統的一部分，通過對耦合系統的研究可以了解左心輔助裝置的運行狀態及對系統的輔助性能，對心臟泵的設計尤其是控制策略的制定具有重要意義。Choi等[10]通過最小二乘法擬合得到最早的軸流式心臟泵數學模型，該模型為泵壓差關于泵流量和轉速的多元微分方程，被學者廣泛用于心臟輔助裝置控制建模中。在此基礎上Faragallah等[11]通過能量守恒定律推導出改進型泵模型，該模型將機械功率和電功率相互轉化，用電機電流代替轉速作為控制變量，能大幅度簡化實際控制策略的設計難度。
基于以上研究，本文將在課題組已建立的改進型血液循環系統基礎上，通過水力特性和泵的相似性原理提出新型的心臟泵模型，建立優化的基于左心輔助的血液循環系統模型。通過控制泵的轉速，研究不同工況下系統的血液動力學特性，并且采用泵最小流量(心臟舒張期泵流量)變化斜率作為反饋機制，能在提供最大泵輸出量的同時避免抽吸，為心室輔助裝置的設計優化及應用提供理論依據。
1 耦合心臟泵的血液循環系統
左心室是心臟泵血的主要動力源，當左心室功能嚴重衰退、泵血不足時，需左心輔助裝置輔助。隨著旋轉電機制造工藝的進步，適于長期植入人體的左心輔助裝置設計和控制策略研究也不斷完善。本文左心輔助裝置類型為軸流式心臟泵，以并聯方式將軸流泵橋接于左心室和主動脈之間，輔助左心室共同泵血，如圖 1所示。
圖1
左心輔助裝置連接心臟循環示意圖
Figure1.
Cross section of the heart with the cardiovascular-pump model superimposed
1.1 心臟泵模型
用于耦合心血管模型的旋轉泵數學模型最早由Chio等[10]提出，該數學模型是由泵壓差與泵流量以及轉速的關系擬合得到，并未考慮泵的壓差和流量之間的關系，即水力特性關系。為了準確描述心臟泵的水力特性及轉速與壓差、流量之間的關系，本文從泵理想工作點出發，根據泵的水力特性和旋轉泵的相似性原理進行推導。
泵的理想特性曲線為直線(理想狀況下，不考慮任何損失)。當泵的轉速為ω0時，泵的實際揚程H0和理想揚程Ht之間存在以下關系，其中ΔH為水力損失
| $~{{H}_{0}}={{H}_{t}}-\Delta H$ |
由于理想的水力特性關系中，理想揚程Ht和流量Q0之間為線性關系，因此有
| ${{H}_{t}}={{\beta }_{0}}{{Q}_{0}}+b$ |
對于旋轉式心臟泵，在設計工況點附近存在沖擊水力損失，由于該水力損失與流量Q0和設計流量Qopt之間的偏離值ΔQ的平方成正比[12]，同時考慮水力慣性，則綜合得到轉速為ω0時的揚程公式：
| ${{H}_{0}}={{\beta }_{0}}{{Q}_{0}}-k{{({{Q}_{0}}-{{Q}_{opt}})}^{2}}+{{\beta }_{1}}\frac{d{{Q}_{0}}}{dt}+b$ |
根據泵相似理論，某一轉速ω與已知轉速ω0之間的流量和揚程關系如下：
| $\frac{{{Q}_{0}}}{Q}=\frac{{{\omega }_{0}}}{\omega },\frac{{{H}_{0}}}{H}={{(\frac{{{\omega }_{0}}}{\omega })}^{2}}$ |
將式(4)代入式(3)，推廣可得新型旋轉泵數學模型：
| $H=\frac{{{\beta }_{0}}\omega }{{{\omega }_{0}}}\cdot Q-k{{(Q-{{Q}_{opt}}\cdot \frac{\omega }{{{\omega }_{0}}})}^{2}}+b\cdot {{(\frac{\omega }{{{\omega }_{0}}})}^{2}}+\frac{{{\beta }_{1}}\omega }{{{\omega }_{0}}}\cdot \frac{dQ}{dt}$ |
可改寫為：
| $H={{R}_{P}}Q+{{L}_{P}}\frac{dQ}{dt}+{{\beta }_{2}}{{\omega }^{2}}$ |
其中，。
該模型從泵的水力特性曲線出發，同時考慮水力損失，通過控制轉速ω，能夠實現不同工況下對血液循環系統的交互作用。
1.2 耦合模型
運用基爾霍夫定律模擬血液動力學，搭建心血管集總參數電路模型是研究左心輔助裝置與心血管交互作用的有效方法。根據臨床常見情況，假設右心室和肺循環健康正常，則左心室和輔助裝置才是影響血液動力學的主要因素。模型前負荷和肺循環可簡化為左心室的順應性，后負荷可由Windkessel四元件表征[13-14]。在此基礎上，本課題組基于時變電阻二尖瓣模型建立改進型左心血液循環系統，用于模擬健康狀況下血液的動力特性，并仿真模擬二尖瓣關閉期間的反流狀況。旋轉泵和左心血液循環系統耦合后等效電路如圖 2所示，由電阻、電感及電容搭建而成，分別代表血液的黏滯性、慣性以及器官的順應性，用節點電流和回路電壓等效心臟各腔室和血管的血壓和血流，其中H為泵壓差。表 1為等效電路圖系統狀態變量，分別為：左心室壓(left ventricular pressure,LVP)，左心房壓(left atrial pressure,LAP)，動脈壓(arterial pressure,AP)，主動脈壓(aorta pressure,AOP )；模型參數如表 2所示。
圖2
基于左心輔助裝置的左心血液循環系統等效電路圖
Figure2.
Combined cardiovascular and LVAD model
表1
系統狀態變量
Table1.
State variables of the system
表2
模型參數
Table2.
Model parameters
左心室是負責心臟泵血工作的最主要器官，同時也是該電路模型主要動力源，C(t)用以表征左心室壓力-容積變化。本文采用Suga等[15]提出的時變彈性函數的方法對左心室進行建模：
| $E\left( t \right)=\frac{LVP\left( t \right)}{LVV\left( t \right)-{{V}_{0}}}$ |
其中，LVP(t)為左心室壓力，LVV(t)為左心室容積(left ventricular volume,LVV)，V0為左心室在收縮末期血液排空時的初始體積。左心室時變彈性函數又稱時變倒電容，其獨立于心臟的前、后負荷，工程上一般如式(6)和式(7)所示：
| $E\left( t \right)=({{E}_{MAX}}-{{E}_{MIN}}){{E}_{n}}({{t}_{n}})+{{E}_{MIN}}$ |
| ${{E}_{n}}({{t}_{n}})=1.55\left( \frac{{{(\frac{{{t}_{n}}}{0.7})}^{1.9}}}{1+{{(\frac{{{t}_{n}}}{0.7})}^{1.9}}} \right)\left( \frac{1}{1+{{({{t}_{n}}1.17)}^{21.9}}} \right)$ |
為了能夠模擬二尖瓣工作中的返流現象，課題組提出了基于實驗結果的改進型二尖瓣模型[9]。該模型采用時變電阻RM和電感LM代替傳統二極管瓣膜模型，能夠模擬二尖瓣返流現象。如式(8)所示，為RM的表達式，其中Td為心室彈性曲線E(t)達到最大值時間點時的取值，為0.325 s，T1為二尖瓣開始閉合時間點取0.28×Td，T2為二尖瓣完全閉合時間點取0.53×Td，T3為二尖瓣張開時間點取1.3×Td：
| ${{R}_{M}}\left( t \right)=\left\{ \begin{align} &\begin{matrix} 0.004&if~0\le t<T \\ \end{matrix}{{~}_{1}} \\ &\begin{matrix} 1E+11\times {{t}^{12.82}}&if~{{T}_{1}}\le t<{{T}_{2}} \\ \end{matrix}~ \\ &{{\begin{matrix} 15&if~{{T}_{2}}\le t<T \\ \end{matrix}}_{3}} \\ &\begin{matrix} 0.004&f~{{T}_{3}}\le t\le {{t}_{c}} \\ \end{matrix}i \\ \end{align} \right.$ |
為模擬左心室抽 吸，引入時變非線性電阻Rk[16]，表達式如式(9)所示，其中α和X1抽吸函數的權重和發生閾值分別為3.5和1:
| ${{R}_{k}}=\left\{ \begin{matrix} 0&if~{{x}_{1}}\left( t \right)>\overline{{{x}_{1}}} \\ \alpha ({{x}_{1}}\left( t \right)-\overline{{{x}_{1}}})~&if~{{x}_{1}}\left( t \right)\le \overline{{{x}_{1}}} \\ \end{matrix} \right.$ |
根據基爾霍夫定律，耦合左心輔助裝置的心血管循環系統可列寫為七階狀態方程，如式(10)所示，控制變量為u(t)=ω(t)2：
| $\dot{x}={{A}_{p}}\left( t \right)\times {{x}_{i}}+{{B}_{p\left( t \right)}}\times {{C}_{p}}\left( x \right)+b\cdot u(t)$ |
其中
| $b={{[\begin{matrix} 0&0&0&0&0&0&-{{\beta }_{2}}/{{L}^{*}} \\ \end{matrix}]}^{T}}$ |
| ${B_{p\left( t \right)}} = \left( {\begin{array}{*{20}{c}} 0&{ - \frac{1}{{C\left( t \right)}}}\\ 0&0\\ 0&0\\ 0&{\frac{1}{{{C_A}}}}\\ 0&0\\ 0&0\\ 0&0 \end{array}} \right){C_p} = \left( {\begin{array}{*{20}{c}} {\frac{1}{{{R_M}\left( t \right)}}\left( {{x_2} - {x_1}} \right)}\\ {\frac{1}{{{R_A}D}}({x_1} - {x_4})} \end{array}} \right)$ |
| ${A_{p\left( t \right) = }}\left( {\begin{array}{*{20}{c}} { - \frac{{C\prime \left( t \right)}}{{C\left( t \right)}}}&0&0&0&0&{\frac{1}{{C\left( t \right)}}}&{ - \frac{1}{{C\left( t \right)}}}\\ 0&{ - \frac{1}{{{R_S}{C_R}}}}&{\frac{1}{{{R_S}{C_R}}}}&0&0&{ - \frac{1}{{{C_R}}}}&{\frac{1}{{{C_R}}}}\\ 0&{\frac{1}{{{R_S}{C_S}}}}&{ - \frac{1}{{{R_S}{C_S}}}}&0&{\frac{1}{{{C_S}}}}&0&{ - \frac{1}{{{C_S}}}}\\ 0&0&0&0&{ - \frac{1}{{{C_A}}}}&0&{\frac{1}{{{C_A}}}}\\ 0&0&{ - \frac{1}{{{L_S}}}}&{\frac{1}{{{L_S}}}}&{ - \frac{{{R_C}}}{{{L_S}}}}&0&0\\ { - \frac{1}{{{L_M}}}}&{\frac{1}{{{L_M}}}}&0&0&0&{ - \frac{{{R_M}\left( t \right)}}{{{L_M}}}}&0\\ {\frac{1}{{{L^*}}}}&0&0&{ - \frac{1}{{{L^*}}}}&0&0&{ - \frac{{{R^*}}}{{{L^*}}}} \end{array}} \right)$ |
其中，L*=LI+LO+LP，R*=RI+RO+Rk+Rp，D(ξ)=
2 反饋控制
左心輔助裝置控制器的主要目標是根據患者的需求，輔助輸送血液以滿足患者生理要求的心搏出量，同時還需避免發生抽吸異常現象。抽吸是由于靜脈回流量不足或泵導管與心室壁接觸所產生的異常狀況，會導致嚴重的溶血和組織破壞等嚴重醫療事故。對于病變的心臟，目前還沒有能夠完善解決抽吸的辦法。
為了保持最大的心輸出量同時及時檢測和避免抽吸異常，本文采用根據泵最小流量(心臟舒張期)變化斜率作為反饋的控制策略[17]。如圖 3所示為本文采用的泵轉速反饋控制器。該控制器由三部分組成：第一部分提取泵流量每一個波動周期內的最小值；第二部分計算由最小量組成包絡線的斜率dx7/dt，此處x7為泵搏出量x7最小值所組成的包絡線；第三部分為更新轉速模塊，如式(11)所示：
| $\omega \left( k+1 \right)=\omega \left( k \right)+c\frac{d{{{\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{x}}}_{7}}}{dt}$ |
式中c為增益常量(c值應該小于泵轉速轉換速率)，常量大小可由醫生根據患者狀況進行確定，過小會導致較長的調節過程，過高可能導致泵轉速調節過頭甚至引發抽吸現象。當控制器開始工作時，左心輔助裝置啟動轉速低于抽吸點，此時dx7/dt＞0，根據式(12)可知，轉速保持上升，當接近抽吸點時dx7/dt近似為零，此時泵轉速可維持低于抽吸轉速。
圖3
泵轉速反饋控制器
Figure3.
Pump speed feedback controller
3 數值仿真
3.1 數值方法
為了評估基于左心輔助的血液循環系統，本文采用Matlab軟件對模型進行了仿真研究，模擬心臟泵在不同工況下對系統的血液動力學特性的影響和控制器性能。通過改變心室收縮期壓容比(EMAX)，對血液循環系統進行健康狀態和擴張型心肌病(dilate cardial myopathy,DCM)狀態(因心肌擴張而導致心室收縮無力、射血不足，是心衰后期主要癥狀)模擬仿真[18]，為帶有左心輔助裝置治療的模型提供對照。并對心肌擴張狀態心血管系統附加心臟泵，觀察其在開環加速控制下的血液動力特性以及抽吸現象。最后，課題組通過進行反饋控制，評估該反饋控制器對抽吸現象的預防。
3.2 數值結果
3.2.1 開環左心輔助血液循環動力學特性
如圖 4所示，為健康狀態和心肌擴張狀態以及左心輔助裝置輔助下的左心室壓力-容積環。左心室健康狀態下的壓力-容積環飽滿；而當心肌擴張時，為了維持必要心搏出量，心室容積增大而壓力下降，導致壓力-容積環整體左移。可見，心衰導致左心室壓也下降。臨床上環路面積代表心臟每搏功，自然心臟泵血功能受損，心臟每搏功下降，無法提供足夠的灌注量，因而射血出現問題需要裝置輔助泵血。對心肌擴張患者采用左心輔助裝置輔助治療后，其血液動力學行為發生改變且受轉速影響。對比8 000、9 000、10 000 r/min的壓力容積曲線，可以得出，隨著轉速升高，環面積減小并左移，即左心室壓和心室容積均減小。原因在于隨著心臟泵轉速的增加，心臟泵的輔助程度加大，從而對自然心臟的左心室起到卸載作用，減輕患者左心室負荷，降低左心室的壓力和每搏輸出量。
圖4
左心室壓容關系曲線
Figure4.
Left ventricular pressure-volume loops for health and DCM and LVAD supported models
表 3為左心輔助裝置不同轉速下血液循環系統的主要血液動力學參數值。當轉速為8 000 r/min時系統流量(即總搏出量)為4.26 L/min，主動脈收縮壓為80.5 mm Hg，當轉速為10 000 r/min時搏出量達到7.23 L/min，主動脈收縮壓達到82.6 mm Hg，滿足人體正常生理灌注，可見適當程度輔助有助于心衰患者治療。
表3
基于左心輔助裝置的血液循環系統動力學參數
Table3.
Hemodynamic parameter for DCM with LVAD
圖 5為健康狀態、心肌擴張狀態以及采用左心輔助裝置下的血液動力學參數曲線。維持足量的泵血是心臟最重要的功能，心肌擴張狀態下每搏輸出量不足34 mL，在進行左心輔助治療過程中(設置轉速為9 000 r/min)，平均每搏輸出量上升至82 mL，達到健康狀態的心輸出量要求。健康狀態左心室壓峰值為118 mm Hg，主動脈壓在78~118 mm Hg之間變化，而心肌擴張狀態下左心室壓峰值降至83 mm Hg，主動脈壓也隨之下降，在左心輔助后主動脈壓回升到77~100 mm Hg。
圖5
血液動力特性參數
Figure5.
Hemodynamic parameters
可見，在心肌擴張狀態下，左心室功能下降導致泵血不足，心室壓和每搏輸出量均降低。進行左心輔助后心搏出量和平均主動脈壓均有明顯回升。由于心臟泵幫助左心卸載，使得主動脈量上升但左心室壓降低，可在確保灌注量的同時分擔左心室工作壓力，此過程中，主動脈壓高于左心室壓，主動脈瓣保持關閉狀態。
3.2.2 抽吸及反饋控制仿真
抽吸現象是進行左心輔助過程中嚴重的醫療事故，可能給患者帶來生命危險，其主要原因在于泵轉速過高。將心臟泵初步轉速設為8 000 r/min并逐步提高泵轉速，驗證心肌擴張狀態中開環響應下血液動力特性，并模擬抽吸異常。如圖 6所示，隨著泵轉速逐漸升高，泵凈流量隨著轉速增加而上升，當t＞5 s時，左心室壓持續低于主動脈壓，主動脈保持關閉狀態。在48 s時，轉速為12 800 r/min，泵流量波形開始劇烈振蕩，此后泵平均流量不再隨著轉速的增加而上升。這是由于心臟泵轉速過高，使得一個心動周期內泵和左心室的總搏出量大于靜脈回流量，導致發生了抽吸異常情況。
圖6
左心輔助裝置線性加速下血液循環動力學變化及抽空異常現象
Figure6.
Pump flow and pressure of linearly increasing pump speed,and suction
圖 7為采用反饋控制的泵搏出量和心臟腔室壓力，由于起始收縮期泵最小流量斜率dx7/dt較大，經過反饋控制泵流量迅速上升。隨著轉速上升接近12 800 r/min，泵流量下包絡線接近水平，即其斜率dx7/dt近似為零，轉速穩定下來。結果顯示該反饋控制能夠增加并維持轉速在抽吸點以下，在保持足量灌注的同時避免了抽吸現象。
圖7
反饋控制下血液循環系統流量及壓力變化
Figure7.
Pump flow and pressure with feedback control
4 結論
本文根據旋轉泵的相似性原理和水力動力學建立了一種新型旋轉心臟泵模型，并耦合課題組提出的左心血液循環系統模型，采用基于泵最小流量變化斜率的反饋策略，提出了基于旋轉式左心輔助裝置的血液循環系統的改進型控制模型，對左心輔助裝置與血液循環系統的交互作用進行仿真研究，并對過度輔助造成的抽吸異常現象及可采取的避免措施進行了模擬。本文分別針對健康狀態、心肌擴張狀態及采用心臟泵左心輔助后的血液循環系統的動力學特性進行了仿真研究。結果表明，進行左心輔助后，心肌擴張模型患者心搏出量回升至健康水平(人體靜息時每搏心輸出量約為70 mL)，并且主動脈壓也有一定提升，同時對左心室壓起卸載作用，使得左心室壓下降，證明左心輔助裝置對心衰患者能起到良好的輔助卸載作用，且對左心室輔助程度與心臟泵轉速成正相關。但臨床上衰弱的心臟對左心輔助裝置的高轉速特別敏感，在輔助灌注的過程中常伴隨著抽吸異常的危險。本文模擬實驗中，隨著轉速的不斷增加，若不經過反饋控制，到達抽吸點后將會發生嚴重的抽吸異常；而當我們采用預防抽吸的反饋控制后，左心輔助裝置能在提供足量搏出量的同時很好地避免抽吸情況。
該改進型模型能夠較好地模擬心臟泵對血液循環系統動力學特性的影響，可為復雜的人體生理時變系統控制策略研究提供理論參考。但本文基于假設左心為心臟功能主要載體，未考慮右心室和肺循環系統等，為了系統的簡易性也忽略了其對血液循環系統的影響。另外，控制模型的數值結果的有效性還需要借助體外物理實驗加以驗證。 左心輔助過程中，隨著轉速的上升，可能會使左心室壓持續低于主動脈壓，這會導致主動脈瓣一直保持關閉狀態，不利于心臟后期康復。因此，更加符合生理和臨床需求的多目標控制策略仍需深入研究和改進。