摘要:心脏出现电螺旋波和时空混沌会导致心律失常, 心肌细胞与成纤维细胞电耦合(M-F耦合)导致心肌组织结构改变也会引起心律失常, 如何消除由螺旋波、时空混沌导致的心律失常是科学工作者十分关注的问题. 本文构造了由心肌细胞和成纤维细胞组成的双层复合介质, 采用Luo-Rudy相I心脏模型及被动成纤维细胞模型研究了M-F耦合对螺旋波形成的影响, 以及螺旋波和时空混沌的控制, 提出用提高细胞之间的耦合强度来控制复合介质中的螺旋波和时空混沌, 数值模拟结果表明: M-F耦合对螺旋波动力学有重要影响, 随着成纤维细胞密度的增加, M-F耦合会导致螺旋波漫游和破碎成时空混沌, 甚至会产生从时空混沌(或螺旋波)到无波的相变. 通过提高细胞之间耦合强度来消除复合介质中的螺旋波和时空混沌, 只在大部分情况下有效, 依赖于成纤维细胞所起的作用. 当成纤维细胞起到电流吸收器的作用时, 虽然大部分情况下可以通过提高细胞之间耦合强度来消除螺旋波和时空混沌, 但是可控区较小. 当成纤维细胞起到电流源的作用时, 提高细胞之间耦合强度使之超过临界值基本都可以有效控制螺旋波和时空混沌, 且可控区比前者大为增加, 提高控制效果的关键是提高心肌细胞之间的耦合强度.
关键词: 螺旋波/
时空混沌/
心肌细胞与成纤维细胞电耦合/
控制

English Abstract

全文HTML

--> --> -->

众所周知, 心脏是一种复杂多细胞器官, 由心肌细胞和非心肌细胞组成, 非心肌细胞包括心肌成纤维细胞、平滑肌细胞、常驻干细胞、免疫细胞等^[1], 其中成纤维细胞是数量最多的非心肌细胞. 正常成人心脏中, 心肌细胞、成纤维细胞体积分别占心脏约75%和约25%, 而它们的数量分别占心脏细胞总数的30%—40%和50%以上^[1,2], 例如心肌局部缺血、心脏肥厚、心衰等老化和病理状态还会导致成纤维细胞增殖和心脏纤维化^[3].
成纤维细胞是一种非可激发细胞, 一些心脏疾病与成纤维细胞被激活和成纤维细胞表型转变有关, 成纤维细胞与心肌细胞存在耦合已经在细胞培养物中得到证实, 并发现它们主要通过连接蛋白43(Cx43) 缝隙链接与心肌细胞形成电耦合^[4], 心肌细胞之间也是主要通过Cx43缝隙链接形成电耦合^[1], 而且一个心肌细胞平均与11到12个其他细胞耦合^[5]. 成纤维细胞与心肌细胞电耦合(简称M-F耦合)在心脏电活动中可扮演电桥作用和两种不同角色, 实验发现: 被若干成纤维细胞阻隔的两个心肌细胞可以建立同步的收缩, 这就是成纤维细胞扮演的电桥作用^[4], 当成纤维细胞的静息电位比心肌细胞的静息电位高很多时, 成纤维细胞起到电流源的作用, 在一定条件下可以导致心肌细胞自发激发, 否则纤维细胞起到电流吸收器的作用, 降低心肌细胞的激发性^[6]. 此外M-F耦合还会导致心肌细胞动作电位延长, 因为通常成纤维细胞的静息电位高于心肌细胞, 因此当M-F耦合较弱时, M-F耦合会抬高心肌细胞的静息膜电位, 导致传导速度增大; 对于中等M-F耦合, 由于心肌细胞的静息膜电位抬高达到饱和, 电流流向成纤维细胞, 将导致传导速度减小和心肌细胞膜电位的最大上升速度减少; 强M-F耦合则会使各种效应达到饱和, 在强耦合和高成纤维细胞密度下, 心肌细胞动作电位会延长^[6], 在一定条件下, M-F耦合可导致动作电位交替^[7]和早期后除极化^[8].
成纤维细胞在心脏的电功能失调中扮演的重要角色促使人们从理论和实验角度去研究成纤维细胞的作用^[9-13], 提出了各种成纤维细胞模型^[9], 如被动模型、各种主动模型、心房成纤维细胞模型等. 2016年, Greisas和Zlochiver^[3]用细胞模型研究了M-F耦合对心肌组织传导速度的影响, 提出了一个心肌细胞占据一个格点, n_f个成纤维细胞占据一个格点的模拟方案, 在数值模拟中n_f 值最大值为20, 因为成纤维细胞比心肌细胞小很多. 2009年, Shajahan等^[14]使用数值模拟方法对心脏中螺旋波湍流态(即时空混沌态)的控制问题进行了研究, 他们发现, 在由心肌细胞和成纤维细胞组成的复合介质中, 局部电刺激无法消除螺旋波湍流态, 于是他们提出了网线控制方法, 即用一定大小正方形网格线覆盖在介质上, 将电刺激施加在网线上, 使与网线接触的细胞被电击, 这种控制方法可以有效控制螺旋波和螺旋波湍流态, 显然这种控制方法在临床上不容易实施. 2013年, Nayak等^[15]采用由心肌细胞和成纤维细胞组成的复合介质, 用数值模拟方法研究了螺旋波动力学, 发现在低M-F耦合强度时出现传导失败, 可以通过调节M-F耦合强度和成纤维细胞的静息电位控制螺旋波的旋转速度和稳定性, 观察到螺旋波被成纤维细胞分布不均匀区钉扎, 并证实了用网线控制方法可以有效控制螺旋波, 但是他们没有提出更容易操作的控制方法, 因此需要进一步探索复合介质中螺旋波和时空混沌的控制.
最近, Kofron等^[16]通过实验发现: 通过腺病毒基因在成纤维细胞中过度表达G_q蛋白α-亚基或使Gα_q持续活性, 可以将成纤维细胞的静息电位从–33 mV增加到–20 mV, 并且增加Cx43和Cx45的表达, Cx43表达增加5倍, Cx45表达增加2.5倍, 从而较大地增加M-F耦合强度. 考虑到可以实现调节Cx43表达来调节细胞之间的耦合强度, 成纤维细胞以片状或条状的方式沿平行肌纤维方向分布, 或者分布在两层心肌组织之间的缝隙中^[17], 本文采用Luo和 Rudy^[18]提出的心脏模型(记为L-R相I心脏模型)和文献[7]使用的成纤维细胞被动模型和参数研究了老化心肌组织中螺旋波和时空混沌的控制, 构造一种复合介质, 它由两层老化的心肌组织组成, 每层心肌组织都随机嵌入一定比率的成纤维细胞, 并在两层心肌组织之间夹有成纤维细胞, 且每个心肌细胞可与多个成纤维细胞耦合. 根据成纤维细胞所起的作用, 提出通过提高细胞之间耦合强度(即心肌细胞之间、成纤维细胞之间和心肌细胞与成纤维细胞之间的耦合强度)来控制复合介质中的螺旋波和时空混沌, 由于两个细胞之间的耦合强度与接触面积有关, 因此我们在模拟中假设成纤维细胞之间和心肌细胞与成纤维细胞之间的耦合强度是一样的, 它与心肌细胞之间的耦合强度是不相同的. 通过数值模拟发现, 如果成纤维细胞起到电流源的作用时, 只要将细胞之间耦合强度增加到足够大时就可以有效消除螺旋波和时空混沌, 且成纤维细胞比率越大控制效果越好. 下面先介绍本文模型, 然后介绍数值模拟结果, 最后给出结论.

