在臨床機械通氣過程中選擇合適通氣模式是一個難點。本文采用非線性多室肺模型并優化呼吸氣流模式得到呼吸氣周期的最小呼吸功、最小輔助吸入量、最小呼吸周期彈性勢能和氣流速率變化。使用sigmoidal函數對建立的非線性呼吸功能方程進行平滑處理后轉化成求解非線性邊界條件的兩點邊值問題, 最后使用梯度下降法求解。實驗結果表明肺容積和氣流速率經過優化后具體良好的敏感性和收斂速度。研究結果為開發多變量控制器機械通氣監護重癥患者提供理論基礎。
引用本文: 蔡永銘, 谷凌雁, 陳夫華. 基于非線性多室肺模型優化呼吸氣流模式的研究. 生物醫學工程學雜志, 2015, 32(1): 32-37. doi: 10.7507/1001-5515.20150006 復制
引言
呼吸衰竭通氣不足是危重病急救醫學臨床上常見的問題。一旦確診通氣不足原因后,呼吸衰竭患者通常需要機械通氣支持。但是臨床上在采用機械通氣過程中,由于過量容積或壓力過高,很容易引起急性肺損傷。通過控制算法可以實現部分自動控制容量和壓力。這些算法需要一個參考模型用于辨別出臨床上合理的呼吸模式。然而,目前文獻[1-2]已經提出來的肺呼吸模型通常假設肺的同質性。早期呼吸控制優化研究都是使用簡單的同質肺模型處理呼吸的頻率問題。Otis等[3]通過使用最小工作率準則預測呼吸頻率。這是涉及靜態優化問題,假定氣流模式是一個固定的正弦函數。Yamashiro等[4]研究了采用固定的吸氣和呼氣相呼吸周期預測呼吸氣流模式。這些研究結果被Hamalainen等[5]進行擴展,考慮了兩級分層模型來控制呼吸,上層模型負責對下層模型的氣體模型總體控制變量進行優化,下層模型用最小化高層模型的標準來選擇呼吸參數的方法確定氣流模型。
雖然尋找最佳的呼吸模式問題已經被Sundaresan等[6]提出來,但是在這些研究中呼吸模型都是以固定呼吸參數及單個腔體肺模型來預測控制。然而,人體的肺,特別是患病的肺是異構的。從解剖角度看肺包括許多的肺泡子結構腔室,在這些子結構里面的氣體交換的能力是不同的。仿真的肺模型應該考慮這種肺異質性。此外,氣流的阻力和肺的順應性也不是恒定的,隨著肺容量而變化,肺的順應性也跟肺容量相關。實現機械通氣的自適應控制,采用比單室肺的容積或壓力控制通氣更詳細的控制結構模型是必要的,尤其是要考慮非線性氣體阻力、順應性和肺異質性。氣道中的流量和壓力并不是簡單的波形圖,雖然目前臨床上普遍是這樣描述,但是事實上在這些波形中嵌入了大量的信息。本文將提供一個嚴格的數學框架用于表達非線性多室肺力學模型系統并優化呼吸氣流模式。
本文在陳夫華等[7]前期研究結果基礎上,用非線性多室模型優化肺力學系統呼吸氣流模式。這里的優化是指微積分方程的最優解,而不是在生理學意義上的最優呼吸模式。首先擴展Chellaboina等[8]中的線性多室肺癌模型成非線性模型。然后擴展Hamalainen等[9]吸氣和呼氣的呼吸循環功能函數,采用典型微積分方法推導最佳氣流模式。包括了生理學上的最小化呼吸工作和幫助患者吸氣增加肺容積,最小化彈性勢能和呼氣階段改變氣流速度變化等方面。最后,用求非線性兩點邊界值解確定吸氣和呼氣循環過程的最佳氣流模式。
1 非線性多室肺模型
1.1 單室肺模型
為了使肺模型更加符合實際肺組織結構和工作過程,這里擴展線性多室肺模型成非線性模型并模擬整個吸氣過程,同時假設肺的順應性是關于肺容量的非線性函數[10]。
單室肺模型是通過一個單室腔體表示肺生理解剖結構,如圖 1所示。用非線性函數 C(x)表示單室肺的順應性,R表示氣流的阻力。當時間t=0時,驅動壓力Pin(t)表示加載在父層的氣道開口,Pin(t)可能是由呼吸肌肉或機械通氣作用下產生的。該壓力發生在0≤t≤Tin之間,也就是呼吸周期的吸氣階段,Tin是吸氣時間。當時間t=Tin時,該壓力將被釋放,緊接著被呼氣壓力替代,此時外部壓力等于大氣壓。在Tin≤t≤Tin+Tex時,氣體壓力為Pex(t),Tex是呼氣階段持續的時間。

