一种左心室压力个性化估计模型参数子集选择方法

柳 军， 郝丽玲， 何光宇， 徐礼胜,

(1. 东北大学 医学与生物信息工程学院， 辽宁 沈阳 110169； 2. 中国医科大学 生物医学工程系， 辽宁 沈阳 110122；3. 沈阳东软智能医疗科技研究院有限公司， 辽宁 沈阳 110167)

心脏中左心室压力波形(left ventricle pressure waveform, LVPW)直接反映了左心室的血流动力学特征，在临床压力容积环指标中是不可或缺的[1]，对于主动脉瓣狭窄等临床病症有重要的诊断和治疗意义[2].

目前LVPW的测量主要是有创测量和无创模型估计两种方法.有创测量是在临床中通过导管介入手术实现，但是病人需要承受一定的创伤和风险，并且医疗成本较高.无创测量的方法主要是通过各种特定背景下的心血管模型估计[3-4]，同时可以适当结合患者的外周无创测量信息.模型估计方法由于无需手术、患者没有创伤成为目前的一个研究方向[5].而针对每个不同病人，模型参数个性化成为近年来的研究热点[6].

通常一个模型中有很多参数，在临床应用中不可能对模型中所有参数进行个性化取值，因为有的模型参数可能无法测量或者临床中的医疗成本过高，所以应当挑选出模型中对LVPW有较大影响的参数进行个性化取值.灵敏度分析(sensitivity analysis, SA)方法近年来逐渐应用于心血管模型参数个性化选择中[7].

本文以一个体循环模型估计LVPW作为实例，采用自适应稀疏多项式混沌展开(adaptive sparse polynomial chaos expansion,asPCE)算法构建出原始模型的元模型(meta model)，采用索贝尔法(Sobol method)进行灵敏度指标计算，最后遴选出对LVPW有较大影响的模型参数作为参数子集.

1 体循环模型

本文以一个体循环模型作为估计LVPW的例子，在最近的研究中也有类似的模型用于估计LVPW[8]，模型的具体原理公式的推导在文献[9]中有详细的介绍，本文不再赘述.总体结构上模型由一个左心室单纤维模型作为驱动，外周血管由Windkessel模型构成，主动脉瓣和二尖瓣由二极管模型模拟.体循环模型构成原理图如图1所示.

图1 体循环模型构成原理图

模型的输出为LVPW以及4个波形特征，在本文中4个波形特征分别定义为一个心动周期内LVPW的平均压pMB，收缩压pSB，从二尖瓣关闭时刻到收缩压pSB的时间间隔tr，从收缩压pSB时刻到二尖瓣打开时刻的时间间隔td，具体波形特征在LVPW中的含义如图2所示.

图2 左心室压力波形特征

本文体循环模型中涉及的模型参数分为三类：生理类参数，Windkessel外周模型参数，单纤维模型心肌纤维参数.表1～表3分别给出模型中各个参数的单位、描述及其范围.

表1 体循环模型中的生理类参数Table 1 Physiological parameters of systemic circulation model

表2 体循环模型中的外周参数Table 2 Peripheral parameters of systemic circulation model

表3 体循环模型中的心肌纤维参数Table 3 Myocardial fiber parameters of systemic circulation model

2 模型参数灵敏度计算

2.1 灵敏度指标

本文灵敏度分析采用Sobol法进行计算.该方法是基于方差分解[10]的思想，用模型输出的方差来衡量输出不确定性，输出的方差可以分解成偏方差的组合.这些偏方差可以表示每个输入参数或者多个输入参数相互影响的不确定性对输出的影响.

假设模型输出输入关系为

Y=f(X) .

(1)

其中：Y是模型输出量;X是模型输入参数，假设是M维的向量，可以表示为X=(X1X2…XM).本文计算的灵敏度指标包括主灵敏度指标和总灵敏度指标，主灵敏度指标定义如下：

(2)

其中：Si是主灵敏度指标；V表示求方差运算；E表示求期望运算.主灵敏度指标提供的信息表明应该准确确定哪个输入参数来获得最大的方差降幅.总灵敏度指标定义如下：

(3)

其中：ST,i是总灵敏度指标；-i表示去除Xi的参数集合.总灵敏度指标用来决定哪个参数是可以在其数值范围内固定而对最后的结果几乎没有影响.

