跳转至

第 6 章:可激发系统中的波传播 (Wave Propagation in Excitable Systems)

章首引言把第 5 章的兴奋性 (excitability) 概念从单点 (space-clamped) 推进到空间延展 (spatially extended) 介质。沿神经轴突的电流问题比第 4 章的树突网络复杂得多;第 5 章已说明离子电流的电压依赖可以导致兴奋性动作电位 (action potential)。本章要指出的是,当可激发膜 (excitable membrane) 被嵌入一条非线性电缆方程 (nonlinear cable equation) 时,它可以产生沿轴突以固定速度传播的兴奋波——这一性质是 Hodgkin–Huxley 方程如此重要的原因之一。除了解释乌贼巨轴突中动作电位的产生之外,Hodgkin 与 Huxley 还展示了该动作电位能以可计算的速度沿轴突传播。然而,神经轴突只是众多空间延展可激发系统的一个例子——骨骼肌和心肌组织、视网膜、脑皮层、以及多种单细胞中都能观察到电活动或化学活动的传播波。本章系统描述这些波现象,先讨论一维电缆上电活动的传播,再简要讨论高维可激发系统中的波。第 5 章末尾也指出:兴奋介质的空间延展版本(一维雷管引线、二维干草地)允许波传播;干草地比雷管多一个恢复 (recovery) 特性——烧过的草地经几个月后会重新长出燃料,可再次燃烧。本章正是要给出这种"恢复"特性在数学上如何被形式化。

6.1 波传播概述 (Brief Overview of Wave Propagation)

6.1 节把兴奋性波按"行波 (traveling wave)"的形态做一次分类。生物系统波传播文献极多——Murray (2002)、Britton (1986)、Grindrod (1991) 等专著之外还有大量期刊文献。激发性系统 (excitable system) 中的波源自细胞自身的兴奋性——这种波的一个模型是神经轴突上动作电位的传播,另一个是草原上火焰的传播。当基础动力学是振荡 (oscillatory) 而非兴奋 (excitable)、且许多振荡单元通过扩散耦合时,所产生的行为是振荡波与周期波列 (periodic wave trains);本章只关注激发性介质中的波,把耦合振子理论留到第 12 章和第 18 章。从一开始就要强调:行波是指偏微分方程在无穷域上的一个解(显然是虚构的对象)以常速、保持固定形状地传播。在激发性系统中有两种最重要的行波:(A) 行波前沿 (traveling front)——前方 \(v\) 处于某低稳态值,后方处于某高稳态值,整个解像一个移动的"平台",在两个稳态之间做过渡;(B) 行波脉冲 (traveling pulse)——前后 \(v\) 处于同一个稳态值,中间是一个"鼓包"。这两种波可以用第 5 章的 Hodgkin–Huxley 快-慢相平面来解释:当恢复变量 (recovery variable) 固定在稳态时,快-慢相平面有两个稳定稳态 \(v = v_r\)\(v = v_e\)(即双稳 (bistable)),适当条件下存在 \(v = v_r\) 在前、\(v = v_e\) 在后的行波前沿——它像拉链一样把介质从静息"拉到"激发;然而若恢复变量 \(n\) 允许变化,则解最终会被迫返回静息,行波解就成为行波脉冲前沿与脉冲的主要区别在于:前沿情形下没有恢复(或恢复是静态的),而脉冲情形下恢复扮演关键的动力学角色。最简单的生物波传播模型之一是 Fisher 方程——它在种群生物学和生态学中广泛使用,但在生理学上关联不大,故本章不讨论(见练习 14 和 Fife, 1979)。再复杂一层是双稳方程 (bistable equation)——它有两个稳定静息点,等价于不带恢复的 FitzHugh–Nagumo 方程;该方程有行波前沿但无行波脉冲。加入恢复变量后则得到空间分布的 FitzHugh–Nagumo 方程——该方程有行波脉冲(以及其它波型)。FitzHugh–Nagumo 方程中的波传播在高于一维时仍未被完全理解。最高复杂度是 Hodgkin–Huxley 型的空间分布模型——这些方程组抵抗分析处理。

6.2 节讨论行波前沿,并以双稳方程为基础。3 个子节分别处理:连续的双稳方程 (6.2.1)、生理学上重要的髓鞘化电缆 (6.2.2)、以及离散双稳方程 (6.2.3)。

6.2.1 双稳方程 (The Bistable Equation)

6.2.1 节给出双稳方程的定义并用相平面方法研究行波前沿的存在性与唯一性。双稳方程是电缆方程 (4.18) 的一个具体形式:

\[ \frac{\partial V}{\partial t} = \frac{\partial^2 V}{\partial x^2} + f(V) \tag{6.1} \]

其中 \(f(V)\)\(0\)\(\alpha\)\(1\) 三个零点处为零,\(0 < \alpha < 1\)\(V=0\)\(V=1\) 是常微分方程 \(\mathrm{d}V/\mathrm{d}t = f(V)\) 的稳定稳态解。变量 \(V\) 被标度使 \(0\)\(1\) 成为 \(f(V)\) 的零点;在标准无量纲形式下 \(f'(0)=-1\)(参见 (4.13) 中关于被动电缆电阻在静息时斜率为 1 的约定),但这一限制常被忽略。\(f(V)\) 的两个常用例子是三次多项式

\[ f(V) = aV(V-1)(\alpha - V), \quad 0 < \alpha < 1 \tag{6.2} \]

分段线性函数(其中 \(H\) 是 Heaviside 函数)

\[ f(V) = -V + H(V - \alpha), \quad 0 < \alpha < 1 \tag{6.3} \]

(Mckean, 1970)。分段线性函数既不连续也没有三个零点,但在双稳方程的行波解研究中仍然有用——因为它是一个分析上可解的模型、保留了许多重要的定性特征。行波解是指方程 (6.1) 的一个平移不变解——它在两个稳定静息点(\(f(V)\) 的零点)之间做过渡,并以常速传播。即寻找 (6.1) 的如下形式的解

\[ V(x,t) = U(x + ct) = U(\xi) \tag{6.4} \]

其中 \(c\) 是待定波速。行波坐标 \(\xi\) 的特点是固定值以固定速度 \(c\) 移动;用 \(\xi = x + ct\) 意味着 \(c>0\) 对应从右向左传播的波(也可取 \(x-ct\))。把 (6.4) 代入 (6.1) 立即得到行波必须满足的常微分方程

\[ U_{\xi\xi} - c U_{\xi} + f(U) = 0 \tag{6.5} \]

它比原偏微分方程更易分析;为提供静息点之间的过渡,需要 \(f(U(\xi)) \to 0\)\(\xi \to \pm \infty\)。把 (6.5) 写成两个一阶方程更方便:

\[ U_{\xi} = W, \qquad W_{\xi} = cW - f(U) \tag{6.6,6.7} \]

为了求双稳方程的行波前沿解,要在 \((U, W)\) 相平面上寻找连接两个静息点 \((U, W) = (0,0)\)\((U, W) = (1,0)\) 的轨道——这种连接不同稳态的轨道叫做异宿轨道 (heteroclinic trajectory),由 \(\xi\) 参数化:\(\xi \to -\infty\) 时趋近 \((0,0)\)\(\xi \to +\infty\) 时趋近 \((1,0)\)(见图 6.2A 的虚线)。\(U=0\)\(U=1\) 处的稳态都是鞍点 (saddle)\(U=\alpha\) 处的稳态则两个特征值的实部同号\(c>0\) 时为负,\(c<0\) 时为正),即要么是节点要么是螺线 (spiral)。既然 \(U=0\)\(U=1\) 是鞍点,问题归结为:能否选择参数 \(c\) 使得从 \(U=0\) 离开(\(\xi=-\infty\))的轨道与 \(U=1\)\(\xi=+\infty\))的鞍点连接?这种数学过程叫做打靶法 (shooting),图 6.2A 给出了若干尝试轨道。首先可以确定 \(c\) 的符号。假设存在单调增\(U_\xi > 0\))的连接轨道,把 (6.5) 乘以 \(U_\xi\) 并从 \(\xi=-\infty\) 积分到 \(\xi=\infty\),得

\[ c \int_{-\infty}^{+\infty} W^2\,\mathrm{d}\xi = \int_0^1 f(u)\,\mathrm{d}u \tag{6.8} \]

——也就是说:若行波解存在,则 \(c\) 的符号与 \(f(u)\)\(u=0\)\(u=1\) 之间的面积同号。若该面积为正,则行波解把状态变量 \(U\)\(U=0\) 推到 \(U=1\)\(U=1\) 称为支配态 (dominant);在 (6.2) 和 (6.3) 两个特例中,当 \(\alpha < 1/2\)\(U=1\) 是支配态。

接下来考察 \(c=0\)\(c\) 很大时轨道的极限行为。\(c=0\) 时把 (6.5) 乘以 \(U_\xi\) 积分得

\[ \frac{W^2}{2} + \int_0^U f(u)\,\mathrm{d}u = 0 \tag{6.9} \]

——若该轨道在某 \(W\) 处达到 \(U=1\),则

\[ \frac{W^2}{2} + \int_0^1 f(u)\,\mathrm{d}u = 0 \tag{6.10} \]

