第 1 章：多物种波与实际应用（Multi-Species Waves and Practical Applications）

1.1 直觉预期（Intuitive Expectations）

本章把注意力从 Volume I 中讨论过的单组分反应扩散系统，推广到多个物种（细胞、反应物、种群、细菌等）相互作用的反应扩散-趋化机制。典型的一维系统写成 $\partial_t u = f(u) + D\nabla^2 u$，其中 $u$ 是反应物向量，$f$ 是非线性反应动力学，$D$ 是常扩散系数矩阵。在这一节里，作者先不给严格的数学分析，而是先做"直觉预热"，勾勒出后续各章将要看到的多物种波现象的丰富图谱。

从单组分系统在三个均匀稳态 $u_1, u_2, u_3$ 的简单情形出发，作者回忆 Volume I §13.5 的结论：如果初始在 $u = u_1$ 处全空间一致，然后突然把 $x < 0$ 区域改到 $u_3$，那么扩散会引发一个从 $u_3$ 到 $u_1$ 的行波前沿，把整个空间带到 $u_3$。对于多组分系统，作者认为类似的连接两稳态的行波前沿应当存在；尽管数学上可能存在解的全谱，但这里只关心非负解。捕食者-被捕食者模型中会出现一类"追逐-逃避波"（pursuit-evasion wave），在 §1.2 详述；而 §1.5 则把模型对接到 Belousov–Zhabotinskii (BZ) 化学反应的行波前沿，并与实验对照。实践中还有两类竞争波的实例：转基因生物的扩散，以及英国灰松与红松的分布变化。

在 Volume I 第 3 章的回顾中，单种群或单反应物不会产生稳定的极限环，但多反应物反应动力学可以出现稳定的极限环，并通过 Hopf 分岔从稳态分叉出来。假设参数 $\gamma > \gamma_c$ 使整个系统在空间上稳定振荡，如果在小区域 $0 < |x| \le \varepsilon$ 内短时扰动，使该处的振荡相位与外部不同，这就形成了一个局部"起搏点"（pacemaker）。扩散试图抹平这种差异，但每次起搏点重新进入新相位都相当于重新点燃一个行波前沿，因此反应扩散系统会承载一串行波波列（wave train）穿过介质，这种结构在 §1.7 详细讨论。当涉及三方程以上甚至单时滞方程时，混沌振荡是可能的；如果把一个小的混沌振荡区放到规则振荡的背景里，可能产生相当复杂的波现象。即便只在一维空间，解的行为也可能相当"巴洛克"；混沌行为可以不存在混沌起搏点，详见 §1.9 的图 1.23。

进入二维空间，一个以异于周围频率振荡的小圆形区域会形成从中心向外传播的同心圆波列——目标图案（target patterns）。Zaikin 与 Zhabotinskii（1970）首次在 BZ 反应中实验观察到目标图案（图 1.1(a)）。把图像再推进一步：让一个"起搏点"沿一个小圆环持续移动，每经过一个位置就引发一道新波向外传播，所有这些波的包络就形成螺旋波（spiral wave），核心是极限环起搏器。螺旋波在 BZ 反应中也有实验记录（图 1.1(b)；Winfree 1974，Müller et al. 1985，Agladze 与 Krinskii 1982），Kuramoto 与 Koga (1981) 和 Agladze 与 Krinskii (1982) 还演示了混沌波图案的诞生（图 1.23）。三维情形下，拓扑结构更加惊人：每一个"二维"螺旋的局部本身也是螺旋；Winfree (1974)、Welsh et al. (1983) 给出三维波的实验照片，Winfree 与 Strogatz (1984)、Winfree (2000) 讨论了拓扑性质；Mimura 一派对球面波做了大量解析与数值工作（Yagisita et al. 1998）。目标图案和螺旋波在生物学中很常见，尤其是螺旋波，在心脏学与神经生物学中有重要的实际意义。黏菌 Dictyostelium discoideum 是一个著名的生物学例子（图 1.1(c) 与图 1.18，Newell 1983）。

反应动力学只有一个稳定稳态但具有阈值行为的情形，与 Volume I §3.8、§7.5 的 FitzHugh–Nagumo (FHN) 模型同源。假设整个空间处在稳定稳态，扰动一个小区域使局部触发阈值行为，虽然扰动最终会衰减，但在衰减之前会经历一个大的相空间偏移，暂时看起来像两个不同稳态共存，扩散从而可以引发一个行波前沿。阈值机制给行波脉冲（travelling pulse wave）提供了基础，§1.6 详细讨论。横向耦合的波可以用基本可激发模型系统来分析（Gáspár et al. 1991），他们展示出相互作用的圆波如何产生螺旋波，以及复杂平面波图案如何演化；Petrov et al. (1994) 考察了带三次自催化的反应扩散系统的波反射和波分裂；Pascual (1993) 数值地证明沿空间梯度扩散的某些标准捕食者-被捕食者模型可以在固定空间点呈现时间混沌，并提出准周期路径的证据；Sherratt et al. (1995) 研究了一维捕食者-被捕食系统的漫画模型，证明入侵波的尾流中可出现混沌；Merkin et al. (1996) 用带三次反应项的两物种模型研究了波诱导混沌；Epstein 与 Showalter (1996) 综述了非线性化学动力学中振荡、波图案与混沌的复杂性。Maini (1995) 编辑的论文集展示了时空波现象在生物医学领域的普遍性和多样性，包括创伤愈合、肿瘤生长、胚胎学、群体中的个体运动、细胞-细胞相互作用等。对于趋化机制，例如 Dictyostelium 黏菌（cf. Volume I §11.4），Keller 与 Segel (1971)、Keller 与 Odell (1975) 证明在一定参数范围内存在行波；Tyson et al. (1998, 1999) 描述了留下假稳态空间图案的更复杂细菌波，第 5 章会详细讨论。作者最后总结：多物种反应扩散-趋化机制中的空间波现象比单物种情形丰富得多；如果再允许混沌起搏点、时滞动力学等，现象谱会更宽；尽管已有大量研究，许多实际波问题尚未被研究并将带来新的时空现象。

1.2 捕食者-被捕食者系统中的追逐与逃避波（Waves of Pursuit and Evasion in Predator–Prey Systems）

当捕食者与被捕食者在空间上分布时，种群数量会随时间-空间变化：捕食者向猎物移动、猎物逃离捕食者。行波带（travelling bands）已在海洋浮游生物（Wyatt 1973）、动物迁徙、真菌与植被（Lefever 与 Lejeune 1997，Lejeune 与 Tlidi 1999）等现象中观察到；在原始生物入侵营养源的运动中也相当常见，第 5 章将详细讨论细菌在趋化信号下产生的复杂空间波和空间现象。本节主要作为分析技术的示例，研究一个带扩散的简单捕食者-被捕食者系统，并展示行波前沿解如何出现。模型是 Lotka–Volterra 系统的修改版（参见 Volume I §3.1）：猎物有 logistic 增长，捕食者和猎物都通过扩散分散。Dunbar (1983, 1984) 详细研究过这个模型：

\[ \partial_t U = AU\left(1 - \frac{U}{K}\right) - BUV + D_1\nabla^2 U, \quad \partial_t V = CUV - DV + D_2\nabla^2 V, \]

其中 $U$ 是猎物、$V$ 是捕食者，$A, B, C, D$ 和 $K$（猎物环境容纳量）均为正常数，$D_1, D_2$ 是扩散系数。无量纲化（设 $u = U/K$, $v = BV/A$, $t^* = At$, $x^* = x\sqrt{A/D_2}$, $D^* = D_1/D_2$, $a = CK/A$, $b = D/CK$）后，丢弃星号得到一维方程

\[ \partial_t u = u(1 - u - v) + D \partial_{xx} u, \quad \partial_t v = av(u - b) + \partial_{xx} v, \]

我们只关心非负解。

空间无关系统的分析是 Volume I 第 3 章相平面分析的直接应用：三个稳态为 $(0,0)$、$(1,0)$（无捕食者、猎物处于环境容纳量）、$(b, 1-b)$（两物种共存，假设 $b < 1$）。$(0,0)$ 和 $(1,0)$ 都不稳定，$(b, 1-b)$ 在 $4a \le b/(1-b)$ 时为稳定结点，$4a > b/(1-b)$ 时为稳定螺线。事实上，方程 $(1.3)$（去掉空间项）有 Lyapunov 函数

\[ L(u,v) = a\left(u - b - b\ln\frac{u}{b}\right) + \left(v - 1 + b - (1-b)\ln\frac{v}{1-b}\right), \]

$L(b, 1-b) = 0$，其它正象限 $(u,v)$ 处 $L$ 为正且 $dL/dt < 0$，故 $(b, 1-b)$ 是全局稳定的（Jordan 与 Smith 1999 是 Lyapunov 函数的可读介绍）。回想 Volume I §3.1 中，最简 Lotka–Volterra 系统（即去掉猎物饱和项的 $(1.2)$）的非零共存稳态只是中性稳定，没有实际用处；修改后的系统 $(1.2)$ 更符合实际。

通过设 $u(x,t) = U(z)$, $v(x,t) = V(z)$, $z = x + ct$（$c > 0$ 是待定波速，向左传播），将 $(1.3)$ 化为常微分方程组

\[ cU' = U(1 - U - V) + DU'', \quad cV' = aV(U - b) + V'', \]

其中撇号表示对 $z$ 求导。这是四维相空间问题。作者考虑一个简化情形：猎物的扩散远小于捕食者的扩散，取 $D = D_1/D_2 = 0$（即对应一个浮游-草食动物系统中只有草食动物能动的情形）。Dunbar (1984) 证明 $D = 0$ 与 $D \ne 0$ 系统的定性行为基本相同。$D = 0$ 时的三阶一阶常微系统为