考虑到随着心脏的老化, 成纤维细胞会不断增殖, 使纤维细胞数量逐渐增加, 而高成纤维细胞密度又会导致单层复合心肌组织容易出现传导障碍, 因此本文使用由上下两层二维复合心肌组织组成的介质, 以避免出现单层传导障碍, 更贴近实际心脏动力学行为. 每一层心肌组织大小为8.4 cm × 8.4 cm, 离散成300 × 300格点, 一个心肌细胞占据一个格点, 而成纤维细胞可以$n$个占据一个格点, 因为通常心肌细胞所占体积是成纤维细胞的几倍^[7]. 每一层二维心肌组织中随机嵌入一定比率的成纤维细胞, 每一层复合心肌组织中, 成纤维细胞占格点数与总格点数之比称为成纤维细胞密度$\rho $.
两层复合心肌组织之间对应点可能直接耦合, 也可能通过成纤维细胞间接耦合, 看起来好像两层复合心肌组织之间夹有一层成纤维细胞层(中间层), 该中间层的成纤维细胞分布按如下方法随机分布在300 × 300正方点阵上, 并且$n$个成纤维细胞占据一个格点, 并与上下两层复合心肌组织层对应点上的细胞形成电耦合, 在中间层中, 相邻格点之间的成纤维细胞无耦合. 下面给出中间层细胞的分布: 假设分别在两复合心肌组织层上随机覆盖占比为50%的成纤维细胞, 在对应格点上, 如果两复合心肌组织层都覆盖有成纤维细胞, 则中间层在该点有成纤维细胞占据, 而且与两层复合心肌组织层对应点的细胞都有电耦合; 如果对应点只有一层复合心肌组织层覆盖有成纤维细胞, 另一层没有, 则中间层在该对应点上也有成纤维细胞, 只是该成纤维细胞只与被覆盖心肌组织层对应点的细胞有电耦合; 如果对应点上下两层皆没有覆盖的成纤维细胞, 则中间层在该点没有成纤维细胞, 这时在该对应点上下两复合心肌组织层的细胞直接耦合.
心肌细胞的动力学方程采用L-R相I模型^[18], 成纤维细胞采用文献[7]使用的被动模型, 每个格点上心肌细胞和成纤维细胞的动力学方程如下^[18,7]:

$\begin{split}\frac{{\partial {V_{\rm{m}}}}}{{\partial t}} = \;&- \frac{{{I_{{\rm{ion}}}}}}{{{C_{\rm{m}}}}} + \frac{n}{{{C_{{\rm{mz}}}}}}\sum\limits_{k = 1}^{{\rm{num}}1} {{G_{{\rm{gap}}}}\left( {V_{f,k}^{} - {V_{\rm{m}}}} \right)}\\&+ \frac{D}{{{\varDelta ^2}}}\sum\limits_{k = 1}^{{\rm{num}}2} {\left( {V_{{\rm{m}},k}^{} - {V_{\rm{m}}}} \right)}, \end{split}\tag{1a} $

$\begin{split}\frac{{\partial {V_{\rm{f}}}}}{{\partial t}} =\;& - \frac{{{I_{\rm{f}}}}}{{{C_{\rm{f}}}}} + \frac{1}{{{C_{\rm{f}}}}}\sum\limits_{k = 1}^{{\rm{num}}3} {{G_{{\rm{gap}}}}\left( {V_{{\rm{f}},k}^{} - {V_{\rm{f}}}} \right)}\\&+ \frac{1}{{{C_{\rm{f}}}}}\sum\limits_{k = 1}^{{\rm{num}}4} {{G_{{\rm{gap}}}}\left( {V_{{\rm{m}},k}^{} - {V_{\rm{f}}}} \right)} , \end{split}\tag{1b}$

${I_{\rm{f}}} = {G_{\rm{f}}}\left( {{V_{\rm{f}}} - {E_{\rm{f}}}} \right), \tag{1c}$

${I_{{\rm{ion}}}} = {I_{{\rm{Na}}}} + {I_{{\rm{Si}}}} + {I_{\rm{K}}} + {I_{{\rm{k1}}}} + {I_{{\rm{kp}}}} + {I_{\rm{b}}}, \tag{1d}$