这只能在 \(\int_0^1 f(u)\,\mathrm{d}u \leq 0\) 时发生,与假设矛盾,因此 \(c=0\) 的轨道不能达到 \(U=1\)。另一方面,当 \(c\) 很大时:在 \((U, W)\) 相平面上,从 \(U=0\) 离开的不稳定轨道的斜率是 \(\lambda^2 - c\lambda + f'(0) = 0\)正根,该正根总大于 \(c\)(练习 1)。设 \(K\) 是使 \(f(u)/u \leq K\) 对所有 \(0 < u \leq 1\) 成立的最小正数;\(\sigma\) 是任意固定正数。在直线 \(W = \sigma U\) 上,轨道的斜率满足

\[ \frac{\mathrm{d}W}{\mathrm{d}U} = \frac{c - f(U)}{W} \geq \frac{c - K}{\sigma} \tag{6.11} \]

——把 \(c\) 取得足够大可使 \(c - K/\sigma > \sigma\),所以一旦轨道高于 \(W = \sigma U\) 这条线,它就停留在该线上方。因此对足够大的 \(c\),从 \(U=0\) 出发的轨道从一开始就在该曲线上方,并停留在 \(W = \sigma U\) 上方,从而穿过静息点 \((U, W) = (1, 0)\)。把两个极限连起来:\(c=0\) 时的轨道通过在 \(U<1\) 处穿过 \(W=0\) 轴而错失 \(U=1\)\(c\) 很大时的轨道则通过停留在 \(U=1\) 上方而错失 \(U=1\);由解对参数的连续依赖性,存在至少一条轨道恰好命中 \((U, W) = (1, 0)\)。该 \(c\) 的值是唯一的——这是因为由 (6.11) 知 \(\mathrm{d}W/\mathrm{d}U\) 关于 \(c\) 单调递增;故对 \(c > c_0\) 不可能再存在连接轨道,\(c < c_0\) 也不可能。对大多数 \(f(V)\) 必须数值地计算行波前沿的传播速度;然而在 (6.2) 和 (6.3) 两个特例中可显式计算。分段线性情形 (6.3) 下

\[ c = \frac{1-2\alpha}{\sqrt{\alpha - \alpha^2}} \tag{6.12} \]

(见练习 4)。三次多项式情形 (6.13) 下 \(f(u) = -A^2 (u - u_0)(u - u_1)(u - u_2)\),其中零点排序为 \(u_0 < u_1 < u_2\)。寻找连接最小零点 \(u_0\) 与最大零点 \(u_2\) 的异宿连接,猜测 \(W = -B (U - u_0)(U - u_2)\)(式 (6.14))。代入 (6.5) 得

\[ B^2 (2U - u_0 - u_2) - cB - A^2 (U - u_1) = 0 \tag{6.15} \]

——作为 \(U\) 的线性函数,要使其恒为零,必须取 \(B = A/\sqrt{2}\)

\[ c = \frac{A}{\sqrt{2}} (u_2 - 2u_1 + u_0) \tag{6.16} \]

由 (6.14) 得行波解的显式形式

\[ U(\xi) = \frac{u_0 + u_2}{2} + \frac{u_2 - u_0}{2} \tanh\!\left(\frac{A}{\sqrt{2}}\,\frac{u_2 - u_0}{2}\,\xi\right) \tag{6.17} \]

——\(u_1\) 无关。当 \(u_0 = 0\)\(u_1 = \alpha\)\(u_2 = 1\) 时,波速化简为

\[ c = \frac{A}{\sqrt{2}} (1 - 2\alpha) \tag{6.18} \]

——即波速是 \(\alpha\) 的递减函数、且在 \(\alpha = 1/2\) 处改变方向;波廓为

\[ U(\xi) = \frac{1}{2}\left[1 + \tanh\!\left(\frac{A}{2\sqrt{2}}\,\xi\right)\right] \tag{6.19} \]

图 6.2B 给出该行波廓。一旦 (6.1) 的解已知,可简单地用物理参数写为

\[ V(x,t) = U\!\left(\frac{x}{\lambda_m} + c\,\frac{t}{\tau_m}\right) \tag{6.20} \]

其中 \(\lambda_m\)\(\tau_m\) 是第 4 章中的电缆空间常数和时间常数;行波的物理速度

\[ s = \frac{c \lambda_m}{\tau_m} = \sqrt{\frac{c}{2 C_m}}\,\sqrt{\frac{d}{R_m R_c}} \tag{6.21} \]

——可见波速依赖于电容、膜电阻、胞质电阻和轴突直径;速度对离子通道电导的依赖包含(但隐藏在)\(c\) 中。经验上动作电位沿轴突速度的良好估计

\[ s = \sqrt{\frac{d}{10^{-6}\,m}}\ \mathrm{mm/ms} \tag{6.22} \]

——用 \(d = 500\ \mu\mathrm{m}\)(乌贼轴突)代入得 \(s = 22.4\ \mathrm{mm/ms}\),与实测 \(s = 21.2\ \mathrm{mm/ms}\) 吻合。标度论证还可以揭示波速对其它参数的依赖。例如假设一种药物按比例封锁膜上所有离子通道、\(\rho\) 是剩余工作通道的比例,则传播速度降低 \(\sqrt{\rho}\) 倍——这是因为减少通道数后的双稳方程 \(V'' - sV' + \rho f(V) = 0\)(式 (6.23))可以通过 \(V(\xi) = U(\sqrt{\rho}\,\xi)\)\(s = c\sqrt{\rho}\) 回到原方程 (6.5)。表 6.1 列出若干可激发组织中钠通道密度 (channels/μm²):哺乳动物迷走神经(非髓鞘)110、Ranvier 节 2100、骨骼肌 205–560;其它动物如乌贼巨轴突 166–533、青蛙缝匠肌 280、电鳗电器官 550、雀鳝嗅神经 35、龙虾步行足神经 90。

阈值与稳定性 (Thresholds and Stability):双稳方程的解还满足一个比较性质 (comparison property)——任意两个解 \(u_1(x,t)\)\(u_2(x,t)\),若在某时刻 \(t=t_0\) 满足 \(u_1(x,t_0) \leq u_2(x,t_0)\),则 \(u_1(x,t) \leq u_2(x,t)\) 对所有 \(t \geq t_0\) 成立。利用比较论证可证明若干额外事实(Aronson and Weinberger, 1975):双稳方程表现出阈值行为 (threshold behavior)——若初始数据足够小,则解当 \(t \to \infty\) 时趋于 \(0\);但存在紧支撑、介于 0 和 1 之间的初始函数,使解当 \(t \to \infty\) 时趋于 \(1\)根据比较定理,任何更大的初始函数也趋于 \(1\)——这样的初始数据称为超阈值 (superthreshold)。此外,行波解在非常强的意义下稳定(Fife, 1979; Fife and McLeod, 1977)——从任何在 \(x \to -\infty\) 时介于 0 和 \(\alpha\) 之间、\(x \to +\infty\) 时介于 \(\alpha\) 和 1 之间的初始数据出发,解在大时间极限下趋近于行波解的某个相移

6.2.2 髓鞘化 (Myelination)

6.2.2 节讨论髓鞘 (myelin) 对神经纤维传播的影响,并把髓鞘化纤维建模为离散电缆。大多数神经纤维被一层叫髓鞘的脂质包覆,中间有周期性的Ranvier 节暴露(图 6.3)。髓鞘由单个施旺细胞 (Schwann cell) 构成,它被多次(约 100 次)包裹在轴突膜周围。这种包裹使有效膜电阻提高约 100 倍、膜电容降低约 100 倍——粗略数据是:\(R_m\) 对细胞膜是 \(10^3\ \Omega\,\mathrm{cm}^2\),对髓鞘是 \(10^5\ \Omega\,\mathrm{cm}^2\)\(C_m\) 对细胞膜是 \(10^{-6}\ \mu\mathrm{F/cm}^2\),对髓鞘化纤维是 \(10^{-8}\ \mu\mathrm{F/cm}^2\)。髓鞘段长度通常为 1–2 mm(接近 \(100\,d\)\(d\) 为纤维直径),Ranvier 节宽度约 \(1\ \mu\mathrm{m}\)髓鞘化纤维的传播速度比非髓鞘化纤维快——这是因为在髓鞘化段跨膜离子电流和电容电流都很小,轴突可当作一个简单电阻。动作电位不沿髓鞘化纤维连续传播,而是从一个节跳到下一个节——这种节-到-节传播叫做跳跃式 (saltatory) 传播(拉丁语 saltare,"跳跃"之意)。髓鞘损伤会损害中枢神经系统的神经冲动传递,这种病理状态叫做多发性硬化 (multiple sclerosis)——一种自身免疫病,常发于 18–40 岁青年、女性略多于男性。病理上是对大脑和脊髓白质的免疫反应导致多处脱髓鞘,但底层神经轴突和胞体通常不被破坏;髓鞘丧失减慢或停止动作电位的传递,导致肌肉疲劳、无力、或"沉重感"等症状

为建立髓鞘化纤维中电活动的模型,假设在髓鞘段忽略电容电流和跨膜离子电流,故沿髓鞘的轴向电流 (6.24)

\[ I_e = -\frac{1}{r_e}\,\frac{\partial V_e}{\partial x}, \qquad I_i = -\frac{1}{r_i}\,\frac{\partial V_i}{\partial x} \tag{6.24} \]

是常数(沿用第 4 章记号);同时假设 \(V\) 在每个 Ranvier 节内不变(节是等电位的 (isopotential)),\(V_n\) 是第 \(n\) 个节的电位。第 \(n\)\(n+1\) 节之间的轴向电流是

\[ I_e = -\frac{1}{L r_e}(V_{e,n+1} - V_{e,n}), \qquad I_i = -\frac{1}{L r_i}(V_{i,n+1} - V_{i,n}) \tag{6.25} \]

其中 \(L\) 是节间髓鞘长度。一个节上的总跨膜电流为 (6.26)

\[ \pi \mu p \!\left[C_m\,\frac{\partial V_n}{\partial t} + I_{\text{ion}}\right] = I_{i,n} - I_{i,n+1} = \frac{1}{L(r_i + r_e)}(V_{n+1} - 2V_n + V_{n-1}) \tag{6.26} \]

其中 \(\mu\) 是节长度。引入无量纲时间 \(\tau = t/(C_m R_m) = t/\tau_m\)引入无量纲空间)将 (6.26) 改写为

\[ \frac{\mathrm{d}V_n}{\mathrm{d}\tau} = f(V_n) + D(V_{n+1} - 2V_n + V_{n-1}) \tag{6.27} \]

其中 \(D = R_m / [\mu L p(r_i + r_e)]\)耦合系数 (coupling coefficient)。该方程称为离散电缆方程 (discrete cable equation)

6.2.3 离散双稳方程 (The Discrete Bistable Equation)

6.2.3 节研究离散双稳方程——方程组 (6.27) 中 \(f(V)\) 取双稳形式 (6.2) 或 (6.3)。离散双稳方程的研究连续版本 (6.1) 困难得多——虽然形式上是连续双稳方程的有限差分近似,但两者的解有显著不同的行为。证明离散系统行波解的存在性是高度非平凡的(Zinner, 1992)。然而若行波解存在,它满足特殊关系 \(V_{n+1}(\tau) = V_n(\tau - \tau_d)\):第 \(n+1\) 节的时间历程与第 \(n\) 节相同,仅延迟 \(\tau_d\);若 \(V_n(\tau) = V(\tau)\),则由 (6.27) 得 \(V(\tau)\) 必须满足时滞微分方程

\[ \frac{\mathrm{d}V}{\mathrm{d}\tau} = D\bigl(V(\tau + \tau_d) - 2V(\tau) + V(\tau - \tau_d)\bigr) + f(V(\tau)) \tag{6.28} \]

\(V(\tau)\) 足够光滑、\(\tau_d\) 足够小,把 \(V(\tau + \tau_d)\) 用 Taylor 级数

\[ V(\tau + \tau_d) = \sum_{n=0}^{\infty} V^{(n)}(\tau) \tau_d^n / n! \]

展开,则 (6.28) 近似为

\[ D\!\left(\tau_d^2 V_{\tau\tau} + \frac{\tau_d^4}{12} V_{\tau\tau\tau\tau}\right) - V_\tau + f(V) = 0 \tag{6.29} \]

——忽略 \(\tau_d^6\) 及以上项。现在设 \(\tau_d\) 很小,首阶方程是

\[ D \tau_d^2 V_{\tau\tau} - V_\tau + f(V) = 0 \tag{6.30} \]

——它有解 \(V_0(\tau) = U(c\tau)\),其中 \(D \tau_d^2 = 1/c^2\)\(U\) 是连续双稳方程 (6.5) 的行波前沿解,\(c\) 是无量纲波速。波速 \(s\) 是节间距 \(L + \mu\) 除以时间延迟 \(\tau_m \tau_d\),即

\[ s = \frac{L + \mu}{\tau_m \tau_d} = \frac{(L+\mu) c}{\tau_m \sqrt{D}} \tag{6.31} \]

对髓鞘化神经纤维 \(D = R_m / [\mu L p(r_i + r_e)]\);忽略细胞外电阻,速度的首阶近似是

\[ s = \frac{L+\mu}{\sqrt{\mu L}}\,\sqrt{\frac{c}{2 C_m}}\,\sqrt{\frac{d}{R_m R_c}} \tag{6.32} \]

——与非髓鞘化纤维相比速度增加 \(\frac{L+\mu}{\sqrt{\mu L}}\) 倍;用 \(L = 100 d\)\(\mu = 1\ \mu\mathrm{m}\) 代入,该因子为 \(10 \sqrt{d / 10^{-6} m}\)——这是很大的加速。经验上髓鞘化纤维相对非髓鞘化纤维的加速约 \(6 \sqrt{d / 10^{-6} m}\) 倍。

高阶近似:用标准正则摄动 (regular perturbation) 可得到传播速度的高阶近似。设 \(\epsilon = 1/D\),寻找 (6.29) 的解

\[ V(\tau) = V_0(\tau) + \epsilon V_1(\tau) + \cdots, \qquad \tau_d^2 = \epsilon\!\left(\frac{1}{c^2} + \epsilon^2 \tau_1 + \cdots\right) \tag{6.33,6.34} \]

把 (6.29) 按 \(\epsilon\) 展开并令系数为零,第一条方程是 (6.30),第二条是

\[ \mathcal{L}[V_1] = \frac{1}{c^2} V_1'' - V_1' + f'(V_0) V_1 = -\frac{V_0^{(4)}}{12 c^4} - \tau_1 V_0'' \tag{6.35} \]

——其中 \(\mathcal{L}[\cdot]\) 表示一个线性微分算子。目标是求 (6.35) 在无穷域上平方可积的解,使 \(V_0\) 是"好的"近似。由 \(\mathcal{L}[V_0'(\tau)] = 0\)(对 (6.30) 关于 \(\tau\) 求导一次)知 \(\mathcal{L}\) 在该空间不可逆;由 Fredholm 择一性定理 (Fredholm alternative)(Keener, 1998)知 (6.35) 有解当且仅当右端与伴随算子 \(\mathcal{L}^*\) 的零空间正交。伴随算子是

\[ \mathcal{L}^*[V] = \frac{1}{c^2} V'' + V' + f'(V_0) V \tag{6.36} \]

其零空间(\(\mathcal{L}^*[V] = 0\) 的解)由下式给出:

\[ V^*(\tau) = \exp(-c^2 \tau) \cdot V_0'(\tau) \tag{6.37} \]

由此得可解性条件

\[ \tau_1 \int_{-\infty}^{+\infty} \exp(-c^2 \tau) V_0'(\tau) V_0''(\tau)\,\mathrm{d}\tau = -\frac{1}{12 c^4} \int_{-\infty}^{+\infty} \exp(-c^2 \tau) V_0'(\tau) V_0^{(4)}(\tau)\,\mathrm{d}\tau \tag{6.38} \]

——\(\tau_1\) 可通过两个积分(解析或数值地)算出,传播速度确定为

\[ s = \frac{(L + \mu)\,c}{\tau_m \sqrt{D}} \left(1 - \frac{\tau_1 c^2}{2 D} + O\!\left(\frac{c^2}{D}\right)^{\!2}\right) \tag{6.39} \]

——这个练习在数值分析观点下也很有意思:它显示了数值离散化对传播速度的影响;该方法可用于其它带行波解的方程的数值格式(练习 20)。

传播失败 (Propagation Failure):离散方程与连续方程最显著的差异是——离散系统存在耦合阈值,低于该阈值则传播失败,而连续模型在所有耦合强度下都允许传播。由 (6.21) 容易看出,对连续电缆方程物理参数的连续变化导致传播速度的连续变化,除非直径为 0 或电阻、电容无穷大,否则速度不能被驱动到 0。离散系统是这样:耦合过弱时(即有效节间电阻过大),从激发节流向未激发节的电流可能太小而不能超过未激发节的阈值,传播不能继续。为研究传播失败,求离散方程 (6.27) 的驻立 (time-independent, i.e., \(\mathrm{d}V_n/\mathrm{d}\tau = 0\)) 解——其动机来自最大值原理比较论证:若离散双稳方程的两组初始数据在初始时有序,对应解对所有时间保持有序;由此若离散双稳方程有一个单调增驻立前沿解,则可能再有行波前沿解。离散双稳方程的驻立前沿解是满足有限差分方程

\[ 0 = D(V_{n+1} - 2V_n + V_{n-1}) + f(V_n) \tag{6.40} \]

的序列 \(\{V_n\}\),其中 \(V_n \to 1\)\(n \to +\infty\)\(V_n \to 0\)\(n \to -\infty\)。可以证明(Keener, 1987)——对任何双稳 \(f\)存在\(D^* > 0\) 使 \(D \leq D^*\) 时离散双稳方程驻立解(传播失败)。为简单理解该耦合阈值的行为,在分段线性动力学 (6.3) 下求 (6.40):由于该情形下离散方程线性,齐次解可写成某数 \(\lambda\) 的幂的线性组合

\[ V_n = A \lambda^n + B \lambda^{-n} \tag{6.41} \]

其中 \(\lambda\) 是特征多项式

\[ \lambda^2 - \left(2 + \frac{1}{D}\right)\lambda + 1 = 0 \tag{6.42} \]

的根;由此 (6.43) \(D = \lambda / (\lambda - 1)^2\)。特征方程有两个正根,一个大于 1、一个小于 1;记小于 1 的那个为 \(\lambda\)。把 \(\pm\infty\) 的条件纳入后,解可写为

\[ V_n = \begin{cases} 1 + A \lambda^n, & n \geq 0 \\ B \lambda^{-n}, & n < 0 \end{cases} \tag{6.44} \]

该表达式对 \(n = -1, 0\) 也必须满足分段线性离散双稳方程——从而

\[ D(V_1 - 2V_0 + V_{-1}) = V_0 - 1, \qquad D(V_0 - 2V_{-1} + V_{-2}) = V_{-1} \tag{6.45,6.46} \]

——其中假设 \(V_n \geq \alpha\) 对所有 \(n \geq 0\)\(V_n < \alpha\) 对所有 \(n < 0\)。把 (6.43) 代入 \(D\) 并解出 \(A\)\(B\),得 \(B = A + 1 = 1/(1+\lambda)\)。这对所有 \(n\) 都是解只要 \(V_0 \geq \alpha\);由于 \(V_0 = B = 1/(1+\lambda)\),需

\[ \frac{1}{1+\lambda} \geq \alpha \quad\Longleftrightarrow\quad \lambda \leq \frac{1-\alpha}{\alpha} \]

\(\lambda < 1\)\(D\)\(\lambda\) 的增函数,故 \(\lambda \leq (1-\alpha)/\alpha\) 当且仅当

\[ D \leq D\!\left(\frac{1-\alpha}{\alpha}\right) = \frac{\alpha(1-\alpha)}{(2\alpha - 1)^2} = D^* \tag{6.47} \]

——即当耦合足够小 (\(D \leq D^*\)) 时存在驻立解、传播不能发生\(\alpha\) 衡量介质的兴奋性 (excitability):当介质弱兴奋\(\alpha\) 接近 \(1/2\))时 \(D^*\) 很小的电阻就足以阻断传播;反之当 \(\alpha\) 小(介质强兴奋)时电阻阈值大难以阻断传播。

6.3 节把研究对象从行波前沿(连接不同稳态)推进到行波脉冲(起止于同一稳态)。行波脉冲(常称孤立脉冲,solitary pulse) 是行波坐标中起止于同一稳态的解——这种轨道叫做同宿轨道 (homoclinic orbit)。求行波脉冲有三种主要方法:(1) 用分段线性函数近似非线性函数,得到超越方程的精确解;(2) 利用不同时间尺度摄动,得到近似解析表达;(3) 用数值模拟直接求解。下面依次用三种子节说明。

6.3.1 FitzHugh–Nagumo 方程 (The FitzHugh–Nagumo Equations)

6.3.1 节是 6.3 节的核心——在 FitzHugh–Nagumo 方程上系统地构造行波脉冲。方程组为

\[ \epsilon\,\frac{\partial v}{\partial t} = \epsilon^2\,\frac{\partial^2 v}{\partial x^2} + f(v,w) \tag{6.48} \]
\[ \frac{\partial w}{\partial t} = g(v,w) \tag{6.49} \]

其中 \(\epsilon\) 是小正数。无损一般性,空间已被标度使扩散系数\(\epsilon^2\)——这意味着物理扩散系数的大小,仅是为使波前沿在新坐标下显得陡峭,便于把波作为整体研究。\(v\) 与空间扩散耦合\(w\) 与空间耦合——这是因为 \(v\) 代表膜电位、\(w\) 代表慢离子电流或门控变量。为研究行波,引入行波坐标 \(\xi = x - ct\)\(c > 0\) 是待定波速(注意:这里的 \(\xi\) 不同于 6.2 节中用的 \(x + ct\),故 \(c>0\) 在此对应从左到右的波)。把 (6.48)–(6.49) 变为常微分方程

\[ \epsilon^2 v_{\xi\xi} + c \epsilon v_\xi + f(v,w) = 0 \tag{6.50} \]
\[ c w_\xi + g(v,w) = 0 \tag{6.51} \]

分段线性模型 (Rinzel and Keller, 1973):先考虑最简单的分段线性动力学 (6.52)–(6.53)

\[ f(v,w) = H(v - \alpha) - v - w, \qquad g(v,w) = v \tag{6.52,6.53} \]

——由于动力学分段线性,精确解可以分段构造。寻找图 6.4 所示形式的解,沿 \(\xi\) 轴的波位置由 \(v(0) = v(\xi_1) = \alpha\) 固定;\(\xi_1\) 待定。记 I、II、III 分别为 \(\xi < 0\)\(0 < \xi < \xi_1\)\(\xi_1 < \xi\) 三个区域;每个区域内微分方程线性、可精确求解。三段区域解在 \(\xi = 0\)\(\xi = \xi_1\)\(v\)\(w\) 连续、\(v\) 的导数连续六个条件连接。在 I 和 III 区 \(v < \alpha\),微分方程为

\[ \epsilon^2 v_{\xi\xi} + c \epsilon v_\xi - v - w = 0 \tag{6.54} \]
\[ c w_\xi + v = 0 \tag{6.55} \]

——设 \(v = A \mathrm{e}^{\lambda \xi}\)\(w = B \mathrm{e}^{\lambda \xi}\),则 \((A, B)\) 必须满足系数矩阵 \((6.56)\);有非平凡解当且仅当行列式为零 (6.57),即 \(\lambda\)特征多项式

\[ \epsilon^2 p(\lambda) = \epsilon^2 \lambda^3 + \epsilon c \lambda^2 - \lambda + 1/c = 0 \]

(式 (6.58))的根。该多项式恰有一个负根 \(\lambda_1\),另两个根 \(\lambda_2\)\(\lambda_3\) 的实部为正。在 II 区微分方程为

\[ \epsilon^2 v_{\xi\xi} + c \epsilon v_\xi + 1 - v - w = 0 \tag{6.59} \]
\[ c w_\xi + v = 0 \tag{6.60} \]

——非齐次解为 \(w = 1\)\(v = 0\),齐次解为 \(e^{\lambda_i \xi}\) 形式的指数和。要求解当 \(\xi \to \pm \infty\) 时趋于 0,行波脉冲可写为 (6.61):在 \(\xi \geq \xi_1\) 区域用 \(\mathrm{e}^{\lambda_1 \xi}\) 单一指数、在 \(0 \leq \xi \leq \xi_1\) 区域用三个指数之和、在 \(\xi \leq 0\) 区域用 \(\mathrm{e}^{\lambda_2 \xi}\)\(\mathrm{e}^{\lambda_3 \xi}\) 之和(\(v = -c w_\xi\))。要求 \(w(\xi)\)\(v(\xi)\)\(v_\xi(\xi)\)\(\xi = 0\)\(\xi_1\) 处连续且 \(v(0) = v(\xi_1) = \alpha\)六个未知常数两个未知参数 \(c\)\(\xi_1\) 须由六个连续性条件两个约束确定。消去 \(A\)\(B_i\)\(C_i\) 后,留两个约束

\[ \mathrm{e}^{\lambda_1 \xi_1} + \epsilon^2 p'(\lambda_1) \alpha - 1 = 0 \tag{6.62} \]
\[ \frac{\mathrm{e}^{-\lambda_2 \xi_1}}{p'(\lambda_2)} + \frac{\mathrm{e}^{-\lambda_3 \xi_1}}{p'(\lambda_3)} + \frac{1}{p'(\lambda_1)} + \epsilon^2 \alpha = 0 \tag{6.63} \]

——两个未知数 \(c\)\(\xi_1\) 两个方程。更方便的做法是把 \(c\) 当作已知、\(\alpha\) 当作未知,求 \(\alpha\) 关于 \(c\) 的函数。设 \(s = \mathrm{e}^{\lambda_1 \xi_1}\),(6.63) 变为

\[ h(s) = 2 - s + \frac{p'(\lambda_1)}{p'(\lambda_2)} \mathrm{e}^{-\lambda_2 \ln s / \lambda_1} + \frac{p'(\lambda_1)}{p'(\lambda_3)} \mathrm{e}^{-\lambda_3 \ln s / \lambda_1} = 0 \tag{6.64} \]

——求 \(0 < s < 1\)\(h(s) = 0\) 的解。注意到 \(h(0) = 2\)\(h(1) = 0\)\(h'(1) = 0\)\(h''(1) = p'(\lambda_1)/\lambda_1^2 - 2\):第一个关系来自 \(\lambda_2\)\(\lambda_3\) 的实部与 \(\lambda_1\) 异号\(s \to 0\) 时指数项指数趋于 \(-\infty\) 而消失;第二个关系 \(h(1) = 0\) 来自 \(1/p'(\lambda_1) + 1/p'(\lambda_2) + 1/p'(\lambda_3) = 0\) (练习 9);后两个关系类似,留作练习(练习 9, 10)。若 \(h''(1) < 0\),则 \(s = 1\)\(h(s)\)局部极大值\(s\) 略小于 1 时 \(h(s) < 0\);由 \(h(0) > 0\)\(0 < s < 1\)必有根。当 \(\lambda_2\)\(\lambda_3\) 为实数时 \(h(s)\)\(0 < s < 1\)至多一个拐点——这是因为 \(h''(s) = 0\) 可写为 \(\mathrm{e}^{(\lambda_2 - \lambda_3) \xi_1} = c\)至多一个根。因此若 \(h''(1) < 0\)恰有一个根;若 \(h''(1) > 0\)无根。当 \(\lambda_2\)\(\lambda_3\) 为复数时唯一性不保证,但 \(h''(1) < 0\) 保证至少一个根。对定义多项式 (6.58) 关于 \(\lambda\) 求导,可得 \(h''(1) < 0\) 等价于 \(\epsilon^2 \lambda_1^2 + 2 c \epsilon \lambda_1 - 1 < 0\) ;由特征多项式 \(\epsilon^2 \lambda_1^2 - 1 = -c \epsilon \lambda_1 + \epsilon^2/(\lambda_1 c)\) ,知 \(h''(1) < 0\)\(\lambda_1 < -1/(c\sqrt{\epsilon})\)。由于 \(p(\lambda)\)\(\lambda_1\) 单调增,\(\lambda_1 < -1/(c\sqrt{\epsilon})\) \(p(-1/(c\sqrt{\epsilon})) > 0\),即

\[ c^2 > \epsilon \tag{6.65} \]

——只要 \(c > \sqrt{\epsilon}\)\(0 < s < 1\)保证有根\(s\) 一旦已知,\(\alpha\) 可由 (6.62) 得出

\[ \alpha = \frac{1 - s}{\epsilon^2 p'(\lambda_1)} \tag{6.66} \]

图 6.5 给出 (6.64) 的数值结果——波速 \(c\) 关于 \(\alpha\) 画三条 \(\epsilon\) 曲线,虚线是 \(\epsilon \to 0\) 时 (6.12) 的渐近极限关键现象是:对每个 \(\alpha\) 和足够小的 \(\epsilon\)两个行波脉冲,\(\alpha\) 没有行波脉冲。图 6.6A 是快行波脉冲\(c = 2.66\))、图 6.6B 是慢行波脉冲\(c = 0.34\)),都是 \(\alpha = 0.1\)\(\epsilon = 0.1\)\(v(\xi)\) 实线 \(w(\xi)\) 虚线。慢脉冲的振幅远小于快脉冲;一般而言快脉冲稳定(Jones, 1984; Yanagida, 1985),慢脉冲不稳定(Maginu, 1985)。注意这些波解的构造不要求 \(\epsilon\) 很小

奇异摄动理论 (Singular Perturbation Theory):利用 \(\epsilon\) 的小性可进一步提取行波脉冲信息(Keener, 1980a;另一方法见 Rauch and Smoller, 1978)。预期解在可能的情况下停留在 \(f(v,w) = 0\) 的零测线上、两个外支之间有急剧的过渡(与第 5 章图 5.15 的无扩散相图类似)。将 \(\epsilon\) 置零于 (6.48) 得外方程

\[ w_t = g(v,w), \qquad f(v,w) = 0 \tag{6.67} \]

由于 \(f(v,w)=0\) 假设有三个关于 \(v\) 的解,其中上、下两支稳定(参见 5.2 节与图 5.13),外方程 (6.67) 化为

\[ \frac{\partial w}{\partial t} = G_\pm(w) \tag{6.68} \]

——\(v = V_+(w)\) 的区域是激发区 (excited region)\(v = V_-(w)\) 的区域是恢复区 (recovering region);外方程在扩散不大时有效,但预期存在扩散很大界面 (interface) 区域,(6.68) 失效。为研究扩散很大时发生什么,重新标度空间和时间。设 \(y(t)\) 是波前沿位置,取 \(\tau = t\)\(\xi = [x - y(t)]/\epsilon\),则原方程组 (6.48)–(6.49) 变为

\[ v_{\xi\xi} + y'(\tau) v_\xi + f(v,w) = \epsilon\,\frac{\partial v}{\partial \tau} \tag{6.69} \]
\[ -y'(\tau) w_\xi = \epsilon\!\left[g(v,w) - \frac{\partial w}{\partial \tau}\right] \tag{6.70} \]

\(\epsilon = 0\)约化内方程

\[ v_{\xi\xi} + y'(\tau) v_\xi + f(v,w) = 0 \tag{6.71} \]
\[ y'(\tau) w_\xi = 0 \tag{6.72} \]

——内方程虽是偏微分方程,但 \(\tau\) 只作为参数出现,故 (6.71)–(6.72) 可如同常微分方程求解(因为行波在动坐标 \((\xi, \tau)\) 中是静止的);由 (6.72) 知 \(w\) \(\xi\) 无关(但不一定与 \(\tau\) 无关)。内方程在两个外动力学区域之间提供过渡层,要求匹配条件 \(f(v,w) \to 0\)\(\xi \to \pm\infty\)。这里用 \(y(t)\) 而非 \(ct\) 标记波前沿,故 \(y'(\tau)\)瞬时波速。把 (6.71) 视为双稳方程:对固定 \(w\),若 \(f(v,w)=0\) 有三个根、其中两个是 \(\mathrm{d}v/\mathrm{d}t = f(v,w)\) 的稳定解,则存在\(c = c(w)\) 使

\[ v'' + c(w) v' + f(v,w) = 0 \tag{6.73} \]

有连接 \(f(v,w) = 0\) 的两个稳定根的异宿轨道——这对应一个以速度 \(c\) 移动的过渡层。注意由于 \(f(v,w)=0\) 的根是 \(w\) 的函数,\(c\) 也是 \(w\) 的函数。特指地,定义 \(c(w)\) 为使 (6.73) 有满足 \(v \to V_-(w)\)\(\xi \to +\infty\))、\(v \to V_+(w)\)\(\xi \to -\infty\))的解的唯一参数值。若 \(c(w) > 0\),该过渡是从左到右上跳 (upjump);若 \(c(w) < 0\),则是下跳 (downjump)现在可以给出行波传播的一般图景:在大多空间,外动力学 (6.68) 成立;在两类外动力学的任一过渡处,\(v\) 的急剧过渡通过保持 \(w\) 连续而实现——若 \(v = V_-(w)\) 在右、\(v = V_+(w)\) 在左,则过渡层以速度 \(y'(t) = c(w)\) 移动;若 \(v = V_+(w)\) 在右、\(v = V_-(w)\) 在左,则以速度 \(y'(t) = -c(w)\) 移动,\(w\) 是过渡层内恢复变量的值。当过渡层通过空间某点时,外动力学会从一个解支切换到另一个解支

这种奇异摄动描述可用于更细致地考察行波脉冲特例。孤立脉冲的相图绘于图 6.7——一个行波脉冲由单个激发前沿后跟单个恢复背沿组成。假设右方远处介质处于静息,介质必存在静息点 \(G_-(w_+) = 0\)从左到右的波以 \(v = V_-(w_+)\) 在右、\(v = V_+(w_+)\) 在左传播,速度 \(y'(t) = c(w_+)\)必有 \(c(w_+) > 0\)。沿用推导 (6.8) 的过程可得

\[ c(w) = \frac{\int_{V_-(w)}^{V_+(w)} f(v,w)\,\mathrm{d}v}{\int_{-\infty}^{+\infty} v_\xi^2\,\mathrm{d}\xi} \tag{6.74} \]

——故 \(c(w_+) > 0\) 当且仅当

\[ \int_{V_-(w_+)}^{V_+(w_+)} f(v,w_+)\,\mathrm{d}v > 0 \]

(式 (6.75))。若 (6.75) 不成立,介质兴奋性不足以维持传播脉冲。若 \(f(v,w)\)广义 FitzHugh–Nagumo 形式,\(c(w)\)\((W_*, W^*)\) 区间有唯一零点\(W_*\)\(W^*\) 由 5.2 节定义)。激发前沿左方紧邻区域处于激发态,满足上支 \(v = V_+(w)\) 上的外动力学;因(假设)\(G_+(w) > 0\),该状态至多能保持有限时间,之后外动力学迫使另一过渡层出现。该第二个过渡层在右边的激发区与左边的恢复区之间提供过渡,以速度 \(y'(t) = -c(w)\) 移动(\(w\) 是过渡层内恢复变量的值)——负号是因为第二个过渡层必须是下跳。为形成稳态传播脉冲,上跳与下跳的速度必须相同。故下跳处的 \(w\)\(w_-\) 必须使 \(c(w_-) = -c(w_+)\)。可能没有 \(w_-\) 满足该方程——若如此,下跳必须发生在膝点 (knee),此时波称为相位波 (phase wave),因为下跳的时序完全由外动力学的相位决定,与扩散过程无关可以用标准论证证明这种波以大于某个最小速度的任意速度传播。相位波的动力学与双稳方程不同——下跳必须是鞍点与鞍-节点之间的异宿连接(膝点处双稳方程三个稳态解中的两个合并为一个),证明存在性的方法与 Fisher 方程类似——非线性 \(f(v,w)\) 在 (6.73) 中有两个简单零点(而非双稳的三个);在相位波问题中 \(f(v,w)\) 的一个零点是简单的,而是二次型——典型例子是 \(f(v,w) = v^2(1-v)\)进一步展开,说这样的波存在(见练习 16)。

总结:由奇异摄动理论——行波脉冲前方的 \(w\) 值由稳态 \(w_+\) 给出,升前沿的波速由 (6.73) 取 \(w = w_+\) 决定;前沿把 \(v\)\(v = V_-(w_+)\)(前沿前方)切到 \(v = V_+(w_+)\)(前沿后方)。波背\(w = w_-\) 出现,\(w_-\)\(c(w_-) = -c(w_+)\) 决定;波背把 \(v\)\(v = V_+(w_-)\) 切到 \(v = V_-(w_-)\)。脉冲的激发相持续时间

\[ T_e = \int_{w_+}^{w_-} \frac{\mathrm{d}w}{G_+(w)} \tag{6.76} \]

绝对不应期

\[ T_{ar} = \int_{w_-}^{w_0} \frac{\mathrm{d}w}{G_-(w)} \tag{6.77} \]

其中 \(w_0\)\(c(w) = 0\)\(w\) 值(练习 11)。该近似解称为奇异解,因为在 \(\epsilon \to 0\) 极限下解的导数变为无穷大(奇异)

6.3.2 Hodgkin–Huxley 方程 (The Hodgkin–Huxley Equations)

6.3.2 节指出行波脉冲对 HH 方程必须数值计算。最直接的方法是模拟长一维空间域上的偏微分方程;也可用打靶法 (shooting)——事实上 Hodgkin 与 Huxley 在 1952 年的论文中就用打靶法展示了 HH 方程支持行波解。打靶法也是给出行波存在性严格证明的方法(Hastings, 1975; Carpenter, 1977)。打靶论证如下。把 HH 方程写为

\[ \tau_m\,\frac{\partial v}{\partial t} = \lambda_m^2\,\frac{\partial^2 v}{\partial x^2} + f(v, m, n, h) \tag{6.78} \]
\[ \frac{\mathrm{d}w}{\mathrm{d}t} = \alpha_w(v)(1-w) - \beta_w(v) w, \quad w = n, m, h \tag{6.79} \]

寻找 \((x,t)\) 的解为平移变量 \(\xi = x/c + t\) 的函数,化为常微分方程组

\[ \frac{\lambda_m^2}{c^2}\,\frac{\mathrm{d}^2 v}{\mathrm{d}\xi^2} + f(v,m,n,h) - \tau_m\,\frac{\mathrm{d}v}{\mathrm{d}\xi} = 0 \tag{6.80} \]
\[ \frac{\mathrm{d}w}{\mathrm{d}\xi} = \alpha_w(v)(1-w) - \beta_w(v) w, \quad w = n, m, h \tag{6.81} \]

把 (6.80)–(6.81) 在静息解 \(v = 0\) 附近线性化——有四个负特征值一个正特征值忽略 \(g_K\)\(g_{Na}\) 变化可得不稳定流形的合理近似

\[ v(t) = v_0 \mathrm{e}^{\mu t} \tag{6.82} \]

其中

\[ \mu = \frac{1}{2}\!\left[\frac{\tau_m}{c^2 \lambda_m^2} + \frac{c}{\lambda_m}\sqrt{\frac{\tau_m^2}{c^2 \lambda_m^2} + 4}\right] \]

为实施打靶——选一个 \(c\) 值和接近静息点但在不稳定流形 (6.82) 上的初值,数值积分直到(很可能)膜电位变得很大;可能大正也可能大负。若能找到 \(c\) 既能产生大正又能产生大负,则用二分法逐步逼近当 \(\xi \to \infty\) 时返回静息点的同宿轨道。对 HH 方程找到行波脉冲 \(c = 3.24\ \lambda_m\ \mathrm{ms}^{-1}\);用乌贼轴突典型值(表 4.1,\(\lambda_m = 0.65\ \mathrm{cm}\))得 \(c = 21\ \mathrm{mm/ms}\)接近 Hodgkin 与 Huxley 实验测得\(21.2\ \mathrm{mm/ms}\)。Hodgkin 与 Huxley 估计乌贼轴突的空间常数 \(\lambda_m = 0.58\ \mathrm{cm}\),由此算得 \(c = 18.8\ \mathrm{mm/ms}\);他们的计算速度与实验数据非常吻合,因此他们的仅基于离子电导测量的模型被用来精确预测轴突的宏观行为定量模型能如此成功地被应用是非常罕见的。表 6.2 列出几类可激发组织的传播速度:髓鞘化神经纤维——大直径 16–20 μm 100–120 m/s、中直径 10–12 μm 60–70 m/s、小直径 4–6 μm 30–50 m/s;非髓鞘化神经纤维——中直径 3–5 μm 15–20 m/s;骨骼肌纤维 6 m/s;心脏——Purkinje 纤维 1.0 m/s、心肌 0.5 m/s;平滑肌 0.05 m/s。

6.4 节研究周期波列 (periodic wave trains)。可激发系统有两个特征——兴奋性 (excitability)不应期 (refractoriness):在系统对超阈值刺激产生大偏离响应后,存在一个不应期(该期内不能再激发),随后是恢复期(兴奋性逐渐恢复)。兴奋性恢复后再激发下一个波;下一个波的速度强烈依赖于上一次激发后恢复时间的长短——一般(但不总是)恢复时间越长,新波越快可猜想神经轴突除支持单个行波脉冲外,还支持周期动作电位波列恢复是单调过程,可预期周期波列的传播比行波脉冲——后续动作电位发生时介质未完全恢复,故 Na\(^+\) 升支比行波脉冲。速度与周期的关系叫做色散曲线 (dispersion curve)

数值计算 HH 方程的色散曲线至少有两种方法。最直接的是构造一个环 (ring)——一个一维周期边界的域,在环上启动一个单向传播的脉冲,数值求解直到解在时间上周期。然后用该波形作为稍不同长度的环的初值,再做一次计算。该方法简单,但主要缺点要求周期解稳定——色散曲线常有不稳定周期解的区域,该方法不能找到它们。当然,只有稳定解物理上可实现,所以对实验者而言这一缺点不严重。第二种方法是在行波坐标 (6.80)–(6.81)数值延拓 (continuation) 找周期解(如 AUTO 等自动延拓程序)——这种方法依赖于稳定性,整条色散曲线可以算出。可激发系统的色散曲线有典型形状(图 6.8、6.9):两条分支(快波与慢波),在绝对不应期膝点/拐角 (knee/corner) 处相汇;短于该周期周期解。快支解是典型动作电位,通常(不总是)稳定;慢支解是小振幅振荡,不稳定

6.4.1 分段线性 FitzHugh–Nagumo 方程 (Piecewise-Linear FitzHugh–Nagumo Equations)

6.4.1 节在分段线性 FitzHugh–Nagumo 方程上显式构造周期波列、从而解析地得色散曲线。图 6.9 的色散曲线是对分段线性函数 (6.52)–(6.53) 的 FitzHugh–Nagumo 系统 (6.48)–(6.49) 算出的;该计算行波脉冲的类似(Rinzel and Keller, 1973)。由于该系统分段线性,可把它的解写成三个指数之和

\[ w(\xi) = \sum_{i=1}^{3} A_i \mathrm{e}^{\lambda_i \xi} \]

(式 (6.83))在区间 \(0 \leq \xi < \xi_1\) 上,以及

\[ w(\xi) = 1 + \sum_{i=1}^{3} B_i \mathrm{e}^{\lambda_i \xi} \]

(式 (6.84))在 \(\xi_1 \leq \xi \leq \xi_2\) 上,其中 \(v = -c w_\xi\)。假设 \(v > \alpha\)\(\xi_1 \leq \xi \leq \xi_2\) 上。\(\lambda_i\)\(i=1,2,3\))是特征多项式 (6.58) 的根。要求 \(w(\xi)\)\(v(\xi)\)\(v'(\xi)\)\(\xi = \xi_1\) 处连续,\(w(0) = w(\xi_2)\)\(v(0) = v(\xi_2)\)\(v'(0) = v'(\xi_2)\)周期性),\(v(0) = v(\xi_1) = \alpha\)总共八个方程,未知数 \(A_1, \ldots, A_3\)\(B_1, \ldots, B_3\)\(\xi_1\)\(\xi_2\)\(c\)九个消去系数后留 \(\xi_1\)\(\xi_2\)\(c\) 三个未知数的两个方程

\[ \sum_{i=1}^{3} \frac{\mathrm{e}^{\lambda_i (P - \xi_1)} - 1}{p'(\lambda_i) (\mathrm{e}^{\lambda_i P} - 1)} + \epsilon^2 \alpha = 0 \tag{6.85} \]
\[ \sum_{i=1}^{3} \frac{\mathrm{e}^{\lambda_i P} - \mathrm{e}^{\lambda_i \xi_1}}{p'(\lambda_i) (\mathrm{e}^{\lambda_i P} - 1)} + \epsilon^2 \alpha = 0 \tag{6.86} \]

其中 \(P = \xi_2/c\)注意:由于两个方程三个未知数,(6.85)–(6.86) 定义了一族周期波(由周期或波速参数化);该波族周期与速度的关系就是色散曲线。图 6.9 给出 \(\alpha = 0.1\) 时若干 \(\epsilon\) 取值下的色散曲线;改变 \(\alpha\) 对该图定性影响很小。虚线显示 \(\epsilon \to 0\)上支(快波)的极限行为。该图的关键是——存在快波与慢波大波长极限下周期波趋近图 6.5 表示的孤立行波脉冲(练习 12)。事实上周期解看起来很像是均匀间隔的(截断的)孤立脉冲的周期重复。分段线性 FitzHugh–Nagumo 系统的色散曲线典型于可激发介质的色散曲线——快、慢支在拐角相汇;拐角位置一般依赖于介质的兴奋性(这里即参数 \(\alpha\))和时间尺度比 \(\epsilon\)

6.4.2 奇异摄动理论 (Singular Perturbation Theory)

6.4.2 节用奇异摄动理论一般 FitzHugh–Nagumo 系统色散曲线快支\(\epsilon \to 0\) 极限)。周期波由一系列上跳与下跳组成、其间是外动力学区域;图 6.10 绘出周期波列的相图。为周期,若 \(w_+\) 是上跳处恢复变量的值、速度为 \(c(w_+)\),则下跳处恢复变量的值 \(w_-\) 必须满足 \(c(w_+) = -c(w_-)\)。在激发支上花费的时间为

\[ T_e = \int_{w_+}^{w_-} \frac{\mathrm{d}w}{G_+(w)} \tag{6.87} \]

恢复支上花费的时间为

\[ T_r = \int_{w_-}^{w_+} \frac{\mathrm{d}w}{G_-(w)} \tag{6.88} \]

色散曲线就是速度 \(c(w_+)\) 与周期

\[ T = T_e + T_r \tag{6.89} \]

的关系,以 \(w_+\) 参数化。该近似色散曲线(数值地)算得,即图 6.9 的虚线色散曲线的慢支也可用摄动方法求得——由于速度是 \(\epsilon\) 阶小量,正则摄动适用;具体展开超出本书范围,感兴趣的读者可参 Dockery and Keener (1989)。一般地,慢周期解不稳定(Maginu, 1985),因此物理意义比快解。行波解的稳定性理论超出本书范围

6.4.3 运动学 (Kinematics)

6.4.3 节提出运动学理论 (kinematic theory)不是所有波都是周期的——可以存在不规则间隔不同速度的动作电位波列。运动学理论是一种不追踪脉冲结构细节追踪单个动作电位进度的理论(Rinzel and Maginu, 1984)。最简单的运动学理论局部解读色散曲线。已知稳定周期波列的速度关于周期的函数 \(c = C(T)\)。设波列由动作电位组成,第 \(n\) 个动作电位在 \(t_n(x)\) 时刻到达位置 \(x\)记录到达时间

\[ \frac{\mathrm{d}t_n}{\mathrm{d}x} = \frac{1}{c} \tag{6.90} \]

——配以 \(c = C(t_n(x) - t_{n-1}(x))\),意识到 \(t_n(x) - t_{n-1}(x)\) 是位置 \(x\) 处介质感受到的波列瞬时周期更精细的运动学理论可由 FitzHugh–Nagumo 方程的奇异解推导。推导时假设恢复总通过相位波实现,恢复值为 \(W_*\)。设第 \(n\) 个动作电位的前沿速度为 \(c(w_n)\),对应过渡层内恢复值 \(w_n\)记录到下一动作电位的时间

\[ t_{n+1}(x) - t_n(x) = T_e(x) + T_r(x) = \int_{w_n}^{W_*} \frac{\mathrm{d}w}{G_+(w)} + \int_{W_*}^{w_{n+1}} \frac{\mathrm{d}w}{G_-(w)} \tag{6.91,6.92} \]

对 (6.92) 关于 \(x\) 求导,得 \(w_{n+1}(x)\) 的微分方程

\[ \frac{1}{G_-(w_{n+1})}\,\frac{\mathrm{d}w_{n+1}}{\mathrm{d}x} = \frac{1}{G_+(w_n)}\,\frac{\mathrm{d}w_n}{\mathrm{d}x} + \frac{1}{c(w_{n+1})} - \frac{1}{c(w_n)} \tag{6.93} \]

——由此可追踪 \(w_{n+1}\) 关于 \(x\) 的函数(已知 \(w_n(x)\)),重建\(n+1\) 个动作电位前沿的速度与到达时间。该公式对 FitzHugh–Nagumo 方程有用,但更广泛地可推广:因前沿速度与前沿内 \(w\)一一对应,可把这些函数用前沿速度表示为

\[ t_{n+1}(x) - t_n(x) = A(c_n) + t_r(c_{n+1}) \tag{6.94} \]

其中 \(A(c_n) = T_e\) 是第 \(n\) 个上冲后的动作电位持续时间 (APD)\(t_r(c_{n+1}) = T_r\) 是第 \(n+1\) 个上冲之前的恢复时间。对该守恒律关于 \(x\) 求导,得第 \(n+1\) 个前沿速度关于第 \(n\) 个前沿速度的微分方程

\[ t_r'(c_{n+1}) \frac{\mathrm{d}c_{n+1}}{\mathrm{d}x} = \frac{1}{c_{n+1}} - \frac{1}{c_n} - A'(c_n) \frac{\mathrm{d}c_n}{\mathrm{d}x} \tag{6.95} \]

——该公式的优势是 \(A\)\(t_r\) 可由其它途径(如实验数据)得到。一般认为动作电位持续时间函数地依赖于上一动作电位的速度,下一动作电位的速度函数地依赖于上一动作电位结束以来的恢复时间——这些假设是合理的。故模型 (6.95) 适用范围超过 FitzHugh–Nagumo 上下文。:这一思路已用于 Beeler–Reuter 动力学**(Courtemanche et al., 1996)。

6.5 节把一维电缆分析推广到高维介质。不是所有可激发介质都能被看作一维电缆——也不是所有传播活动都是一维的。一维描述不够的组织包括骨骼肌和心肌组织、视网膜、脑皮层等。理解这些组织中的通讯与信号传递需要的数学分析比一维电缆更复杂。研究二维、三维波传播时,一个诱惑是把一维电缆方程扩展到高维——把一阶空间导数换为空间梯度、二阶导数换为Laplacian 算子。本章讨论的所有模型确实都是这种类型;然而这种替换并不总是合适。某些细胞(如 Xenopus 卵母细胞 (frog eggs)大到足以在单个细胞内维持化学活动波。这罕见——大多数正常生理情况下的波都服务于细胞间通讯。对单细胞内的化学波,合理的第一猜测是空间耦合化学扩散实现。若局部化学动力学由 \(\partial u/\partial t = k f\) 描述,则带空间耦合的动力学为

\[ \frac{\partial u}{\partial t} = \nabla \cdot (D \nabla u) + k f \tag{6.96} \]

——\(D\) 是化学物种的标量扩散系数,\(\nabla\)三维梯度算子。时间常数 \(k\)(单位为 \(\mathrm{time}^{-1}\))的引入使 \(f\)\(u\) 同量纲对许多细胞类型——细胞间通讯通过间隙连接 (gap junctions)紧邻细胞间实现,故扩散空间均匀。这时,第一猜测(或希望)是所描述现象的长度常数远大于典型细胞大小,可用均匀化 (homogenization)有效扩散系数 \(D_e\)(如第 8 章所述)。那么 (6.96) 配以 \(D_e\) 就是合理的模型;注意没有先验理由相信细胞耦合是各向同性的或 \(D_e\)标量长度常数远大于细胞大小、均匀化近似合理,那么我们就只能痛苦地研究离散耦合细胞间的通讯——第 7 章给出了这种研究的描述。对电活性细胞(如心细胞、肌细胞),情况更复杂——信号是跨膜电位,确定它需同时知道胞内与胞外电位。因此空间耦合能由跨膜电位的 Laplacian 表示。第 12 章讨论心肌组织中的波时会处理这一问题。更复杂的是神经网络——轴突可延伸长距离、细胞可与许多(远多于最近邻)其他细胞连接,故空间耦合非局部 (nonlocal)总之——(6.96) 是研究高维波的一个有趣的起点,但绝不明显所有研究结果都适用于细胞介质。

6.5.1 传播前沿 (Propagating Fronts)

6.5.1 节处理高维介质中的平面波弯曲波

平面波 (Plane Waves):高维介质中最简单的波是平面波。设典范问题 \(U'' + c_0 U' + f(U) = 0\)(式 (6.97))是双稳的、有行波前沿解 \(U(\xi)\)\(c_0\) 唯一、依赖于 \(f\)。为求 (6.96) 的平面波解,假设 \(u\)单变量 \(\xi = \mathbf{n} \cdot \mathbf{x} - ct\) 的函数,\(\mathbf{n}\)沿波前沿传播方向的单位向量。在行波坐标 \(\xi\) 中,时间导数 \(\mathrm{d}/\mathrm{d}t\) 替换为 \(-c \,\mathrm{d}/\mathrm{d}\xi\),空间梯度算子 \(\nabla\) 替换为 \(\mathbf{n} \,\mathrm{d}/\mathrm{d}\xi\),故控制方程化为

\[ (\mathbf{n} \cdot D \mathbf{n}) u'' + c u' + k f(u) = 0 \tag{6.98} \]

把 (6.98) 与典范方程 (6.97) 比较,简单重新标度得

\[ u(\mathbf{x}, t) = U\!\left(\mathbf{n} \cdot \mathbf{x} - c t\,\ell(\mathbf{n})\right) \tag{6.99} \]

其中 \(c = c_0 k \ell(\mathbf{n})\)方向依赖的速度、 \(\ell(\mathbf{n}) = \sqrt{\mathbf{n} \cdot D \mathbf{n} / k}\)方向依赖空间常数

曲率波 (Waves with Curvature):二维或三维介质中的波前沿预期是平面波——它们通常特定位置被启动故可以是圆形;介质可能结构不均匀或有复杂几何,所有这些都使波前沿带曲率曲率对可激发介质中波前沿的传播起重要作用——直观地说,向内收缩的圆形波前沿由于不同部分都在激发同一点,故正前方的激发应平面波更快向外扩展的圆形波前沿则平面波(因为激发努力被分散)。尽管曲率效应在很多情境中广为人知,本章关心的是其定量描述。6.5.1 节推导程函-曲率方程 (eikonal-curvature equation)——目的是显示曲率对波前沿速度的贡献。程函-曲率方程在许多生物学情境中被使用:可激发 Belousov–Zhabotinsky 试剂(Foerster et al., 1989; Keener, 1986; Keener and Tyson, 1986; Tyson and Keener, 1988; Ohta et al., 1989),Xenopus 卵母细胞中的 Ca\(^{2+}\) 波(Lechleiter et al., 1991b; Sneyd and Atri, 1993; Jafri and Keizer, 1995),心肌组织研究(Keener, 1991a; Colli-Franzone et al., 1990, 1993;也见第 12 章)。程函-曲率方程在其它科学领域也有长久历史——晶体生长(Burton et al., 1951)、火焰前沿传播(Frankel and Sivashinsky, 1987, 1988)。程函-曲率方程的推导使用奇异摄动理论的标准数学论证。关键观察是:(6.96) 内部隐藏了双稳方程 (6.97);要探索的思路是——在某个待定动坐标系中,(6.96) 可由 (6.97) 良好近似。目标是把 (6.96) 改写为形如 (6.97)动坐标系。在三维中需三个空间坐标——其中一个局部正交于波前沿,另两个是描述波前沿面的坐标。假设 \(u\) 近似与波前沿坐标无关对变量做标度使对第一个坐标的导数最重要、其它导数次要。由此可知道坐标系必须如何运动才能保持自己是一个波前沿坐标系——这一坐标系运动的定律就是程函-曲率方程首先,引入一般(尚待定)动坐标系

\[ \mathbf{x} = \mathbf{X}(\boldsymbol{\xi}, \tau), \qquad t = \tau \tag{6.100} \]

链式法则

\[ \frac{\partial}{\partial \xi_i} = \frac{\partial X^j}{\partial \xi^i}\,\frac{\partial}{\partial x^j}, \qquad \frac{\partial}{\partial \tau} = \frac{\partial}{\partial t} + \frac{\partial X^j}{\partial \tau}\,\frac{\partial}{\partial x^j} \tag{6.101} \]

求和约定重复指标从 1 到 3 求和,除非另有说明)于是

\[ \frac{\partial}{\partial x^i} = \alpha_{ij}\,\frac{\partial}{\partial \xi_j}, \qquad \frac{\partial}{\partial t} = \frac{\partial}{\partial \tau} - \frac{\partial X^j}{\partial \tau}\,\alpha_{jk}\,\frac{\partial}{\partial \xi_k} \tag{6.102} \]

其中 \(\alpha_{ij}\)\(\partial X^j/\partial \xi^i\)(坐标变换 (6.100) 的 Jacobian)的矩阵的元。 \(\xi_1\) 视为正交于 \(u\)等值面的坐标,\(\xi_2\)\(\xi_3\)动等值面的坐标。定义切向量

\[ \mathbf{r}_i = \partial \mathbf{X}/\partial \xi^i \]

\(i=1,2,3\))和法向量 \(\mathbf{n}_i = \mathbf{r}_j \times \mathbf{r}_k\)\(i \neq j, k\)\(j < k\))。不失一般性\(\mathbf{r}_1 = \sigma (\mathbf{r}_2 \times \mathbf{r}_3)\) ,故 \(\mathbf{r}_1\) 始终正交于 \(u\) 的等值面,\(\sigma\)任意的标度因子。\(\mathbf{r}_2\)\(\mathbf{r}_3\) 切于动等值面,但未必正交——虽然可强制 \(\mathbf{r}_2\)\(\mathbf{r}_3\) 正交,在动曲面上构造这种坐标系通常非常困难优选一个对坐标系无额外限制的运动方程——因为运动应与描述曲面的坐标系无关可显式计算 \(\alpha_{ij}\)——由 Cramer 法则(练习 23)得

\[ \alpha_{ij} = \frac{(n_j)_i}{\mathbf{r}_j \cdot \mathbf{n}_j} \quad\text{(no summation)} \tag{6.103} \]

其中 \((n_j)_i\) 是第 \(j\) 个法向量 \(\mathbf{n}_j\) 的第 \(i\) 分量。现在写出完整的变量变换。计算(把 \(\alpha_{ij}\) 视为 \(x\) 的函数)得

\[ \frac{\partial u}{\partial t} = \frac{\partial u}{\partial \tau} - \frac{\partial X^j}{\partial \tau}\,\alpha_{jk}\,\frac{\partial u}{\partial \xi_k} \tag{6.104} \]
\[ \nabla^2 u = \alpha_{ip} \alpha_{iq}\,\frac{\partial^2 u}{\partial \xi_p \partial \xi_q} + \frac{\partial \alpha_{ip}}{\partial x_i}\,\frac{\partial u}{\partial \xi_p} \tag{6.105} \]

把 (6.96) 用新变量改写(\(D\)常数标量情形)得

\[ 0 = D \alpha_{ip} \alpha_{iq}\,\frac{\partial^2 u}{\partial \xi_p \partial \xi_q} + D \frac{\partial \alpha_{ip}}{\partial x_i}\,\frac{\partial u}{\partial \xi_p} - \left(u_\tau - \frac{\partial X^j}{\partial \tau}\,\alpha_{jk}\,\frac{\partial u}{\partial \xi_k}\right) + k f(u) \tag{6.106} \]

现在调用两个重要假设——\(\xi_1\) 方向上的空间变化尺度短于 \(\xi_2\)\(\xi_3\) 方向上的空间变化尺度。量化为:存在小参数 \(\epsilon\)\(\alpha_{j1} = O(1)\)\(\alpha_{jk} = O(\epsilon)\) 对所有 \(j\)\(k \neq 1\)此外假设 \(\epsilon\)首阶 \(u\) \(\xi_2\)\(\xi_3\)\(\tau\) 无关。因此 (6.106) 中不是所有项都同等重要——计入 \(\alpha_{ij}\)\(\epsilon\) 依赖后,(6.106) 化为

\[ D |\alpha|^2 \frac{\partial^2 u}{\partial \xi_1^2} + \left(D \nabla \cdot \alpha + \frac{\partial X}{\partial \tau} \cdot \alpha\right) \frac{\partial u}{\partial \xi_1} + k f(u) = O(\epsilon) \tag{6.107} \]

其中 \(\alpha\)分量为 \(\alpha_{j1}\) 的向量,因而(由 (6.103))正比于法向量 \(\mathbf{n}_1\)。(6.107) 左端所有项都比 \(\epsilon\) 。这里看到一个形如双稳方程 (6.97) 的方程。若 (6.107) 的系数为常数,可令

\[ \frac{\partial \mathbf{X}}{\partial \tau} \cdot \alpha + D \nabla \cdot \alpha = k c_0 \tag{6.108} \]

同时要求 \(D |\alpha|^2 = k\)(6.108) 告诉我们坐标系应如何运动;由于 \(\alpha\) 正比于 \(\mathbf{n}_1\)\(D |\alpha|^2 = k\) 决定了正交于波前沿的坐标尺度——即波前沿的厚度实际上 (6.107) 的系数是常数,这两个要求过度约束完整坐标变换 \(\mathbf{X}(\boldsymbol{\xi}, \tau)\)。为克服该困难,假设波前沿与坐标系都是空间上慢变的,把 (6.108) 解读为只确定坐标系中线的运动(前沿梯度最大处)、确定整个坐标系的运动。(6.108) 正是所求——描述动作电位前沿运动的方程,称为程函-曲率方程对数值模拟几乎无用——数值算法可靠地模拟该方程极其困难代替地,引入一个函数 \(S(\mathbf{x}, t)\) 作为前沿的指示函数(视 \(S\)冲击位置确定器):\(S(\mathbf{x}, t) > 0\) 时介质激活\(S(\mathbf{x}, t) < 0\) 时介质静息。取 \(\alpha\) 沿向前的波前沿运动方向,即 \(\alpha = -\sqrt{k/D} \nabla S/|\nabla S|\)(用 \(|\alpha| = \sqrt{k/D}\));\(S(\mathbf{x}, t)\)零等值面表示波前沿位置,故 \(S\) 沿前沿为常数,由 \(0 = \nabla S \cdot \mathbf{X}_t + S_t\)\(\mathbf{X}_t \cdot \alpha = \sqrt{k/D}\,S_t/|\nabla S|\)。故

\[ S_t = |\nabla S| c_0 \sqrt{Dk} + D |\nabla S| \nabla \cdot \frac{\nabla S}{|\nabla S|} \tag{6.109} \]

指示函数 \(S(\mathbf{x}, t)\) 决定界面运动的方法称为水平集方法 (level set method)(Osher and Sethian, 1988)——强大易实现。(6.109) 称为程函-曲率方程——每个项都有物理解释。忽略扩散项则得程函方程 (eikonal equation)

\[ \frac{\partial S}{\partial t} = |\nabla S|\,c_0 \sqrt{Dk} \tag{6.110} \]

\(R\)\(S(\mathbf{x}, t)\) 的等值面、\(\mathbf{n}\) 是该面某点的单位法向量,则 (6.110) 蕴含 \(R\) 在法方向的速度 \(R_t \cdot \mathbf{n}\) 满足

\[ R_t \cdot \mathbf{n} = c_0 \sqrt{Dk} \tag{6.111} \]

——即前沿沿法方向 \(\mathbf{n}\)速度 \(c = c_0 \sqrt{Dk}\) 运动。(6.111) 是几何惠更斯 (Huygens) 构造的基础,但 (6.110) 或 (6.111) 的数值积分充满困难——特别是尖点 (cusp) 奇异性,指示函数 \(S(\mathbf{x}, t)\)有限时间(通常很快)变得病态。(6.109) 右端第二项曲率修正 (curvature correction)——因为项 \(\nabla \cdot (\nabla S/|\nabla S|)\) 在三维空间是平均曲率的两倍、在二维空间是曲率(练习 24)。程函-曲率方程可写为

\[ R_t \cdot \mathbf{n} = c_0 \sqrt{Dk} - D \kappa \tag{6.112} \]

\[ \tau R_t \cdot \mathbf{n} = c_0 \ell - \ell^2 \kappa \tag{6.113} \]

其中 \(\kappa\)曲率(二维)或平均曲率的两倍(三维)、\(\ell = \sqrt{D/k}\)空间常数\(\tau = 1/k\)时间常数尽管曲率修正通常只是法速度的修正,物理稳定性原因使它重要——防止奇异性形成。曲率修正的符号使带涟漪的前沿逐渐平滑平面波Belousov–Zhabotinsky 试剂实验验证了速度与曲率的该关系——例如 Foerster et al. (1988) 在旋转螺旋波不同位置和两个螺旋波的交点(从而不同符号的曲率)测得速度和曲率,发现法速度与曲率的关系良好地斜率为快速反应物种扩散系数直线**近似。

6.5.2 空间模式与螺旋波 (Spatial Patterns and Spiral Waves)

6.5.2 节研究高维介质中螺旋波 (spiral waves) 的产生与性质。已知波前沿如何在可激发介质中传播,接下来要确定可能产生的空间模式。最常见的模式由单个源产生并向外扩展;若介质足够大可同时存在多个波前沿,则称为靶模式 (target patterns)——靶模式需要周期源,故不能均匀非振荡介质中存在。第二种空间模式是螺旋波 (spiral waves)——螺旋波需要周期源来维持,因它们通常自持 (self-sustained)由于螺旋波自持,它们通常出现在病理状态——也就是说,对一个依赖信号忠实传播的系统,被自持模式接管通常是好事。因此心脏上的螺旋波致命、皮层上的螺旋波可能导致癫痫发作、视网膜或视觉皮层上的螺旋波可能导致幻觉。螺旋波在可激发介质中的著名例子是 Belousov–Zhabotinsky 反应(Winfree, 1972, 1974),在这里讨论。螺旋波的数学讨论集中于带兴奋性动力学和扩散耦合的微分方程组的周期解具体例子是二维带扩散耦合的 FitzHugh–Nagumo 方程

\[ \epsilon \frac{\partial v}{\partial t} = \epsilon^2 \nabla^2 v + f(v,w) \tag{6.114} \]
\[ \frac{\partial w}{\partial t} = g(v,w) \tag{6.115} \]

首阶奇异摄动分析(\(\epsilon = 0\)建议把域分为两个——其中外动力学

\[ \frac{\partial w}{\partial t} = G_\pm(w) \tag{6.116} \]

成立(沿用 6.3.1 节记号)。\(\partial w/\partial t = G_+(w)\) 的区域是激发区\(\partial w/\partial t = G_-(w)\) 的区域是恢复区。分隔两区的是移动界面 (interface)——其中 \(v\) 急剧变化(空间与时间常数 \(\epsilon\)),故扩散重要\(w\) 保持基本常数。在空间任一点解应关于时间周期——故在大半径处(波前沿近于平面)解应位于色散曲线上

关于界面应如何运动的第一猜测是——界面近于平面、故有与平面波相同的速度

\[ \mathbf{R}_t = c(w) \mathbf{n} \tag{6.117} \]

其中 \(\mathbf{R}\)界面的位置向量、\(\mathbf{n}\)\(\mathbf{R}\)单位法向量\(c(w)\)作为 \(w\) 函数的平面波速度。为看出程函方程的蕴涵,假设螺旋界面是曲线 \(\mathbf{R}\)

\[ X(r,t) = r \cos(\theta(r) - \omega t), \qquad Y(r,t) = r \sin(\theta(r) - \omega t) \tag{6.118} \]

界面是\(r\) 参数化的曲线,故切向量\((X_r, Y_r)\)。由此计算

\[ \mathbf{R}_t = \begin{pmatrix} -\omega r \sin(\theta - \omega t) \\ \omega r \cos(\theta - \omega t) \end{pmatrix} \]

(式 (6.119))以及

\[ \sqrt{1 + r^2 \theta'^2} \mathbf{n} = \begin{pmatrix} -\sin(\theta - \omega t) - r \theta' \cos(\theta - \omega t) \\ \cos(\theta - \omega t) - r \theta' \sin(\theta - \omega t) \end{pmatrix} \]

(式 (6.120))故程函方程变为

\[ c(w) \sqrt{1 + r^2 \theta'^2} = \omega r \tag{6.121} \]

积分得

\[ \theta(r) = \rho(r) - \tan(\rho(r)), \qquad \rho(r) = \sqrt{r^2/r_0^2 - 1} \tag{6.122} \]

其中 \(r_0 = c/\omega\),故界面为

\[ X = r_0 \cos s + r_0 \rho(r) \sin s, \qquad Y = r_0 \sin s - r_0 \rho(r) \cos s \tag{6.123} \]

\(s = \rho(r) - \omega t\)。该界面是半径 \(r_0\) 的圆的渐开线 (involute)渐开线从圆上展开的线段端点的轨迹)。该解作为螺旋解有显著困难——最显著的是它仅在 \(r \geq r_0\) 时存在\(r_0\) 任意但正,故该螺旋绕一个有限大小的"洞"旋转。旋转频率与色散曲线一致的要求决定。螺旋有波长 \(2\pi r_0\)周期 \(2\pi/\omega\),故 \(c = r_0 \omega\)。由于色散曲线通常有膝点——短于一定波长周期波,存在使该关系能成立半径下界。数值研究表明螺旋没有这样的内核半径下界、并表明无中心洞的介质有唯一的螺旋频率。不幸地仅用程函方程不能揭示唯一频率的选定方式,故需另辟蹊径——程函-曲率方程

应用程函-曲率方程旋转螺旋波:假设波前沿可写为 (6.118) 的形式,故曲率

\[ \kappa = \frac{X' Y'' - Y' X''}{(X'^2 + Y'^2)^{3/2}} = \frac{\psi'}{(1 + \psi^2)^{3/2}} + \frac{\psi}{r (1 + \psi^2)^{1/2}} \tag{6.124} \]

其中 \(\psi = r \theta'(r)\) 称为形函数程函-曲率方程 (6.113) 变为

\[ r \frac{\mathrm{d}\psi}{\mathrm{d}r} = (1 + \psi^2)^{1/2} \left[\frac{r c(w)}{\epsilon} \sqrt{1 + \psi^2} - \frac{\omega r^2}{\epsilon} - \psi\right] \tag{6.125} \]

假设 \(w\) 沿螺旋前沿为常数,(6.125) 可由\(r = \infty\) 打靶数值求解。\((r, \psi)\) 平面若干轨道示于图 6.11。(6.125) 的轨道是刚性的 (stiff)——即对 \(r\),曲线 \(c(1 + \psi^2)^{1/2} - \omega r = 0\)强吸引子刚性可由变量 \(\phi = r \psi / \sqrt{1 + \psi^2}\) 改写 (6.125) 清楚地看出

\[ \frac{\epsilon}{r}\,\frac{\mathrm{d}\phi}{\mathrm{d}r} = c - \omega \sqrt{r^2 - \phi^2} \tag{6.126} \]

——\(\epsilon\) 导数项。\(r = \infty\) 积分,(6.125) 的轨道在原点附近通过向上向下爆破趋近原点。原点鞍点——若参数恰当地选择,轨道趋近 \(\psi = 0\)。故存在 \(\omega\)\(c\) 唯一关系 \(\omega / \epsilon = F(c/\epsilon)\) 使轨道一直到达 \(r = 0\)\(\psi = 0\)注意变量重标 \(r \to \alpha r\)\(c \to c/\alpha\)\(\omega \to \omega/\alpha^2\) 保持 (6.125) 不变,故使轨道趋近原点鞍点\(\omega\)\(c\) 关系

\[ \frac{\omega}{\epsilon \alpha^2} = F\!\left(\frac{c}{\epsilon \alpha}\right) \tag{6.127} \]

——由此 \(F(c/\alpha) = (1/\alpha^2) F(c)\),从而 \(F(x) = m_* x^2\)(某常数 \(m_*\)),即

\[ \omega = \frac{c^2 m_*}{\epsilon} \tag{6.128} \]

数值地确定 \(m_* = 0.330958\)。图 6.12 给出该螺旋前沿的一个例子(实线),与同一参数下的渐开线螺旋 (6.123)(虚线比较;该图 \(c/\epsilon = 3.0\)到此我们得到一族可能螺旋轨道——它们有正确的大 \(r\) 渐近行为\(r \to 0\)趋近 \(\psi = 0\)。该族由速度 \(c\) 参数化。为决定该族中哪一个是正确的螺旋前沿,还需螺旋是周期波——即必须满足色散关系这两个要求唯一地决定螺旋性质为看出为何如此,图 6.13 绘出临界曲线 (6.128) 和近似色散曲线 (6.89);该图使用分段线性动力学 \(f(v,w) = H(v - \alpha) - v - w\)\(g(v,w) = v - \gamma w\)\(\alpha = 0.1\)\(\gamma = 0\)\(\epsilon = 0.05\)

更多关于螺旋的讨论 (More About Spirals):本节对高维波的讨论仅是冰山一角(或螺旋之巅),仍有许多未解有趣问题。生理系统中观察到的螺旋确实共享某些定性相似性,但细节当然不同这里讨论的 FitzHugh–Nagumo 方程仅展示了一般可激发系统的定性行为,故很少定量关联。其它生理系统可能由其它动力学支配——例如,皮层传播性抑制 (spreading cortical depression) 的模型已由 Tuckwell and Miura (1978; Miura, 1981) 提出,数值模拟已显示旋转螺旋;然而这些方程的详细数学研究尚未给出类似地Ca\(^{2+}\) 螺旋波已在某些细节上数值研究(第 7 章),但分析上未被很好理解。对 FitzHugh–Nagumo 系统的解析计算基于奇异摄动理论、数学上不严格。事实上目前还没有对可激发介质中螺旋波存在性严格证明虽然这里给出的近似解已知渐近有效的,螺旋核心的结构不正确这一问题已由 Pelce and Sun (1991) 和 Keener (1992) 处理——针对单扩散变量的 FitzHugh–Nagumo 模型,以及 Keener (1994) 和 Kessler and Kupferman (1996) 处理——针对双扩散变量的 FitzHugh–Nagumo 模型,适用于化学反应系统第二个关切是螺旋的稳定性——这也是主题,在这里处理感兴趣的读者可参 Winfree (1991)、Jahnke and Winfree (1991)、Barkley (1994)、Karma (1993, 1994)、Panfilov and Hogeweg (1995)、Kessler and Kupferman (1996) 的工作。由于可激发介质的分析研究太困难已寻更简单的模型——有限状态自动机 (finite-state automata) 相当流行有限状态自动机把状态空间划分为几个离散值(如 \(v = 0\)\(1\))、空间域划分为离散格子、时间划分为离散步——再制定格子状态如何随时间变化的规则有限状态自动机极其容易编程和可视化——它们对理解可激发介质的行为有用具迷惑性、可给出不易察觉的"错误"答案有限状态自动机的文献极其丰富(例如 Moe et al., 1964; Smith and Cohen, 1984; Gerhardt et al., 1990)。二维区域上的螺旋波三维中的显然推广是卷波 (scroll wave)(Winfree, 1973, 1991; Keener and Tyson, 1992)。卷波被数值地观察到(Jahnke et al., 1988; Lugosi and Winfree, 1988)、在三维 BZ 试剂中(Gomatam and Grindrod, 1987)、在心脏组织中(Chen et al., 1988)——虽然实验上极其困难可视化数值模拟可能有趣拓扑启动卷波——闭合卷环、打结卷波、链合卷环对等。卷波的数学理论也处于起步阶段闭合卷波拓扑属性已由 Winfree and Strogatz (1983a,b,c; 1984) 完成一般渐近理论(用于其演化)已由 Keener (1988) 提出、并圆卷环和螺旋卷做过数值实验对比这里没有足够空间讨论卷波理论——卷波在后面心脏波及节律章节中还会简略提及**。

本章个人批注

本章是全书的物理核心——它把第 5 章的"兴奋性"和"动作电位"从单点时间动力学推进到空间延展介质的行波传播。对我而言这一章的两个东西最触动

第一是"行波"作为数学对象的"实在性":作者在 6.1 节用一段就建立了行波 (traveling wave) = PDE 在无穷域上以常速保形传播的解这个定义;接下来 6.2 节把双稳方程 (6.1) 的行波前沿问题等价为 ODE (6.5) 的异宿轨道问题——这一降维是所有后续分析的起点。从 (6.5) 出发,作者用 (6.8) 的积分恒等式给 \(c\)符号一个面积解释;再用 \(c=0\) 极限和 \(c\) 大极限的夹逼存在性;用 \(c\) 单调性证唯一性。这一单调+连续+夹逼的方法对我自己的工作(如凝血动力学中的空间模型)直接可用——我以前多次在证明游走解唯一性时用类似的"两个极限夹一族"论证,作者给了一个教科书级的样板。

第二是奇异摄动作为"图像":6.3.1 节用 \(\epsilon\) 的小性把 FitzHugh–Nagumo 方程 (6.48)–(6.49) 分解为外方程 (6.68)(激发区/恢复区的慢动力学)和内方程 (6.71)(过渡层上的快动力学)。这一内外分离给行波脉冲一个几乎是工程图的图像:脉冲 = 一个上跳(升前沿)后跟一段激发相(持续时间 \(T_e\) 由 (6.76) 决定),再后跟一个下跳(波背)。速度的"匹配条件" \(c(w_+) = -c(w_-)\) 给我留下了深刻印象——它把"脉冲的稳态传播"等价一个代数方程:找到 \(w_-\) 使上升速度与下降速度互为相反数。若该方程无解,则下跳必须发生在膝点——波变成相位波 (phase wave),动力学发生质变(从双稳鞍-鞍-节点的异宿连接)。这一正则-奇异转换图像对我写神经组织中的电传导模型很有借鉴——可激发系统的两类解("扩散型" vs "相位型")在数学上是有根本不同的。作者明确说分析不严格(核心结构错误,目前无严格存在性证明),并指明 Keener (1992) 和 Kessler and Kupferman (1996) 的工作——这坦诚地承认理论的局限值得学习。

6.4 节色散曲线运动学重构也很有启发。作者把波列追踪问题等价为一个关于前沿速度 \(c_n\)递推 ODE (6.95)——这绕开完整 PDE 的解。这一思路与我所知的心脏电生理 APD 恢复模型惊人地相似:实际生理学家解完整 PDE,而是用APD-恢复曲线传导速度曲线作为经验数据驱动一个关于脉冲到达时间的差分/ODE 系统。作者明确说 (6.95) 的"优势是 \(A\)\(t_r\) 可由实验数据给出"——这恰恰是临床/工程取向与纯数学取向的桥接

6.5.2 节的螺旋波数学是本章有挑战性的——程函方程给出渐开线螺线但仅在 \(r \geq r_0\) 有效程函-曲率方程给出唯一的关系 \(\omega = c^2 m_*/\epsilon\)\(m_* = 0.330958\)数值确定的常数),且与色散曲线的交点唯一确定螺旋频率速度。这一渐近分析虽然不严格,但给出螺旋频率介质参数显式依赖。作者点明严格存在性至今证——这提醒我奇异摄动给出有用工程级预测,但证明还需新工具(如水平集方法 (6.109)消除尖点的技术)。最后作者提到卷波——3D 螺旋——及其拓扑属性(闭合卷环、打结卷波、链合卷环对)——这给我打开了一扇我以前充分注意的纯数学/拓扑**世界。

跨章节连接:6.2.2 的髓鞘化 → 第 4 章电缆理论;6.3.2 的 HH 行波脉冲 → 第 5 章 HH 模型;6.5 的高维波 → 第 12 章心肌波、第 7 章 Ca\(^{2+}\) 螺旋(已在 6.5.2 末尾预告)、第 18 章耦合振子(6.1 节预告);螺旋波在病理中的意义 → 第 12 章心脏致命螺旋。本章承接第 5 章 HH 方程导出行波,导出色散曲线作为周期性动作电位的速度-周期关系;第 7、12 章的Ca\(^{2+}\) 与心脏应用奠基

与上下章的衔接(一段话)

本章是第 5 章兴奋性的空间延展——前两章(4 章、5 章)一个奠定电缆理论、一个奠定兴奋性的单点动力学,本章把两者结合可激发介质的行波传播。第 5 章末尾预告的"一维雷管引线、二维干草地"在本章得到形式化:一维的双稳方程FitzHugh–Nagumo 方程给出行波前沿行波脉冲;二维的程函-曲率方程给出靶模式螺旋波。因此本章在全书结构居于"理论→应用"的转折点——前 5 章建立动力学与数值方法(细胞体积、离子通道、HH 方程、FitzHugh 简化、随机过程),本章首次把兴奋性放在空间域上系统处理。向下看,本章的许多结果会被引用到后续章节:6.2.2 的髓鞘化模型在第 4 章基础上自然延伸;6.3 的FitzHugh–Nagumo 行波为第 7 章的Ca\(^{2+}\) 行波提供方法论;6.4 的色散曲线为第 12 章的心肌波的周期响应分析直接相关;6.5.2 的螺旋波则在第 12 章心脏和第 18 章耦合振子中成为病理学的数学语言——毫不夸张地说,本章是从"单细胞"到"组织/器官"的数学桥梁