2.2 元模型构建和灵敏度指标计算

2.2.1 元模型构成及灵敏度计算

各个参数的索贝尔灵敏度指标最初是由蒙特卡洛统计学方法[11]求出，由于蒙特卡洛方法需要进行大量的模型样本运算，计算量巨大，后来有学者[12]提出了采用多项式混沌展开来替代原始模型，根据展开式的系数来计算索贝尔灵敏度指标.本文采用自适应稀疏多项式混沌展开算法进行元模型的构建[13]，然后根据多项式系数计算主灵敏度和总灵敏度两个灵敏度指标.

模型的输出可以表示为输入参数的一系列截断后的有限项正交多项式的组合，本文采用勒让德正交多项式，输入参数需要转换到[-1,+1]的范围内，以满足勒让德正交多项式的要求.模型的输出输入关系为

(4)

(5)

(6)

多变量互相阶数和定义如下：

(7)

(8)

多项式的单变量阶数和与多变量互相阶数和可以由设计者设定，考虑到多项式的高阶项在多项式中的影响很小，本文设定如下：

A{3,3}={α|Degree(α)≤3 and Order (α)≤3} .

(9)

本文多项式组合的系数是通过最小化模型真值和多项式计算值差的平方的方法求得.在求得多项式系数后，利用多项式的系数计算索贝尔灵敏度指标.元模型的方差通过下式求出：

(10)

其中:Hα是一个规范化因子，其具体取值取决于多项式的类型.本文采用勒让德正交多项式，其计算公式如下：

(11)

最后可以计算得到主灵敏度指标：

(12)

其中Ai是包括多项式中那些αi是正值而其他阶次是零的集合.总灵敏度指标的计算如下：

(13)

其中AT,i是包括多项式中所有αi的阶次是正值的集合.

2.2.2 asPCE算法构建元模型

在构建元模型和求解多项式系数的过程中，本文采用asPCE的算法.该算法的优点是可以减少模型运算样本数量.在构建元模型的过程中，有选择地挑选多项式组合中对输出结果有较大影响的项，作为最终元模型多项式组合的组成部分.在算法的迭代过程中，将重复进行两步重要多项式的挑选操作：前向选择和后向选择.在前向选择中，从零开始逐渐增加多项式的最高阶次和互相阶次，挑选出可以显著提高元模型准确度的部分，构成前向选择的多项式组合；在构建完前向选择的多项式组合后，再进行后向选择，在后向选择中，去掉那些前向选择中留下的似乎重要但是并没有显著提高元模型准确度的部分.多项式的最高阶次和互相阶次从零开始逐渐增大到设定值，在此过程中，算法需要的模型运算样本点数量将逐渐增加，前向选择和后向选择的操作在不断地迭代进行，直到最后确定元模型所包括的多项式组合.

(14)

(15)

其中fPCEi(→X(i))表示除去样本X(i)所构建的元模型.

2.3 模型样本运算

由于asPCE算法中，求解多项式系数需要用到线性回归方法，因此需要提前进行模型样本的运算.在通用的多项式混沌展开(general polynomial chaos expansion,gPCE)方法中，至少需要NP次模型运算样本.而在asPCE算法中，模型运算样本的数量可以大大减少，本文按照gPCE方法进行了NP=4 060次模型运算，但是实际asPCE算法迭代中只是用到了其中一部分.

在进行模型运算之前，模型的输入参数需要进行随机取值.本文采用Sobol序列准随机数生成方法[14]，该方法以2为底数来形成更精细的单位间隔，并且可以生成输入参数的基于均匀分布的随机序列.

由于模型的输入参数是高维随机序列生成，在模型运算中，模型输出值有很多情况是不符合实际人体的生理数据，因此，要把这些不正确的模型输出值排除在最后元模型的构建过程中.本文根据健康人群和各种心血管病人群的特点以及文献[15]的相关类似设置，设置了以下的模型输出过滤规则：

1) 模型输出的计算是基于连续多个心动周期逐渐收敛稳定，因此设定相邻两个周期的2范数阈值为0.05，当超过15个心动周期并且相邻2个心动周期波形的2范数小于0.05，就认为模型输出波形达到收敛;2) 左心室压力波形的平均值范围为1～10 kPa;3) 左心室压力波形的收缩压范围为5～24 kPa;4) 左心室压力波形在舒张期内的压力范围为0.3～3 kPa;5) 左心室压力波形中收缩期时长与舒张期时长的比值范围为0.5～1.2.