吸氣階段的狀態方程表示為
$ \begin{align} & {{R}_{\text{in}}}{{x}^{'}}\left(t \right)+\left(1/-{{C}_{\text{in}}}\left(x \right)\right)x\left(t \right)-{{P}_{\text{in}}}\left(t \right), \\ & x\left(0 \right)=x_{0}^{\operatorname{in}}, 0\le t\le {{T}_{\text{in}}}, \\ \end{align} $ |
其中x(t)∈,t≥0,表示肺容量。Rin∈,表示吸氣階段氣體受到的阻力。Cin:→+是一個非線性函數,表示在吸氣階段肺的順應性。x0in∈ +,表示吸氣開始時候的肺容量作為系統初始條件。式(1)實際是一個簡單的壓力平衡方程,應用在肺室模型時候,順應性和氣道阻力與肺容量成比例。這里假設呼氣動作是由于肺室彈性舒展的被動結果。
呼氣階段的狀態方程表示為
$ \begin{align} & {{R}_{\text{ex}}}{{x}^{'}}\left(t \right)+\left(1/-{{C}_{\text{ex}}}\left(x \right)\right)x\left(t \right)-{{P}_{\text{ex}}}\left(t \right), \\ & x\left({{T}_{\text{in}}} \right)=x_{\text{0}}^{\text{ex}}, {{T}_{\text{in}}}\le t\le {{T}_{\text{in}}}+{{T}_{\text{ex}}}, \\ \end{align} $ |
其中x(t)的意義跟吸氣狀態方程一樣,Rex∈ +,表示呼氣階段氣體受到的阻力。Cex: → +是一個非線性函數,表示在呼氣階段肺的順應性。x0ex∈ +,表示呼氣開始時候的肺容量。
1.2 多室肺模型
為了更加接近真實的肺結構,模型中假設支氣管樹成二分叉結構,也就是在每一級氣道單元重復分叉成兩個子的氣道單元,圖 2表示了4個子室的模型肺。這里定義2n(n≥0)室肺模型,由2n個小單元肺室組成,這些肺單元由n層壓力源連接到氣道單元。因為每個氣道分叉成兩個后繼級子氣道,所以整個肺部是由2n個子室組成。

假設xi(i=1, 2, …, 2n)表示第i個肺子室的容量,ciin(i=1, 2, …, 2n)表示對應于第i個肺子室在吸氣階段的順應性。Rj, iin(i=1, 2, …, 2j)表示多肺室模型位于第j層壓力源的第i個肺子室在吸氣階段的氣體阻力。
在單室模型中,假設氣體壓力Pin(t),t≥0,表示吸氣階段由吸氣肌肉或機械通氣作用下產生的壓力。現在可將吸氣階段的狀態方程改為
$ \begin{align} & R_{n}^{\text{in}}{{x}^{'}}_{i}\left(t \right)+\frac{1}{c_{i}^{\text{in}}\left({{x}_{i}}\left(t \right)\right)}{{x}_{i}}\left(t \right)+\sum\limits_{j=0}^{n-1}{R_{j, {{k}_{j}}}^{\text{in}}}\sum\limits_{l=\left({{k}_{j}}-1 \right){{2}^{n-j}}+1}^{{{k}_{j}}{{2}^{n-j}}}{{{x}^{'}}_{l}\left(t \right)={{P}_{\text{in}}}\left(t \right)}, \\ & {{x}_{i}}\left(0 \right)=x_{{{i}_{o}}}^{\text{in}}, 0\le t\le {{T}_{\text{in}}}, i=1, 2, \cdots, {{2}^{n}}, \\ \end{align} $ |
其中ciin(xi),i=1, 2, …,2n是關于xi(i=1, 2, …,2n)的非線性方程。
$ c_{i}^{\text{in}}\left({{x}_{i}} \right)=\left\{ \begin{align} & a_{{{i}_{1}}}^{\text{in}}+b_{{{i}_{1}}}^{\text{in}}{{x}_{i}}, \ \ \ \ \text{if}\ \ \ \text{0}\le {{x}_{i}}\le x_{{{i}_{1}}}^{\text{in}} \\ & a_{{{i}_{2}}}^{\text{in}}, \ \ \ \ \ \ \ \ \ \ \ \text{if}\ \ \ x_{{{i}_{1}}}^{\text{in}}\le {{x}_{i}}\le x_{{{i}_{2}}}^{\text{in}} \\ & a_{{{i}_{3}}}^{\text{in}}+b_{{{i}_{3}}}^{\text{in}}{{x}_{i}}, \ \ \text{if}\ \ \ \ x_{{{i}_{2}}}^{\text{in}}\le {{x}_{i}}\le {{V}_{{{T}_{i}}}} \\ \end{align} \right., $ |
這里aimin(m=1, 2, 3)且bimin(m=1, 3)是模型參數,bi1in>0且bi3in<0, ximin(m=1, 2), 表示容積變化在xi1in和xi2in范圍內時,順應性為固定常數ai2in,VTi表示潮氣量。
$ {{k}_{j}}=\left[\frac{{{k}_{j+1}}-1}{2} \right]+1, j=0, \cdots, n-1, {{k}_{n}}=i, $ |
其中表示一個不大于k的整數,kj表示第i個子室氣體經過第j級氣道編號。
對于呼氣狀態方程,如在單室模型中假設呼氣過程是被動及外部施加壓力的結果,pex(t),t≥0。類似吸氣過程的處理方法,可以得到呼氣階段的狀態方程:
$ R_{n}^{\text{ex}}{{x}^{'}}_{i}\left(t \right)+\sum\limits_{j=0}^{n-1}{R_{j, {{k}_{j}}}^{\text{ex}}}\sum\limits_{l=\left({{k}_{j}}-1 \right){{2}^{n-j}}+1}^{{{k}_{j}}{{2}^{n-j}}}{{{x}^{'}}_{l}\left(t \right)+}\frac{1}{c_{i}^{\text{ex}}\left({{x}_{i}}\left(t \right)\right)}{{x}_{i}}\left(t \right)={{P}_{\text{ex}}}\left(t \right), {{x}_{i}}\left({{T}_{\text{in}}} \right)=x_{{{i}_{0}}}^{\text{ex}}, {{T}_{\text{in}}}\le t\le {{T}_{\text{in}}}+{{T}_{\text{ex}}}, i=1, 2, \cdots, {{2}^{n}}, $ |
其中ciex(xi)由如下定義:
$ c_{i}^{\text{ex}}\left({{x}_{i}} \right)=\left\{ \begin{align} & a_{{{i}_{0}}}^{\text{ex}}+b_{{{i}_{1}}}^{\text{ex}}{{x}_{i}}, \ \ \ \ \text{if}\ \ \ \text{0}\le {{x}_{i}}\le x_{{{i}_{1}}}^{\text{ex}} \\ & a_{{{i}_{2}}}^{\text{ex}}, \ \ \ \ \ \ \ \ \ \ \ \text{if}\ \ \ x_{{{i}_{1}}}^{\text{ex}}\le {{x}_{i}}\le x_{{{i}_{2}}}^{\text{ex}} \\ & a_{{{i}_{3}}}^{\text{ex}}+b_{{{i}_{3}}}^{\text{ex}}{{x}_{i}}, \ \ \text{if}\ \ \ \ x_{{{i}_{2}}}^{\text{ex}}\le {{x}_{i}}\le {{V}_{{{T}_{i}}}} \\ \end{align} \right. $ |
下面對非線性順應性函數通過sigmoidal函數進行平滑處理。對于呼氣狀態方程,ciin(xi)可以近似為
$ c_{i}^{\text{in}}\left({{x}_{i}} \right)\approx a_{{{i}_{2}}}^{\text{in}}\left(s_{a, b}^{\left(\beta \right)}\left({{x}_{i}} \right)-s_{c, d}^{\left(\beta \right)}\left({{x}_{i}} \right)\right), i=1, 2, \cdots, {{2}^{n}}, $ |
這里a=-ai1in/bi1in, b=(ai2in+bi1in)+a, c=-ai3in/bi3in, d=ai2in/bi3in+c, sa, b(β)=1/(b-a)ln(σb-β(xi))1/β, σb-β=1/(1+e-β(xi-a)),β>0是一個近似參數。圖 3(a)表示平滑逼近分段線性吸氣階段順應性函數ciin(xi)。圖 3(b)表示相應的呼氣階段順應性函數ciex(xi)。其中虛線表示吸氣和呼氣順應性函數與肺容量的關系,實線表示經過平滑處理后的效果。

(a)吸氣過程;(b)呼氣過程
Figure3. Smoothed compliance function of inspiration and expiration(a) inspiration; (b) expiration
最后,采用向量矩陣靜態空間形式改寫式(3)和式(6)。定義靜態向量x=[x1, x2, …, x2n]T, xi表示第i個肺子室的容量。式(3)可改寫為
$ \begin{align} & {{R}_{\text{in}}}{{x}^{'}}_{t}\left(t \right)+{{c}_{\text{in}}}\left(x\left(t \right)\right)x\left(t \right)={{p}_{\text{in}}}\left(t \right)e, \\ & x\left(0 \right)=x_{0}^{\text{in}}, 0\le t\le {{T}_{\text{in}}}, \\ \end{align} $ |
這里e=[1, …, 1]T表示向量2n, cin(x)表示對角矩陣函數,cin(x)和Rin定義如下:
$ {{c}_{\text{in}}}\left(x \right)=\text{diag}\left[\frac{1}{c_{1}^{\text{in}}\left({{x}_{1}} \right)}, \cdots, \frac{1}{c_{{{2}^{n}}}^{\text{in}}\left({{x}_{{{2}^{n}}}} \right)} \right], $ |
$ {{R}_{\text{in}}}=\sum\limits_{j=0}^{n}{\sum\limits_{k=1}^{{{2}^{j}}}{R_{j, k}^{\text{in}}}}{{Z}_{j, k}}Z_{j, k}^{T} $ |
其中Zj, k∈2n, 第l個Zj, k為1,而其它則為0,l=(k-1)2n-j+1, (k-1)2n-j+2, …, k2n-j, k=1, …, 2j, j=0, 1, …, n。
類似地式(6)也可改寫為
$ \begin{align} & {{R}_{\text{ex}}}{{x}^{'}}_{t}\left(t \right)+{{c}_{\text{ex}}}\left(x\left(t \right)\right)x\left(t \right)={{p}_{\text{ex}}}\left(t \right)e, \\ & x\left({{T}_{\text{in}}} \right)=x_{0}^{\text{ex}}, {{T}_{\text{in}}}\le t\le {{T}_{\text{in}}}+{{T}_{\text{ex}}}, \\ \end{align} $ |
這里cex(x)和Rex定義如下:
$ {{c}_{\text{ex}}}\left(x \right)=\text{diag}\left[\frac{1}{c_{1}^{\text{ex}}\left({{x}_{1}} \right)}, \cdots, \frac{1}{c_{{{2}^{n}}}^{\text{ex}}\left({{x}_{{{2}^{n}}}} \right)} \right], $ |
$ {{R}_{\text{ex}}}=\sum\limits_{j=0}^{n}{\sum\limits_{k=1}^{{{2}^{j}}}{R_{j, k}^{\text{ex}}}}{{Z}_{j, k}}Z_{j, k}^{T} $ |
因為Rin和Rex都是正定的,所以Rin和Rex都是可以轉置的矩陣。
2 優化確定吸氣和呼氣的氣流
通過第1節建立的模型,其中式(9)和式(12)描述了呼吸周期性動態系統,下面將對這個模型優化預測呼吸氣流模式。優化標準是以最小化呼吸肌做工(氧消耗)以及肺容量的變化率為目標。呼吸過程的氧消耗主要是由于呼吸肌需要克服阻力以及伸展肺和胸壁。一般呼吸系統力學模型是通過定義PV曲線來描述這個過程,其中P是肺內氣體壓力,V是肺容積[11-15]。在肺部,氣體交換效率跟氣體容積變化率有關,因為快速變化的肺容量可能引發不適甚至無效的肌肉控制收縮[16-18]。此外,氣體容積過高可能會導致肺過度膨脹以致損傷肺部組織,過度的通氣也使呼吸肌肉容易產生疲勞[19-22]。
2.1 吸氣階段
假設呼吸過程是從一個給定的初始狀態x0in,然后伴隨著以吸氣末的狀態作為呼氣的初始狀態。一個吸氣過程后接著呼氣過程稱為一個呼吸周期。此外,假設每個呼吸周期后跟另一個呼吸周期,前面呼吸周期結束狀況是后者的初始狀態。由于呼吸過程的周期性,所以就只需要專注研究一個呼吸周期。下面是通過對吸氣模式采用微積分方程的方法給出的最優解x*(t), 0≤t≤Tin。
根據式(9)得到目標函數,函數前項表示了肺容積變化率的加速度,后項表示了PV曲線選擇最佳呼吸潮氣量。
$ \begin{align} & \text{Min}{{\psi }_{\text{in}}}\left(x \right)=\int\limits_{0}^{{{T}_{\text{in}}}}{\left[{{x}^{''T}}\left(t \right){{x}^{''}}\left(t \right)+{{a}_{1}}{{p}_{\text{in}}}\left(t \right){{e}^{T}}{{x}^{'}}\left(t \right)\right]}\text{d}t, {{a}_{1}} \\ & \ge 0 \\ & \text{s}\text{.t}\text{.}\ \ {{x}_{0}}={{V}_{0}}, {{x}^{'}}\left(0 \right)=0, x\left({{T}_{\text{in}}} \right)={{V}_{0}}+{{V}_{T}}, {{x}^{'}}\left({{T}_{\text{in}}} \right)=0 \\ \end{align} $ |
其中V0∈2n是呼氣末的容積,VT∈2n是潮氣量。如果a1>0,那么有
$ \begin{align} & {{x}^{*}}\left(t \right)={{d}_{1}}+{{d}_{2}}t+\exp \left(\sqrt{{{a}_{1}}}R_{\text{in}}^{1/2}t \right){{d}_{3}}+\\ & \ \ \ \ \ \ \ \ \ \ \ \ \ \ \exp \left(-\sqrt{{{a}_{1}}}R_{\text{in}}^{1/2}t \right){{d}_{4}}, t\ge 0 \\ \end{align} $ |
如果a1=0,那么有
$ {{x}^{*}}\left(t \right)={{d}_{1}}+{{d}_{2}}t+{{d}_{3}}{{t}^{2}}+{{d}_{4}}{{t}^{4}}, t\ge 0, $ |
其中Rin1/2表示Rin的正定平方根,d1、d2、d3和d4∈2n是由常向量的邊界條件決定。如果a1=0,可得d1=V0,d2=0,d3=(3/Tin2)VT及d4=-(2/Tin3)VT。因此,x′(t)=d2+2d3t+3d4t2=(6t/(Tin2))VT(1-(t/Tin))≥≥0,0≤t≤Tin, 蘊含著有解V0i≤xi*(t)≤V0i+VTi, i=1, 2, …, 2n。當a1>0時,可以得到類似的結果。
由第1節肺模型及優化算法可以得到呼吸氣流的最佳模式和相應的驅動壓力。這些模式的輸入信號可以作為最佳呼吸氣流模式的呼吸機驅動壓力參數的輸入。吸氣過程Ψin(x)是關于肺容量變化率和呼吸肌克服阻力的加權平方和。
2.2 呼氣階段
對于呼氣模式,根據式(12)給出的最優解x*(t)(Tin≤t≤Tin+Tex),得到目標函數:
$ \begin{align} & \text{Min}{{\psi }_{\text{ex}}}\left(x \right)=\int\limits_{{{T}_{_{\text{in}}}}}^{{{T}_{\text{in}}}+{{T}_{\text{ex}}}}{\left[{{x}^{''T}}\left(t \right){{x}^{''}}\left(t \right)+{{a}_{2}}p_{\text{ex}}^{2}\left(t \right){{e}^{T}}x \right]}\text{d}t, \\ & \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {{a}_{2}}\ge 0 \\ & \text{s}\text{.t}\text{.}\ \ x\left({{T}_{\text{in}}} \right)={{V}_{0}}+{{V}_{T}}, {{x}^{'}}\left({{T}_{\text{in}}} \right)=0, x\left({{T}_{\text{in}}}+{{T}_{\text{ex}}} \right)={{V}_{0}}, {{x}^{'}}\left({{T}_{\text{in}}}+{{T}_{\text{ex}}} \right)=0 \\ \end{align} $ |
if a2>0,那么最優解滿足:
$ \begin{align} & {{x}^{\left(4 \right)}}\left(t \right)-{{a}_{2}}R_{\text{ex}}^{2}{{x}^{\left(2 \right)}}\left(t \right)+{{a}_{2}}C_{\text{ex}}^{2}\left(x \right)x\left(t \right)+{{a}_{2}}\left[{{C}_{\text{ex}}}\left(x \right)\right.\\ & {{R}_{\text{ex}}}{{x}^{'}}\left(t \right)-{{R}_{\text{ex}}}{{C}_{\text{ex}}}\left(x \right){{x}^{'}}\left(t \right)+X\left(t \right){{C}^{'}}_{\text{ex}}\left(x \right){{R}_{\text{wex}}}{{x}^{'}}\left(t \right)-\\ & \left.{{R}_{\text{ex}}}{{C}^{'}}_{\text{ex}}\left(x \right)X\left(t \right){{x}^{'}}\left(t \right)+X\left(t \right){{C}^{'}}_{\text{ex}}\left(x \right){{C}_{\text{ex}}}\left(x \right)x\left(t \right)\right]=0 \\ \end{align} $ |
其中X(t)=diag[xi(t)]和C′ex(x)=diag[], i=1, …, 2n。如果a2=0,則x*(t)=d1+d2t+d3t2+d4t3(t≥0), d1、d2、d3和d4∈2n是常向量由邊界條件決定。
3 實驗結果
為得到優化肺容量曲線,可以將模型轉化為求解兩點非線性邊界值問題。求解此類問題的數值方法常見有打靶法和梯度下降法。這里使用梯度下降法在MATLAB下,利用求解數值微積分方程式(15)和(16)得到x*(t)(t≥0)優化肺容量曲線。
在模擬實驗中,首先考慮雙室肺模型,設置xi1in=0.3 L,xi2in=0.48 L,xi1ex=0.23 L,xi1ex=0.43 L,i=1, 2。其它參數設置如表 1所示。

假設支氣管樹是二分叉結構。氣道阻力隨著分支級別而變化,典型的數據可參見文獻[23]。此外,呼氣阻力比吸氣阻力高一個影響因素2~3。在這里,假設系數為2.5。又假設Tin為2 s、Tex為3 s,因此兩個權重參數a1和a2是不同的。在自主呼吸時,一般權重參數設置為a1=2.0 L/s3 cm H2O和a2=0.1 L2/s4 cm H2O[24]。
圖 4說明了當a1和a2改變時容積和氣流速率的優化敏感性。從該圖中可以看出,吸氣氣流速率是對稱的,吸氣氣流速率下降的最大值是關于a1增加的函數。此外,呼氣氣流速率的非對稱說明了當a2增加時最低值下降得更快。具體來說,如果設置的權重參數a2=0,得到的呼氣氣流變化曲線是一個拋物線弧。圖 4的氣流模式顯示了在自主呼吸時候的典型特性,也就是吸氣氣流速率曲線相對平坦和呼氣氣流速率波形有個初始槽是不對稱的,非常接近真實的氣流速率。

機械通氣是利用在氣管口施加周期性變化的正氣壓來實現氣體進、出肺部。患者在完全沒有主動呼吸的情況下,通常我們可以觀察到肺容量的吸氣末和呼氣末迅速收斂到穩定值。但是這并不能夠簡單的把肺看作是一個穩定的動態系統。從肺解剖結構可以看到它是一個樹形結構,并且不斷重復分支成小的氣道,最終形成一個個肺泡的氣體交換單位[25]。這里通過設計多室肺的數學模型來展示肺的呼吸過程。通過非線性函數擴展后可以實現比普通的線性二分叉結構更復雜的模型。
4 結語
本文研究了一個非線性多室模型模擬肺力學系統的最佳呼吸氣流模式。考慮肺呼吸肌氧消耗、肺容量和肺彈性勢能等條件,使用經典微積分推導出最佳氣體容量的優化曲線和優化準則。本文的結果可以作為開發多變量控制器機械通氣監護重癥患者的基礎。下一步的工作是設計多變量和自適應控制模型,并在這些模型的架構下實現完全自動化的機械通氣重癥監護,確保患者的安全通氣。
引言
呼吸衰竭通氣不足是危重病急救醫學臨床上常見的問題。一旦確診通氣不足原因后,呼吸衰竭患者通常需要機械通氣支持。但是臨床上在采用機械通氣過程中,由于過量容積或壓力過高,很容易引起急性肺損傷。通過控制算法可以實現部分自動控制容量和壓力。這些算法需要一個參考模型用于辨別出臨床上合理的呼吸模式。然而,目前文獻[1-2]已經提出來的肺呼吸模型通常假設肺的同質性。早期呼吸控制優化研究都是使用簡單的同質肺模型處理呼吸的頻率問題。Otis等[3]通過使用最小工作率準則預測呼吸頻率。這是涉及靜態優化問題,假定氣流模式是一個固定的正弦函數。Yamashiro等[4]研究了采用固定的吸氣和呼氣相呼吸周期預測呼吸氣流模式。這些研究結果被Hamalainen等[5]進行擴展,考慮了兩級分層模型來控制呼吸,上層模型負責對下層模型的氣體模型總體控制變量進行優化,下層模型用最小化高層模型的標準來選擇呼吸參數的方法確定氣流模型。
雖然尋找最佳的呼吸模式問題已經被Sundaresan等[6]提出來,但是在這些研究中呼吸模型都是以固定呼吸參數及單個腔體肺模型來預測控制。然而,人體的肺,特別是患病的肺是異構的。從解剖角度看肺包括許多的肺泡子結構腔室,在這些子結構里面的氣體交換的能力是不同的。仿真的肺模型應該考慮這種肺異質性。此外,氣流的阻力和肺的順應性也不是恒定的,隨著肺容量而變化,肺的順應性也跟肺容量相關。實現機械通氣的自適應控制,采用比單室肺的容積或壓力控制通氣更詳細的控制結構模型是必要的,尤其是要考慮非線性氣體阻力、順應性和肺異質性。氣道中的流量和壓力并不是簡單的波形圖,雖然目前臨床上普遍是這樣描述,但是事實上在這些波形中嵌入了大量的信息。本文將提供一個嚴格的數學框架用于表達非線性多室肺力學模型系統并優化呼吸氣流模式。
本文在陳夫華等[7]前期研究結果基礎上,用非線性多室模型優化肺力學系統呼吸氣流模式。這里的優化是指微積分方程的最優解,而不是在生理學意義上的最優呼吸模式。首先擴展Chellaboina等[8]中的線性多室肺癌模型成非線性模型。然后擴展Hamalainen等[9]吸氣和呼氣的呼吸循環功能函數,采用典型微積分方法推導最佳氣流模式。包括了生理學上的最小化呼吸工作和幫助患者吸氣增加肺容積,最小化彈性勢能和呼氣階段改變氣流速度變化等方面。最后,用求非線性兩點邊界值解確定吸氣和呼氣循環過程的最佳氣流模式。
1 非線性多室肺模型
1.1 單室肺模型
為了使肺模型更加符合實際肺組織結構和工作過程,這里擴展線性多室肺模型成非線性模型并模擬整個吸氣過程,同時假設肺的順應性是關于肺容量的非線性函數[10]。
單室肺模型是通過一個單室腔體表示肺生理解剖結構,如圖 1所示。用非線性函數 C(x)表示單室肺的順應性,R表示氣流的阻力。當時間t=0時,驅動壓力Pin(t)表示加載在父層的氣道開口,Pin(t)可能是由呼吸肌肉或機械通氣作用下產生的。該壓力發生在0≤t≤Tin之間,也就是呼吸周期的吸氣階段,Tin是吸氣時間。當時間t=Tin時,該壓力將被釋放,緊接著被呼氣壓力替代,此時外部壓力等于大氣壓。在Tin≤t≤Tin+Tex時,氣體壓力為Pex(t),Tex是呼氣階段持續的時間。

吸氣階段的狀態方程表示為
$ \begin{align} & {{R}_{\text{in}}}{{x}^{'}}\left(t \right)+\left(1/-{{C}_{\text{in}}}\left(x \right)\right)x\left(t \right)-{{P}_{\text{in}}}\left(t \right), \\ & x\left(0 \right)=x_{0}^{\operatorname{in}}, 0\le t\le {{T}_{\text{in}}}, \\ \end{align} $ |
其中x(t)∈,t≥0,表示肺容量。Rin∈,表示吸氣階段氣體受到的阻力。Cin:→+是一個非線性函數,表示在吸氣階段肺的順應性。x0in∈ +,表示吸氣開始時候的肺容量作為系統初始條件。式(1)實際是一個簡單的壓力平衡方程,應用在肺室模型時候,順應性和氣道阻力與肺容量成比例。這里假設呼氣動作是由于肺室彈性舒展的被動結果。
呼氣階段的狀態方程表示為
$ \begin{align} & {{R}_{\text{ex}}}{{x}^{'}}\left(t \right)+\left(1/-{{C}_{\text{ex}}}\left(x \right)\right)x\left(t \right)-{{P}_{\text{ex}}}\left(t \right), \\ & x\left({{T}_{\text{in}}} \right)=x_{\text{0}}^{\text{ex}}, {{T}_{\text{in}}}\le t\le {{T}_{\text{in}}}+{{T}_{\text{ex}}}, \\ \end{align} $ |
其中x(t)的意義跟吸氣狀態方程一樣,Rex∈ +,表示呼氣階段氣體受到的阻力。Cex: → +是一個非線性函數,表示在呼氣階段肺的順應性。x0ex∈ +,表示呼氣開始時候的肺容量。
1.2 多室肺模型
為了更加接近真實的肺結構,模型中假設支氣管樹成二分叉結構,也就是在每一級氣道單元重復分叉成兩個子的氣道單元,圖 2表示了4個子室的模型肺。這里定義2n(n≥0)室肺模型,由2n個小單元肺室組成,這些肺單元由n層壓力源連接到氣道單元。因為每個氣道分叉成兩個后繼級子氣道,所以整個肺部是由2n個子室組成。

假設xi(i=1, 2, …, 2n)表示第i個肺子室的容量,ciin(i=1, 2, …, 2n)表示對應于第i個肺子室在吸氣階段的順應性。Rj, iin(i=1, 2, …, 2j)表示多肺室模型位于第j層壓力源的第i個肺子室在吸氣階段的氣體阻力。
在單室模型中,假設氣體壓力Pin(t),t≥0,表示吸氣階段由吸氣肌肉或機械通氣作用下產生的壓力。現在可將吸氣階段的狀態方程改為
$ \begin{align} & R_{n}^{\text{in}}{{x}^{'}}_{i}\left(t \right)+\frac{1}{c_{i}^{\text{in}}\left({{x}_{i}}\left(t \right)\right)}{{x}_{i}}\left(t \right)+\sum\limits_{j=0}^{n-1}{R_{j, {{k}_{j}}}^{\text{in}}}\sum\limits_{l=\left({{k}_{j}}-1 \right){{2}^{n-j}}+1}^{{{k}_{j}}{{2}^{n-j}}}{{{x}^{'}}_{l}\left(t \right)={{P}_{\text{in}}}\left(t \right)}, \\ & {{x}_{i}}\left(0 \right)=x_{{{i}_{o}}}^{\text{in}}, 0\le t\le {{T}_{\text{in}}}, i=1, 2, \cdots, {{2}^{n}}, \\ \end{align} $ |
其中ciin(xi),i=1, 2, …,2n是關于xi(i=1, 2, …,2n)的非線性方程。
$ c_{i}^{\text{in}}\left({{x}_{i}} \right)=\left\{ \begin{align} & a_{{{i}_{1}}}^{\text{in}}+b_{{{i}_{1}}}^{\text{in}}{{x}_{i}}, \ \ \ \ \text{if}\ \ \ \text{0}\le {{x}_{i}}\le x_{{{i}_{1}}}^{\text{in}} \\ & a_{{{i}_{2}}}^{\text{in}}, \ \ \ \ \ \ \ \ \ \ \ \text{if}\ \ \ x_{{{i}_{1}}}^{\text{in}}\le {{x}_{i}}\le x_{{{i}_{2}}}^{\text{in}} \\ & a_{{{i}_{3}}}^{\text{in}}+b_{{{i}_{3}}}^{\text{in}}{{x}_{i}}, \ \ \text{if}\ \ \ \ x_{{{i}_{2}}}^{\text{in}}\le {{x}_{i}}\le {{V}_{{{T}_{i}}}} \\ \end{align} \right., $ |
這里aimin(m=1, 2, 3)且bimin(m=1, 3)是模型參數,bi1in>0且bi3in<0, ximin(m=1, 2), 表示容積變化在xi1in和xi2in范圍內時,順應性為固定常數ai2in,VTi表示潮氣量。
$ {{k}_{j}}=\left[\frac{{{k}_{j+1}}-1}{2} \right]+1, j=0, \cdots, n-1, {{k}_{n}}=i, $ |
其中表示一個不大于k的整數,kj表示第i個子室氣體經過第j級氣道編號。
對于呼氣狀態方程,如在單室模型中假設呼氣過程是被動及外部施加壓力的結果,pex(t),t≥0。類似吸氣過程的處理方法,可以得到呼氣階段的狀態方程:
$ R_{n}^{\text{ex}}{{x}^{'}}_{i}\left(t \right)+\sum\limits_{j=0}^{n-1}{R_{j, {{k}_{j}}}^{\text{ex}}}\sum\limits_{l=\left({{k}_{j}}-1 \right){{2}^{n-j}}+1}^{{{k}_{j}}{{2}^{n-j}}}{{{x}^{'}}_{l}\left(t \right)+}\frac{1}{c_{i}^{\text{ex}}\left({{x}_{i}}\left(t \right)\right)}{{x}_{i}}\left(t \right)={{P}_{\text{ex}}}\left(t \right), {{x}_{i}}\left({{T}_{\text{in}}} \right)=x_{{{i}_{0}}}^{\text{ex}}, {{T}_{\text{in}}}\le t\le {{T}_{\text{in}}}+{{T}_{\text{ex}}}, i=1, 2, \cdots, {{2}^{n}}, $ |
其中ciex(xi)由如下定義:
$ c_{i}^{\text{ex}}\left({{x}_{i}} \right)=\left\{ \begin{align} & a_{{{i}_{0}}}^{\text{ex}}+b_{{{i}_{1}}}^{\text{ex}}{{x}_{i}}, \ \ \ \ \text{if}\ \ \ \text{0}\le {{x}_{i}}\le x_{{{i}_{1}}}^{\text{ex}} \\ & a_{{{i}_{2}}}^{\text{ex}}, \ \ \ \ \ \ \ \ \ \ \ \text{if}\ \ \ x_{{{i}_{1}}}^{\text{ex}}\le {{x}_{i}}\le x_{{{i}_{2}}}^{\text{ex}} \\ & a_{{{i}_{3}}}^{\text{ex}}+b_{{{i}_{3}}}^{\text{ex}}{{x}_{i}}, \ \ \text{if}\ \ \ \ x_{{{i}_{2}}}^{\text{ex}}\le {{x}_{i}}\le {{V}_{{{T}_{i}}}} \\ \end{align} \right. $ |
下面對非線性順應性函數通過sigmoidal函數進行平滑處理。對于呼氣狀態方程,ciin(xi)可以近似為
$ c_{i}^{\text{in}}\left({{x}_{i}} \right)\approx a_{{{i}_{2}}}^{\text{in}}\left(s_{a, b}^{\left(\beta \right)}\left({{x}_{i}} \right)-s_{c, d}^{\left(\beta \right)}\left({{x}_{i}} \right)\right), i=1, 2, \cdots, {{2}^{n}}, $ |
這里a=-ai1in/bi1in, b=(ai2in+bi1in)+a, c=-ai3in/bi3in, d=ai2in/bi3in+c, sa, b(β)=1/(b-a)ln(σb-β(xi))1/β, σb-β=1/(1+e-β(xi-a)),β>0是一個近似參數。圖 3(a)表示平滑逼近分段線性吸氣階段順應性函數ciin(xi)。圖 3(b)表示相應的呼氣階段順應性函數ciex(xi)。其中虛線表示吸氣和呼氣順應性函數與肺容量的關系,實線表示經過平滑處理后的效果。

(a)吸氣過程;(b)呼氣過程
Figure3. Smoothed compliance function of inspiration and expiration(a) inspiration; (b) expiration
最后,采用向量矩陣靜態空間形式改寫式(3)和式(6)。定義靜態向量x=[x1, x2, …, x2n]T, xi表示第i個肺子室的容量。式(3)可改寫為
$ \begin{align} & {{R}_{\text{in}}}{{x}^{'}}_{t}\left(t \right)+{{c}_{\text{in}}}\left(x\left(t \right)\right)x\left(t \right)={{p}_{\text{in}}}\left(t \right)e, \\ & x\left(0 \right)=x_{0}^{\text{in}}, 0\le t\le {{T}_{\text{in}}}, \\ \end{align} $ |
這里e=[1, …, 1]T表示向量2n, cin(x)表示對角矩陣函數,cin(x)和Rin定義如下:
$ {{c}_{\text{in}}}\left(x \right)=\text{diag}\left[\frac{1}{c_{1}^{\text{in}}\left({{x}_{1}} \right)}, \cdots, \frac{1}{c_{{{2}^{n}}}^{\text{in}}\left({{x}_{{{2}^{n}}}} \right)} \right], $ |
$ {{R}_{\text{in}}}=\sum\limits_{j=0}^{n}{\sum\limits_{k=1}^{{{2}^{j}}}{R_{j, k}^{\text{in}}}}{{Z}_{j, k}}Z_{j, k}^{T} $ |
其中Zj, k∈2n, 第l個Zj, k為1,而其它則為0,l=(k-1)2n-j+1, (k-1)2n-j+2, …, k2n-j, k=1, …, 2j, j=0, 1, …, n。
類似地式(6)也可改寫為
$ \begin{align} & {{R}_{\text{ex}}}{{x}^{'}}_{t}\left(t \right)+{{c}_{\text{ex}}}\left(x\left(t \right)\right)x\left(t \right)={{p}_{\text{ex}}}\left(t \right)e, \\ & x\left({{T}_{\text{in}}} \right)=x_{0}^{\text{ex}}, {{T}_{\text{in}}}\le t\le {{T}_{\text{in}}}+{{T}_{\text{ex}}}, \\ \end{align} $ |
這里cex(x)和Rex定義如下:
$ {{c}_{\text{ex}}}\left(x \right)=\text{diag}\left[\frac{1}{c_{1}^{\text{ex}}\left({{x}_{1}} \right)}, \cdots, \frac{1}{c_{{{2}^{n}}}^{\text{ex}}\left({{x}_{{{2}^{n}}}} \right)} \right], $ |
$ {{R}_{\text{ex}}}=\sum\limits_{j=0}^{n}{\sum\limits_{k=1}^{{{2}^{j}}}{R_{j, k}^{\text{ex}}}}{{Z}_{j, k}}Z_{j, k}^{T} $ |
因為Rin和Rex都是正定的,所以Rin和Rex都是可以轉置的矩陣。
2 優化確定吸氣和呼氣的氣流
通過第1節建立的模型,其中式(9)和式(12)描述了呼吸周期性動態系統,下面將對這個模型優化預測呼吸氣流模式。優化標準是以最小化呼吸肌做工(氧消耗)以及肺容量的變化率為目標。呼吸過程的氧消耗主要是由于呼吸肌需要克服阻力以及伸展肺和胸壁。一般呼吸系統力學模型是通過定義PV曲線來描述這個過程,其中P是肺內氣體壓力,V是肺容積[11-15]。在肺部,氣體交換效率跟氣體容積變化率有關,因為快速變化的肺容量可能引發不適甚至無效的肌肉控制收縮[16-18]。此外,氣體容積過高可能會導致肺過度膨脹以致損傷肺部組織,過度的通氣也使呼吸肌肉容易產生疲勞[19-22]。
2.1 吸氣階段
假設呼吸過程是從一個給定的初始狀態x0in,然后伴隨著以吸氣末的狀態作為呼氣的初始狀態。一個吸氣過程后接著呼氣過程稱為一個呼吸周期。此外,假設每個呼吸周期后跟另一個呼吸周期,前面呼吸周期結束狀況是后者的初始狀態。由于呼吸過程的周期性,所以就只需要專注研究一個呼吸周期。下面是通過對吸氣模式采用微積分方程的方法給出的最優解x*(t), 0≤t≤Tin。
根據式(9)得到目標函數,函數前項表示了肺容積變化率的加速度,后項表示了PV曲線選擇最佳呼吸潮氣量。
$ \begin{align} & \text{Min}{{\psi }_{\text{in}}}\left(x \right)=\int\limits_{0}^{{{T}_{\text{in}}}}{\left[{{x}^{''T}}\left(t \right){{x}^{''}}\left(t \right)+{{a}_{1}}{{p}_{\text{in}}}\left(t \right){{e}^{T}}{{x}^{'}}\left(t \right)\right]}\text{d}t, {{a}_{1}} \\ & \ge 0 \\ & \text{s}\text{.t}\text{.}\ \ {{x}_{0}}={{V}_{0}}, {{x}^{'}}\left(0 \right)=0, x\left({{T}_{\text{in}}} \right)={{V}_{0}}+{{V}_{T}}, {{x}^{'}}\left({{T}_{\text{in}}} \right)=0 \\ \end{align} $ |
其中V0∈2n是呼氣末的容積,VT∈2n是潮氣量。如果a1>0,那么有
$ \begin{align} & {{x}^{*}}\left(t \right)={{d}_{1}}+{{d}_{2}}t+\exp \left(\sqrt{{{a}_{1}}}R_{\text{in}}^{1/2}t \right){{d}_{3}}+\\ & \ \ \ \ \ \ \ \ \ \ \ \ \ \ \exp \left(-\sqrt{{{a}_{1}}}R_{\text{in}}^{1/2}t \right){{d}_{4}}, t\ge 0 \\ \end{align} $ |
如果a1=0,那么有
$ {{x}^{*}}\left(t \right)={{d}_{1}}+{{d}_{2}}t+{{d}_{3}}{{t}^{2}}+{{d}_{4}}{{t}^{4}}, t\ge 0, $ |
其中Rin1/2表示Rin的正定平方根,d1、d2、d3和d4∈2n是由常向量的邊界條件決定。如果a1=0,可得d1=V0,d2=0,d3=(3/Tin2)VT及d4=-(2/Tin3)VT。因此,x′(t)=d2+2d3t+3d4t2=(6t/(Tin2))VT(1-(t/Tin))≥≥0,0≤t≤Tin, 蘊含著有解V0i≤xi*(t)≤V0i+VTi, i=1, 2, …, 2n。當a1>0時,可以得到類似的結果。
由第1節肺模型及優化算法可以得到呼吸氣流的最佳模式和相應的驅動壓力。這些模式的輸入信號可以作為最佳呼吸氣流模式的呼吸機驅動壓力參數的輸入。吸氣過程Ψin(x)是關于肺容量變化率和呼吸肌克服阻力的加權平方和。
2.2 呼氣階段
對于呼氣模式,根據式(12)給出的最優解x*(t)(Tin≤t≤Tin+Tex),得到目標函數:
$ \begin{align} & \text{Min}{{\psi }_{\text{ex}}}\left(x \right)=\int\limits_{{{T}_{_{\text{in}}}}}^{{{T}_{\text{in}}}+{{T}_{\text{ex}}}}{\left[{{x}^{''T}}\left(t \right){{x}^{''}}\left(t \right)+{{a}_{2}}p_{\text{ex}}^{2}\left(t \right){{e}^{T}}x \right]}\text{d}t, \\ & \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {{a}_{2}}\ge 0 \\ & \text{s}\text{.t}\text{.}\ \ x\left({{T}_{\text{in}}} \right)={{V}_{0}}+{{V}_{T}}, {{x}^{'}}\left({{T}_{\text{in}}} \right)=0, x\left({{T}_{\text{in}}}+{{T}_{\text{ex}}} \right)={{V}_{0}}, {{x}^{'}}\left({{T}_{\text{in}}}+{{T}_{\text{ex}}} \right)=0 \\ \end{align} $ |
if a2>0,那么最優解滿足:
$ \begin{align} & {{x}^{\left(4 \right)}}\left(t \right)-{{a}_{2}}R_{\text{ex}}^{2}{{x}^{\left(2 \right)}}\left(t \right)+{{a}_{2}}C_{\text{ex}}^{2}\left(x \right)x\left(t \right)+{{a}_{2}}\left[{{C}_{\text{ex}}}\left(x \right)\right.\\ & {{R}_{\text{ex}}}{{x}^{'}}\left(t \right)-{{R}_{\text{ex}}}{{C}_{\text{ex}}}\left(x \right){{x}^{'}}\left(t \right)+X\left(t \right){{C}^{'}}_{\text{ex}}\left(x \right){{R}_{\text{wex}}}{{x}^{'}}\left(t \right)-\\ & \left.{{R}_{\text{ex}}}{{C}^{'}}_{\text{ex}}\left(x \right)X\left(t \right){{x}^{'}}\left(t \right)+X\left(t \right){{C}^{'}}_{\text{ex}}\left(x \right){{C}_{\text{ex}}}\left(x \right)x\left(t \right)\right]=0 \\ \end{align} $ |
其中X(t)=diag[xi(t)]和C′ex(x)=diag[], i=1, …, 2n。如果a2=0,則x*(t)=d1+d2t+d3t2+d4t3(t≥0), d1、d2、d3和d4∈2n是常向量由邊界條件決定。
3 實驗結果
為得到優化肺容量曲線,可以將模型轉化為求解兩點非線性邊界值問題。求解此類問題的數值方法常見有打靶法和梯度下降法。這里使用梯度下降法在MATLAB下,利用求解數值微積分方程式(15)和(16)得到x*(t)(t≥0)優化肺容量曲線。
在模擬實驗中,首先考慮雙室肺模型,設置xi1in=0.3 L,xi2in=0.48 L,xi1ex=0.23 L,xi1ex=0.43 L,i=1, 2。其它參數設置如表 1所示。

假設支氣管樹是二分叉結構。氣道阻力隨著分支級別而變化,典型的數據可參見文獻[23]。此外,呼氣阻力比吸氣阻力高一個影響因素2~3。在這里,假設系數為2.5。又假設Tin為2 s、Tex為3 s,因此兩個權重參數a1和a2是不同的。在自主呼吸時,一般權重參數設置為a1=2.0 L/s3 cm H2O和a2=0.1 L2/s4 cm H2O[24]。
圖 4說明了當a1和a2改變時容積和氣流速率的優化敏感性。從該圖中可以看出,吸氣氣流速率是對稱的,吸氣氣流速率下降的最大值是關于a1增加的函數。此外,呼氣氣流速率的非對稱說明了當a2增加時最低值下降得更快。具體來說,如果設置的權重參數a2=0,得到的呼氣氣流變化曲線是一個拋物線弧。圖 4的氣流模式顯示了在自主呼吸時候的典型特性,也就是吸氣氣流速率曲線相對平坦和呼氣氣流速率波形有個初始槽是不對稱的,非常接近真實的氣流速率。

機械通氣是利用在氣管口施加周期性變化的正氣壓來實現氣體進、出肺部。患者在完全沒有主動呼吸的情況下,通常我們可以觀察到肺容量的吸氣末和呼氣末迅速收斂到穩定值。但是這并不能夠簡單的把肺看作是一個穩定的動態系統。從肺解剖結構可以看到它是一個樹形結構,并且不斷重復分支成小的氣道,最終形成一個個肺泡的氣體交換單位[25]。這里通過設計多室肺的數學模型來展示肺的呼吸過程。通過非線性函數擴展后可以實現比普通的線性二分叉結構更復雜的模型。
4 結語
本文研究了一個非線性多室模型模擬肺力學系統的最佳呼吸氣流模式。考慮肺呼吸肌氧消耗、肺容量和肺彈性勢能等條件,使用經典微積分推導出最佳氣體容量的優化曲線和優化準則。本文的結果可以作為開發多變量控制器機械通氣監護重癥患者的基礎。下一步的工作是設計多變量和自適應控制模型,并在這些模型的架構下實現完全自動化的機械通氣重癥監護,確保患者的安全通氣。