式中, ${V_{\rm{m}}}$和${V_{\rm{f}}}$分别为心肌细胞和成纤维细胞膜电压, 单位为mV; C_m为心肌细胞的膜电容, ${C}_{\rm{m}}= $$ 1 \;{\rm{\text{μ} F/cm}}^{2}$, C_mz为整个心肌细胞膜电容, ${C}_{\rm{mz}}= $$ 125\; {\rm{pF}}$; C_f为成纤维细胞的膜电容, $ {C}_{\rm{f}}=25\; {\rm{pF}}$; t为时间, 单位为ms; ${I_{{\rm{ion}}}}$为心肌细胞总膜电流, 包括钠离子电流${I_{{\rm{Na}}}}$、钙离子电流${I_{{\rm{Si}}}}$、外行钾离子电流${I_{\rm{K}}}$、不含时外行钾离子电流${I_{{\rm{k1}}}}$、饱和钾离子电流${I_{{\rm{kp}}}}$和背景电流${I_{\rm{b}}}$, 单位为${\rm{\text{μ}F/c}}{{\rm{m}}^2}$; ${I_{\rm{f}}}$为成纤维细胞总膜电流, 单位为nA; $n$是占据一个格点的成纤维细胞个数, 也是一个心肌细胞与相邻一个格点上的成纤维细胞耦合的数量, 称为细胞耦合数; ${G_{{\rm{gap}}}}$为成纤维细胞之间和成纤维细胞与心肌细胞之间的耦合强度, 单位为nS; ${G_{\rm{f}}}$为成纤维细胞的膜电导率, 单位为nS; ${E_{\rm{f}}}$为成纤维细胞的静息电位; D既为心肌细胞的扩散系数, 也是心肌细胞之间的耦合强度, $D = 0.001\, \, {\rm{c}}{{\rm{m}}^2}{\rm{/ms}}$. 在方程(1)的右边求和中, ${\rm{num}}X$意义为: 一个心肌细胞的5个近邻格点中有num1个格点被成纤维细胞占据, 有num2个格点被心肌细胞占据, 而一个成纤维细胞的近邻格点中有num3个格点被成纤维细胞占据, 有num4个格点被心肌细胞占据. 除特别指出外, 心肌细胞的各离子电流、门变量、离子浓度的表达式和使用参数均与文献[18]一致.
数值求解方程(1)时, 固定取钾离子最大电导率$\overline {{G_{\rm{K}}}} = 0.705\, \, {\rm{mS/c}}{{\rm{m}}^2}$, 钙离子最大电导率$\overline {{G_{{\rm{si}}}}} = $$ 0.02\, \, {\rm{mS/c}}{{\rm{m}}^2}$, 成纤维细胞静息电位$ {E}_{\rm{f}}=-50 \;{\rm{mV}}$, 空间步长$ \varDelta =0.028 \;{\rm{cm}}$, 时间步长$ \Delta t=0.02\;{\rm{ms}}$, 使用无流边界条件, 时间求导使用一阶欧拉向前差分法, 空间二阶导数采用中心差分方法, 方程(1a)右式最后一项就是离散后的扩散项.
为了描述控制效果, 引入平均膜电位差如下:

$\overline V = \frac{1}{{2 \times {{300}^2}}}\mathop \sum \limits_{i,j = 1}^{300} \left( {\left| {V_{i,j}^1 - {V_{\rm{r}}}} \right| + \left| {V_{i,j}^{\rm{2}} - {V_{\rm{r}}}} \right|} \right), $

式中, $V_{i, j}^1$, $V_{i, j}^2$分别表示第一、二层复合介质在格点$(i, j)$上的细胞膜电位值; 固定取值如下: 当格点上为心肌细胞时${V_{\rm{r}}} \approx - 80.0\;{\rm{ mV}}$, 若该格点上为成纤维细胞则使用${V_{\rm{r}}} = - 55\;{\rm{ mV}}$. 当螺旋波和时空混沌因控制作用消失时, 所有细胞处于静息态, $\overline V $基本不再随时间变化, $\overline V \text{-} t$曲线变成水平直线, 平均膜电位差最初成水平直线对应的时刻就是螺旋波和时空混沌被控制所用时间, 此时$\overline V $不为0, 原因是: 在不同成纤维细胞密度和控制参数下的静息膜电位${V_{\rm{r}}}$其实也会发生变化, 不同环境下细胞的静息电位也不相同, 而在(2)式中${V_{\rm{r}}}$取固定值, 这就导致$\overline V $不为0. 如果在控制作用下大部分波因传导障碍消失, 剩余的波会演化成螺旋波和时空混沌, 这时$\overline V $达到一个较小值后会增大, 不会维持在一个较小的值, 此后$\overline V $规则变化对应系统处于螺旋波态, $\overline V $不规则变化则对应系统处于时空混沌态, 总之${V_{\rm{r}}}$取其他值不影响判断螺旋波和时空混沌是否被控制.