图3 asPCE算法流程与灵敏度计算

3 灵敏度指标结果及元模型与真实模型输出对比

本文模型有27个输入参数和4个左心室压力波形特征输出指标，经过以上主灵敏度和总灵敏度指标的计算，所有输入参数的灵敏度指标计算结果如图4所示.从最后的主灵敏度和总灵敏度结果可以看出，4个左心室压力波形特征分别对应的重要的输入参数有所区别.取4个波形特征所对应的所有重要参数的并集作为最后模型输出左心室压力波形的重要参数子集，共计12个参数，并在表4给出了每个参数对应的4个波形特征的主灵敏度最大值和总灵敏度最大值.

为了检验元模型与真实模型的逼近程度，本文进行了左心室压力波形特征的对比，在输入参数相同的情况下，比较从真实模型和元模型输出的左心室压力波形的4个特征，具体对比结果如图5所示.其中所有相关性图中，r=0.99,P<0.001.并且4个特征分别对应的两组数据没有显著差异(P>0.05).从上面的LVPW的波形特征相关性曲线和一致性结果图中，可以清楚地看到元模型与真实模型具有很好的相关性和一致性.

表4 对应左心室压力波形特征的重要输入参数Table 4 Important parameters corresponding to the features of LVPW

图4 灵敏度分析结果

4 模型输出LVPW对比及灵敏度计算方法对比

为了验证当模型中不重要参数取值固定的左心室压力波形情况，进行了左心室压力波形和波形特征的对比.把模型的参数分成两个部分，重要参数和非重要参数.重要参数取值相同(随机取值)，非重要参数一种情况随机取值，另外一种情况固定取值(经验值).在这两种情况下，对真实模型输出的左心室压力波形和波形特征进行对比，结果如图6和图7所示，其中所有的相关性曲线中，r=0.99,P<0.001.LVPW和波形特征在两种情况下没有显著差异(P>0.05).本文进行了3 431次实验，随机取其中一例波形对比图如图6所示.经过灵敏度选择后，重要参数子集的数量比全部参数数量显著减少，从对比结果图中可以看出，在重要参数取值相同，非重要参数取值固定与取值随机两种情况下，最后输出的左心室压力波形和波形特征都具有很好的一致性.经过灵敏度选择后，重要参数作为重点关注的参数子集，非重要参数取值固定对LVPW波形的估计几乎没有重大影响.

图5 真实模型和元模型输出的波形特征相关性曲线和对应的Bland-Altman图

图6 不同参数取值情况下模型输出LVPW对比

传统的灵敏度计算经常采用蒙特卡洛方法，其劣势在于需要巨量的模型计算样本，计算量大，特别对于参数数量较多的模型，计算量呈几何数量级增长.采用asPCE算法构建原始模型的元模型进行灵敏度指标的计算，一方面可以对不了解原始模型内部原理的研究者提供一个简单的多项式组合的元模型，另一方面，基于此算法进行灵敏度指标的计算，可以显著降低模型计算样本数量的需求，对模型参数较多的情况，效果更加明显.本文模型有27个参数，多项式的最高阶次为3，假定蒙特卡洛样本数量10 000的情况下，蒙特卡洛算法，gPCE算法和本文提出的asPCE算法所需要的模型运算样本数量依次是290 000,4 060和2 582.

图7 不同参数取值情况下模型输出波形特征的相关性分析和Bland-Altman分析

5 结 语

当临床所用心血管模型参数很多时，LVPW估计的个性化和参数优化计算工作将会很复杂.本文提出基于灵敏度分析的方法对模型中的参数进行重要参数子集选择，最终选出对LVPW有较大影响的参数子集，这样的参数子集可为临床患者的LVPW个性化估计提供参考，临床中可以重点关注这些灵敏度指标较高的参数.参数子集中的参数数量减少，可以明显降低模型参数优化计算的复杂度.实验结果表明，基于灵敏度选择的参数子集的LVPW估计，与模型全参数估计具有很高的一致性.