\[ U' = \frac{U(1-U-V)}{c}, \quad V' = W, \quad W' = cW - aV(U-b). \]

在 $(U, V, W)$ 相空间中有两个不稳定稳态 $(0,0,0)$、$(1,0,0)$ 和一个稳定稳态 $(b, 1-b, 0)$，只在 $b < 1$ 时关心。仿照 Volume I §13.2 中 Fisher–Kolmogoroff 方程的分析经验，从 $(1,0,0)$ 到 $(b, 1-b, 0)$ 或从 $(0,0,0)$ 到 $(b, 1-b, 0)$ 都可能存在行波解。本节只关心边界条件 $(1.7)$：$U(-\infty) = 1, V(-\infty) = 0$ 与 $U(\infty) = b, V(\infty) = 1-b$。

关于 $(1, 0, 0)$ 的线性化给出特征值 $\lambda_1 = -1/c$ 与 $\lambda_{2,3} = \frac{1}{2}(c \pm \sqrt{c^2 - 4a(1-b)})$。存在不稳定流形（关联 $\lambda_2, \lambda_3$ 为正）。若 $c^2 < 4a(1-b)$，$(1, 0, 0)$ 以振荡方式不稳定；非负 $U, V$ 的行波前沿解存在的唯一可能是

\[ c \ge [4a(1-b)]^{1/2}, \quad b < 1. \]

在此条件下，一个具有波速下界的现实解可能存在，$u = 1, v = 0$ 作为 $z \to -\infty$ 时的状态。这让人想起 Volume I 第 13 章的行波前沿解。

但这里的解与 Volume I 第 13 章在定性上可以不同。在 $(b, 1-b, 0)$ 附近线性化 $(1.6)$，特征值由三次特征多项式 $p(\lambda) \equiv \lambda^3 - \lambda^2(c - b)/c - \lambda b - ab(1-b)/c = 0$ 给出。其极值点 $\lambda_{M,m}$ 与 $a$ 无关。$a = 0$ 时根为 $\lambda = 0$ 和两个其它根（见原文 $(1.11)$ 上下的式子）。当 $a$ 从零增大，效应是从 $p(\lambda; a=0)$ 曲线上处处减去 $ab(1-b)/c$，由于极值与 $a$ 无关，曲线情形如图 1.2 所示。对 $0 < a < a^*$ 有两个负实根和一个正实根；$a = a^*$ 时两个负实根相等；$a > a^*$ 时两个负实根变成带负实部的复数根（连续性论证保证 $a$ 刚大于 $a^*$ 时即如此）。$a^*$ 可解析确定；Routh–Hurwitz 判据（Volume I 附录 B）也能给出相同结论但直观性较差。

临界 $a^*$ 的存在意味着：$a > a^*$ 时满足 $(1.7)$ 的行波前沿解 $(U, V)$ 振荡式地趋近稳态 $(b, 1-b)$；$a < a^*$ 时单调趋近。图 1.3 示意两种解的形态。完全捕食者-被捕食者系统 $(1.3)$（双方均扩散）也产生可呈振荡的行波前沿（Dunbar 1983, 1984）；这些波的存在性证明需要在相平面上仔细分析存在一条正象限内的连接相关奇点的轨迹。这些波常被称为"追逐-逃避波"，尽管图 1.3 所示解中很难看出猎物的明显逃避——除了安静繁殖外，猎物基本上就是等着被吃。

对流式捕食者-被捕食者追逐与逃避模型。 完全不同的"追逐-逃避"捕食者-被捕食系统是：猎物试图逃避捕食者，而捕食者只在双方相遇时才去抓猎物。这会引发一种本质不同的空间相互作用。作为示例，作者用一维形式简述一个可能模型。设猎物 $u$ 与捕食者 $v$ 可分别以速度 $c_1, c_2$ 移动，扩散对种群分散不起主要作用，每个种群按自身动力学演化、达到自身稳态。考虑图 1.4(a)：两个种群不相互作用，由于没有扩散分散，任意空间位置的种群都生长或衰减到该种群自己的稳态。这种动力学情形如图 1.4(a)，两个种群以未受干扰的速度 $c_1, c_2$ 移动，没有空间扩散，因此带的宽度 $w_u, w_v$ 在 $u, v$ 趋向稳态时保持不变。现在假设捕食者追上猎物时，猎物通过以正比于捕食者梯度的额外速度远离捕食者，类似地捕食者以正比于猎物梯度的额外速度向猎物集中移动。在最基本的非平凡层面，可以用包含对流效应的守恒律（参见 Volume I 第 11 章）建模为

\[ u_t - [(c_1 + h_1 v_x) u]_x = f(u, v), \quad v_t - [(c_2 - h_2 u_x) v]_x = g(v, u), \]

其中 $f, g$ 表示种群动力学，$h_1, h_2$ 是与捕食者追击、猎物逃避相关的正常数。这些是 $u, v$ 的守恒律，所以左端必须是散度形式。

交互项 $f$ 和 $g$ 可以是任何捕食者-被捕食者情形。典型地，$f(u, 0)$ 表示没有捕食者时猎物的动力学（种群简单生长或衰减到非零稳态）；捕食者的作用是降低猎物稳态大小，所以 $f(u, 0) > f(u, v > 0)$。类似地，$g(v, u \ne 0)$ 产生的稳态大于 $g(v, 0)$ 产生的稳态。

要理解对流项的物理意义，假设 $(1.12)$ 中 $h_1 = 0$，则 $u_t - c_1 u_x = f(u, v)$，这正是以速度 $c_1$ 移动的行波坐标系下猎物的动力学。如果取 $z = x + c_1 t$ 作为独立变量，方程简化为 $u_t = f(u, v)$。如果 $c_2 = c_1$，则 $h_2 = 0$ 时捕食者方程变为 $v_t = g(v, u)$。于是我们看到种群变化的行波直到稳态（图 1.4(a)），之后变成恒定形状的行波（"礼帽波"）。

现在考虑 $h_1, h_2$ 都为正、$c_1 \ne c_2$ 的更复杂情形。参考图 1.4(b) 的重叠区域，$(1.12)$ 中 $h_1 v_x$ 项（$v_x > 0$ 故为正）使猎物的波在局部上向左移动加速；$-h_2 u_x$（$u_x < 0$ 故为正）增大捕食者的局部对流速度。相互作用的复杂本质取决于解的形式、参数相对大小以及动力学。由于对流项（以及动力学）的非线性，$u, v$ 可能存在不连续跳跃的激波解（见 Murray 1968, 1970, 1973；反应扩散例子见 Volume I 第 13 章 §13.5）。

将系统 $(1.12), (1.13)$ 改写：展开左端，得

\[ u_t - [(c_1 + h_1 v_x)] u_x = f(u, v) + h_1 u v_{xx}, \]

\[ v_t - [(c_2 - h_2 u_x)] v_x = g(v, u) - h_2 v u_{xx}. \]

这个形式中 $h_1, h_2$ 项表示交叉扩散，一正一负。交叉扩散（Volume I §11.2 定义）出现在多物种模型中，扩散矩阵非严格对角时：一个物种的方程中出现另一个物种的扩散型项。例如在 $u$ 方程中，$h_1 u v_{xx}$ 像是 $v$ 的扩散项，扩散系数为 $h_1 u$。这正是一个在建模问题中自然产生交叉扩散的例子，并不多见。

系统 $(1.12)$–$(1.14)$ 的数学分析是个挑战，相关研究还不充分。Hasimoto (1974)、Yoshikawa 与 Yamaguti (1974) 考察了 $h_1 = h_2 = 0$ 的情形，Murray 与 Cohen (1983) 研究了 $h_1, h_2$ 非零的系统。Hasimoto 在 $f(u, v) = l_1 u v, g(u, v) = l_2 u v$（$l_1, l_2$ 为常数）情形下得到解析解，展示了 blow-up 在某些情形下可发生。对一般系统 $(1.12)$–$(1.14)$ 类型，新的解行为是值得期待的。

二维情形下的对流追逐与逃避问题具有生态学意义且尤其具有挑战性，目前尚未被研究。本书第一版曾假设：建模狼-驼鹿捕食者-被捕食者情形时，引入狼的领地边界效应，看领地间的"无人区"是否能为猎物提供部分安全庇护。直观理由是狼有较少的倾向去窜入相邻领地。已有证据表明驼鹿沿狼领地边界移动。沿这个方向的研究已经完成，将在第 14 章详细讨论。

对流耦合动力学还能产生相关波现象，例如固定柱中的生化离子交换。单反应动力学耦合对流的过程，Goldstein 与 Murray (1959) 详细研究过，光滑初始数据可以演化出有趣的激波解。这里发展的数学技术与上述问题直接相关。当多个离子交换同时发生，就是色谱法，它是生化分析中一种强大的技术。

1.3 英国灰松空间扩散的竞争模型（Competition Model for the Spatial Spread of the Grey Squirrel in Britain）

背景与生态学事实。 20 世纪初左右，北美灰松（Sciurus carolinensis）在英国多地被释放，最重要的释放地在东南部。此后灰松成功扩散到英国大部分地区，最北到达苏格兰低地；与此同时，本地红松（Sciurus vulgaris）从这些地方消失。Lloyd (1983) 注意到，灰松的涌入通常与红松的衰退和随后的消失相吻合，二者在分布上的重叠期仅持续几年。英国的松鼠分布记录似乎表明灰松对红松有明确的负面效应（Williamson 1996）。MacKinnon (1978) 给出了三条可能假设（竞争、环境变化、红松疾病如"松鼠流感"）中竞争最可能成立的理由（Reynolds 1985）；这些假设并不互斥。

灰松引入之前，红松在没有种间竞争的环境下进化，选择倾向于中等繁殖力和低数量损耗。灰松则在强种间竞争背景下（与美国红松和狐松）进化，选择倾向于过度繁殖。两种松鼠一年可繁殖两次，但红松体型较小，每胎通常不超过 2–3 只幼崽，而灰松常一胎 4–5 只（Barkalow 1967）。在北美，红松和灰松占据几乎不重叠的生态位：灰松偏好混合硬木林，红松偏好北方针叶林。但在英国，没有灰松的环境下，本地红松必须进化到既适应硬木林也适应针叶林。Holm (1987) 的工作也支持以下假设：灰松在落叶林地带可能具有竞争优势，本地红松大部分已被灰松取代。北美灰松体型大、健壮，体重约为红松的两倍。在分开的栖息地中，两种松鼠表现出相似的社会组织、觅食和活动生态学；但在同一栖息地中，对资源的利用应更为相似，因此如此相近的两个物种无法在共享同一资源的情况下共存。

总之，似乎合理地假设：两种松鼠的相互作用（很可能主要是通过对资源的间接竞争，但也包括直接相互作用如追逐）有利于灰松把红松从英国落叶林大部分地区赶走。Okubo et al. (1989) 调查了灰松对红松的这种取代，并基于以上事实提出和研究了竞争模型。本节沿用他们的工作。他们还用模型模拟了灰松被随机引入红松区域后如何扩散殖民，并将模型结果与可获得的数据进行了比较。

竞争模型系统。 记 $S_1(X, T)$、$S_2(X, T)$ 分别为灰松和红松在位置 $X$、时间 $T$ 的种群密度。假设两者竞争相同的食物资源，一个可能的模型是带扩散的 Lotka–Volterra 竞争系统（参见 Volume I 第 5 章）：

\[ \partial_T S_1 = D_1 \nabla^2 S_1 + a_1 S_1(1 - b_1 S_1 - c_1 S_2), \]

\[ \partial_T S_2 = D_2 \nabla^2 S_2 + a_2 S_2(1 - b_2 S_2 - c_2 S_1), \]

其中对 $i=1,2$，$a_i$ 是净出生率，$1/b_i$ 是环境容纳量，$c_i$ 是竞争系数，$D_i$ 是扩散系数，均为非负。交互（动力学）项就是带竞争的 logistic 增长。基于上述讨论，假设灰松胜过红松：

\[ b_2 > c_1, \quad c_2 > b_1. \]

无量纲化（设 $\theta_i = b_i S_i$, $t = a_1 T$, $x = \sqrt{a_1/D_1} X$, $\gamma_1 = c_1/b_2$, $\gamma_2 = c_2/b_1$, $\kappa = D_2/D_1$, $\alpha = a_1/a_2$），$(1.15)$ 化为

\[ \partial_t \theta_1 = \nabla^2 \theta_1 + \theta_1(1 - \theta_1 - \gamma_1 \theta_2), \]

\[ \partial_t \theta_2 = \kappa \nabla^2 \theta_1 + \alpha \theta_2(1 - \theta_2 - \gamma_2 \theta_1). \]

由 $(1.16)$ 可得 $\gamma_1 < 1, \gamma_2 > 1$。Volume I 第 5 章已详细分析过无扩散时该竞争模型的相平面：三个齐次稳态 $(0,0)$（不稳定结点）、$(1,0)$（稳定结点）、$(0,1)$（鞍点）。加上扩散后，按通常方法，存在一条从 $(0,1)$ 到 $(1,0)$ 的解轨迹和连接这些临界点的行波。这对应灰松（$\theta_1$）竞争排除红松（$\theta_2$）的生态学情形。

一维情形下，设行波解 $\theta_i = \theta_i(z)$, $z = x - ct$, $c > 0$，$\theta_1(z), \theta_2(z)$ 表示以速度 $c$ 沿 $x$ 正方向运动的恒定形状行波。代入 $(1.18)$ 得

\[ \frac{d^2 \theta_1}{dz^2} + c \frac{d\theta_1}{dz} + \theta_1(1 - \theta_1 - \gamma_1 \theta_2) = 0, \]

\[ \kappa \frac{d^2 \theta_2}{dz^2} + c \frac{d\theta_2}{dz} + \alpha \theta_2(1 - \theta_2 - \gamma_2 \theta_1) = 0, \]

边界条件为 $\theta_1 = 1, \theta_2 = 0$（在 $z = -\infty$），$\theta_1 = 0, \theta_2 = 1$（在 $z = \infty$），即波以速度 $c$ 传播，灰松把红松取代。

Hosono (1988) 在 $(1.19)$ 和 $(1.22)$ 条件下、参数值在某些范围时，研究了 $(1.18)$ 系统的行波存在性。一般情况下 $(1.18)$ 的常微系统无法解析求解。但在特殊情形 $\kappa = \alpha = 1$, $\gamma_1 + \gamma_2 = 2$ 下，可以得到一些解析结果：将两方程相加得到

\[ \frac{d^2 \theta}{dz^2} + c \frac{d\theta}{dz} + \theta(1 - \theta) = 0, \quad \theta = \theta_1 + \theta_2, \]

即著名的 Fisher–Kolmogoroff 方程（Volume I 第 13 章），已知它在 $\pm \infty$ 处有合适的边界条件下的行波解。然而此处的边界条件（$\theta = 1$ 在 $z = \pm \infty$）与经典 Fisher–Kolmogoroff 方程不同；这意味着对所有 $z$ 有 $\theta = 1 \Rightarrow \theta_1 + \theta_2 = 1$。代入 $(1.21)$ 第一式得

\[ \frac{d^2 \theta_1}{dz^2} + c \frac{d\theta_1}{dz} + (1 - \gamma_1) \theta_1 (1 - \theta_1) = 0, \]

又是 $\theta_1$ 的 Fisher–Kolmogoroff 方程，边界条件为 $(1.22)$。由波速结果得到灰松行波前沿速度

\[ c \ge c_{\min} = 2(1 - \gamma_1)^{1/2}, \quad \gamma_1 < 1. \]

类似地，对红松

\[ c \ge c_{\min} = 2(\gamma_2 - 1)^{1/2}, \quad \gamma_2 > 1. \]

由于 $\gamma_1 + \gamma_2 = 2$（且 $\kappa = \alpha = 1$），两个最小波速相等。回到有量纲量，可得最小波速

\[ C_{\min} = 2 \left[ a_1 D_1 \left( 1 - \frac{c_1}{b_2} \right) \right]^{1/2}. \]

参数估计。 现在必须把分析结果与英国真实的竞争情形对接。行波波速依赖于模型 $(1.15)$ 中的参数，所以需要估计这些参数，以便将理论波速与可获得的数据比较。这正是真实建模的关键方面。

考虑本征净增长率 $a_1, a_2$。Okubo et al. (1989) 用了一个由 Williamson 与 Brown (1986) 详细描述的修改 Leslie 矩阵。原则上估计应该是种群密度为零处的值，但人口学数据通常指接近平衡密度时的种群。估计本征净增长率时考虑三个成分：性别比、出生率、死亡率。性别比取一比一。但出生率和死亡率难以确定，取决于窝仔数与频度、它们对年龄的依赖性、年龄分布、数据采集的时间和地点、食物资源水平和预期寿命；Okubo et al. (1989) 的论文展示了相关细节。在对众多（有时互相矛盾）来源的仔细分析后，他们估计灰松的本征出生率为 $a_1 = 0.82$/年，稳定年龄分布为接近 3 只幼崽对 1 只成年；红松的 $a_2 = 0.61$/年，年龄分布为略多于 2 只幼崽对 1 只成年。

估计环境容纳量 $1/b_1, 1/b_2$ 也涉及对已有文献的类似详细考察（Okubo et al. 1989 同样做了）。他们建议灰松和红松的环境容纳量分别为 $1/b_1 = 10$/公顷和 $1/b_2 = 0.75$/公顷。

很遗憾，竞争系数 $c_1, c_2$ 没有定量信息。但模型中只需要比值 $c_1/b_2 = \gamma_1$ 和 $c_2/b_1 = \gamma_2$ 来估计行波最小波速。灰松的传播速度只依赖于 $\gamma_1$ 的估计；记得 $0 < \gamma_1 < 1$。由于 $\gamma_1$ 以 $(1 - \gamma_1)^{1/2}$ 形式出现在最小波速中，$\gamma_1$ 不太大时（实际上不超过 0.6 左右时）波速对它不敏感。我们预期竞争系数 $c_1$（红松对灰松）应有较小值。再加上 $b_2^{-1}$ 较小的事实，可以合理地假设 $\gamma_1$ 接近零，于是从 $(1.30)$ 灰松行波最小速度近似为 $2(D_1 a_1)^{1/2}$。Okubo et al. (1989) 在数值模拟中对不同的 $\gamma$ 用了多个值，因为上面的分析只对一些特殊值有效。

接下来考虑扩散系数 $D_1, D_2$。它们是波传播中的关键参数，但极难估计。（同样的扩散估计问题在后面讨论狂犬病在狐狸种群中的空间扩散、细菌图案和脑肿瘤细胞时还会出现。）直接观察分散既困难又通常是短期的。报道的运动数值变化范围很大。还要考虑林地之间的迁移。

基于个体运动，灰松的一维扩散系数最大值估计为 1.25 平方公里/年，二维扩散系数为 0.63 平方公里/年。然而这可能不对应灰松在林地之间的运动。如果灰松的年度分散主要发生在林地之间而非林地内部，这些扩散系数的值可能太小而无法作为代表。Okubo et al. (1989) 推测：灰松的扩散系数合理预期为 10–20 平方公里/年级别，而非 1 平方公里/年级别。他们给出了一个启发式论证（这里转述）来支持这些大得多的扩散系数。

考虑一片林地，每片面积 $A$ 公顷，四个相邻林地，彼此间距 $l$ 公里。假设一片林地充满灰松，环境容纳量为 10/公顷。该片林地容纳 $n = 10A$ 只灰松。设本征增长率 $a_1 = 0.82$/年，下一年林地将有 $22.7A$ 只动物（$e^{0.82} = 2.27$），其中 $12.7A$ 只必须分散。假设动物分散到最近的相邻林地，$12.7/4A = 3.175A$ 只将到达相邻林地。这片林地将在 $\tau = 1.4$ 年内被灰松充满（$10A = 3.175 A e^{0.82 \tau}$），之后会再次分散。换言之，灰松平均每 1.4 年向最近邻林地分散一次。于是灰松的二维扩散系数估计为

\[ D_1 = l^2 / (4 \times 1.4) \text{ 平方公里/年} = l^2/5.6 \text{ 平方公里/年}. \]

表 1.1 给出由 $(1.31)$ 计算出的 $D_1$ 关于 $l$ 的函数。Williamson 与 Brown (1986) 估计灰松的分散速度为 7.7 公里/年。取此值，由表中得 $D_1 = 17.9$ 平方公里/年。所以，相邻林地平均间距 10 公里是合理的，由此得到的行波最小速度与数据吻合。

理论扩散速度与数据的比较。 英国灰松扩散信息最好的来源之一是 Reynolds (1985) 对 1960–1981 年东英吉利的详细研究。东英吉利的灰松殖民相对较新。1959 年没有发现灰松，红松在 1959 年和 1971 年的调查中仍基本遍布诺福克郡。然而到 1971 年，灰松也已记录在诺福克郡约一半的区域。

Reynolds (1985) 用 $5 \times 5$ 平方公里网格构建了一系列地图，展示 1960–1981 年灰松和红松的年度分布。Williamson 与 Brown (1986) 基于这些地图计算了 1965–1981 年期间灰松的扩散速度，得到 5 至 10 公里/年；灰松的平均扩散速度计算为 7.7 公里/年（即上面提到的值），用 Fisher–Kolmogoroff 最小波速即可估计灰松扩散系数约为 $D_1 = 17.9$ 平方公里/年。所以，对上面给出的扩散系数的启发式估计有了一定的数据支持。

无量纲模型系统 $(1.18)$ 的解必须用数值方法求出（$\gamma$ 值不满足 $\gamma_1 + \gamma_2 = 2$ 时）。一维情况下，波的形态如边界条件和方程形式所预期，$\gamma$ 的取值无关大局。灰松有一个与典型 Fisher–Kolmogoroff 波相似的推进波，红松则有一个几乎对称镜像的撤退波；Okubo et al. (1989) 给出了数值结果。当然殖民是二维的，没有解析的推进行波前沿解，这是一个非常困难的问题。在灰松径向对称分布的特例中，入侵波速度小于对应的一维情形（因为 Laplacian 中的 $(1/r)(\partial \theta_1/\partial r)$ 项和 $\theta_2$ 的对应项）。但数值解相对容易得到。作者从以红松为主、零星灰松斑块的初始分布出发，这些小灰松区域向外扩展、与其它灰松区域合并，最终把红松完全赶走。图 1.5 展示了一组典型参数下的数值解。

基于扩散和增长的种群扩散基础模型（如 Fisher–Kolmogoroff 模型）从一个初始种子开始径向扩散，最终因为 $(1/r)(\partial \theta_1/\partial r)$ 项随 $r \to \infty$ 趋于零而变成有效的一维波。本文讨论的模型也有相同表现，只是竞争推进波较慢，这并不奇怪，因为灰松的有效出生率小于简单的 logistic 增长。Okubo et al. (1989) 中有众多灰松扩散、红松撤退的地图（可追溯到 1930 年）。图 1.5 所示行为是所见主要格局的合理代表。所用的参数值基于 Reynolds (1985) 的详细调查。参数和竞争进程不可避免地随气候、树木密度和类型而变化。看来灰松取代红松的广泛特征可以由这个简单竞争模型捕捉到，它是 Volume I 第 5 章中讨论的竞争排除原则的实际例子。

1.4 转基因生物的扩散（Spread of Genetically Engineered Organisms）

利用重组 DNA 技术改造植物（和动物）以执行特殊农业功能的应用正在迅速增加。然而，对生态系统甚至气候系统可能因释放这类转基因生物而受到破坏的担忧也在增加。研究转基因生物在自然环境中的时空动力学显然变得越来越重要。科学家尚未就转基因生物的风险或控制达成共识。对植物而言，初步时间尺度比转基因树木要短。例如果树已被改造以杀死落在叶子上的害虫。也在研究用树木清理污染土地。杂交植物当然已被广泛使用很长时间，但没有直接的基因设计输入。对于更具争议的转基因（以及克隆）动物，存在其他严重风险。它们的使用与用于人体移植的动物发育相关，带来不同的流行病学问题。无论 Luddites 抗议者们说什么或做什么，对植物和动物（包括人类）的基因操控都将持续存在。

关于释放工程生物的主要关注之一是：在不同生态情景和管理计划下，它们可能扩散多远、多快。对释放这类生物相关风险的无偏评估应导致有效控制疫情暴发的策略。估计扩散速率和分析可能控制策略的可靠定量信息仍然很少。Cruywagen et al. (1996) 沿这些方向做了一些初步工作。

转基因微生物特别适合数学分析，因为它们不断繁殖、缺乏复杂行为、其种群动力学可被简单模型良好描述。一种这样的微生物是 Pseudomonas syringae（ice-minus 细菌），它能通过占据作物叶面、排除引起冻害的 Pseudomonas syringae 菌株来减少对作物的冻害（Lindow 1987）。

本节发展一个模型，对转基因生物在空间异质环境中的时空扩散给出定量结果；沿用 Cruywagen et al. (1996) 的工作。我们从释放地点的工程种群的扩散率和增长率以及竞争种的相应参数中，获得关于工程种群暴发风险的信息。环境的性质在生物扩散中起关键作用。我们特别关注：能否通过地理屏障（如水、不同作物或贫瘠土地）保证控制。

基本模型从上一节灰松空间扩散用的两个竞争和扩散的物种系统出发。如上所述，该模型为外部引入的灰松入侵、本地红松被驱赶到灭绝区域的竞争扩散提供了解释。

大多数入侵模型把入侵视为在齐次环境中传播的行波。然而由于环境变化（自然的或人为的），这几乎从来不是这样。空间异质性不仅是最明显的自然特征之一，也可能是影响种群动力学的更重要因素。

Shigesada et al. (1986)（另见 Shigesada 与 Kawasaki 1997 的书）对 Fisher–Kolmogoroff 方程（描述具有 logistic 增长和分散的单一物种）在异质无界栖息地中的传播行波做了首次分析。这里再次使用带扩散的 Lotka–Volterra 竞争模型来描述自然微生物和竞争工程微生物的种群动力学。然而，我们通过假设一个由好坏斑块周期交替构成的周期性变化区域，对模型加以修改以适应空间异质环境。好斑块代表有利于微生物释放的有利区域，坏斑块则模拟抑制微生物扩散的不利屏障。我们特别关注工程种群的入侵和控制条件。

虽然本节讨论的动机是确定转基因生物扩散的条件，但模型和分析也适用于其它外来物种的引入（在这些情形下，控制或有时是有意传播是目标）。

记 $E(x, t)$ 为工程微生物、$N(x, t)$ 为未改造的微生物。这里只考虑一维情形。用经典 Lotka–Volterra 动力学描述工程微生物与自然微生物的竞争，并允许关键模型参数随空间变化以反映栖息地异质性。所以系统动力学为

\[ \partial_t E = \partial_x\left(D(x) \partial_x E\right) + r_E E [G(x) - a_E E - b_E N], \]

\[ \partial_t N = \partial_x\left(d(x) \partial_x N\right) + r_N N [g(x) - a_N N - b_N E], \]

其中 $D(x), d(x)$ 是空间依赖的扩散系数，$r_E, r_N$ 是生物的本征增长率。这些量被标度使得承载容量函数 $G(x), g(x)$ 的最大值为 1。正参数 $a_E, a_N$ 度量种内竞争效应，$b_E, b_N$ 是种间竞争系数。

本节通过让扩散和承载容量 $D(x), d(x), G(x), g(x)$ 为空间周期函数来建模环境异质性。设 $l$ 是环境变化的周期，则定义 $D(x) = D(x+l)$，$d(x), G(x), g(x)$ 类似。初始假设没有工程微生物，即 $E(x, 0) \equiv 0$，所以自然微生物 $N(x, 0)$ 满足

\[ \partial_x\left(d(x) \partial_x N\right) + r_N N(g(x) - a_N N) = 0. \]

工程生物随后在释放地点（取为原点）引入。$E(x, t)$ 的初始分布由下式表示：

\[ E(x, 0) = \begin{cases} H > 0, & |x| \le x_c, \\ 0, & |x| > x_c, \end{cases} \]

其中 $H$ 是正常数。

为引入有利和不利斑块的概念，假设环境由两种齐次斑块组成：长度为 $l_1$ 的斑块 1（有利斑块）和长度为 $l_2$ 的斑块 2（不利斑块），沿 $x$ 轴交替排列，$l = l_1 + l_2$。不利斑块中生物的扩散和承载容量小于有利斑块。这可能因为不利斑块是恶劣环境，限制了种群或干扰其分散。对应地，$D(x), d(x), G(x), g(x)$ 是 $x$ 的周期函数。在斑块 1 中（$ml < x < ml + l_1$, $m = 0, 1, 2, \ldots$）：

\[ D(x) = D_1 > 0, \quad d(x) = d_1 > 0, \quad G(x) = 1, \quad g(x) = 1, \]

在斑块 2 中（$ml - l_2 < x < ml$）：

\[ D(x) = D_2 > 0, \quad d(x) = d_2 > 0, \quad G(x) = G_2, \quad g(x) = g_2. \]

由于斑块 1 有利，$D_1 \ge D_2, d_1 \ge d_2, 1 \ge G_2, 1 \ge g_2$。图 1.6 展示了工程微生物扩散系数在空间中如何变化的一个例子。

在斑块边界 $x = x_i$（$x_{2m} = ml, x_{2m+1} = ml + l_1$, $m = 0, \pm 1, \pm 2, \ldots$），种群密度和通量连续，所以

\[ \lim_{x \to x_i^+} E(x, t) = \lim_{x \to x_i^-} E(x, t), \quad \lim_{x \to x_i^+} N(x, t) = \lim_{x \to x_i^-} N(x, t), \]

\[ \lim_{x \to x_i^+} D(x) \partial_x E(x, t) = \lim_{x \to x_i^-} D(x) \partial_x E(x, t), \quad \lim_{x \to x_i^+} d(x) \partial_x N(x, t) = \lim_{x \to x_i^-} d(x) \partial_x N(x, t). \]

现在数学问题已定义。我们要回答的关键问题是：(i) 工程生物在稀有时什么条件下能成功入侵？(ii) 如果入侵成功，工程物种是会把自然种群驱赶到入侵者主导状态，还是会达到共存状态？这里沿用 Shigesada et al. (1986) 的方法，考虑无穷区域上的问题，假设扩散和承载容量在不同斑块类型间变化。我们聚焦系统对从非常少的工程生物开始的入侵的稳定性。数学上，这意味着可以对稳态解做时空扰动的线性分析。

无量纲化。 通过引入 $e = a_E E$, $n = a_N N$, $t^* = r_E t$, $x^* = x\sqrt{r_E/D_1}$, $d^*(x) = d(x)/D_1$, $D^*(x) = D(x)/D_1$, $r = r_N/r_E$, $\gamma_e = b_E/a_N$, $\gamma_N = b_N/a_E$, $l^* = l\sqrt{r_E/D_1}$, $l_1^* = l_1\sqrt{r_E/D_1}$, $l_2^* = l_2\sqrt{r_E/D_1}$ 进行无量纲化。略去星号，无量纲模型方程为

\[ \partial_t e = \partial_x\left(D(x) \partial_x e\right) + e [G(x) - e - \gamma_e n], \]

\[ \partial_t n = \partial_x\left(d(x) \partial_x n\right) + r n [g(x) - n - \gamma_n e], \]

其中

\[ D(x) = \begin{cases} 1, & ml < x < ml + l_1, \\ D_2, & ml - l_2 < x < ml, \end{cases} \quad d(x) = \begin{cases} d_1, & ml < x < ml + l_1, \\ d_2, & ml - l_2 < x < ml, \end{cases} \]

函数 $G(x), g(x)$ 与 $(1.37)$–$(1.39)$ 中相同。在斑块边界 $x = x_i$（$x_i = ml$ 对 $i = 2m$，$x_i = ml + l_1$ 对 $i = 2m+1$）处，无量纲条件为

\[ \lim_{x \to x_i^+} e(x, t) = \lim_{x \to x_i^-} e(x, t), \quad \lim_{x \to x_i^+} n(x, t) = \lim_{x \to x_i^-} n(x, t), \]

\[ \lim_{x \to x_i^+} D(x) \partial_x e(x, t) = \lim_{x \to x_i^-} D(x) \partial_x e(x, t), \quad \lim_{x \to x_i^+} d(x) \partial_x n(x, t) = \lim_{x \to x_i^-} d(x) \partial_x n(x, t). \]

无斑块性与控制条件。 当整个区域有利时，不利斑块长度为 $l_2 = 0$。因此 $D(x) = 1, d(x) = d_1, G(x) = 1, g(x) = 1$ 处处成立。这回到上一节考虑过的带扩散的 Lotka–Volterra 竞争模型。

初始稳态简化为 $e_1 = 0, n_1 = 1$（自然主导稳态）。另外两个相关稳态是：入侵者主导稳态 $e_2 = 1, n_2 = 0$，以及共存稳态

\[ e_3 = \frac{\gamma_n - 1}{\gamma_n \gamma_e - 1}, \quad n_3 = \frac{\gamma_e - 1}{\gamma_n \gamma_e - 1}. \]

共存稳态只有在其为正时才有意义，这要求要么 $\gamma_e < 1$ 和 $\gamma_n < 1$（即两个物种的种间竞争都弱），要么 $\gamma_e > 1$ 和 $\gamma_n > 1$（即种间竞争都强）。平凡零稳态在此无意义。在这些特定竞争交互下，Turing 意义下（即零通量边界条件）没有其它稳态解。

如红松-灰松竞争所示，存在连接自然主导稳态 $(e_1, n_1)$ 与共存稳态 $(e_2, n_2)$ 或入侵者主导稳态 $(e_3, n_3)$ 的行波解。这些解对应微生物入侵波——或将自然物种驱赶到入侵者主导、或到新的（但较低的）稳态。

对初始自然主导稳态 $(e_1, n_1)$ 的通常线性稳定性分析（寻找形如 $e^{ikx + \lambda t}$ 的解）给出弥散关系（留作练习）：

\[ \lambda(k^2) = \frac{1}{2}\left[-b(k^2) \pm \sqrt{b^2(k^2) - 4c(k^2)}\right], \]

\[ b(k^2) = k^2(d_1 + 1) + (\gamma_e - 1) + r, \quad c(k^2) = d_1 k^4 + [\gamma_e - 1 + r] k^2 + r(\gamma_e - 1). \]

如果存在 $k^2$ 使 $\lambda(k^2) > 0$，则自然主导稳态线性不稳定。由弥散关系可知，若 $\gamma_e > 1$，则 $b(k^2), c(k^2)$ 始终为正，初始稳态总是线性稳定。但如果 $\gamma_e < 1$，则存在使稳态不稳定的 $k$ 值，工程物种 $e$ 的入侵将成功。

对其它稳态做类似线性化（另一个练习），可得稳定性：入侵者主导稳态 $(e_2, n_2)$ 在 $\gamma_n > 1$ 时稳定，$\gamma_n < 1$ 时不稳定；共存稳态 $(e_3, n_3)$ 在 $\gamma_n < 1$ 和 $\gamma_e < 1$ 时稳定，$\gamma_n > 1$ 和 $\gamma_e > 1$ 时不稳定。若 $\gamma_e < 1 < \gamma_n$ 或 $\gamma_n < 1 < \gamma_e$，则 $(1.47)$ 中 $e_3$ 或 $n_3$ 为负，共存稳态不再相关。平凡零稳态总是线性不稳定（因为在没有本地物种时自然菌株、工程菌株或两者都会入侵）。

由于上一节已经考虑过这种情形下的入侵行波，这里不再重复。总结：若自然主导稳态不稳定，则若 $\gamma_e < 1$，连接自然主导稳态与入侵者主导稳态的行波出现，但仅当 $\gamma_n > 1$ 时才如此。另一方面，连接自然主导稳态与共存稳态的行波仅在 $\gamma_n < 1$ 时出现。$\gamma_e < 1$ 和 $\gamma_n < 1$ 时的数值解与图 1.5 相似。

$\gamma_e < 1$（自然主导稳态不稳定）的需求，按 $(1.41)$ 的原始有量纲变量，意味着自然生物 $n$ 对工程物种 $e$ 的种间竞争效应弱于自然物种的种内竞争。若 $\gamma_n > 1$，自然物种被驱赶到入侵者主导；按原始有量纲参数，这发生在工程物种密度增加比自然物种的种群增长率更抑制其自身增长率时。若 $\gamma_n < 1$，情形相反。

另一方面，若自然主导稳态稳定（$\gamma_e > 1$），则当 $\gamma_n > 1$ 时入侵者主导稳态可同时稳定。这也是共存稳态不稳定的条件。在此情形下，自然主导稳态的稳定性依赖于初始条件 $(1.36)$。若 $H$ 代表对 $e = 0$ 的小扰动，则自然主导稳态仍为最终稳态解。然而 Cruywagen et al. (1996) 通过数值实验发现：对于对应于非常大的工程生物初始释放的极大扰动，行波解出现，入侵者主导稳态成为最终解。所以，对所有初始释放策略，只有当 $\gamma_e > 1$ 和 $\gamma_n < 1$ 时才能保证控制。

若考虑整个区域不利（取 $l_1 = 0$ 而非 $l_2 = 0$），则得到类似结果。然而此情形下非零稳态不同：自然主导稳态为 $e_1 = 0, n_1 = g_2$，入侵者主导稳态为 $e_2 = G_2, n_2 = 0$，共存稳态为

\[ e_3 = \frac{\gamma_n G_2 - g_2}{\gamma_e \gamma_n - 1}, \quad n_3 = \frac{\gamma_e g_2 - G_2}{\gamma_e \gamma_n - 1}. \]

稳定性条件现在由 $\gamma_e, \gamma_n$ 是否分别大于 $G_2/g_2$ 和 $\gamma_n < g_2/G_2$ 决定：共存稳态稳定，所有其它稳态不稳定。

空间变化扩散。 继续对各稳态做线性稳定性分析，但现在考虑扩散的空间变化，即斑块性影响扩散函数。这里只研究工程物种入侵能力如何受空间变化扩散系数的影响。

入侵条件。 初始自然主导稳态为 $(e_1, n_1)$，其中 $e_1 = 0$；$n_1$ 依赖于 $g(x)$，可能是 $n_1 = 1$ 或 $x$ 的周期函数（周期与斑块长度相关）。

为简化，首先假设 $g(x) = 1$，则 $n_1$ 与 $x$ 无关。Cruywagen et al. (1996) 考虑了 $g(x)$ 为 $x$ 周期函数时更为复杂的问题。

为确定初始自然主导稳态的稳定性，对 $(e_1, n_1) = (0, 1)$ 线性化得

\[ \partial_t e = \partial_x\left(D(x) \partial_x e\right) + e [G(x) - \gamma_e], \]

\[ \partial_t n = \partial_x\left(d(x) \partial_x n\right) + r [-n - \gamma_n e], \]

这里 $e, n$ 表示对稳态 $(n_1, e_1)$ 的小扰动（$|e| \ll 1, |n| \ll 1$）。

这里只通过工程物种的方程 $(1.51)$ 即可确定系统稳定性，因为它与 $n$ 无关。这把线性稳定性问题化简为分析

\[ \partial_t e = \partial_{xx} e + e (1 - \gamma_e) \text{ in Patch 1}, \quad \partial_t e = D_2 \partial_{xx} e + e(G_2 - \gamma_e) \text{ in Patch 2}. \]

代入 $e(x, t) = e^{-\lambda t} f(x)$ 给出 Hill 方程形式

\[ \frac{d}{dx}\left(D(x) \frac{df}{dx}\right) + [G(x) - \gamma_e + \lambda] f = 0, \]

其中 $G(x) - \gamma_e$ 和 $D(x)$ 都是周期 $l$ 的周期函数。Hill 方程周期系数理论告诉我们：存在一个严格递增的实特征值 $\lambda$ 的无穷序列，$-\infty < \lambda_0 < \tilde{\lambda}_1 \le \tilde{\lambda}_2 < \lambda_1 \le \lambda_2 \le \tilde{\lambda}_3 \le \tilde{\lambda}_4 < \cdots$，对应 $(1.55)$ 给出非零解。$\lambda = \lambda_i$ 时解周期为 $l$，$\lambda = \tilde{\lambda}_i$ 时解周期为 $2l$。$\lambda = \lambda_0$ 对应的解无零点，在空间意义下全局不稳定（即 $f \to \infty$ 当 $|x| \to \infty$）；Shigesada et al. (1986) 详细讨论了这点。详细理论见 Coddington 与 Levinson (1972)。

所以偏微分方程系统 $(1.53), (1.54)$ 的自然主导稳态稳定性由 $\lambda_0$ 的符号决定。$\lambda_0 < 0$ 时 $(1.55)$ 的平凡解 $e_0 = 0$ 动力学不稳定；$\lambda_0 > 0$ 时动力学稳定。Cruywagen et al. (1996) 得到了 $\lambda_0$ 的一个界，从而给出我们需要的控制条件，下面推导。

定义函数

\[ Q(x) = G(x) - \frac{l_1 + G_2 l_2}{l}, \]

可将 $(1.55)$ 写成

\[ \frac{d}{dx}\left(D(x) \frac{df}{dx}\right) + \left[Q(x) + \frac{l_1 + G_2 l_2}{l} - \gamma_e + \lambda\right] f = 0. \]

作为准备，先推导与方程

\[ \frac{d}{dx}\left(D(x) \frac{du}{dx}\right) + [\sigma + Q(x)] u = 0 \]

相关的一个结果，其中 $D(x), Q(x)$ 是周期 $l$ 的周期函数，$D(x) > 0$。由 Hill 方程理论知，对应最小特征值 $\sigma = \sigma_0$ 的周期 $l$ 的周期解 $u(x) = u_0(x)$ 无零点。可以假设 $u_0(x) > 0$ 对所有 $x$，然后定义积分因子 $h(x) = d[\ln u_0(x)]/dx$。于是 $h(x)$ 是周期 $l$ 的周期函数，满足

\[ \frac{d}{dx}[D(x) h(x)] + D(x) h^2(x) = -\sigma_0 - Q(x). \]

在一个周期长度 $l$ 上积分，由于 $D(x), h(x), Q(x)$ 都周期，得 $\int_{\zeta l}^{(\zeta+1)l} D(x) h^2(x) dx = -l \sigma_0$（对实数 $\zeta$）。于是如果 $h^2(x)$ 上的积分为零则 $\sigma_0 = 0$；否则因为 $D(x) > 0$，$\sigma_0 < 0$。又因为 $\int_{\zeta l}^{(\zeta+1)l} Q(x) dx = 0$（对任意实数 $\zeta$），比较 $(1.58)$ 和 $(1.59)$，用上面刚推出的结果，得

\[ \frac{l_1 + G_2 l_2}{l} - \gamma_e + \lambda_0 < 0. \]

于是得到 $\lambda_0 < 0$（即系统 $(1.53), (1.54)$ 不稳定）的充分条件：

\[ (1 - \gamma_e) l_1 \ge (\gamma_e - G_2) l_2. \]

现在有三种相关情形要讨论。记 $G_2 < 1$。

当 $\gamma_e > 1 > G_2$ 时，自然主导稳态在两个斑块孤立考虑时都稳定（参考上面关于有利或不利斑块稳定性条件的详细讨论；注意这里 $g_2 = 1$）。虽然看起来合理，但仅由 $(1.65)$ 不能得出自然主导稳态在全问题下稳定，因为 $(1.65)$ 只是不稳定的充分条件。

另一方面，当 $1 > G_2 > \gamma_e$ 时，自然主导稳态在两个斑块孤立考虑时都不稳定。不仅如此，正如所料，从 $(1.65)$ 可得自然主导稳态在所考虑的全区域问题上也不稳定。所以，如果工程微生物在不利斑块（斑块 2）的承载容量（用 $G_2$ 表示）超过其因自然生物种间竞争效应（用 $\gamma_e$ 表示）所受的损失，则工程生物总会入侵。

然而，如果 $1 > \gamma_e > G_2$，自然主导稳态在有利斑块不稳定但在不利斑块孤立考虑时稳定。哪个斑块主导自然主导稳态的实际稳定性依赖于它们的相对大小，可从不等式 $(1.65)$ 看出。增大有利斑块长度 $l_1$ 和/或减小不利斑块长度 $l_2$，自然主导稳态将变得不稳定，从而入侵确实发生。所以条件 $(1.65)$ 不是不稳定的必要条件，因此并未给出保证自然主导稳态稳定性的精确条件。

我们从为两种斑块推导 $(1.53), (1.54)$ 的可分离解开始。因为预期周期解，所以用 Fourier 级数展开求解。

在斑块 1 中经过一些代数得

\[ e(x, t) = \sum_{i=0}^{\infty} A_i e^{-\gamma_i t} \cos\left[\frac{|x - l_1/2 - ml|}{1 - \gamma_e + \lambda_i}\right], \]

在斑块 2 中

\[ e(x, t) = \sum_{i=0}^{\infty} B_i e^{-\lambda_i t} \cos\left[\sqrt{\frac{|x + l_2/2 - (m+1)l|}{(G_2 - \gamma_e + \lambda_i)/D_2}}\right], \]

其中 $A_i, B_i$ 为常数。应用连续性条件 $(1.45), (1.46)$，下述等式必须成立：

\[ \sqrt{1 - \gamma_e + \lambda_i} \tan\left[\frac{l_1}{2} \sqrt{1 - \gamma_e + \lambda_i}\right] = -D_2 \sqrt{\frac{G_2 - \gamma_e + \lambda_i}{D_2}} \tan\left[\frac{l_2}{2} \sqrt{\frac{G_2 - \gamma_e + \lambda_i}{D_2}}\right], \]

对 $i = 0, 1, 2, \ldots$。如果平方根内表达式为负，则用恒等式 $\tan iz = i \tanh z, \arctan iz = i \operatorname{arctanh} z$。

我们关心满足上述等式的最小特征值 $\lambda = \lambda_0$ 的符号。容易证明 $\lambda_0 < 0$ 当且仅当 $1 - \gamma_e + \lambda_0$ 和 $G_2 - \gamma_e + \lambda_0$ 在平方根内有相反符号。由于 $G_2 < 1$，这只在 $\gamma_e < 1$ 时才发生。所以，由于在空间均匀系数问题中就是这种情形，自然主导稳态不稳定的必要条件（即允许工程微生物入侵）是自然物种对工程物种的竞争效应 $b_E$ 小于自然物种的种内竞争效应 $a_N$；参见 $(1.41)$ 中的无量纲形式。

若 $G_2 \ge \gamma_e$，则 $\lambda_0 < 0$，入侵必将成功，与其它参数和斑块大小无关，如前所述。然而若 $G_2 < \gamma_e < 1$，则根据各参数值不同，$\lambda_0$ 可正可负。下面更详细地考虑这种情形，即自然主导稳态在有利斑块不稳定但在不利斑块稳定的情况。如上面所看到的，斑块的相对大小此时变得重要。

在临界值 $\lambda_0 = 0$ 处，下式成立：

\[ \sqrt{1 - \gamma_e} \tan\left[\frac{l_1}{2} \sqrt{1 - \gamma_e}\right] = D_2 \sqrt{\frac{\gamma_e - G_2}{D_2}} \tanh\left[\frac{l_2}{2} \sqrt{\frac{\gamma_e - G_2}{D_2}}\right], \]

由此确定斑块 1 的临界长度 $l_1^*$ 为

\[ l_1^* = \frac{2}{\sqrt{1 - \gamma_e}} \arctan\left[\sqrt{\frac{D_2(\gamma_e - G_2)}{1 - \gamma_e^3}} \tanh\left(\frac{l_2}{2} \sqrt{\frac{\gamma_e - G_2}{D_2}}\right)\right]. \]

对 $l_1 < l_1^*$，自然主导稳态将稳定（$\lambda_0 > 0$）；对 $l_1 > l_1^*$，$\lambda_0 < 0$，自然主导稳态不稳定。所以如本节前面所示，如果有利斑块相对不利斑块足够大，入侵就会成功。

注意当 $l_2 \to \infty$ 时，边界曲线渐近到

\[ \lim_{l_2 \to \infty} l_1(l_2) = l_1^c = \frac{2}{\sqrt{1 - \gamma_e}} \arctan\left[\sqrt{\frac{D_2(\gamma_e - G_2)}{1 - \gamma_e}}\right]. \]

所以，如果 $l_1 \ge l_1^c$，无论不利斑块大小，入侵总会成功。又因为

\[ l_1^c \le l_1^m = \frac{2 \arctan \infty}{\sqrt{1 - \gamma_e}} = \frac{\pi}{\sqrt{1 - \gamma_e}}, \]

如果 $l_1 \ge l_1^m$，无论 $l_2, G_2(<\gamma_e), D_2$ 取何值，入侵都会成功。当 $G_2 < \gamma_e$ 时，$l_1$ 和 $l_2$ 平面上的稳定区域如图 1.7(a) 所示。

类似地可画出 $\gamma_e$ 对 $l_2$ 的稳定曲线。已经证明若 $\gamma_e < G_2$ 入侵总成功，与斑块 2 长度 $l_2$ 无关。然而当 $\gamma_e$ 增大超过 $G_2$ 时，稳定曲线从某临界值 $\gamma_e = \gamma_e^c$ 处的无穷远处出现。曲线的渐近 $\gamma_e^c$ 由下述非线性关系得到：

\[ D_2 \frac{\gamma_e^c - G_2}{1 - \gamma_e^c} = \tan^2\left[\frac{l_1}{2} \sqrt{1 - \gamma_e^c}\right]. \]

此稳定区域如图 1.7(b) 所示。注意当 $l_1$ 增大趋向 $l_1^c$ 时，稳定曲线会出现在越来越大的 $\gamma_e^c$ 值处，而当 $l_1 \ge l_1^c$ 时稳定曲线根本不会出现。

$G_2$ 对 $l_2$ 的稳定图与 $\gamma_e$ 对 $l_2$ 的图性质类似。$G_2 > \gamma_e$ 时入侵成功，然而 $G_2$ 值减小超过 $\gamma_e$ 时稳定曲线从无穷远处在渐近

\[ G_2^c = \gamma_e + \frac{\gamma_e - 1}{D_2} \tan^2\left[\frac{l_1}{2} \sqrt{1 - \gamma_e}\right] \]

处出现。工程物种在不利斑块（斑块 2）的扩散率 $D_2$ 在确定自然主导稳态稳定性时也起重要作用（当 $l_1 < l_1^*$ 时）。对足够小的 $D_2$ 值，无论 $l_2$ 取何值入侵都成功。从生物学角度看，这意味着不利斑块中扩散太小，对有利斑块影响最小。然而当 $D_2$ 增大时，稳定曲线从临界渐近值 $D_2 = D_2^c$ 处的无穷远处出现，

\[ D_2^c = \frac{1 - \gamma_e}{\gamma_e - G_2} \tan^2\left[\frac{l_1}{2} \sqrt{1 - \gamma_e}\right]. \]

另一方面，当 $D_2 \to \infty$ 时，稳定曲线趋近渐近线 $l_2^c$。由于从 $(1.70)$，

\[ l_2 = 2 \sqrt{\frac{D_2}{\gamma_e - G_2}} \operatorname{arctanh}\left[\sqrt{\frac{1 - \gamma_e}{D_2(\gamma_e - G_2)}} \tan\left(\frac{l_1}{2} \sqrt{1 - \gamma_e}\right)\right], \]

应用 L'Hôpital 法则得

\[ \lim_{D_2 \to \infty} l_2(D_2) = l_2^c = \frac{2\sqrt{1 - \gamma_e}}{\gamma_e - G_2} \tan\left[\frac{l_1}{2} \sqrt{1 - \gamma_e}\right]. \]

所以 $l_2 < l_2^c$ 时工程微生物对任何扩散率都入侵。图 1.7(c) 展示 $D_2$ 对 $l_2$ 的稳定曲线。

可以得出结论：对 $\gamma_e < 1$，控制的必要条件是 $G_2 < \gamma_e$，以及

\[ l_1 < l_1^c, \quad l_2 > l_2^c, \quad \gamma_e > \gamma_e^c, \quad G_2 < G_2^c, \quad D_2 > D_2^c. \]

这些不等式表明：通过减小 $l_1$ 或 $G_2$，或增大 $\gamma_e, D_2$ 或 $l_2$ 可以保证控制。回想前面，$\gamma_e < 1$ 意味着自然主导稳态在有利斑块不稳定，而 $G_2 < \gamma_e$ 意味着在不利斑块稳定。然而保证自然主导稳态稳定的最简策略仍然是 $\gamma_e > 1$，如上所述。

需要强调的是：上述线性稳定性分析是针对整个斑块化区域的。即使在整区域上入侵成功，依赖于初始条件，入侵可能只是局部的，事实上被局限在某个斑块中。当 $1 > \gamma_e > G_2$ 时可能就是这种情况。虽然自然主导稳态在完整问题上可以稳定，但在有利斑块中局部不稳定。由该斑块中 $e$ 的小非零初始扰动引起的局部入侵在某些情况下可能会发生。Cruywagen et al. (1996) 考虑过这样的例子。

至此得到了自然主导稳态 $(e_1, n_1)$ 的稳定性条件。其它三个可能稳态是零稳态 $(e_0, n_0)$、入侵者主导稳态 $(e_2, n_2)$、共存稳态 $(e_3, n_3)$。可以确定在什么条件下解会演化为这些稳态。理所当然地，这些情形下预计会出现连接不稳定自然主导稳态与非自然主导稳态的行波解。Cruywagen et al. (1996) 用类似方法考察了这些其它稳态，证明 $\gamma_e > 1$ 和 $\gamma_n < 1$ 是对任何初始微生物释放密度都安全的唯一控制策略。比较这个结果与前面得到的，变化扩散并不影响稳定性条件。

Cruywagen et al. (1996) 还通过把 $G(x), g(x)$ 视为空间周期函数考虑了变化承载容量的效应，这意味着三个非零稳态都有空间周期解。确定其稳定性很困难，所以用摄动技术得到近似的稳定性条件。

虽然 Cruywagen et al. (1996) 的模型实际上过于简单而不现实，但它确定了异质环境影响竞争系统入侵的几个关键场景。尽管模型相对简单，详细分析却很复杂，涉及大量簿记。Cruywagen et al. (1996) 给出了相当复杂的总结表，预测转基因微生物与野生型菌株竞争时的共存、灭绝和入侵，分别按各种参数以及扩散和可变承载容量（含斑块大小）的形式列出。首先，尽管对合适的栖息地岛屿周围较大的不利"护城河"，入侵的可能性较小，但一旦内部岛屿变得足够大，周围恶劣区域无论多大都不能阻止入侵。一般地，缩减本地物种占优势的斑块类型的相对大小，会增加外来物种入侵的机会。也许最有趣的场景是：入侵在某个斑块内局部成功，但在全局失败。正是在这种情形下结果不直观，数学模型特别有用。

Lewis et al. (1996) 对包含对流输运以及静止和可动类别区分的模型做了进一步工作；另见 Lewis 与 Schmitz (1996)。该工作考虑了微生物进入不同的可动和不可动部分，如植物根、各种宿主、地下水或风。Lewis (1997) 的论文和 Tilman 与 Kareiva (1997) 编辑的论文集对入侵物种的空间扩散尤其相关。最后，对所考虑基本模型的任何研究和扩展都应与现场研究相结合，以估计种间竞争参数、确保现实的现场预测（Kareiva 1990）。

1.5 Belousov–Zhabotinskii 反应中的行波前沿（Travelling Fronts in the Belousov–Zhabotinskii Reaction）

研究 Belousov–Zhabotinskii (BZ) 反应中化学波的一个原因是：它只是众多展示类似波和图案形成现象的反应之一。这些反应被用作生物图案形成系统的典范，并因此引发了大量胚胎发育实验，极大增强了我们对心脏和身体其它器官复杂波行为的理解。所以，尽管研究这些化学波本身就有趣，我们在这里讨论它们，是为了用它们来加深理解生物图案形成机制这个教学目的。

图 1.1(a) 中的波是 BZ 反应中由局部起搏点产生的化学浓度行波带。本节推导并分析一个描述这种波传播前沿的模型。远离中心处波本质上是平面的，所以我们考虑一维问题，部分沿用 Murray (1976) 的分析。研究这个具体问题的理由是假设波前沿速度主要取决于关键化学物质：溴酸（HBrO₂）和溴离子（Br⁻）的浓度（分别记为 $x, y$）。参考 Volume I 第 8 章 §8.1 关于模型反应动力学的细节。本节可以独立阅读，只需从下面的 $(1.80)$ 反应图开始。假设这些反应物以扩散系数 $D$ 扩散。我们相信波前沿由反应过程 I 主导，即 (i) 把溴离子浓度减小到小值、(ii) 把溴酸增大到最大浓度、(iii) 铈离子催化剂处于 Ce³⁺ 状态的反应序列。由于 Volume I §8.1 中 Ce⁴⁺ 的浓度记为 $z$，最后的假设意味着 $z = 0$。从 Volume I 的 $(8.2)$ 出发（无铈反应且 $z = 0$）的简化反应序列为

\[ A + Y \xrightarrow{k_1} X + P, \quad X + Y \xrightarrow{k_2} 2P, \quad A + X \xrightarrow{k_3} 2X, \quad 2X \xrightarrow{k_4} P + A, \]

其中 $X, Y$ 分别代表溴酸和溴离子，$k$ 为速率常数。$P$（化合物 HBrO）不出现在分析中，$A$（BrO₃⁻）的浓度是常数。

对这个图应用质量作用定律（参见 Volume I 第 6 章），用小写字母表示浓度，加入 $X, Y$ 的扩散项，得到

\[ \partial_t x = k_1 a y - k_2 x y + k_3 a x - k_4 x^2 + D \partial_{ss} x, \]

\[ \partial_t y = -k_1 a y - k_2 x y + D \partial_{ss} y, \]

其中 $s$ 是空间变量。适当的无量纲化（$u = k_4 x/(k_3 a)$, $v = k_2 y/(k_3 a r)$, $s^* = \sqrt{k_3 a/D} \, s$, $t^* = k_3 a t$, $L = k_1 k_4/(k_2 k_3)$, $M = k_1/k_3$, $b = k_2/k_4$）后，$(1.81)$ 化为（略去星号）

\[ \partial_t u = L r v + u(1 - u - r v) + \partial_{ss} u, \]

\[ \partial_t v = -M v - b u v + \partial_{ss} v. \]

用 Volume I 第 8 章 $(8.4)$ 给出的各速率常数和参数的估计值，可得 $L \approx M = O(10^{-4})$, $b = O(1)$，参数 $r$ 可在约 5 到 50 之间变化。

无量纲化 $(1.82)$ 后，现实稳态为 $u = v = 0$ 和 $u = 1, v = 0$，可预期 $u, v$ 都是 $O(1)$ 有界。所以，作为一阶近似，由于 $(1.83)$ 中 $L \ll 1, M \ll 1$，可以忽略这些项，从而得到 BZ 反应行波前沿的模型

\[ \partial_t u = u(1 - u - r v) + \partial_{ss} u, \quad \partial_t v = -b u v + \partial_{ss} v, \]

其中 $r, b$ 是 $O(1)$ 的正参数。注意这个模型近似引入了一个新稳态 $(0, S)$，$S > 0$ 可取任意值。这是因为这只是一个前沿的模型而非整个波脉冲，在波两侧 $v \to 0$。

现在寻找 $(1.85)$ 的行波前沿解，波从高溴酸浓度区移动到低溴酸浓度区，同时降低溴离子水平。用 $(1.84)$ 寻找边界条件为 $u(-\infty, t) = 0, v(-\infty, t) = 1, u(\infty, t) = 1, v(\infty, t) = 0$ 的行波，波向左移动。

寻找行波解前应注意一些特殊情形，它们把问题化为 Fisher–Kolmogoroff 方程。设 $v = (1 - b)/r (1 - u)$（$b \ne 1, r \ne 0$），系统 $(1.85)$ 化为

\[ \partial_t u = b u (1 - u) + \partial_{ss} u, \]

即 Volume I 的 Fisher–Kolmogoroff 方程 $(13.4)$，它有从 $u = 0$ 到 $u = 1$ 的单调行波前沿解，速度为 $c \ge 2\sqrt{b}$。由于只关心非负的 $u, v$，必须在 $(1.87)$ 中 $b < 1$。取初始条件 $u(s, 0) \sim O(\exp[-\beta s])$ 当 $s \to \infty$，由 Volume I 第 13 章得行波前沿的渐近速度为

\[ c = \begin{cases} \beta + b/\beta, & 0 < \beta \le \sqrt{b}, \\ 2\sqrt{b}, & \beta > \sqrt{b}. \end{cases} \]

由 Fisher–Kolmogoroff 波加 $(1.87)$ 给出的行波前沿解并不具有实际意义，除非 $1 - b = r$（因为要求 $u = 1$ 时 $v = 0$、$u = 0$ 时 $v = 1$）。取合适初始条件 $u(s, 0) = 0$（$s < s_1$）、$h(s)$（$s_1 < s < s_2$）、$1$（$s > s_2$），其中 $h(s)$ 是正的单调连续函数，$h(s_1) = 0, h(s_2) = 1$，则波速为 $c = 2\sqrt{b} = 2\sqrt{1 - r}$（由 $(1.88)$），$v = 1 - u$。必须有 $0 < r \le 1$。

可以进一步利用 Fisher–Kolmogoroff 方程的结果和抛物方程的最大值原理。设 $u^f(s, t)$ 是

\[ \partial_t u^f = u^f(1 - u^f) + \partial_{ss} u^f, \quad u^f(-\infty, t) = 0, \quad u^f(\infty, t) = 1 \]

以初始条件 $(1.89)$ 开始的唯一解。渐近行波前沿解速度为 $c = 2$。写 $w(s, t) = u(s, t) - u^f(s, t)$，让 $u(s, t)$ 与 $u^f$ 有同样的初始条件 $(1.89)$。从 $(1.85)$ 的 $u$ 方程减去 $(1.90)$ 并用 $w$ 的定义，最后得

\[ w_{ss} - w_t + [1 - (u + u^f)] w = r u v. \]

我们把解限制在 $0 \le u \le 1$ 内；由于 $0 \le u^f \le 1$，有 $[1 - (u + u^f)] \le 1$，所以不能立即用通常的最大值原理。设 $W = w \exp[-K t]$，$K > 0$ 是有限常数，最后一个方程变为

\[ W_{ss} - W_t + [1 - (u + u^f) - K] W = r u v e^{-K t} \ge 0. \]

选 $K > 1$，则 $[1 - (u + u^f) - K] < 0$，最大值原理可对 $W$ 方程使用。它说明 $W$（从而 $w$）在 $t = 0$ 或 $s = \infty$ 处取最大值。但 $w_{\max} = (u - u^f)_{\max} = 0$ 在 $t = 0$ 和 $s = \pm \infty$ 处，所以

\[ u(s, t) \le u^f(s, t) \text{ 对所有 } s, t > 0. \]

这说明 $(1.85)$ 的 $u$ 解在所有点上都小于或等于从初始条件 $(1.89)$ 演化的 Fisher–Kolmogoroff 解 $u^f$。所以，如果 $(1.85)$ 的解有边界条件 $(1.86)$ 和等价于 $(1.89)$ 的初始条件的行波解，则波速 $c$ 必须受 Fisher–Kolmogoroff 波速的界限制，对所有 $r, b$ 值有 $c(r, b) \le 2$。直观上我们预期 $(1.85)$ 任何这样的行波解都有 $c \le 2$，因为 $u v \ge 0$，第一式中 $-r u v$ 项像是除动力学 $u(1 - u)$ 之外的汇项。这抑制 $u$ 在任意点的增长，与 Fisher–Kolmogoroff 波解相比，我们预期 $u$ 和它的速度被 Fisher–Kolmogoroff 解从上方界定。

波速 $c$ 关于 $r, b$ 的各种极限值可以从方程系统 $(1.85)$ 导出。推导时必须小心，因为存在非一致极限情形；这些将在适当地方指出。

若 $b = 0$，$(1.85)$ 中 $v$ 方程变为基本扩散方程 $v_t = v_{ss}$，不能有波解。这意味着 $(1.85)$ 中 $u$ 的第一个方程也不能有波解，因为波解要求 $u, v$ 有相同的传播速度。这说明极限 $c(b \to 0, r) = 0$（$r > 0$）。若 $b \to \infty$，$(1.85)$ 给出 $v = 0$（排除平凡解 $u = 0$），此时 $c(b \to \infty, r) = 2$ 对所有 $r \ge 0$。若 $r = 0$，$u, v$ 方程不耦合，$u$ 方程为基本 Fisher–Kolmogoroff 方程 $(1.90)$，以初始条件 $(1.89)$ 时有波速 $c = 2$ 的波前解；这意味着相关 $v$ 解也有速度 2。这给出极限情形 $c(b, r \to 0) = 2$（$b > 0$）。若 $r \to \infty$，则 $u = 0$ 或 $v = 0$，两者都意味着没有波解，所以 $c(b, r \to \infty) = 0$。注意这个情形的非一致极限：$r \to \infty$ 时 $v \ne 0$ 与 $v = 0$ 然后 $r \to \infty$ 的情形不同。在后者中，$u$ 由 Fisher–Kolmogoroff 方程决定，$r$ 无关紧要。然而如上所述，我们关心的是 $u, v$ 都不恒为零的行波。综合起来：

\[ c(0, r) = 0, r > 0; \quad c(\infty, r) = 2, r \ge 0; \]

\[ c(b, 0) = 2, b > 0; \quad c(b, \infty) = 0, b \ge 0. \]

第一式并未给出小 $b$ 时的完整图景（见下文）。

对系统 $(1.85)$ 的行波前沿问题，用行波变换 $u(s, t) = f(z), v(s, t) = g(z), z = s + ct$ 和边界条件 $(1.86)$，化为

\[ f'' - c f' + f(1 - f - r g) = 0, \quad g'' - c g' - b f g = 0, \]

\[ f(\infty) = g(-\infty) = 1, \quad f(-\infty) = g(\infty) = 0. \]

利用各种关于 $(1.92)$ 的 $f \ge 0, g \ge 0$ 单调解的界和估计技术，Murray (1976) 得到 $c$ 关于 $r, b$ 的一般界：

\[ \left[\left(\frac{r^2 + 2b}{3}\right)^{1/2} - r\right] [2(b + 2r)]^{-1/2} \le c \le 2. \]

系统 $(1.85)$ 带初始和边界条件 $(1.86), (1.89)$ 被数值求解（Murray 1976），部分结果见图 1.8。注意在图 1.8(b) 中由 $b = 0, c^2 = 4b$ 和 $c = 2$ 围成的区域内非负解不存在。极限曲线 $c^2 = 4b$ 由特殊情形 $v = 1 - u$ 和 $b = 1 - r$（$r < 1$）得到。Manoranjan 与 Mitchell (1983) 对模型系统 $(1.85)$ 进行了更全面的数值研究。

现在回到实验情形。从图 1.8(a)，如果固定 $b$，增大 $r$（按 $(1.82)$ 的无量纲化等价于增加上游溴离子浓度）的效果是展平 $v$ 曲线。即波前沿变得不那么陡。另一方面，固定 $r$ 增大 $b$ 时，前沿变陡。虽然不精确，但我们可以从计算出的波前沿解的宽度 $\omega$ 估计实际波前沿宽度。有量纲下这是 $\omega_D$，由 $(1.82)$ 得 $\omega_D = \sqrt{D/(k_3 a)} \, \omega \approx 4.5 \times 10^{-4} \omega$ cm，其中取 $D \approx 2 \times 10^{-5}$ cm²/s（合理小分子的典型值）和 $k_3 a \approx 10^2$ s⁻¹（从第 8 章 $(8.4)$ 的参数值得到）。从图 1.8(a) 得 $\omega$ 约为 10，所以 $\omega_D$ 约为 $10^{-3}$ cm。这与实验观察到的量级一致；前沿非常薄。

从图 1.8(b) 得到的另一个实际预测是：$b$ 大于约 2 时，对固定 $r$ 波速相对独立于 $b$。对 $b$ 大到约 50 的计算证实了这一点。实验上也观察到了。

从无量纲化 $(1.82)$，有量纲波速 $c_D$ 为 $c_D = \sqrt{k_3 a D} \, c(r, b)$，其中 $r$ 是上游溴离子浓度的度量，$b = k_2/k_4$。从 Volume I 第 8 章 $(8.4)$ 的参数估计，$b \approx 1$。赋值 $r$ 不容易，实验上 5–50 是合理的。取 $r$ 为 $O(10)$、$b$ 约为 1，从图 1.8(b) 得无量纲波速为 $O(10^{-1})$；$c$ 的精确值可从模型系统计算。用上面的 $D, k_3 a$ 值，$c_D$ 为 $O(4.5 \times 10^{-3}$ cm/s$)$ 或 $O(2.7 \times 10^{-1}$ cm/min$)$，又是在实验观察范围内。鉴于与实验的合理定量比较以及本节从模型推导出的结果，我们认为 BZ 波前沿的传播速度主要由前缘决定而非后缘。

最后，关于反应扩散波前沿与简单扩散的传播速度比较：让波前沿移动 1 cm 的时间为 $O(10/2.7$ min$)$，即约 4 分钟，而扩散时间 $O(1$ cm²/D$)$，即 $O(5 \times 10^4$ s$)$ 约 850 分钟。所以，作为通过化学浓度变化传递信息的手段，可以说反应扩散波比纯扩散快几个数量级（距离不是很小时）。后文将详细讨论胚胎学背景下的图案形成问题，其中关心的距离是细胞直径量级，所以扩散再次是传递信息的相关机制。然而我们将看到，在胚胎学背景下，扩散不是唯一可能的机制。

模型系统 $(1.85)$ 已被多个作者研究。Gibbs (1980) 例如证明了行波的存在性和单调性。Ortoleva 与 Schmidt (1985) 把这类行波现象有趣地表述为 Stefan 问题，并给出大 $r, b$ 极限下对应 Murray 模型 $(1.85)$ 的奇异摄动分析。

Showalter 和他的同事们对行波前沿领域做出了重要贡献，主要与实验有关，许多涉及 BZ 反应。例如 Merkin et al. (1996) 考察了基于带化学反馈的简单三次自催化模型的波致混沌。他们的模型涉及用反应速率 $k_1 a b^2$ 表示的三次自催化反馈 $A + 2B \to 3B$，耦合 $B \to C$ 速率 $k_2 b$ 的衰变步。设置包括凝胶反应区，一侧是保持恒定 $a_0$ 的储库可交换 $a$，另一侧是保持 $b_0$ 的储库可交换 $b$。反应扩散系统为

\[ \partial_t a = D_A \nabla^2 a + k_f(a_0 - a) - k_1 a b^2, \]

\[ \partial_t b = D_B \nabla^2 b + k_f(b_0 - b) + k_1 a b^2 - k_2 b, \]

其中 $k_f$ 是与反应物流入流出相关的常数参数，$D$ 是各自的扩散系数。虽然这些方程看起来与许多已研究过的方程没太大不同（除反应细节外），但它们展示出有趣的解谱。Merkin et al. (1996) 通过线性、行波和数值（一维和二维）分析表明，存在无量纲参数 $\mu = k_1 a_0^2/k_f$ 和 $\phi = (k_f + k_2)/k_f$ 的范围，使行波前沿、行波脉冲、Hopf 分岔以及有趣的混沌行为都可出现，后者是行波后被诱导的。

1.6 可激发介质中的波（Waves in Excitable Media）

可激发行为的最广泛研究系统之一是神经细胞通过电信号的神经通信。Volume I 第 7 章讨论了重要的 Hodgkin–Huxley 模型，并推导了它的数学漫画 FitzHugh–Nagumo (FHN) 方程。本节首先通过一个例子考察 FHN 模型的时空行为，并展示只在超过某个阈值扰动时才能传播的行波脉冲的存在。这里说的"脉冲"是指从稳态出发再回到稳态的行波，像水上的孤立波（参见图 1.10）。考虑动力学为 $u_t = f(u, v), v_t = g(u, v)$ 的模型，针对 $f, g$ 的特定类。所描述的方法相当一般，可应用于一大类可激发介质定性模型，它们的零斜线与图 1.9(a) 定性类似。本节可以独立阅读，只需把方程系统 $(1.94)$ 视为可激发介质的具体模型例子。

Volume I §7.5 中讨论过无空间变化（即空间钳制）情形下的 FHN 方程，展示了如图 7.12、§3.8 所述的时间阈值行为。无施加电流（$I_a = 0$）的 FHN 系统，但允许跨膜电位的空间"扩散"，并为本章记号稍作修改，得

\[ \partial_t u = f(u) - v + D \partial_{xx} u, \quad \partial_t v = b u - \gamma v, \quad f(u) = u(a - u)(u - 1). \]

这里 $u$ 直接与膜电位相关（Volume I §7.5 中的 $V$），$v$ 起多种变量的作用（贡献于钠、钾和其它离子的膜电流）。"扩散"系数 $D$ 与轴突中的轴向电流相关，参考 §7.5 中的电流守恒方程 $(7.38)$，电位的空间变化 $V$ 在右端贡献 $(d/(4 r_i)) V_{xx}$，其中 $r_i$ 是电阻率，$d$ 是轴突直径。参数 $0 < a < 1$, $b$, $\gamma$ 都为正。$(u, v)$ 平面上的动力学零斜线如图 1.9(b) 所示。

本节要展示具有可激发动力学的反应扩散系统如何产生行波解。除了用 FHN 模型模拟神经动作电位的传播外，还有几个重要的生理学应用。一个重要应用是肌肉组织（特别是心肌）中产生的波：在二维和三维情形下，这些可激发波与心房扑动和纤颤问题密切相关（参见 Winfree 1983a, b）。另一个例子是大脑皮层的回响性皮层抑郁波（Shibata 与 Bureš 1974）。二维和三维可激发波也可以在 BZ 反应和其它反应中出现。下面会回到这些应用。系统 $(1.94)$ 已被详细研究过，下面只是来自长长不断增加的参考文献的一小部分。Rinzel (1981) 的综述专门讨论了神经生物学中的模型。Rinzel 与 Keller (1973) 考虑了 $(1.94)$ 的分段线性漫画（$\gamma = 0$），得到了行波脉冲和周期波列的解析结果；这里的分析方法具有更广的适用性，常常是非线性问题解析的唯一方法。McKean (1970) 研究了把 $f(u)$ 替换为分段线性近似 $f(u) = H(u - a) - u$（$H$ 是 Heaviside 函数：$H(x) = 0$ 若 $x < 0$，$H(x) = 1$ 若 $x > 0$）的漫画形式；Feroe (1982) 研究了该漫画多脉冲解的稳定性。Ikeda et al. (1986) 考虑了 Hodgkin–Huxley 系统，证明了某些慢波解的不稳定性。Rinzel 与 Terman (1982) 研究了 $(1.94)$ 中 $b, \gamma$ 使得 $v = b u/\gamma$ 与 $u$ 零斜线相交给出三个稳态的情形。Keener (1980)、Zykov (1988)、Tyson 与 Keener (1988) 对可激发介质中的波作了综述讨论；Carpenter (1979) 讨论了可激发膜中的周期暴发现象。Keener 与 Sneyd (1998) 的书也讨论了这一现象。

$(1.94)$ 的行波解中 $u, v$ 只是行波坐标 $z = x - ct$ 的函数，满足 $(1.94)$ 的行波坐标形式

\[ D u'' + c u' + f(u) - v = 0, \quad c v' + b u - \gamma v = 0, \quad z = x - ct, \]

其中撇号表示对 $z$ 求导，波速 $c$ 待定。对应孤立脉冲的边界条件为

\[ u \to 0, \quad u' \to 0, \quad v \to 0 \text{ 当 } |z| \to \infty, \]

脉冲典型形式如图 1.10(a) 所示。$(u, v)$ 平面上的对应相轨迹示意于图 1.10(b)。

初始条件在行波脉冲存在性中起关键作用。可以直观地看出原因。假设空间域内 $(u, v)$ 初始在零静止态，并对小区域内的 $u$ 作局部升高、$v = 0$ 的扰动，如图 1.11(a)。如果扰动最大值 $u$ 小于图 1.10(b) 中的阈值 $u_A$（也是图 1.11(c)），则动力学使 $u$ 回到原点，空间扰动简单地消失。另一方面，如果扰动大于阈值 $u_A$，则动力学在 $u$ 和 $v$ 中引发大的偏移，如图 1.11(c) 中 $0 B C D 0$ 所示。波被引发时，尾缘在相平面上由 $C D$ 表示。虽然直观上看前缘应在 $0 B$，$C D$ 的位置不那么明显。下面考虑行波脉冲的这一重要方面。

分析上更容易看出在做什么，如果我们考虑 $(1.94)$ 中 $b, \gamma$ 小的情形，写 $b = \varepsilon L, \gamma = \varepsilon M, 0 < \varepsilon \ll 1$，则 $(1.94)$ 变为

\[ u_t = D u_{xx} + f(u) - v, \quad v_t = \varepsilon(L u - M v). \]

现在参考图 1.10(a)，考虑前缘 $0 A B$。在极限情形 $\varepsilon \to 0$ 中，第二个方程说 $v \approx$ 常数，从图 1.10(a), 1.10(b) 看这个常数为零。$(1.97)$ 中 $u$ 方程变为

\[ u_t = D u_{xx} + f(u), \quad f(u) = u(a - u)(u - 1), \]

其中 $f(u)$ 作为 $u$ 的函数在图 1.10(b) 中示意。它有三个稳态 $u = 0, a, 1$。在无扩散时 $(1.98)$ 意味着 $u = 0$ 和 $u = 1$ 线性稳定，$u = a$ 不稳定。于是可以有连接 $u = 0$ 与 $u = 1$ 的行波解。$(1.98)$ 是 Volume I §13.5 中研究的特例，具体是方程 $(13.73)$，它有精确解析解 $(13.78)$，唯一波速由 $(13.77)$ 给出。所以这里的波解有

\[ u = u(z), \quad z = x - ct; \quad c = \sqrt{\frac{D}{2}} (1 - 2a), \]

波速只在 $a < 1/2$ 时为正。这与 Volume I $(13.70)$ 的符号决定给出的条件一致；对 $(1.98)$，

\[ c \gtrless 0 \text{ 若 } \int_0^1 f(u) du \gtrless 0. \]

参考图 1.10(b)，由 $0 A$ 和曲线 $v = f(u)$ 围成的面积小于 $A B$ 和曲线 $v = f(u)$ 围成的面积，所以 $c > 0$；积分得 $c > 0$ 对所有 $a < 1/2$。

我们通过忽略 $(1.97)$ 中 $\varepsilon$ 项推导出波脉冲前沿方程组。这给出对应相平面中 $0 A B$ 的脉冲贡献。沿 $BC$，$v$ 变化。由 $(1.97)$ 知 $v$ 的变化需要长时间，实际上是 $O(1/\varepsilon)$，因为 $v_t = O(\varepsilon)$。要得到解的这一部分，需要做奇异摄动分析（参见例如 Keener 1980），其要点是给出一个慢过渡期，$u$ 变化不大而 $v$ 变化较大。这对应图 1.10(a), 1.10(b) 中标为 $BC$ 的脉冲部分。

一个关键问题立即出现：下一次快过渡发生在哪里，换言之 $C$ 在相轨迹上的位置。记住我们正在研究以不变形状传播的脉冲解存在性。为此，尾缘（即从 $C$ 经 $P$ 到 $D$ 的波前沿解）的速度必须与前缘 $0 A B$ 的速度相同（图 1.10(a) 和图 1.11(c)）。在该轨迹这部分 $v \approx v_C$，从 $(1.97)$ 得尾缘波前沿方程为

\[ u_t = D u_{xx} + f(u) - v_C. \]

该方程的行波前沿解必须有 $u = u(z), z = x - ct; u(-\infty) = u_D, u(\infty) = u_C$。Volume I §13.5 的分析再次给出其唯一波速为

\[ c = \sqrt{\frac{D}{2}} (u_C - 2 u_P + u_D). \]

由 $(1.98)$ 中 $f(u)$ 的表达式，$f(u) = v_C$ 的根 $u_C, u_D, u_P$ 由 $v_C$ 决定。$(1.101)$ 中波速 $c$ 就是 $c(v_C)$，是 $v_C$ 的函数。现在通过要求 $c(v_C)$ 等于之前算出的脉冲前缘波速 $c = \sqrt{D/2}(1 - 2a)$（由 $(1.99)$）来确定 $v_C$。原则上可以这样做，因为 $v_C$ 的表达式是一个多项式的解。

要完成波脉冲的解析确定，现在必须考虑解和相轨迹 $D 0$ 的部分（图 1.10, 1.11(c)）。对 $BC$ 部分，在该阶段 $v$ 在 $O(1/\varepsilon)$ 时间内又变化 $O(1)$。这称为现象的不应期。图 1.12 展示系统 $(1.94)$ 的一个计算例子，其中三次 $f(u)$ 用分段线性表达式 $f(u) = H(u - a) - u$ 近似。

阈值波对更一般的可激发介质模型也存在。为突出解析概念，考虑一个反应中一个是快的两物种系统。为方便分析，考虑反应扩散系统

\[ \varepsilon u_t = \varepsilon^2 D_1 u_{xx} + f(u, v), \quad v_t = \varepsilon^2 D_2 v_{xx} + g(u, v), \]

其中 $0 < \varepsilon \ll 1$，动力学 $f, g$ 的零斜线如图 1.13(a) 所示，并利用 $\varepsilon$ 小的优势。$f(u, v) = 0$ 的关键定性形状是三次曲线。这种形式在涉及激活和抑制的反应中很典型（参见 Volume I §6.6, §6.7）。

$(1.102)$ 在无扩散时可激发。Volume I §3.8 中阈值机制的描述相当模糊。本节给出一个更精确的机制定义：一种机制如果大小足够的刺激能引发在介质中传播的行波脉冲，则它是可激发的。

对 $0 < \varepsilon \ll 1$，$(1.102)$ 的 $O(1)$ 形式是 $f(u, v) = 0$，可解出 $u$ 作为 $v$ 的多值函数。从图 1.13(a) 看，对所有给定的 $v_m < v < v_M$，$f(u, v) = 0$ 对 $u$ 有三个解：它们是 $v = $ 常数线与零斜线 $f(u, v) = 0$ 的交点。与 FHN 系统讨论类似，我们有连接 $S$ 与 $A$ 的波前型解和从 $B$ 到 $C$ 的尾随波前，中间为慢过渡。$u$ 从 $S$ 处的值变到 $A$ 处的值是快的，时间为 $O(\varepsilon)$。从 $A$ 经过 $B$ 到 $C$ 则需相对长时间 $O(1)$，而 $B C$ 又在 $O(\varepsilon)$ 内覆盖。

脉冲解的解析研究相当复杂，详细分析已由 Keener (1980) 给出。这里用前缘分析作为示例说明：即考虑图 1.13(a), 1.13(b) 中从 $S$ 到 $A$ 的过渡。过渡发生得快，时间 $O(\varepsilon)$，空间上是厚度 $O(\varepsilon)$ 的陡前沿（图 1.13(b)）。这些尺度由对 $(1.102)$ 的奇异摄动估计指出。这提示我们通过变换引入新的自变量

\[ \tau = t/\varepsilon, \quad \xi = (x - x_T)/(\varepsilon \sqrt{D_1}), \]

其中 $x_T$ 是过渡前缘的位置（这里不需要）。引入 $x_T$ 只是为了使波脉冲的前缘在 $\xi$-平面的原点 $\xi = 0$。代入 $(1.102)$，令 $\varepsilon \to 0$，$\tau, \xi$ 固定（按通常奇异摄动方法，见 Murray 1984 或 Kevorkian 与 Cole 1996）得 $O(1)$ 系统

\[ u_\tau = u_{\xi\xi} + f(u, v), \quad v_\tau = 0. \]

考虑线 $S A$，第二个方程就是 $v = v_S$，于是需要解

\[ u_\tau = u_{\xi\xi} + f(u, v_S). \]

这正是 $u$ 的标量方程，其中 $f$ 有三个稳态 $u_S, u_D, u_A$，定性与上一章研究过的相同。本质上与方程 $(13.62)$ 相同，在 Volume I §13.5 详细讨论过。在无扩散时，$S, D, A$ 处的稳态分别为稳定、不稳定、稳定。Volume I §13.5 已展示如何存在以唯一波速传播的单调行波解，它能连接 $u_S$ 和 $u_A$（图 1.13(b)）。

完整解需要确定波速和脉冲的其它部分，即 $AB, BC, CS$，确保它们都连续地连接起来。这是一个有趣的奇异摄动分析。Keener (1980) 做了此分析，也给出了数值解和二维阈值波的分析。

另一种实际感兴趣的阈值波发生在 $f(u, v) = 0, g(u, v) = 0$ 的零斜线如图 1.14(a) 相交时。此时有三个稳态。用 $(1.102)$ 的标度，我们有从 $A$ 到 $D$ 的陡前沿和较慢的 $DC$ 过渡，本质上同上。我们可用 Volume I §13.5 的方法得到波速符号。此时波有尾巴，因为 $C$ 是线性稳定稳态。图 1.14(b) 是这种阈值前沿波的典型例子。也存在从一致稳态 $C$ 扰动从而实现向 $A$ 的过渡波。Rinzel 与 Terman (1982) 在 FHN 背景下研究过这些波。

阈值波存在于各种各样的真实世界系统中——任何能展示阈值动力学的系统。例如 Britton 与 Murray (1979) 在一类底物抑制振荡器中研究过（参见 Britton 1986 的书）。Showalter 课题组的一篇文章 (Steinbock et al. 1996) 描述了可激发系统中化学波逻辑门的新概念（与该课题组大多数文章一样，思想和分析都有原创和启发性的实验支持）。

1.7 具有振荡动力学的反应扩散系统中的行波波列（Travelling Wave Trains in R-D Systems with Oscillatory Kinetics）

一般反应扩散系统带极限环动力学的波列解已被广泛研究；Kopell 与 Howard (1973) 和 Howard 与 Kopell (1977) 的数学论文是开创性的。Field 与 Burger (1985) 编辑的书中几篇综述文章与本节相关；其它参考文献将在适当处给出。

所关心的一般演化系统是 $(1.1)$，对当前目的我们限制在一个空间维度并为代数简单将扩散系数结合到新的标度空间变量 $x \to x/\sqrt{D_1}$ 中。方程系统为

\[ \partial_t u = f(u) + \partial_{xx} u. \]

假设空间齐次系统 $du/dt = f(u; \gamma)$（$\gamma$ 是分岔参数）有 $\gamma < \gamma_c$ 时的稳定稳态，并通过 Hopf 分岔（参见 Strogatz 1994），在 $\gamma > \gamma_c$ 时演化为稳定极限环解；即对 $\gamma = \gamma_c + \varepsilon$（$0 < \varepsilon \ll 1$），存在稳定的小振幅极限环解。

行平面波列解形式为 $u(x, t) = U(z), z = \sigma t - k x$，其中 $U$ 是 $z$（"相位"）的 $2\pi$-周期函数。$\sigma > 0$ 是频率，$k$ 是波数；波长 $w = 2\pi/k$。波以速度 $c = \sigma/k$ 传播。这种形式只是 Volume I 第 13 章和上面使用的一般行波形的轻微变体，可通过重新标度时间简化。代入 $(1.108)$ 到 $(1.106)$ 得 $U$ 的常微方程

\[ k^2 U'' - \sigma U' + f(U) = 0, \]

其中撇号表示对 $z$ 求导。我们要找 $\sigma, k$ 使最后一个方程对 $U$ 有 $2\pi$-周期解。

不直接考虑一般情形（参见 Kopell 与 Howard 1973 及下文评论），作为示例对 $\lambda$–$\omega$ 模型系统（Volume I §7.4，方程 $(7.30)$）进行代数上更简单的分析是启发性的。后面将把它与生物学实际情形下出现的一般反应扩散系统联系起来。这两个反应物（$u, v$）的模型机制为

\[ \partial_t \begin{pmatrix} u \\ v \end{pmatrix} = \begin{pmatrix} \lambda(r) & -\omega(r) \\ \omega(r) & \lambda(r) \end{pmatrix} \begin{pmatrix} u \\ v \end{pmatrix} + \partial_{xx} \begin{pmatrix} u \\ v \end{pmatrix}, \quad r^2 = u^2 + v^2. \]

$\omega(r), \lambda(r)$ 是 $r$ 的实函数。如果对某 $r_0 > 0$, $r_0$ 是 $\lambda(r)$ 的孤立零点，$\lambda'(r_0) < 0, \omega(r_0) \ne 0$，则空间齐次系统（即 $\partial_{xx} = 0$）有极限环解（见 Volume I §7.4，$(1.113)$）。将变量从 $(u, v)$ 变到极坐标 $(r, \theta)$（$\theta$ 是相位），定义 $u = r \cos \theta, v = r \sin \theta$，$(1.110)$ 变为

\[ r_t = r \lambda(r) + r_{xx} - r \theta_x^2, \quad \theta_t = \omega(r) + \frac{1}{r^2} (r^2 \theta_x)_x. \]

如果存在 $r_0 > 0$ 且 $\lambda'(r_0) < 0$，动力学的不稳定极限环解立即由 $r = r_0, \theta = \theta_0 + \omega(r_0) t$ 给出，$\theta_0$ 是任意相位。代入 $(1.111)$ 得极限环解

\[ u = r_0 \cos[\omega(r_0) t + \theta_0], \quad v = r_0 \sin[\omega(r_0) t + \theta_0], \]

频率为 $\omega(r_0)$，振幅为 $r_0$。

在极坐标形式下寻找 $(1.108)$ 类型的行平面波解，$r = \alpha, \theta = \sigma t - k x$。代入 $(1.112)$ 得这些是行波解的必要充分条件

\[ \sigma = \omega(\alpha), \quad k^2 = \lambda(\alpha). \]

所以，以 $\alpha$ 为便利参数，$(1.110)$ 有一个单参数族的行波波列解

\[ u = \alpha \cos[\omega(\alpha) t - x \lambda^{1/2}(\alpha)], \quad v = \alpha \sin[\omega(\alpha) t - x \lambda^{1/2}(\alpha)]. \]

波速为 $c = \sigma/k = \omega(\alpha)/\lambda^{1/2}(\alpha)$。如果 $r = \alpha \to r_0$（即 $\lambda$–$\omega$ 动力学有极限环解），平面波的波数趋于零。这提示我们应在极限环附近寻找行平面波列解。Kopell 与 Howard (1973) 展示了一般如何做。这里考虑一个具体简单但非平凡的例子：$\lambda(r), \omega(r)$ 满足 $(1.107)$ 的 Hopf 条件，可以简单地进行解析推导行波波列解；主要沿用 Ermentrout (1981) 的分析。

设 $\omega(r) \equiv 1, \lambda(r) = \gamma - r^2$。$(1.110)$ 中的动力学有 $u = v = 0$ 作为稳态，$\gamma < 0$ 时稳定，$\gamma > 0$ 时不稳定。$\gamma = 0$ 是分岔值 $\gamma_c$，在 $\gamma = 0$ 时 $u = v = 0$ 线性化的特征值为 $\pm i$。这是标准 Hopf 分岔要求（见 Strogatz 1994），所以预期 $\gamma = \gamma_c + \varepsilon$（$0 < \varepsilon \ll 1$）时小正 $\gamma$ 给出小振幅极限环解。$\lambda = 0$ 时 $r = \sqrt{\gamma}$，由上述一般解 $(1.117)$，极限环解为 $u_\gamma(t) = \sqrt{\gamma} \cos t, v_\gamma(t) = \sqrt{\gamma} \sin t$（$\gamma > 0$）；极坐标下 $r_0 = \sqrt{\gamma}, \theta = t + \theta_0$，可取 $\theta_0 = 0$。

现在考虑带 $(1.119)$ 中 $\lambda, \omega$ 的反应扩散系统 $(1.112)$。代入 $r = r_0, \theta = \sigma t - k x$ 形式的行平面波解，按 $(1.116)$ 得 $\sigma = 1, k^2 = \gamma - r_0^2$（$0 < r_0 < \sqrt{\gamma}$），从而小振幅行波波列解为

\[ u = r_0 \cos(t - x \sqrt{\gamma - r_0^2}), \quad v = r_0 \sin(t - x \sqrt{\gamma - r_0^2}). \]

图 1.15 示意这些解，振幅 $r_0 < \sqrt{\gamma}$，波长 $L = 2\pi/\sqrt{\gamma - r_0^2}$。

只有当这些行波波列稳定时，它们才与图 1.1 中的目标图案等相关。这个特定系统的线性稳定性可以做，但一般来说远非平凡。这是一个罕见的例子可以相当容易地进行分析。

扩散对具有周期行为的反应动力学的效应是产生行波周期波列解。上述 $\lambda$–$\omega$ 例中的具体非线性 $\lambda(r) = \gamma - r^2$ 是 Hopf 分岔问题的典型。似乎可能的是，反应扩散机制——其中反应动力学本身通过 Hopf 分岔呈现周期极限环行为——也会产生周期波列解。为证明这一点，只须证明具有这种性质的一般反应扩散系统在 Hopf 分岔附近类似于 $\lambda$–$\omega$ 系统。

考虑两物种系统 $u_t = F(u, v; \gamma) + D \nabla^2 u, v_t = G(u, v; \gamma) + D \nabla^2 v$。为代数简单，设 $(1.123)$ 有 $u = v = 0$ 处的稳态，无扩散（$D = 0$）系统在分岔值 $\gamma_c$ 处呈现到极限环的 Hopf 分岔。把 $u, v$ 看作零稳态的扰动，写 $T = (u, v)^T$，$M = (F_u, F_v; G_u, G_v)|_{u=v=0}$，$P = \text{diag}(D, D)$，则 $(1.123)$ 的线性化形式为 $T_t = M T + P \nabla^2 T$，完整系统 $(1.123)$ 可写成 $T_t = M T + P \nabla^2 T + H$，其中 $H$ 包含 $u, v$ 的非线性贡献。

由于动力学在 $\gamma = \gamma_c$ 处经历 Hopf 分岔，$M$ 的特征值 $\sigma$ 满足 $\text{Re } \sigma(\gamma) < 0$ 对 $\gamma < \gamma_c$，$\text{Re } \sigma(\gamma_c) = 0$，$\text{Im } \sigma(\gamma_c) \ne 0$，$\text{Re } \sigma(\gamma) > 0$ 对 $\gamma > \gamma_c$。所以在 $\gamma = \gamma_c$，$\text{Tr } M = 0, \det M > 0 \Rightarrow \sigma(\gamma_c) = \pm i (\det M)^{1/2}$。

引入常数矩阵 $N$ 和非常数矩阵 $W$，设 $T = N W$，从 $(1.125)$ 得 $W_t = N^{-1} M N W + N^{-1} P N \nabla^2 W + N^{-1} H$。选 $N$ 使

\[ N^{-1} M N = \begin{pmatrix} 0 & -k \\ k & 0 \end{pmatrix} \text{ 在 } \gamma = \gamma_c, \Rightarrow k^2 = \det M |_{\gamma = \gamma_c}. \]

在变换后的系统 $(1.127)$ 中，线性化矩阵的系数为 $N^{-1} M N = \begin{pmatrix} \alpha(\gamma) & -\beta(\gamma) \\ \beta(\gamma) & \delta(\gamma) \end{pmatrix}$，其中 $\alpha(\gamma_c) = 0 = \delta(\gamma_c), \beta(\gamma_c) \ne 0$。也就是说，具有 $\gamma_c$ 处稳态 Hopf 分岔的一般系统 $(1.123)$ 在分岔 $\gamma_c$ 附近可以变换为 $\lambda$–$\omega$ 形式（cf. $(1.110)$）。这个结果有一定的重要性，因为对 $\lambda$–$\omega$ 系统有效的分析在许多情形下可推广到实际反应扩散系统。$\lambda$–$\omega$ 系统对 $u, v$ 的扩散系数不同的情况已由 Duffy et al. (1980) 讨论，他们研究了其对螺旋波的影响，这是下一节的主题。

可以预期，如果我们有与振荡动力学相关的行波波列，那么具有周期倍增和混沌的动力学会呈现更复杂的波现象。Merkin et al. (1996) 简要描述的混沌类型不属于这种周期倍增类型。

1.8 螺旋波（Spiral Waves）

旋转螺旋波自然地出现在各种各样的生物、生理和化学背景中。被广泛研究的一类是在 BZ 反应中产生的那些。相对而言，它比生理学中出现的系统简单得多——生理学中我们不知道涉及的详细机制，而 BZ 机制则是已知的。许多人在这种反应中做过螺旋解的实验工作，例如 Winfree (1974)——他是早期研究和后续应用于心脏问题的主要人物，Krinskii et al. (1986)，Müller et al. (1985, 1986, 1987)。后者的新颖实验技术使用光吸收，定量地突出实际浓度水平。图 1.16 以及图 1.19、图 1.20 展示了 BZ 反应中实验观察到的螺旋波；另见图 1.1(b)。虽然这些图中的螺旋是对称的，但这绝不是唯一的图案形式；Winfree (1974) 特别地，以及 Müller et al. (1986) 展示了复杂螺旋图案的引人注目的例子。大量工作投入到螺旋波的数学研究中，特别是 FKN 模型系统的扩散版本。Keener 与 Tyson (1986) 分析了带一般可激发动力学的可激发反应扩散系统中的螺旋波。他们将技术应用于带扩散的 FKN 模型，结果与实验很好吻合。在不同的背景下，见图 1.18 中非对称以及对称螺旋的其它例子。Keener 与 Sneyd (1998) 的书中给出了大量螺旋波例子。Keener (1986) 给出了螺旋波的几何理论一般讨论，Zykov (1988) 在他关于可激发介质中波动过程的书中，Grindrod (1996) 的书也讨论了。

Showalter 和他的同事们对化学螺旋波做了许多新颖和有开创性的工作。例如 Amemiya et al. (1996) 用 BZ 反应的 Field–Körös–Noyes (FKN) 模型系统（类似于 Volume I 第 8 章中研究的）展示可激发动力学，研究三维螺旋波并进行相关实验支持他们的分析；见该处其它参考文献。

螺旋波还有许多其它重要出现。大脑组织可以展示"扩散性抑郁"的电化学波，它扩散通过大脑皮层。这些波的特征是神经元膜去极化和神经活动减少。Shibata 与 Bureš (1972, 1974) 实验研究了这一现象，展示了螺旋波的存在，它们在大鼠皮层脑组织中绕病灶旋转。图 1.17(a) 示意他们观察到的波行为。Keener 与 Sneyd (1998) 一般性地讨论波运动，特别是 Hodgkin–Huxley 方程及其漫画系统 FitzHugh–Nagumo 方程中发现的波运动类型。他们也描述心脏节律性和波传播以及钙波；讨论的一些波现象与本书中涉及的有相当不同，比如波曲率效应。Grindrod et al. (1991) 给出了基于 eikonal 方法的通用程序，可包含曲率效应，特别是曲面上的曲率效应。

正在不断发现反应扩散模型和螺旋波的新现象和新的应用。Winfree 的论文和他文中的参考文献是一个好的出发点。例如 Winfree et al. (1996) 把行波与心肌和神经中"可激发"的运动方面联系起来。他们得到了复杂的周期行波，类似于从涡环辐射的卷轴，涡环是组织中心。Winfree (1994b) 用一般可激发反应扩散系统，表明这些涡线的一般构型是一个混沌缠结。他用的系统为

\[ \partial_t u = \nabla^2 u + \left[\frac{u - u^3}{3} - v\right]/\varepsilon, \quad \partial_t v = \nabla^2 v + \varepsilon \left[u + \beta - \frac{v}{2}\right]. \]

当心肌纤维的协调收缩中断导致死亡时，原因通常是纤颤。在纤颤心脏中，小区域基本上独立收缩。心脏看起来像一把蠕动的虫——是一团颤动的组织。如果这种中断持续几分钟以上，通常会导致死亡。Krinskii (1978) 和 Krinskii et al. (1986) 例如讨论了心律失常数学模型中的螺旋波。Winfree (1983a, b) 考虑了应用于突发性心脏死亡的可能性。他提出纤颤的前兆是电脉冲旋转波的出现。图 1.17(b) 示意 Allessie et al. (1977) 在兔心组织中诱导的这种波。这些作者 (Allessie et al. 1973, 1976, Smeets et al. 1986) 也对心肌中旋转波传播进行了广泛的实验研究。

Winfree (1994a, 1995) 提出了一个有趣的假设：突发性心脏死亡可能涉及电活动的三维旋转子（螺旋型波），当心脏厚度超过某个临界值时它突然变得不稳定；见该处对这一领域其它文章的参考引用。他对过去 20 年肌肉组织中的复杂波现象进行了广泛研究。例如，电方面、激活前缘、各向异性等在心脏生理中很重要，Winfree (1997) 在反应扩散背景下详细讨论了它们。他也讨论了心脏中电位扩散、电湍流、激活前缘曲率和各向异性的作用，以及它们在心脏衰竭中的可能作用。他的工作提出了对严重心脏衰竭（如心室纤颤）补救治疗的可能方案。通过他一系列关于心脏肌肉组织活动建模的论文，Winfree 极大地增强了我们对突发性心脏死亡的理解，并以重大方式改变了之前持有的（医学）观念。

黏菌 Dictyostelium discoideum 中出现的信号图案螺旋同样引人注目（图 1.18）。Tyson et al. (1989a, b) 提出了基于实验动机动力学方案的模型。重要的是（一般来说也是），虽然图 1.1(b) 和图 1.18 之间的相似性很显著，但不应假设 BZ 反应的模型就是黏菌图案的合适模型——机制非常不同。虽然产生正确类型的图案是成功建模的重要和必要方面，但理解基本机制是最终目标。

相互作用的种群情形中大尺度螺旋的可能存在似乎尚未从实际应用的角度被考虑，但鉴于模型的反应扩散特征，它们在理论上的确存在。

从数学观点看，螺旋波是什么？在 BZ 反应的情形，例如它是反应物浓度的旋转、周期、空间结构；见图 1.17, 1.20。固定时刻的快照显示典型螺旋图案。过程的电影显示整个螺旋图案像旋转的钟表弹簧那样运动。图 1.19 显示这样的快照和按固定时间间隔取图案的叠加。陡波前沿是恒定浓度的等值线，即等浓度线。

现在考虑一个绕其中心旋转的螺旋波。如果你站在介质中的固定位置，似乎局部有周期波列通过，因为每次螺旋旋转一次，就有一个波前移过你。

Volume I 第 9 章中，反应物状态或浓度可以描述为相位的函数。讨论螺旋波时显然适合用极坐标 $r, \theta$。简单的旋转螺旋由相位的周期函数描述，$\phi = \omega t \pm m \theta + \psi(r)$，其中 $\omega$ 是频率，$m$ 是螺旋臂数，$\psi(r)$ 是描述螺旋类型的函数。$m \theta$ 项中的 $\pm$ 决定旋转方向。图 1.20 展示 1 臂和 3 臂螺旋的例子，包括后者的实验例子。例如设 $\phi = 0$ 考察稳态情形；从 $(1.131)$ 得到一个简单的螺旋几何描述：例如 1 臂螺旋由 $\theta = \psi(r)$ 给出。具体 $\psi(r)$ 为 $\theta = a r$ 或 $\theta = a \ln r$（$a > 0$）；分别是阿基米德螺线和对数螺线。对中心核的螺旋，对应形式为 $\theta = a(r - r_0)$ 或 $\theta = a \ln(r - r_0)$。

反应物 $u$ 的螺旋构型的数学描述可表达为 $u(r, \theta, t) = F(\phi)$，其中 $F(\phi)$ 是相位 $\phi$（由 $(1.131)$ 给出）的 $2\pi$-周期函数。固定 $t$ 得到螺旋快照，其形式对 $\psi(r)$ 加一定约束；$a r$ 和 $a \ln r$ 在 $(1.132)$ 中只是两个简单情形。混合类型例如有 $\psi(r) = a r + b \ln r$（$a, b$ 常数）。$(1.134)$ 中 $\phi$ 如 $(1.131)$，固定 $r, t$ 绕中心一圈有 $m$-重对称性，$m$ 是臂数；$m = 3$ 的例子见图 1.20(b) 和 Agladze 与 Krinskii (1982) 实验得到的图 1.20(d)。固定 $r, \theta$（即停在固定点），则看到一系列波前依次通过。如果波前在 $t = t_0$（$\phi = \phi_0$）通过，下一个波在 $t = t_0 + 2\pi/\omega$ 通过（此时 $\phi = \phi_0 + 2\pi$）。

如果看螺旋快照并沿射线从中心向外移动，直观上存在一个与螺旋相关的波长；它随着从中心向外移动而变化。如果一波前在 $r_1$，下一波前（向外）在 $r_2$，可定义波长 $\lambda$ 为 $\lambda = r_2 - r_1$，$\theta(r_2) = \theta(r_1) + 2\pi$。从 $(1.131)$，固定 $t$，沿曲线 $\phi = $ 常数，$\phi_\theta + \phi_r (dr/d\theta)|_{\phi = \text{const}} = 0$，所以如果具体取 $(1.131)$ 中 $-m \theta$，则

\[ \left. \frac{dr}{d\theta} \right|_{\phi = \text{const}} = -\frac{\phi_\theta}{\phi_r} = \frac{m}{\psi'(r)}. \]

波长 $\lambda(r)$ 现由

\[ \lambda(r) = \int_{\theta(r)}^{\theta(r) + 2\pi} \left. \frac{dr}{d\theta} \right|_{\phi = \text{const}} d\theta = \int_{\theta(r)}^{\theta(r) + 2\pi} \frac{m}{\psi'(r(\theta))} d\theta \]

给出，$r$ 作为 $\theta$ 的函数由 $(1.131)$ 在 $t$ 常数、$\phi$ 常数（取为零）时给出。对阿基米德螺线 $r = \theta/a$，$\psi' = a$，波长 $\lambda = m/a$。

螺旋的螺距定义为 $\left. dr/d\theta \right|_{\phi = \text{const}} = m/\psi'(r)$。对阿基米德螺线 $\psi'(r) = a$ 给出常数螺距 $m/a$；对数螺线给出螺距 $m r/a$ 因为 $\psi'(r) = a/r$。对大 $r$，后者的螺距大，即疏松缠绕；对小 $r$，螺距小；即紧密缠绕。

在讨论具体反应扩散系统的解析解之前，应该注意 Krinskii et al. (1986) 和 Tsujikawa et al. (1989) 进行的关于螺旋波产生的数值研究。后者考虑了 FitzHugh–Nagumo 可激发机制 $(1.94)$ 并数值研究了有限空间范围的激发波传播。图 1.21 展示激发行波的时间序列，展示了螺旋波的演化；Krinskii et al. (1986) 得到了类似的演化图，他们详细讨论了螺旋波的演化。该图中的演化图样与 BZ 反应中实验观察到的正在发展的螺旋相似。

1.9 λ–ω 反应扩散系统的螺旋波解（Spiral Wave Solutions of λ–ω Reaction Diffusion Systems）

许多作者研究了一般反应扩散模型的螺旋波解，例如 Cohen et al. (1978)，Duffy et al. (1980)，Kopell 与 Howard (1981)，Mikhailov 与 Krinskii (1983)。Keener 与 Tyson (1986) 关于 BZ 反应的论文和 Tyson et al. (1989a, b) 关于黏菌 Dictyostelium 的论文是具体例子。分析通常很复杂，多用渐近方法。$\lambda$–$\omega$ 系统如 §1.7 所见展示波列解，被用作模型系统，因为其分析相对代数简单。$\lambda$–$\omega$ 系统的螺旋波解已由 Greenberg (1981)、Hagan (1982)、Kuramoto 与 Koga (1981)、Koga (1982) 等人研究。参考文献相当广泛；其它相关参考文献见这些论文和 Keener 与 Sneyd (1998) 的书。本节为 $\lambda$–$\omega$ 系统发展一些解。

两反应物的 $\lambda$–$\omega$ 反应扩散机制为

\[ \partial_t \begin{pmatrix} u \\ v \end{pmatrix} = \begin{pmatrix} \lambda(A) & -\omega(A) \\ \omega(A) & \lambda(A) \end{pmatrix} \begin{pmatrix} u \\ v \end{pmatrix} + D \nabla^2 \begin{pmatrix} u \\ v \end{pmatrix}, \quad A^2 = u^2 + v^2, \]

$\omega(A), \lambda(A)$ 是 $A$ 的实函数。（从 $(1.110)$ 改变记号以便用 $r$ 作通常极坐标。）假设动力学支撑极限环振荡；这对 $\lambda, \omega$ 加上通常约束：$A_0$ 是 $\lambda(A)$ 对某 $A_0 > 0$ 的孤立零点，$\lambda'(A_0) < 0$，$\omega(A_0) \ne 0$，则空间齐次系统（$D = 0$）有稳定极限环解 $u^2 + v^2 = A_0$，频率 $\omega(A_0)$（见 Volume I §7.4）。

设 $w = u + i v$，$(1.135)$ 变为单复数方程 $w_t = (\lambda + i \omega) w + D \nabla^2 w$。该方程的形式提示设 $w = A \exp[i \phi]$，其中 $A$ 是 $w$ 的振幅，$\phi$ 是其相位。代入 $(1.136)$，实虚部分离得 $A, \phi$ 的方程组

\[ A_t = A \lambda(A) - D A |\nabla \phi|^2 + D \nabla^2 A, \]

\[ \phi_t = \omega(A) + 2 A^{-1} D (\nabla A \cdot \nabla \phi) + D \nabla^2 \phi, \]

这是 $\lambda$–$\omega$ 系统的极坐标形式。极坐标 $r, \theta$ 是讨论螺旋波时的适当坐标。受 $(1.131)$ 和上节讨论的启发，我们寻找形式为

\[ A = A(r), \quad \phi = \omega t + m \theta + \psi(r) \]

的解，其中 $\omega$ 是未知频率，$m$ 是螺旋臂数。代入 $(1.138)$ 得 $A, \psi$ 的常微方程：

\[ D A'' + \frac{D}{r} A' + A \left[ \lambda(A) - D \psi'^2 - D \frac{m^2}{r^2} \right] = 0, \]

\[ D \psi'' + D \left[ \frac{1}{r} + 2 \frac{A'}{A} \right] \psi' = \omega - \omega(A). \]

第二个方程乘以 $r A^2$，积分得

\[ \psi'(r) = \frac{1}{D r A^2(r)} \int_0^r s A^2(s) [\omega - \omega(A(s))] ds. \]

$(1.140), (1.141)$ 的形式对分析很方便，是许多关于 $\lambda$–$\omega$ 系统螺旋波论文的基础，使用渐近方法、不动点定理、相空间分析等。

在分析 $(1.140)$ 之前，必须决定合适的边界条件。我们要求解在原点正则、$r \to \infty$ 时有界。前者结合 $A, \psi'$ 的方程形式，要求 $A(0) = 0, \psi'(0) = 0$。

如果 $A \to A_\infty$ 当 $r \to \infty$，从 $(1.141)$ 得

\[ \psi'(r) \sim \frac{1}{D r A_\infty^2} \int_0^r s A_\infty^2 [\omega - \omega(A_\infty)] ds = \frac{[\omega - \omega(A_\infty)] r}{2D} \text{ 当 } r \to \infty, \]

所以 $\psi'$ 有界当且仅当 $\omega = \omega(A_\infty)$。$(1.140)$ 的第一个方程确定 $\psi'(\infty) = [\lambda(A_\infty)/D]^{1/2}$。我们于是有弥散关系

\[ \psi'(\infty) = \left[\frac{\lambda(A_\infty)}{D}\right]^{1/2}, \quad \omega = \omega(A_\infty), \]

它展示无穷远处的振幅如何决定频率 $\omega$。

在 $r = 0$ 附近，设 $A(r) \sim r^c \sum_{n=0}^\infty a_n r^n$ 当 $r \to 0$（$a_0 \ne 0$），代入 $(1.140)$ 第一个方程并按 $r$ 的幂次匹配。最低阶 $r^{c-2}$ 的系数设为零给出 $c(c-1) + c - m^2 = 0 \Rightarrow c = \pm m$。为使 $A(r)$ 在 $r \to 0$ 时非奇异，必须选 $c = m$，于是 $A(r) \sim a_0 r^m$ 当 $r \to 0$，$a_0$ 是待定非零常数。数学问题变为确定 $a_0, \omega$ 使 $A(r), \psi'(r)$ 在 $r \to \infty$ 时有界。从 $(1.137), (1.139)$ 和最后方程，$r = 0$ 附近 $u, v$ 的行为为

\[ \begin{pmatrix} u \\ v \end{pmatrix} \propto \begin{pmatrix} r^m \cos[\omega t + m \theta + \psi(0)] \\ r^m \sin[\omega t + m \theta + \psi(0)] \end{pmatrix}. \]

Koga (1982) 对 $\lambda$–$\omega$ 系统的相位奇异和多臂螺旋做了分析的和数值的研究，其中 $\lambda(A) = 1 - A^2, \omega(A) = -\beta A^2$（$\beta > 0$）。图 1.22 展示了他对 1 臂和 2 臂螺旋的计算结果。

寻找解的基本出发点是 $u, v$ 的函数形式假设

\[ \begin{pmatrix} u \\ v \end{pmatrix} = \begin{pmatrix} A(r) \cos[\omega t + m \theta + \psi(r)] \\ A(r) \sin[\omega t + m \theta + \psi(r)] \end{pmatrix}. \]

$A(r)$ 为常数、$\psi(r) \propto \ln r$ 时，这些表示旋转螺旋波。Cohen et al. (1978) 证明对一类 $\lambda(A), \omega(A)$，系统 $(1.136)$ 有形如 $(1.146)$ 的旋转螺旋波，满足渐近到阿基米德和对数螺线的边界条件，即 $\psi \sim c r$ 和 $\psi \sim c \ln r$ 当 $r \to \infty$。Duffy et al. (1980) 展示了如何将 $u, v$ 扩散系数不等的、带极限环动力学的一般反应扩散系统化为 Cohen et al. (1978) 分析的情形。

Kuramoto 与 Koga (1981) 数值研究了具体 $\lambda$–$\omega$ 系统 $\lambda(A) = \varepsilon - a A^2, \omega(A) = c - b A^2$（$\varepsilon > 0, a > 0$）。代入这些 $(1.136)$ 变为 $w_t = (\varepsilon + i c) w - (a + i b) |w|^2 w + D \nabla^2 w$。通过设 $w \to w e^{i c t}$（代数上等价于 $c = 0$）可消去 $c$ 项，然后按

\[ w \to \sqrt{a/\varepsilon} \, w, \quad t \to \varepsilon t, \quad r \to \sqrt{\varepsilon/D} \, r \]

重新标度 $w, t$ 和空间坐标，得更简单的形式

\[ w_t = w - (1 + i \beta) |w|^2 w + \nabla^2 w, \]

其中 $\beta = b/a$。该方程的空间无关形式有极限环解 $w = \exp(-i \beta t)$。

Kuramoto 与 Koga (1981) 建议 Volume I 第 9 章讨论的"无相位"点或黑洞开始出现，导致混沌不稳定。比较 $(1.147)$ 与 $(1.136)$，有 $\lambda = 1 - A^2, \omega = -\beta A^2$，所以 $\beta$ 度量局部极限环频率对振幅 $A$ 的依赖强度。由于 $A$ 随空间坐标 $r$ 变化，我们有一种耦合的、适当同步的振子阵列情形。随着 $|\beta|$ 增大，振子的变化增大。稳定旋转波需要一定的同步性，增大局部"振子"的变化会破坏同步性，导致无相位点和混沌。化学反应扩散机制中波前的混沌或湍流已由 Kuramoto (1980) 详细考虑；见该处对空间混沌这一有趣问题的其它参考文献。

作为本节结尾，让我们看螺旋波的 1 维类比：一个脉冲从核心（位于原点）周期性地、交替地从核心两侧发出。如果脉冲对称发出，那么就是目标图案的类比。考虑 $(1.138)$ 中 $\nabla^2 = \partial_{xx}$，$\lambda(A) = 1 - A^2$，$\omega(A) = q A^2$。设 $x \to x/\sqrt{D}$，$A = A(x)$，$\phi = \omega t + \psi(x)$，得到 $A, \psi$ 的方程

\[ A_{xx} + A (1 - A^2 - \psi_x^2) = 0, \quad \psi_{xx} + 2 \frac{A_x \psi_x}{A} = \omega - q A^2. \]

边界条件为 $A(x) \sim a_0 x$ 当 $x \to 0$，$\psi_x(0) = 0$，$A(x), \psi_x(x)$ 在 $x \to \infty$ 时有界。问题归结为作为 $q$ 函数求 $a_0$ 和 $\omega$ 使得时间相关方程初值问题的解有界。一个这样的解为

\[ A(x) = \left[\frac{\omega}{q}\right]^{1/2} \tanh(x/\sqrt{2}), \]

\[ \psi_x(x) = \left[1 - \frac{\omega}{q}\right]^{1/2} \tanh(x/\sqrt{2}), \quad \omega^2 + \frac{9}{2 q} \omega - \frac{9}{2} = 0, \]

可验证。这些解在原点处周期性地、交替地从两侧产生。

行波特别是螺旋波的稳定性通常很难解析地证明；Feroe (1982) 关于可激发 FHN 波稳定性的论文充分说明了这一点。然而对 $\lambda$–$\omega$ 系统波列解的情形，不经过冗长复杂的分析就能得到一些稳定性结果。螺旋波稳定性的解析确定一般仍远未完成，尽管数值证据表明许多确实稳定。如前所述，最近 Yagisita et al. (1998) 调查了可激发反应扩散系统中球面上的螺旋波。他们证明螺旋尖旋转，并研究了在齐次和非齐次介质中的传播。

存在作为反应扩散系统以外模型机制解的生物波。例如，后面第 6 章将讨论的若干机械化学（mechanochemical）图案和形态生成模型也支撑行波解。通过后留下空间图案的波也非常重要，我们也会在后几章看到。作为本章的总结，也许应该重申波现象在生物学中的重要性。虽然这从本章材料就已清晰，但它们在胚胎发育过程中的组织通信中可能更加重要。稳态空间图案和形态的生成是同等重要的话题，将在后续章节详细讨论。

本章个人批注

读完第 1 章最大的感受是：Murray 把"多物种波"这个相对抽象的话题，通过 BZ 反应、灰松-红松竞争和转基因生物扩散等具体例子，向生态学、流行病学和应用生物学各方向大幅延伸。第 1.1 节的"直觉预期"其实承担了全章的脚手架功能：pacemaker、target pattern、spiral wave、threshold wave、wave train——每一种都在后面某节有具体模型支撑。这种"先图景后数学"的写法对生物背景的读者非常友好，对数学背景的读者则要小心：直觉经常会跑得比严格结果更远。

第 1.2 节的捕食者-被捕食追逐-逃避波特别值得反复读。Dunbar 的模型得到 (1.10) 中 $c \ge [4a(1-b)]^{1/2}$ 的下界，并以 $a^*$ 区分单调与振荡式趋近；这部分数学结论干净漂亮。但 Murray 也清楚交代：这种"追逐-逃避"在解里几乎看不出猎物的"逃避"行为。这是建模的诚实——模型结果不能反推回直觉得出。第 1.2 节后半的对流式追逐-逃避模型 (1.12)–(1.14) 是一个少见的、自然产生交叉扩散的生态学例子，文献里研究并不充分；Hasimoto 的 blow-up 分析、Schmidt-Ortoleva 的渐近展开都是值得追溯的延伸。

第 1.3 节是全章数据最实在的部分。Okubo et al. (1989) 的灰松模型通过仔细的 Leslie 矩阵分析给出 $a_1 = 0.82$/年、$1/b_1 = 10$/公顷等具体数值，再配合 Williamson-Brown 的 7.7 km/年观测速度，估出 $D_1 = 17.9$ km²/年——这个估计又恰好和"相邻林地平均 10 km"这个看似合理的数据一致。这种"模型—参数—数据"循环是这个领域真正能站住脚的标志。但我注意到：作者承认 $c_1, c_2$ 没有定量信息，$\gamma_1$ 的值是基于"应该小"的直觉取的。这是建模的典型妥协。

第 1.4 节把模型推广到周期斑块环境，是数学上最复杂的一节。Cruywagen et al. (1996) 的工作把问题化为 Hill 方程的稳定性分析，得到 $\gamma_e < 1, G_2 < \gamma_e$ 等多个临界条件；这些条件的最终结论是"$\gamma_e > 1$ 和 $\gamma_n < 1$ 是唯一对所有初始释放策略都安全的控制策略"。但作者也提醒：这是线性稳定性结果，初始条件仍可让局部入侵在某些斑块发生。这种"全局稳定但局部不稳定"的情形，作者特意指出"数学模型在这里特别有用"——因为直觉难以把握。

第 1.5 节的 BZ 反应前沿分析是经典工作。Murray (1976) 推导的 $c(r, b) \le 2$ 界，加上数值给出的 c 与 r、b 的关系，使作者能定量地与实验波速（O(2.7 × 10⁻¹ cm/min)）和前沿宽度（O(10⁻³ cm)）对照。这种"无量纲分析 + 量纲还原"的两步法，是教科书级的示范。我对 (1.93) 的"非均匀极限情形"那段讨论印象很深——$r \to \infty$ 时 $v = 0$ 和 $v \ne 0$ 的极限顺序不同，说明对这类反应扩散系统，参数极限不能随意交换。

第 1.6 节对 FHN 模型的阈值波分析是数学上最精巧的一节。Keener (1980) 的奇异摄动分析把脉冲分为 $0AB$（快前沿）、$BC$（慢过渡，$O(1/\varepsilon)$ 时间）、$CD$（后沿）、$D0$（不应期）四段，并要求前后沿波速一致来确定参数 $v_C$。这整个分析框架后来被 Rinzel-Terman (1982)、Feroe (1982) 等人继续推进，是理解神经脉冲和心肌可激发波的基础。我对"激发的精确定义"那段——"a mechanism is excitable if a stimulus of sufficient size can initiate a travelling pulse which will propagate through the medium"——印象很深：它把 Volume I §3.8 的模糊说法变成了可操作的标准。

第 1.7、1.8、1.9 节是螺旋波理论的核心。λ–ω 系统在极坐标下的特殊形式使周期波列解可以解析写出 (1.117)，再把 $A$ 设为常数、$\psi \propto \ln r$ 就得到旋转螺旋波。Duffy et al. (1980) 和 Cohen et al. (1978) 的一般性结果表明，靠近 Hopf 分岔的一般反应扩散系统都局部等价于 λ–ω 系统——这是个有相当广泛性的结论。Kuramoto-Koga 关于 $\beta$ 增大导致混沌的数值实验，以及他们关于"无相位点（phaseless points）"的物理解释，是把同步性破坏与空间混沌联系起来的有趣尝试。

最后提一个问题：本章在最后一节末尾说"biomechanical models in Chapter 6 also sustain travelling wave solutions"——这暗示第 6 章会重新出现波现象，是后续阅读的钩子。

与上下章的衔接（一段话）

第 1 章是全书的"现象学总览"——它承接 Volume I 第 13 章讨论过的单组分 Fisher–Kolmogoroff 行波前沿，把"波"从标量情形推广到多物种反应扩散和趋化机制，并通过八个具体案例（捕食者追逐-逃避、灰松-红松竞争、转基因微生物入侵、BZ 反应波、可激发介质脉冲、λ–ω 振荡波列、螺旋波）展示该推广带来的新现象与新分析技术。Murray 的安排是：先在 §1.1 用"直觉"勾勒目标图案、螺旋波、阈值波、波列等大类，再在 §1.2–§1.9 逐个给出具体模型、分析结果和实验/数值对照。这种"先鸟瞰后特写"的结构为后续章节定下基调：第 2 章转向"spatial pattern formation"——即稳态图案而非行波；第 3–7 章是各种生物图案（动物毛皮、生长域、细菌趋化、机械图样、形态进化约束）的案例；第 8–12 章是血管网络、伤口愈合、脑肿瘤和神经图案；第 13 章是流行病学扩散；第 14 章是领地性。第 1 章既是引子也是"波现象的工具箱"：后面所有需要行波解、稳定性分析、参数估计的章节，都会回到本章的模型和分析方法。