首先采用截断平面波方式在双层复合介质中产生螺旋波, 研究M-F耦合对螺旋波的形成与传播的影响, 将${G_{\rm{f}}}$, ${{\rm{G}}_{{\rm{gap}}}}$, $\rho $和$n$看成可调参数. 图1给出了在$D = 0.001\;{\rm{c}}{{\rm{m}}^{\rm{2}}}{\rm{/ms}}$和不同参数复合介质中通过截断平面波形成螺旋波得到的相图情况, 典型螺旋波和时空混沌斑图如图2所示, 这里只给出了其中一层的斑图, 原因是另一层斑图外观上看是一致的. 计算了两层介质对应点的膜电位差的绝对值的平均值, 图2(a)—(d)对应的平均值分别为6.45, 8.65, 11.40, 5.90, 此平均值不为0的原因有以下三方面: 1)出现很多对应点上细胞不是同类型细胞, 它们的电位有较大的不同; 2)即使对应点上是同一类细胞, 由于所处环境不同, 其膜电位也有差别; 3)当波在第一层传播受到阻碍时, 通过另一层接力得以继续在第一层传播, 导致对应点细胞不是同时被激发. 由于这个平均值与心肌细胞约100 mV的振幅相比并不是很大, 导致外观上看两层斑图是一致, 因此在后面只给出一层上的斑图.
图 1 不同参数下在$\rho \text{-} n$参数平面上的相图(■代表稳定螺旋波; ●代表漫游单螺旋波; ◆代表不稳定螺旋波; ▼代表时空混沌态; ▲代表介质无波)　(a) $ {G}_{\rm{f}}=1\;{\rm{ns}}$, ${G}_{\rm{gap}}= $$ 3\;{\rm{ns}}$; (b) $ {G}_{\rm{f}}=4\;{\rm{ns}}$, $ {G}_{\rm{gap}}=1\;{\rm{ns}}$
Figure1. Phase diagrams on the $\rho \text{-} n$ parameter plane with different parameters: (a) $ {G}_{\rm{f}}=1\;{\rm{ns}}$, $ {G}_{\rm{gap}}=3\;{\rm{ns}}$; (b) ${G}_{\rm{f}}= $$ 4\;{\rm{ns}}$, $ {G}_{\rm{gap}}=1\;{\rm{ns}}$. Black square, stable spiral wave; black circle, meandering of single spiral wave; black rhombus, unstable spiral wave; black down-triangle, spatiotemporal chaos; black up-triangle, no wave can propagate in the medium.

图 2 不同参数下其中一层心肌层的螺旋波斑图　(a) ${G}_{\rm{f}}= $$ 4\;{\rm{ns}}$, $ {G}_{\rm{gap}}=1\;{\rm{ns}}$, $\rho = 20{\text{%}}$, $n = {\rm{1}}$; (b) $ {G}_{\rm{f}}=4\;{\rm{ns}}$, $ {G}_{\rm{gap}}=1\;{\rm{ns}}$, $\rho = 30{\text{%}}$, $n = {\rm{3}}$; (c) $ {G}_{\rm{f}}=4\;{\rm{ns}}$, $ {G}_{\rm{gap}}=1\;{\rm{ns}}$, $\rho = 40{\text{%}}$, $n = {\rm{2}}$; (d) $ {G}_{\rm{f}}= $$ 1\;{\rm{ns}}$, $ {G}_{\rm{gap}}=3\;{\rm{ns}}$, $\rho = {\rm{4}}0{\text{%}}$, $n = $4
Figure2. Pattern of spiral wave in one layer cardiac tissue for different parameters: (a) $ {G}_{\rm{f}}=4\;{\rm{ns}}$, $ {G}_{\rm{gap}}=1\;{\rm{ns}}$, $\rho = 20{\text{%}}$, $n = {\rm{1}}$; (b) $ {G}_{\rm{f}}=4\;{\rm{ns}}$, $ {G}_{\rm{gap}}=1\;{\rm{ns}}$, $\rho = 30{\text{%}}$, $n = {\rm{3}}$; (c) ${G}_{\rm{f}}= 4\;{\rm{ns}}$, $ {G}_{\rm{gap}}=1\;{\rm{ns}}$, $\rho = 40{\text{%}}$, $n = {\rm{2}}$; (d) ${G}_{\rm{f}}= $$ 1 \;{\rm{ns}}$, $ {G}_{\rm{gap}}=3\;{\rm{ns}}$, $\rho = {\rm{4}}0{\text{%}}$, $n = {\rm{4}}$.

