心血管系统的仿真与建模

心血管系统的仿真与建模1．前言心血管循环系统是最复杂的生命系统之一，具有高度的动态特性。当心脏周期性地收缩和舒张时，心室射入主动脉的血流将以波的形式自主动脉根部出发沿动脉管系传播，这种波就是脉搏波。人体心血管系统由心脏、动脉、毛细血管和静脉组成，构成体循环和肺循环两大回路，其中体循环在人体的血液循环系统中占据极其重要的地位。体循环开始于左心室。血液从左心室搏出后，流经主动脉及其派生的若干分支动脉，进入各组织器官。动脉血经组织器官内的毛细血管完成氧气和营养物质的交换后变为静脉血，再由各级静脉汇集到上腔静脉和下腔静脉，回到左心房，从而完成了整个体循环过程。描述心血管系统功能状态的主要生理参数有[1]：血压(收缩压、舒张压、脉压、平均动脉压)、心率、每搏输出量、血管顺应性、血流阻力、心肌收缩能力等，它们与人体状态有着密切的联系：生理参数反映了人体的生理病理状态，而当人体生理病理状态发生改变时，相应的生理参数也会随之改变。本文的主要内容如下：1.单弹性腔仿真模型在Windkessel理论的基础上，建立心血管系统一阶Simulink仿真模型，结合临床数据计算模型参数进行仿真，分析实验结果与实测生理量的符合程度，从而验证模型对于心血管系统的表征能力。2.双弹性腔仿真模型在一阶模型的基础上，引入频率元件L，细化血管顺应性，建立三阶Simlink仿真模型。3.参数性质分析研究模型参数改变时脉搏波的变化趋势，分析各参数对脉搏波的作用。2.单弹性腔基本理论与仿真模型2.1Windkessel模型简介对于图2.1所示的往复泵供水系统，柱塞P在马达M的驱动下做往复运动：当P向前挤压时，吸水阀门2关闭、供水阀门1开启，被挤压的流体经传输管路A，并通过终端阻力R最后流入贮水槽V；当P向后抽吸时，供水阀门1关闭、吸水阀门2开启，贮水槽中的水被吸入缸内，下一个周期又如此重复。图2.1水在往复泵及其传输管路中的流动[2]实际应用中，人们往往在传输管路A中接入一个空气腔K以保证供水阀门关闭期间管路内液体流动的连续性。流体在往复泵与其传输管路中的流动特性给血液在循环系统中的流动以启示。英格兰生理学家StephenHales(1677-1761)指出，由于主动脉的弹性扩张使心脏的周期性射血变成血管中血液的平稳流动。他把主动脉比拟为空气腔，并引入血液流动外周阻力的概念，认为外周阻力主要来自于人体组织中的小血管，由此发展为后来的Windkessel模型。Windkessel模型在心血管系统与往复泵供水系统间建立起如下的比拟[3]：心血管系统往复泵供水系统左心室往复泵H主动脉与大动脉空气泵K大动脉与毛细血管终端阻力R静脉与心房注水槽V主动脉瓣供水阀门1房室瓣（二尖瓣）供水阀门2当心室收缩、房室瓣关闭、主动脉瓣开启时，心室向主动脉射血：从心室射出的血液一部分经过动脉与毛细血管直接进入静脉腔，另一部分则贮留在主动脉与大动脉之中，使其产生弹性扩张(见图2.2a)；在心室舒张、房室瓣开启、主动脉瓣关闭期间，心脏停止向主动脉射血，但在主动脉与大动脉管壁弹性恢复力的作用下，贮留在其中的血液继续向前流动(见图2.2b)。由此可见，主动脉与大动脉管壁的弹性作用使心脏间断的射血转变为动脉管道中连续的血流。a收缩期b舒张期图2.2主动脉与大动脉的弹性腔作用Windkessel模型将主动脉与大动脉比拟为一个弹性腔，对于某一确定时刻，认为腔内的各种流动参数(如血液压力P和流量g等)仅是时间，的函数，而与距离心脏的距离x无关。Windkessel模型只是分析循环系统中血液流动的一种相当粗糙的模型，但是它能反映动脉血流的一些重要特性，至今仍广为采用。2.2单弹性腔模型的原理图2.3动脉弹性腔的血流量图如图2．3所示的单弹性腔模型[4]，𝑞𝑖𝑛表示单位时间内从心脏流入动脉弹性腔的血液体积，𝑞𝑜𝑢𝑡表示单位时间内从动脉弹性腔经微动脉与毛细血管流入静脉腔的血液体积，V为这段动脉血管的容积。生理压力范围内，忽略血液体积因压力而严生的变化，并认为血液在血管中的流动是连续的，得到心室收缩期和舒张期内血流流动的连续性方程[4]𝑞𝑜𝑢𝑡+𝑑𝑉𝑑𝑡=𝑞𝑖𝑛(2.1)𝑞𝑜𝑢𝑡+𝑑𝑉𝑑𝑡=0(2.2)在弹性腔模型假设下，动脉管中的压力处处相等记为p；静脉管中的压力也处处相等，记为𝑝𝑣，令△p=p−𝑝𝑣表示血液由动脉管留至静脉产生压力降，R为外周阻力，由△p𝑞=R有p−𝑝𝑣𝑅=𝑞𝑜𝑢𝑡(2.3)𝐶=p−𝑝𝑣𝑅(2.4)其中C为血管顺应性。根据(2.1)-(2.4)，可得𝐶𝑑𝑝𝑑𝑡+𝑝−𝑝𝑣𝑅=𝑞𝑖𝑛(0≤𝑡≤𝑇𝑠)(2.5)𝐶𝑑𝑝𝑑𝑡+𝑝−𝑝𝑣𝑅=0(𝑇𝑠≤𝑡≤𝑇)(2.6)其中，Ts表示心室收缩期的持续时间，T表示心动周期。(2.5)、(2.6)便是弹性腔模型基础上，动脉压力p在收缩期与舒张期中分别必须满足的两个方程，各参数取值的情况如下：（1）𝑝𝑣表示静脉腔中的血液压力，通常取恒定值𝑝𝑣=9.3×102−17.3×102𝑃𝑎(7-13mmHg)，方便起见，有时也取为零；（2）R表示外周阻力，看座与压力p无关的常量；（3）𝑞𝑖𝑛表示心脏收缩射出的血液流量，随时间t而变；（4）C表示动脉顺应性，可以使与血压无关的常量（线性弹性腔理论），也可以是一个随血压而变化的变量（非线性弹性腔理论）。下面我们将基于非线性弹性腔理论建立心血管系统的Simulink仿真模型。2.3基于单弹性腔的Simulink仿真模型2.3.1人体心脏搏动间歇模型在仿真试验中我们用正弦曲线仿真心脏的搏动间歇流𝑞𝑖𝑛(𝑡)：心脏收缩时，心室射出的血液流量按照正弦规律变化：而当心室舒张时，血流量为零。这样能近似反映出心脏间断的射血过程[5]𝑞𝑖𝑛(𝑡)={𝑞0𝑠𝑖𝑛𝜋𝑡𝑇𝑠(0≤𝑡≤𝑇𝑠)0(𝑇𝑠≤𝑡≤𝑇)(2.9)Simulink仿真模型如图2.4所示。2.3.2人体动脉管的p-V关系对于一阶模型而言，我们近似认为动脉p、V满足如下的指数函数关系[6]𝑉=𝑎𝑒𝑏𝑝+𝑐(2.7)其中：a，b，c均为常数，需要通过拟合法求值。此时动脉顺应性C将表示为𝐶=𝑑𝑉𝑑𝑝=𝑎𝑏𝑒𝑏𝑝(2.8)对于各类动脉，常数a的取值变化范围较大，b值则是比较接近的。我们通过查阅资料，得到b的平均值为𝑏̅=−0.0131。Simulink的关系图为图2.5所示。图2.4心脏搏动间歇模型图2.5p-V关系图2.3.3系统整体框图根据以上的讨论，我们得出心血管系统一阶Simulink的仿真模型如图2.6所示。图2.6心血管系统一阶仿真模型2.3.4参数分析我们通过查阅资料[7]，取b=−0.0131，𝑇𝑠=0.15𝑠，T=0.8s，其他参数取典型值。仿真时间设置为20s，系统达到稳定状态后的基本波形如图2.7所示。图2.7基本仿真波形图对于图2.7的基本仿真波形图。上排曲线表示心脏每搏输出波形（单位：ml/beat⋅s），下排曲线表示心血管系统一腔室模型内的血压变化（单位：mmHg），横坐标为时间轴（单位：s）。经过若干个心脏循环周期，系统逐渐达到稳定的工作状态。对于不同的人体样本，模型中各参数的取值略有不同，而相应的波形也会有变化。2.4本章小结基于Windkessel理论建立起的心血管系统一阶仿真模型，较好地模拟了动脉血压变化，缺点是无法描述脉搏波的波形细节，下一章我们将具体介绍一种双弹性腔模型的建立方法。3.双弹性腔仿真模型的研究由于单弹性腔模型过于简化，无法解释舒张期所出现的一些诸如潮波、重搏波之类的脉搏波波纹，而正是这些波纹在临床上有重要的意义。为了提高分析的精度，根据对脉搏压力曲线和流量曲线进行傅里叶分析后，人们认为应该在模型中引入一些与频率有关的元件，必须考虑血液流动的黏滞惯性。在前人的研究基础上，1967年由美国的戈特温和瓦特共同提出了双弹性腔模型[8]，该模型采用两个串联的弹性腔，以表现血管系统的不同压力，同时还在两个弹性腔体之间加入一个表示血液惯性的环节，使模型输出的脉搏波曲线能很好地反映出舒张期波纹。双弹性腔模型将主动脉及其主要分支视为两个弹性腔，顺应性分别为𝐶1和𝐶2。心室收缩时，从心室泵出的血液进入第一个弹性腔（用来表征主动脉弓及其主要分支的集总顺应性效应）与血液惯性元件三（用来表征主动脉中的集总血液惯性效应），而后进入第二个弹性腔（用来表征腹主动脉及其主要分支的集总顺应性效应)，最后流经集总的外周阻力R(用来表征外周血管床的总阻力)进入静脉腔[9]。下面我们将利用这种双弹性腔理论，建立对应的心血管系统三阶模型。3.1双弹性腔模型的基本方程图3.1循环系统双弹性腔的血流量图如图3.1所示的双弹性腔模型，假设弹性腔I的顺应性为𝐶1，腔内血液压力为𝑝1，弹性腔体积为𝑉1；弹性腔II的顺应性为𝐶2，腔内血液压力为𝑝2，弹性腔体积为𝑉2，两个弹性腔之间连接着一个长为l，截面积为A的血液柱，并用一个集总的流动惯性L来体现血液在主动脉中的惯性效应。假设心室将血液流量𝑞𝑖𝑛函射进弹性腔I和惯性元件L，再进入弹性腔II，最后流经外周阻力R进入静脉腔。动脉中血液的流动是脉动的，其压力和流量都是一些随时间作周期性变化的量，可以看成是由“直流”成分（即平均压力和平均流量）和“交变”成分组成，这样便可以用构成电路相类似的方法引进血流阻力元件、顺应性元件与惯性元件，用这些元件建立起相应的血液流动等效回路[10]。对于血流阻力元件R与惯性元件L，只要将它们依次画在（即串联在）相应的两个压力点𝑝1、𝑝2之间即可。对于反映血液积聚能力的血管顺应性元件C，虽然它是与电容相比拟的，但是却不能像电容一样可以随意串联在回路的任意地方。在建立血液流动的等效回路时，顺应性元件必须与其后继元件并联，而不能串联。这是因为顺应性代表弹性腔室，其压力差为腔室压力与基准压力之差，而基准压力一般都选用人体周围外界的大气压力（可看作零压力），所以顺应性C两端总是一端连接腔室压力另一端连基准压力（零压力）。如果不考虑连接两个弹性腔的血液柱的粘性阻力，取静脉作为等效回路中的地（零压力），不难看到对应的等效回路为图3.2双弹性腔模型的等效电路由这个等效回路，我们可以建立起相应的方程[13]。根据血流量守恒的条件，知𝑞𝑖𝑛=𝑞𝑐1+𝑞𝑙(2.7)𝑞𝑙=𝑞𝑐2+𝑞𝑜𝑢𝑡(2.8)同时，根据定义可知𝑞𝑐1=𝐶1𝑑𝑝1𝑑𝑡(2.9)𝑞𝑐2=𝐶2𝑑𝑝2𝑑𝑡(2.10)𝑞𝑜𝑢𝑡=𝑝2−𝑝𝑣𝑅(2.11)𝑝1−𝑝2=L𝑑𝑞𝑙𝑑𝑡(2.12)式中：𝑝𝑣——静脉腔内血液压力；𝑞𝑙——通过血液惯性原件l的血流量；𝑞𝑐1和𝑞𝑐2——贮存在弹性腔I和II中的血液流量；𝑞𝑜𝑢𝑡——通过外周阻力R流入静脉腔的血流量。由(2.7)-(2.12)，不难得到描述双弹性腔模型的基本方程𝑑𝑞𝑙𝑑𝑡=1𝐿(𝑝1−𝑝2)(2.13)𝑑𝑞1𝑑𝑡=1𝐶1(𝑞𝑖𝑛−𝑞𝑙)(2.14)𝑑𝑞2𝑑𝑡=1𝐶2(𝑞𝑙−𝑝2−𝑝𝑣𝑅)(2.14)在这组基本方程中，静脉压𝑝𝑣和每搏心输出量𝑞𝑖𝑛通常认为是已知的，因而通过这三个方程，可以确定出待求的三个未知量𝑞𝑙，𝑝1与𝑝2。3.2Simunnk仿真模型的建立3.2.1双弹性腔仿真模型的构建对于基本方程中的所有除法运算都用增益表示，所有的微分都用逆向积分运算表示，心输出量𝑞𝑚仍采用一阶仿真中的搏动间歇流模型，并假设从心室搏出的血流中96%流入主动脉[14]，据此可做出心血管系统三阶Simulink仿真模型图3.3双弹性腔仿真模型其中𝑞𝑖𝑛={𝑞0sin𝜋𝑡𝑇𝑆(0≤𝑡≤𝑇𝑆)0(𝑇𝑆≤𝑡≤𝑇)的计算过程如下：假设心脏每搏输出量70mL(即一