微循环血管不断地调节自身结构，以适应组织的功能性需求。微循环结构自调节模型能够仿真这一过程，为生理研究提供辅助，但目前缺少合适的模型参数设置方法，限制了模型的进一步应用。本文提出一种改进的量子粒子群优化算法用于设置模型参数，并在真实的大鼠肠系膜微循环血管网络上进行仿真实验。仿真结果表明，该方法的参数优化能力优于标准粒子群算法、标准量子粒子群算法和相关文献报道的 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 算法解决模型参数的优化问题。结果显示该优化方法优于传统算法，能使模型的仿真结果与动物实验数据更为接近。但模型仿真与动物实验的结果仍存在一定差异，这可能是由于当前模型中每段血管对刺激信号的响应程度是一致的，而实际生理情况下，每段血管的反应可能有所不同。这还需要今后结合大量动物实验加以改进。