由图1可得出: M-F耦合可导致螺旋波漫游和螺旋波破碎, 螺旋波的波臂如在波头附近折断后波头形成孤立波(参见图2(b)), 或螺旋波波臂在波头附近折断后波头形成小螺旋波对, 但是大部分螺旋波的波臂仍完好时, 将这种螺旋波称为不稳定螺旋波, 这就是图1中用符号◆标记的斑图; 螺旋波波臂折断若不局限在波头附近, 在其他地方也出现波臂折断时的态称为时空混沌态, 这种斑图如图2(c)和图2(d)所示, 它们就是图1中用符号▼标记的斑图. M-F耦合还可以使介质无波, 这种斑图在图1中用符号▲标记. 介质无波的原因有两种: 第一种原因是由于介质激发性很低, 平面波无法传播; 第二种原因是$\rho $和$n$很大, 由于模拟中使用的初态使全部细胞处于静息态, 成纤维细胞膜电位的初值为–50 mV比心肌细胞膜电位的初值–84 mV高很多, 成纤维细胞起到电流源的作用, 导致介质大范围出现自激发, 自激发后的心肌细胞由于M-F耦合无法回到静息态, 从而一直处于不应态, 使得介质无法再传播波而导致无波, 图1(b)给出的(ρ, n) = (50%, 5)参数下的结果就属于这种情况.
由图1还可以看出, M-F耦合对螺旋波影响不仅与成纤维细胞密度有关, 还与M-F耦合强度、成纤维细胞膜电导率、细胞耦合数量有关. 在成纤维细胞密度从10%增加到50%的过程中, 不同参数下M-F耦合对螺旋波影响按如下发展变化: 稳定螺旋波→漫游(包括不稳定螺旋波)→时空混沌态→介质无波态. 图1(b)还出现了另一种发展变化: 稳定螺旋波→漫游螺旋波(包括不稳定螺旋波)→介质无波态. 且基本都是在$\rho \leqslant $ 30%时只有螺旋波态, 在$\rho = $ 40%情况下出现时空混沌态, 在$ \rho = $ 50%情况下出现无波态.
通过对比图1(a)和图1(b)可以看出, 两层复合介质的激发性由于${G_{\rm{f}}}$和${G_{{\rm{gap}}}}$的不同而有所不同, 图1(b)对应的这种复合介质激发性更好一些, 因为出现稳定螺旋波情况多了, 出现时空混沌情况少了. 为了解释其中的原因, 在图3给出了在$\rho = $ 30%和不同$n, \;{G_{\rm{f}}}, \;{G_{{\rm{gap}}}}$情况下成纤维细胞静息电位随时间的变化, 这个静息电位就是在系统演化过程中成纤维细胞膜电位达到的最小值. 我们区分如下六种情况: 被成纤维细胞占据的格点周围无成纤维细胞, 被成纤维细胞占据的格点近邻有1个、2个、3个、4个和5个格点有成纤维细胞, 分别用$\sigma = 1, \;2, \;3, \;4, \;5, \;6$表示这六种情况. 由于在低密度下不容易出现$\sigma = 4, \;5, \;6$情况, 所以在图3只给出成纤维细胞密度在30%情况下成纤维细胞静息电位随时间的变化, 图3(a)和图3(b)对应图1(a)情况, 图3(c) 和图3(d)对应图1(b)情况. 从图3可以看出: 随着$\sigma $增加, 成纤维细胞静息电位增加一般都十分明显, 但是$\sigma = 4, \;5$情况下存在例外. 图3(a)和图3(b)显示, 纤维细胞静息电位在[–79.5, –73] mV之间变化, 图3(c)和图3(d)显示, 纤维细胞静息电位在[–67, –55] mV之间变化, 增加细胞耦合数$n$, 也会提高成纤维细胞静息电位. 数值模拟结果表明: 增加成纤维细胞密度$\rho $, 同样也会提高成纤维细胞静息电位的上限, 减少成纤维细胞密度, 则会降低成纤维细胞静息电位的下限.
图 3 不同参数下成纤维细胞静息电位随时间变化　(a) ${G_{\rm{f}}} = 1\;{\rm{ns}}$, ${G_{{\rm{gap}}}} = 3\;{\rm{ns}}$, $n = 2$, $\rho = 30{\text{%}}$; (b) ${G_{\rm{f}}} = 1\;{\rm{ns}}$, ${G_{{\rm{gap}}}} = 3\;{\rm{ns}}$, $n = 5$, $\rho = 30{\text{%}}$; (c) ${G_{\rm{f}}} = 4\;{\rm{ns}}$, ${G_{{\rm{gap}}}} = 1\;{\rm{ns}}$, $n = 2$, $\rho = 30{\text{%}}$; (d) ${G_{\rm{f}}} = 4\;{\rm{ns}}$, ${G_{{\rm{gap}}}} = 1\;{\rm{ns}}$, $n = 5$, $\rho = 30{\text{%}}$
Figure3. Evolution of the resting potential of a fibroblast for different control parameters: (a) ${G_{\rm{f}}} = 1\;{\rm{ns}}$, ${G_{{\rm{gap}}}} = 3\;{\rm{ns}}$, $n = 2$, $\rho = 30{\text{%}}$; (b) ${G_{\rm{f}}} = 1\;{\rm{ns}}$, ${G_{{\rm{gap}}}} = 3\;{\rm{ns}}$, $n = 5$, $\rho = 30{\text{%}}$; (c) ${G_{\rm{f}}} = 4\;{\rm{ns}}$, ${G_{{\rm{gap}}}} = 1\;{\rm{ns}}$, $n = 2$, $\rho = 30{\text{%}}$; (d) ${G_{\rm{f}}} = 4\;{\rm{ns}}$, ${G_{{\rm{gap}}}} = 1\;{\rm{ns}}$, $n = 5$, $\rho = 30{\text{%}}$.

根据图3的结果, 可以解释图1(a)和图1(b)对应复合介质为什么激发性不一样. 在图1(a)情况下, 成纤维细胞膜电导率比较小, 而M-F耦合强度比较大, 当$\rho \leqslant $ 40%时, 成纤维细胞静息电位接近心肌细胞静息电位(约为–80 mV), M-F耦合对心肌细胞的激发没有帮助, 此时成纤维细胞起到电流吸收器的作用, 所以M-F耦合导致心肌细胞激发性显著降低. 在图1(b)情况下, 成纤维细胞的膜电导率较大而M-F耦合强度较小, 当$\rho \leqslant $ 40%时, 大部分情况下成纤维细胞静息电位都在–65 mV以上, 比${E_{\rm{f}}}$低一点但比心肌细胞静息电位高很多, 这时成纤维细胞起到电流源的作用, M-F耦合有利于心肌细胞的激发, 当心肌细胞膜电位高于成纤维细胞时, 纤维细胞又会成为负担, 降低心肌细胞膜电位上升速度, 使得心肌细胞膜电位上升的幅度比无M-F耦合时小, 同时使心肌细胞不容易回到静息态而长时间处于不应态, 这就是增大成纤维细胞密度最终会导致波无法传播的原因之一.
考虑到文献[7]将${G_{\rm{f}}}$限定在[0.1, 4.0] nS范围内变化, ${G_{{\rm{gap}}}}$限定在[0, 20] nS范围内变化, 而实验上可以通过调节连接蛋白Cx43增加细胞之间的耦合强度, 且老化心脏中成纤维细胞数量与心肌细胞数量之比可以达到很高的值, 下面选择图1所示的螺旋波和混沌态为控制的初态, 将这些态从左到右、从下到上依次用(a)—(t)表示, 即$\rho = $ 10%, n = 1对应(a)态, $\rho = $ 10%, n = 2对应(b)态, $\rho = $ 40%, n = 5对应(t)初态等, 研究通过改变$D$和${G_{{\rm{gap}}}}$是否可以达到控制螺旋波和时空混沌的目的, 控制过程中其他参数与图1相同. 简单起见, 规定D的取值范围为[0.001, 0.005] ${\rm{c}}{{\rm{m}}^{\rm{2}}}{\rm{/ms}}$, D的增量为0.001 ${\rm{c}}{{\rm{m}}^{\rm{2}}}{\rm{/ms}}$; ${G_{{\rm{gap}}}}$的取值范围为[1, 15] ${\rm{nS}}$, 其增量为1 nS. 控制持续时间为3000 ms (即3 s), 如果在这个时间内, 平均膜电位差不再随时间变化, 就认为螺旋波或时空混沌已经被控制, 它们从介质中消失了. 图4与图5是分别选择与图1(a)和图1(b)相对应的态为控制初态, 采用上述控制方法在$D \text{-} {G_{{\rm{gap}}}}$参数平面上以相图方式给出的控制结果, 符号△覆盖的区域为慢可控区, 这里的慢控制是指将波态转变为无波态用时远超过200 ms, 原因是: 大部分波在控制下因传导障碍消失, 剩余的波转变成单螺旋波后漫游出系统后消失, 这就增加了控制时间. 如果剩余的波演化成螺旋波对进而转变成靶波消失, 那么控制时间可能会更长. 符号□覆盖的区域为快可控区, 这里的快控制是指将波态转变为无波态用时少于200 ms, 波因传导障碍直接消失.
图 4 与图1(a)对应的各初态下在$D \text{-} {G_{{\rm{gap}}}}$参数平面上的相图, 其中, 第1行到第4行各图参数分别对应$\rho = $ 10%, 20%, 30%, 40%; 第1列到第5列各图参数分别对应$n = 1, \;2, \;3, \;4, \;5$; ■代表不可控, Δ代表慢可控, □代表快可控
Figure4. Phase diagram in the $D \text{-} {G_{{\rm{gap}}}}$ parameter plane for different initial states showed in Fig. 1(a). The parameter $\rho $ of the panels from the first row to the fourth row equals to 10%, 20%, 30%, 40%, respectively. The parameter $n$ of the panels from the first column to the fifth column equals to 1, 2, 3, 4, 5, respectively. Black and hollow squares represent the uncontrolled point and fast control point, respectively. Hollow up-triangle represents the slow control point.

图 5 与图1(b) 对应的各初态下在$D \text{-} {G_{{\rm{gap}}}}$参数平面上的相图, 其中, 第1行到第4行各图参数分别对应$\rho =$ 10%, 20%, 30%, 40%; 第1列到第5列各图参数分别对应$n = 1, \;2, \;3, \;4, \;5$; ■代表不可控, Δ代表慢可控, □代表快可控
Figure5. Phase diagram in the $D \text{-} {G_{{\rm{gap}}}}$ parameter plane for different initial states showed in Fig. 1(b). The parameter $\rho $ of the panels from the first row to the fourth row equals to 10%, 20%, 30%, 40%, respectively. The parameter $n$ of the panels from the first column row to the fifth column equals to 1, 2, 3, 4, 5, respectively. Black and hollow squares represent the uncontrolled point and fast control point, respectively. Hollow up-triangle represents the slow control point.

从图4可以看出: 1)在高成纤维细胞密度$\rho = $ 40%情况下, 当$n = 1, \;2$时, 时空混沌几乎完全不可控(参见图4(p)和图4(q)); 在$n \geqslant 3$时, 混沌态在细胞之间耦合强度($D, {G_{{\rm{gap}}}}$)足够大时才能控制, 且可控区很小; $n = 5$时可控区则存在不连续情况, 表明时空混沌态难控制. 2) 在$\rho = $ 30%情况下, 当$D \leqslant 0.002\;{\rm{c}}{{\rm{m}}^{\rm{2}}}{\rm{/ms}}$时无论控制参数取什么值螺旋波都不可控; 当$D \geqslant 0.003\, \, {\rm{c}}{{\rm{m}}^{\rm{2}}}{\rm{/ms}}$时, 只有在$n \leqslant 2$时, 适当选择控制参数才能控制螺旋波, 当$n \geqslant 3$时, 可控范围大幅度缩小, 甚至出现螺旋波几乎不可控(参见图4(m)), 而且可控区不连续, 分布没有规律, 显示出控制结果带有偶然性, 这是因为此时大部分螺旋波都不是通过传导障碍消失的. 3) 在$\rho \leqslant $ 20%情况下, 可控区存在小范围不连续区情况, 同时也存在较大范围不连续的例外情况(参见图4(j)), 因此一般在$D, {G_{{\rm{gap}}}}$同时大于某个阈值后, 螺旋波得到控制. 4) 在给定$n$情况下, 成纤维细胞密度$\rho $越大螺旋波控制一般越困难, 因为可控区面积变小, 也存在反常情况(参见图4(g)和图4(l)). 5) 在给定$\rho $情况下, 增大$n$时螺旋波控制区既可能增大, 也可能会缩小, 没有规律说明螺旋波控制有偶然性, 也体现在两相邻的初态控制效果有很大不同.
成纤维细胞密度越大螺旋波控制越困难, 以及控制有偶然性的原因是: 图1(a)参数下的成纤维细胞主要起到电流吸收器的作用, 会极大降低介质激发性, 因此不容易通过传导障碍使螺旋波和时空混沌直接消失, 大部分情况是由波通过转变为靶波消失或转变螺旋波后漫游出边界消失, 从而使螺旋波的可控性有偶然性.
从图5可以看出: 1) 在高$\rho = $ 40%的成纤维细胞密度情况下, 增大$n$可以有效增大控制区, 当$n \leqslant 3$时, 一般需要增加心肌细胞之间的耦合强度D才可以有效控制时空混沌; 当$n \geqslant 4$时, 可以不需要提高心肌细胞之间的耦合强度$D$, 只需要使成纤维细胞之间的耦合强度${G_{{\rm{gap}}}}$大于一个阈值就可以有效消除时空混沌. 当$n \geqslant 3$时$D, {G_{{\rm{gap}}}}$一般存在一个阈值, 当$D, {G_{{\rm{gap}}}}$大于阈值, 时空混沌就可以被控制. 这些结果表明: 在成纤维细胞的膜电导率${G_{\rm{f}}}$比较大和控制前细胞间耦合强度${G_{{\rm{gap}}}}$比较小的情况下, 提高细胞耦合数有利于时空混沌的控制. 2) 在$\rho \leqslant $ 30%情况下, 当$D \leqslant 0.002\;{{{\rm{c}}{{\rm{m}}^2}} / {{\rm{ms}}}}$时, 螺旋波很难被控制, 因为只有少数参数下可控; 只有当$D \geqslant 0.003\;{{{\rm{c}}{{\rm{m}}^2}} / {{\rm{ms}}}}$时, 适当选择其他控制参数才能控制螺旋波, 且在一些情况下的可控区存在小范围的不连续情况, 但是在给定$\rho, n$的情况下, $D, {G_{{\rm{gap}}}}$一般存在一个阈值, 当$D, {G_{{\rm{gap}}}}$大于阈值, 螺旋波就可以被控制. 3) 在$n \leqslant 2$情况下, 成纤维细胞密度$\rho $越大, 螺旋波控制越困难, 因为可控区面积一般会变小. 4) 在给定$\rho $情况下, 增大$n$时螺旋波可控区一般会增大, 但是也存在可控区减小情况, 说明螺旋波控制存在偶然性, 不同的初态在相同控制参数下的控制结果也会不同. 对比图4和图5可以看出, 后者控制效果更好, 因为在图1(b)的参数下, 成纤维细胞主要起到电流源的作用, 因此M-F耦合能促进心肌细胞的激发, 容易产生传导障碍, 从而消除螺旋波和时空混沌.
为了了解控制机制, 图6给出了图2(a)、图2(b)、图2(c)三个典型斑图在控制下平均膜电位差随时间的变化, 可以看出, 在适当选取控制参数时, 经过一段时间后膜电位差的平均值不再随时间变化, 表明螺旋波和时空混沌已经被控制, 且控制螺旋波和时空混沌的用时绝大部分不到1 s. 图6(a)显示$\overline V $快速处于稳定的值, 说明不稳定螺旋波被控制. 为了看出控制机制, 图7给出了第一层复合介质在控制过程中介质区域的中间一行($j = 150$)和中间一列($i = 150$)上各细胞的膜电位随时间变化, 可以看出, 螺旋波是遇到传导障碍直接消失的, 因为图7中显示有几条倾斜粗线, 每一条都对应一个波, 粗线没有与系统边界相交, 而是终止在相邻粗线旁, 这表明波遇到传导障碍后消失. 图6(b) 显示$\bar V$先大幅振荡后趋于稳定, 这是因为在控制过程中, 螺旋波遇到传导障碍而破碎, 大部分波消失后, 剩余的波形成多个螺旋波, 最后因螺旋波漫游出边界后消失, 螺旋波消失的快慢取决于形成的波头位置是否处于边界附近, 如图8所示. 图6(c) 显示$\bar V$先小幅振荡后趋于稳定, 这是因为在时空混沌控制过程中, 时空混沌大部分由于传导障碍消失, 少部分未完全消失的波转变成螺旋波对, 进而形成靶波后消失, 时空混沌斑图消失过程如图9所示. 初态通过控制转变成靶波消失有时需要很长的控制时间, 最长可达到2 s, 图6(b)和图6(c)显示的这两种消失方式使控制具有一定的偶然性.
图 6 在不同的初态斑图和控制参数下平均膜电位差随时间的变化　(a)初态为图2(b), $D = 0{\rm{.004}}\;{\rm{c}}{{\rm{m}}^2}/{\rm{ms}}$, ${G_{{\rm{gap}}}} = {\rm{10\; ns}}$; (b)初态为图2(a), $D = 0{\rm{.003}}\;{\rm{c}}{{\rm{m}}^2}/{\rm{ms}}$, ${G_{{\rm{gap}}}} = $$ {\rm{10\; ns}}$; (c)初态为图2(c), $D = 0{\rm{.004}}\;{\rm{c}}{{\rm{m}}^2}/{\rm{ms}}$, ${G_{{\rm{gap}}}} = {\rm{12\; ns}}$
Figure6. Evolution of the average membrane potential difference for different control parameters and initial states: (a) The initial state is shown in Fig. 2(b), $D = 0{\rm{.004}}\;{\rm{c}}{{\rm{m}}^2}/{\rm{ms}}$, ${G_{{\rm{gap}}}} = {\rm{10\; ns}}$; (b) the initial state is shown in Fig. 2(a), $D = 0{\rm{.003}}\;{\rm{c}}{{\rm{m}}^2}/{\rm{ms}}$, ${G_{{\rm{gap}}}} = {\rm{10\; ns}}$; (c) the initial state is shown in Fig. 2(c), $D = 0{\rm{.004}}\;{\rm{c}}{{\rm{m}}^2}/{\rm{ms}}$, ${G_{{\rm{gap}}}} = {\rm{12\; ns}}$.

图 7 在图6(a)参数下第一层复合介质中$j = 150$这一行(a)与$i = 150$这一列(b)各细胞膜电位随时间的变化斑图
Figure7. Spatiotemporal pattern of the membrane potential of the grid points in the 150th row (a) and the 150th column (b) in the first layer composite medium for the parameters corresponding to Fig. 6(a).

图 8 在图6(b)参数下不同时刻的膜电位斑图　(a) $t = $$ {\rm{10\; ms}}$; (b) $t = {\rm{40\; ms}}$; (c) $t = {\rm{8}}0\;{\rm{ms}}$; (d) $t = {\rm{12}}0\;{\rm{ms}}$; (e) $t = $$ {\rm{18}}0\;{\rm{ms}}$; (f) $t = {\rm{22}}0\;{\rm{ms}}$; (g) $t = {\rm{32}}0\;{\rm{ms}}$; (h) $t = {\rm{38}}0\;{\rm{ms}}$; (i) $t = {\rm{44}}0\;{\rm{ms}}$
Figure8. Pattern of membrane potential at different time moments for the parameters corresponding to Fig. 6(b): (a) $t = $$ {\rm{10\; ms}}$; (b) $t = {\rm{40\; ms}}$; (c) $t = {\rm{8}}0\;{\rm{ms}}$; (d) $t = {\rm{12}}0\;{\rm{ms}}$; (e) $t = $$ {\rm{18}}0\;{\rm{ms}}$; (f) $t = {\rm{22}}0\;{\rm{ms}}$; (g) $t = {\rm{32}}0\;{\rm{ms}}$; (h) $t = {\rm{38}}0\;{\rm{ms}}$; (i) $t = {\rm{44}}0\;{\rm{ms}}$.

图 9 在图6(c)参数不同时刻的膜电位斑图　(a) $t = $$ {\rm{1}}0\;{\rm{ms}}$; (b) $t = {\rm{30\; ms}}$; (c) $t = {\rm{8}}0\;{\rm{ms}}$; (d) $t = {\rm{11}}0\;{\rm{ms}}$; (e) $t = $$ {\rm{14}}0\;{\rm{ms}}$; (f) $t = {\rm{17}}0\;{\rm{ms}}$; (g) $t = {\rm{26}}0\;{\rm{ms}}$; (h) $t = {\rm{36}}0\;{\rm{ms}}$; (i) $t = {\rm{48}}0\;{\rm{ms}}$
Figure9. Pattern of membrane potential at different time moments for the parameters corresponding to Fig. 6(c): (a) $t = $$ {\rm{1}}0\;{\rm{ms}}$; (b) $t = {\rm{30\; ms}}$; (c) $t = {\rm{8}}0\;{\rm{ms}}$; (d) $t = {\rm{11}}0\;{\rm{ms}}$; (e) $t = $$ {\rm{14}}0\;{\rm{ms}}$; (f) $t = {\rm{17}}0\;{\rm{ms}}$; (g) $t = {\rm{26}}0\;{\rm{ms}}$; (h) $t = {\rm{36}}0\;{\rm{ms}}$; (i) $t = {\rm{48}}0\;{\rm{ms}}$.

总之, 螺旋波和时空混沌控制可以通过传导障碍直接消失, 当成纤维细胞起到电流源的作用时, 螺旋波和时空混沌大多数以这种方式被控制, 这时控制效果好; 螺旋波和时空混沌控制还可以通过从初态到螺旋波或靶波的转变(简称波转变)消失, 当成纤维细胞起到电流吸收器的作用时, 螺旋波和时空混沌大多数以波转变方式消失. 如果螺旋波和时空混沌需要通过波转变方式消失, 就有可能出现增大细胞之间耦合强度, 螺旋波和时空混沌反而不能被控制的现象, 但是增大细胞之间耦合强度, 还是有利于消除螺旋波和时空混沌.

本文使用L-R相I心脏模型和被动成纤维细胞模型研究了M-F耦合对双层复合介质中螺旋波形成的影响, 发现: 在双层复合介质中M-F耦合对螺旋波形成的影响主要导致螺旋波漫游和破碎, M-F耦合对波传播的影响主要是抑制波的传播, 其中起最主要因素的是成纤维细胞的密度, 在给定$n$时增大$\rho $会导致螺旋波从稳定到漫游变化, 当$\rho = $ 40%时, 形成的螺旋波将破碎成时空混沌, 当$\rho = $ 50%时, 介质不能维持波传播, 表明M-F耦合都会降低介质的整体激发性.
此外还研究了不同初态下在双层复合介质中螺旋波和时空混沌的控制, 发现成纤维细胞的膜电导率对控制效果影响很大. 当成纤维细胞的膜电导率较小, 纤维细胞起到电流吸收器作用时, 时空混沌很难通过提高细胞之间耦合强度来消除. 除个别情况不能通过增加细胞之间耦合强度消除螺旋波外, 一般情况下, 螺旋波的控制只有通过增加心肌细胞之间耦合强度(即提高扩散系数)才能实现控制, 当$D \leqslant 0.002\, \;{\rm{c}}{{\rm{m}}^{\rm{2}}}{\rm{/ms}}$时, 不论如何增加M-F耦合强度都不可能消除螺旋波, 因为M-F耦合极大降低介质激发性, 只有增大心肌细胞之间耦合强度, 才容易出现传导障碍消除螺旋波和时空混沌.
当成纤维细胞的膜电导率较大, 成纤维细胞起到电流源作用时, 螺旋波和时空混沌都可以通过提高细胞之间耦合强度来消除, 给定$\rho, n$的情况下, 控制参数$D, {{{G}}_{{\rm{gap}}}}$一般存在一个阈值, 当$D, {{{G}}_{{\rm{gap}}}}$大于阈值, 螺旋波和时空混沌就可以被控制. 但是一般必须增加心肌细胞之间耦合强度$D$才能消除螺旋波和时空混沌, 当$n$比较大时, 也可以只通过增加M-F耦合强度就能消除时空混沌, 但不能只通过增加M-F耦合强度消除螺旋波.
通过提高细胞之间的耦合强度消除复合介质中螺旋波和时空混沌的机制有两种, 第一种机制是通过传导障碍消失; 第二种机制是先通过传导障碍使大部分波消失, 剩下的波转变为靶波或螺旋波后消失, 由于很多情况下是通过第二种机制消除螺旋波和时空混沌, 所以一些控制情况下的控制结果存在偶然性.
本文研究结果表明: M-F耦合增加了螺旋波和时空混沌的控制难度, 只有在成纤维细胞起到电流源作用的情况下, 增加细胞之间耦合强度的方法才能达到较好控制螺旋波和时空混沌的效果, 希望本文结果能为老化心脏中的心律失常控制提供有用信息.