第 13 章：生物波：单物种模型（Biological Waves: Single-Species Models）

13.1 背景与行波波形（Background and the Travelling Waveform）

本章讨论生物体系中以化学浓度、机械形变、电信号等形式传播的行波现象。受精后早期胚胎发育中可观察到的众多波动事件，是这一现象最直观的写照。脊椎动物卵表面既有化学波也有机械波传播：例如鱼类 Medaka 卵的 Ca²⁺ 波就是从精子入卵点发出的表面横扫波。Belousov–Zhabotinskii 反应中可见的化学浓度波是另一类被广泛记录的视觉冲击强烈的例子。第 11 章关于昆虫扩散以及具有空间效应的种群相互作用模型中也预期了波现象。另一类重要的例子是流行病的传播波——在欧洲蔓延的狂犬病就是既生动又令人忧虑的案例，作者将在本章讨论其模型。趋化性导向的微生物向食物源迁移是又一典型例子，其中以盘基网柄菌 Dictyostelium discoideum 被研究得最为广泛。Winfree（2000）一书中囊括了大量生物波现象；Segel（1980）的分子与细胞生物学数学模型导论也涉及波动问题。Fife（1979）、Britton（1986）、Grindrod（1996）的反应扩散方程专著则从数学角度切入。Zeeman（1977）则从突变论角度考察发育等生物领域中的波现象。

作者强调，生物医学中波现象的广泛存在要求深入研究行波及其建模与分析。本章和第一卷第一章（以及第一卷中其他多章的众多例子）将讨论扩散扮演关键角色的各种波行为。本章研究的波与第 12 章讨论的波截然不同：相关数学文献已极为丰富，本章只能严格地限制议题数量与讨论深度。所涵盖的内容包括该领域已为学界接受的基础理论，以及两个实际问题——一个与昆虫扩散和控制相关，另一个与两栖类卵上的钙波相关。

在发育中的生命系统里，细胞间和细胞内水平几乎持续地进行信息交换。这种信息交流对发育过程及最终形态图式的依次生成是必需的。不同生化浓度的传播波是传递此类生化信息的方式之一。在发育胚胎中，生物化学物质的扩散系数可以非常小：10⁻⁹ 至 10⁻¹¹ cm²·sec⁻¹ 的量级相当常见。如此小的扩散系数意味着若扩散是主要机制，跨越数毫米的宏观距离需要极长的时间。例如对一行标准的一维扩散方程 ∂u/∂t = D∂²u/∂x²，可以用量纲分析、相似解或（11.10）式的经典解估得传播时间量级为 O(L²/D)。当 L 取 1 mm 量级时，由上述扩散系数得到的典型时间是 O(10⁷ 至 10⁹ sec)，对早期胚胎发育中大多数过程而言过于漫长。因此单纯扩散不太可能是信息在显著距离上传导的主要载体。一个可能的例外是蝴蝶翅图案的生成，它发生在蛹期并历时数日。

与单纯扩散形成对照，作者表明当反应动力学与扩散耦合时，可出现化学浓度行波，且其传输生化变化的速度远高于 (13.1) 这类纯扩散方程所支配的过程。这种耦合产生反应扩散方程，在一维标量情形下形如 ∂u/∂t = f(u) + D∂²u/∂x²，其中 u 是浓度，f(u) 代表动力学项，D 是扩散系数（此处取为常数）。

作者继而给出"行波"的精确定义：行波被取作"形状不变、匀速传播"的波，即 u(x, t) = u(x − ct) = u(z)，z = x − ct，其中 c 是恒定速度。在以速度 c 运动的行波坐标系中，波看起来是静止的。沿负 x 方向行进的波则形如 u(x + ct)。在 (x, t) 中寻找形如 (13.3) 的行波解时，∂u/∂t = −c du/dz，∂u/∂x = du/dz，原偏微分方程化为关于 z 的常微分方程。为了物理上合理，u(z) 必须对所有 z 有界且非负。

经典线性抛物方程理论（如 (13.1)）不存在物理上合理的行波解。设 u = u(x − ct) 代入 (13.1) 得 D d²u/dz² + c du/dz = 0，积分得 u(z) = A + Be^{−cz/D}。为使 u 在所有 z 上有界，必须令 B = 0，剩下 u ≡ A 是平凡解而非行波。与之形成鲜明对照的是，非线性反应扩散方程 (13.2) 在适当形式的反应/相互作用项 f(u) 下可以存在行波解。这种解行为是整个反应扩散理论数学领域兴起的主要推动因素之一。

作者指出，尽管多数现实的生物模型涉及多维和多个因变量，仍有多物种系统可合理地化简为一维单物种机制并保留关键特征。因此本章不仅是常见技巧与基础理论的数学教学性阐述，还讨论生态学与发育生物学中两个非常实际的问题——两者都属于建模扮演重要角色的领域。

13.2 Fisher–Kolmogoroff 方程与行波解（Fisher–Kolmogoroff Equation and Propagating Wave Solutions）

非线性反应扩散方程 (13.2) 最经典的简化情形是 ∂u/∂t = ku(1 − u) + D∂²u/∂x²（k、D 为正参数）。Fisher（1937）提出此方程作为种群中优势基因空间扩散随机模型的确定性版本。它也是第 11 章中讨论的 logistic 种群增长模型在种群通过线性扩散扩散时的自然推广。该方程及其行波解已被广泛研究，更一般形式 f(u) 取代 ku(1 − u) 的情形同样被深入讨论。Kolmogoroff、Petrovsky 与 Piscounov（1937）发表的论文是这一领域具有开创性意义的经典论文。Fife（1979）、Britton（1986）、Grindrod（1996）的专著给出该方程的完整讨论与详尽文献。

尽管 (13.4) 现在被称为 Fisher–Kolmogoroff 方程，但化学行波中波现象的发现、研究与分析最早由 Luther（1906）报道。这篇被重新发现的论文已被 Arnold 等（1987）译为英文。Luther 的论文最早在一次会议上报告；论文末尾的讨论（亦收录于 Arnold 等 1988 年的译本）非常有趣。Luther 指出波速是微分方程的简单推论。Showalter 与 Tyson（1987）将 Luther（1906）对化学波的卓越发现与分析置于现代背景之下。Luther 以其所研究反应的相关参数给出了波速，其解析形式与 Kolmogoroff 等（1937）和 Fisher（1937）对 (13.4) 的结果相同。

讨论 (13.4) 时，首先通过 t = kt、x = x√(k/D) 进行重新标度，去掉星号后方程化为 ∂u/∂t = u(1 − u) + ∂²u/∂x²。在空间齐次情形下，平稳态为 u = 0 和 u = 1，分别不稳定和稳定。这提示我们应寻找满足 0 ≤ u ≤ 1 的行波前解。若行波解存在，可写为 u(x, t) = U(z)，z = x − ct，其中 c 是波速。由于 (13.6) 在 x → −x 下不变，c 可正可负；为确定起见取 c ≥ 0。代入得 U 满足 U'' + cU' + U(1 − U) = 0（撇号表示对 z 求导）。典型的波前解在一端（如 z → −∞）处取一个平稳态，在另一端（z → +∞）处取另一个平稳态。因此需对 c 求特征值问题，找寻满足 lim_{z→∞} U(z) = 0、lim_{z→−∞} U(z) = 1 的非负解的存在条件。

在 (U, V) 相平面（U' = V、V' = −cV − U(1 − U)）中考察 (13.8)，相轨迹由 dV/dU = (−cV − U(1 − U))/V 给出。相平面有两个奇点 (0,0) 和 (1,0)。线性稳定性分析给出 (0,0) 的特征值 λ± = (−c ± √(c² − 4))/2：当 c² > 4 时为稳定结点，c² < 4 时为稳定螺旋（注意 c² < 4 时原点附近 U 振荡）；(1,0) 为鞍点。当 c ≥ c_min = 2 时，原点为稳定结点（c = c_min 时为退化结点）。

由相平面轨迹图的连续性论证或直观推理可知，对所有 c ≥ c_min = 2，存在一条从 (1,0) 到 (0,0) 的相轨迹，完全位于 U ≥ 0、U' ≤ 0 象限内，且 0 ≤ U ≤ 1。回到原始量纲方程 (13.4)，波速范围为 c ≥ c_min = 2√(kD)。图 13.1(b) 是典型行波解的示意图。当 c < 2 时也存在行波解，但因 U 在某处取负值而物理上不现实——此时 U 在前缘以围绕 U = 0 的衰减振荡趋于 0。

此时关键的问题是：原始 Fisher–Kolmogoroff 方程 (13.6) 配以什么样的初值 u(x, 0) 会演化出行波解，若解存在其波速 c 是多少。Kolmogoroff 等（1937）证明：若 u(x, 0) 具有紧支集（即在 x₁ < x < x₂ 外为 0），则 (13.6) 的解 u(x, t) 演化为以 c_min = 2 行进的波前解 U(z)，z = x − 2t。对非紧支集的初值，解的行为关键地依赖于 u(x, 0) 在 x → ±∞ 时的行为。

Mollison（1977）提出的简单分析可直观看到波速 c 对初值在无穷远处行为的依赖。考虑波的前缘处 u 很小，可略去 u² 相对于 u 的项，将 (13.6) 线性化为 ∂u/∂t = u + ∂²u/∂x²。设初值渐近行为 u(x, 0) ∼ Ae^{−ax}（x → ∞，a > 0），并寻找线性方程 (13.15) 的行波解 u(x, t) = Ae^{−a(x−ct)}。代入得色散关系 ca = 1 + a²，即 c = a + 1/a。绘制 c 关于 a 的色散曲线可见 c_min = 2 出现在 a = 1 处；对所有其他 a > 0，c > 2。

对大且正的 x，比较 e^{−ax} 与 e^{−x}。若 a < 1，则 e^{−ax} > e^{−x}，波的传播速度由前缘决定并由 (13.18) 给出；若 a > 1，则 e^{−ax} ≤ e^{−x}，波前以 c = 2 传播。于是当初值满足 (13.16) 时，(13.6) 的行波解的渐近波速为 c = a + 1/a（0 < a ≤ 1）或 c = 2（a ≥ 1）。前者由 McKean（1975）证明，后者由 Larson（1978）证明，并由 Manoranjan 与 Mitchell（1983）数值验证。

Fisher–Kolmogoroff 方程在 x 变号下不变，因此存在形式为 u(x, t) = U(x + ct)（c > 0）的行波解，此时 U(−∞) = 0、U(∞) = 1。因此若在 −∞ < x < ∞ 上求解 (13.6)，初值在有限域外为零，则 u(x, t) 将演化为两个波前，一个向左、一个向右传播，均以速度 c = 2 推进。需注意，若 u(x, 0) < 1，则 u(1 − u) 项使解增长直到 u = 1；显然 u(x, t) → 1（t → ∞，对所有 x）。

Fisher–Kolmogoroff 方程的轴对称形式 ∂u/∂t = ∂²u/∂r² + (1/r) ∂u/∂r + u(1 − u) 不存在波前以常速传播的行波解（因 1/r 项使方程在变量 z = r − ct 下不能化为常微分方程）。直观上，给定形如图 13.2 第一幅的 u(r, 0)，u 因 u(1 − u) 项而增长（因 u < 1），同时扩散使波动式地向外扩散。在"波"上 ∂u/∂r < 0，相当于在 (13.20) 右端产生一个等效的负对流或减弱源项 u(1 − u)，其效应是降低出射波的速度。对大 r，(1/r) ∂u/∂r 项可忽略，解渐近趋向速度 c = 2 的行波前解（与一维情形相同）。因此轴对称的波动解可被视为具有"波速" c(r)，它是 r 的函数，对有界远离 r = 0 的 r 而言，c(r) 随 r 单调增加，且 c(r) ∼ 2（r 充分大）。

方程 (13.4) 已成为多种空间扩散模型的基础。Aoki（1987）讨论了基因-文化推进波；Ammerman 与 Cavalli-Sforza（1971, 1983）将该模型直接应用于欧洲早期农业的扩散。

13.3 Fisher–Kolmogoroff 方程波前解的渐近解与稳定性（Asymptotic Solution and Stability of Wavefront Solutions of the Fisher–Kolmogoroff Equation）

(13.6) 的行波前解 U(z) 满足 U'' + cU' + U(1 − U) = 0；c > 2 时存在单调解，U(−∞) = 1、U(∞) = 0。相轨迹由 dV/dU = (−cV − U(1 − U))/V 给出。对一般 c 尚无解析解，但对特定 c > 2 存在一个精确解（见 13.4 节）。方程中有一个小参数 ε = 1/c² ≤ 0.25，提示对 0 < ε ≪ 1 寻找渐近解。Canosa（1973）对 (13.21) 给出了这种渐近解。

由于行波解对坐标原点的任何平移不变（方程在 z → z + 常数下不变），取 z = 0 为 U = 1/2 处。使用标准的奇异摄动技巧：在前缘附近（即 z = 0 附近）引入变量变换以得到关于小参数 ε 的 Taylor 展开。令 U(z) = g(ξ)，ξ = z/c = ε^{1/2} z，代入 (13.21) 得 ε d²g/dξ² + dg/dξ + g(1 − g) = 0，边界条件 g(−∞) = 1、g(∞) = 0，外加 g(0) = 1/2，且 0 < ε ≤ 1/c²_min = 0.25。方程对小 ε 看似标准奇异摄动问题（ε 乘以最高阶导数）；但实际上约化方程（ε = 0）本身就给出均匀有效的一阶近似——原因是非线性项 g(1 − g) 在两个边界上均为零。

将 g(ξ; ε) 展开为 g₀(ξ) + εg₁(ξ) + ...，由边界条件与 g(0; ε) = 1/2 得到 g₀(−∞) = 1、g₀(∞) = 0、g₀(0) = 1/2 以及对 i = 1, 2, ... 的 gi(±∞) = 0、gi(0) = 0。代入 (13.24) 并令 ε 同次幂相等得 O(1) 方程 dg₀/dξ = −g₀(1 − g₀)，其解为 g₀(ξ) = 1/(1 + e^ξ)；O(ε) 方程给出 g₁ 的微分方程。借助 g₀ 表达 g₁ 方程并用条件 (13.26) 积分得 g₁ = −g'₀ ln[4|g'₀|]，其显式形式见 (13.28)。换回原变量 (U, z) 后得均匀有效的渐近解 U(z; ε) = (1 + e^{z/c})⁻¹ + (1/c²) e^{z/c} (1 + e^{z/c})⁻² ln[4e^{z/c} / (1 + e^{z/c})²] + O(1/c⁴)，c ≥ c_min = 2。该渐近解在 c = 2 时精度最差，但与 (13.6) 计算所得 c = 2 的波前解相比，仅 O(1) 项 (1 + e^{z/c})⁻¹ 就在各处仅有百分之几的偏差——这对使用小参数的渐近解而言是令人鼓舞的。

用渐近解 (13.29) 考察波前解的斜率（陡度）与其传播速度的关系。波前梯度处处为负，陡度 s 取为最大梯度 |U'(z)| 的绝对值，即拐点处（U'' = 0）的值。由 (13.23) 与 (13.25)，该点出现在 ξ = 0（即 z = 0）处，梯度为 −U'(0) = s = 1/(4c) + O(1/c⁵)，该式只对 c ≥ 2 成立。结果显示：波速 c 越大，波前越不陡。波宽（严格意义下为 −∞ 到 +∞）的实际量度 L 取为陡度的倒数：L = 1/s = 4c。图 13.3 给出两个不同波速 c₁ > c₂ ≥ 2 下的陡度与宽度。

本节结果可推广到将 logistic 增长替换为适当 f(u) 的单物种种群模型：∂u/∂t = f(u) + ∂²u/∂x²，其中 f(u) 仅有两个零点 u₁ 和 u₂ > u₁。若 f'(u₁) > 0、f'(u₂) < 0，则由类似分析可得单调从 u₁ 到 u₂ 的波前解，波速满足 c ≥ c_min = 2√(f'(u₁))。这与 (13.32) 通过在前缘 u ≈ u₁ 处线性化 f(u) 并与 (13.15) 比较得到的结果一致。

生物模型解的稳定性是重要议题，也常是模型机制可信度的检验。Fisher–Kolmogoroff 方程的行波前解提供了稳定性研究的教学性案例。波前解的传播速度（见 (13.19) 与 (13.16)）敏感地依赖于初值 u(x, 0) 在 |x| → ∞ 处的显式行为，这意味着波前解对远场扰动不稳定。但若 u(x, 0) 具有紧支集（Kolmogoroff 等 1937 使用的初值 (13.14)），则最终的波不依赖于 u(x, 0) 的细节形式，演化的波速为 c = 2。

数值分析中若不预先知晓期望的波速，随机效应（被限制在有限域内）可能干扰结果。实际模型处理的是有限域，因此考察波解对"在有限域外为零的扰动"的稳定性很重要。Canosa（1973）表明在行波运动坐标系下，解对这种有限扰动是稳定的。令 u(x, t) = u(z, t)，z = x − ct，则 (13.6) 化为 u_t = u(1 − u) + cu_z + u_zz。设 u(z, t) = u_c(z) + ωv(z, t)（0 < ω ≪1），代入并保留一阶得 v_t = (1 − 2u_c(z))v + cv_z + v_zz。u_c(z) 对扰动 v(z, t) 稳定若 lim_{t→∞} v(z, t) = 0 或 lim_{t→∞} v(z, t) = du_c(z)/dz（后者表示沿 x 轴的小平移，因为 u_c(z + δz) ≈ u_c(z) + δz du_c(z)/dz）。

设 v(z, t) = g(z) e^{−λt}，代入 (13.35) 得 g'' + cg' + (λ + 1 − 2u_c(z)) g = 0。注意 λ = 0 时 g(z) = du_c(z)/dz 是其解，对应行波沿 z 轴的平移不变性。v(z, t) 仅在有限域内非零意味着 g 在 g(±L) = 0 处有边界条件。令 g(z) = h(z) e^{−cz/2}，特征值问题化为 h'' + (λ − 2u_c(z) − c²/4 + 1) h = 0、h(±L) = 0。由于 c ≥ 2 且 u_c(z) > 0（对 −L ≤ z ≤ L），有 2u_c(z) + c²/4 − 1 ≥ 2u_c(z) > 0。标准理论（Titchmarsh 1946 第 11 章）给出 (13.38) 的所有特征值 λ 均为正实数，因此 v(z, t) → 0（t → ∞），即 u_c(z) 对所有形式为 (13.34) 的小有限域扰动稳定。Larson（1978）等人的研究表明一般问题（行波坐标系外的扰动）分析更复杂。波对有限域扰动的稳定性解释了为什么 Fisher–Kolmogoroff 方程的典型数值模拟得到稳定的、c = 2 的波前解。

13.4 密度依赖扩散-反应扩散模型与若干精确解（Density-Dependent Diffusion-Reaction Diffusion Models and Some Exact Solutions）

第 11 章 11.3 节中已指出，在某些昆虫扩散模型中扩散系数 D 依赖于种群 u，但未包含增长动力学。若考虑更长的时间尺度，应在模型中纳入增长项。在一维情形下，自然的推广是考虑形如 ∂u/∂t = f(u) + ∂/∂x [D(u) ∂u/∂x] 的方程，典型地 D(u) = D₀ uᵐ（D₀、m 为正常数）。此处 f(u) 假设有两个零点（u = 0 和 u = 1）。f ≡ 0 的方程（参见第 11 章）被研究得远多于 f 非零情形。为更具体，考虑 f(u) = kuᵖ(1 − uᵍ)，p、q 为正常数；通过对 t、x 的适当重新标度吸收 k 与 D₀ 后，化为 ∂u/∂t = uᵖ(1 − uᵍ) + ∂/∂x [uᵐ ∂u/∂x]，其中 p、q、m 为正参数。将扩散项展开得 ∂u/∂t = uᵖ(1 − uᵍ) + muᵐ⁻¹ (∂u/∂x)² + uᵐ ∂²u/∂x²，可将非线性扩散视为以"速度" −muᵐ⁻¹ ∂u/∂x 贡献一个等效对流。

(13.40) 的形式看似特殊，但 p、q、m 的选择余地很大，足以定性模仿更复杂、只能数值求解的模型。解析解的用处在于能清楚地看到解对参数的解析依赖，从而定性推断更复杂更现实模型的行为。但常有一些隐患，一重要者在下面会指出。

为与 Fisher–Kolmogoroff 方程的结果对应，先取 m = 0、p = 1，(13.40) 化为 ∂u/∂t = u(1 − uᵍ) + ∂²u/∂x²，q > 0。u = 0、u = 1 是齐次平稳态。寻找行波解 u(x, t) = U(z)，z = x − ct，U(−∞) = 1、U(∞) = 0，c > 0 待定。U(z) 满足的常微分方程为 L(U) = U'' + cU' + U(1 − Uᵍ) = 0。借助 Fisher–Kolmogoroff 渐近波前解第一项的形式 (13.29)，设 U(z) = 1/(1 + a e^{bz})ˢ，a、b、s 为正常数待定。此形式自动满足 (13.42) 的边界条件。由于方程的平移不变性，a 是任意的（可吸收为 b⁻¹ ln a 的平移）。

将 (13.44) 代入 (13.43) 经整理得 L(U) = 0 对所有 z 成立要求 e⁰、e^{bz}、e^{2bz} 各项系数均为零，由此得 2 − sq = 0、1 或 2，即 s = 2/q、1/q 或 sq = 0。后者不可能（s、q 均正）。取 s = 1/q 时由 e^{bz} 系数得 sb(b + c) = 1，由 e^{2bz} 系数得 s(s + 1)b² = 0 ⇒ b = 0（与 b > 0 矛盾）。最后取 s = 2/q 时，e^{bz} 与 e^{2bz} 系数给出 sb(b + c) = 2、s(s + 1)b² = 1，由此 b = 1/√(s(s + 1))、c = 2/(sb) − b。结合 s = 2/q 整理得 s = 2/q、b = q/√(2(q + 2))、c = (q + 4)/√(2(q + 2))。可见波速 c 随 q 增大而增大。波前斜率（拐点处梯度绝对值）由 (13.44) 计算得 S = b / (1 + 1/s)^{s+1} = (1/(2q)) (1 + q/2)^{3/2 + 2/q}，即随 q 增大，c 增大而陡度下降，与 Fisher–Kolmogoroff 情形一致。

q = 1 时 (13.41) 回到 Fisher–Kolmogoroff 方程 (13.6)，由 (13.46) 得 s = 2、b = 1/√6、c = 5/√6。选取 a 使 z = 0 对应 U = 1/2，则 a = √2 − 1，得精确行波解 U(z) = 1 / (1 + (√2 − 1) e^{z/√6})²。该解波速 c = 5/√6 ≈ 2.04，与渐近解 (13.29) 的 O(1) 项相比要陡得多。这一例子突显了精确解的严重问题：它们常常不能给出所有可能的解，事实上可能并未给出最相关的那个（这里波速虽接近 2，但波形差异显著）。要正确分析 (13.43) 的一般形式，必须进行仔细的相平面分析。

另一类精确解可对 (13.40) 中 m = 0、p = q + 1（q > 0）找到，方程化为 ∂u/∂t = u^{q+1}(1 − uᵍ) + ∂²u/∂x²。代入 (13.44) 形式的 U(z) 完全类似地处理，得行波前解存在且波速唯一：U(z) = 1/(1 + a e^{bz})ˢ，s = 1/q、b = q/√(q + 1)、c = 1/√(q + 1)。

更值得注意的有用精确解是 p = q = 1、m = 1 的情形，(13.40) 化为 ∂u/∂t = u(1 − u) + ∂/∂x (u ∂u/∂x)，这是一个具有 logistic 种群增长且密度依赖扩散的非平凡例子。物理上，模型意味着种群在越拥挤的地方扩散到低密度区域的速度越快。该解由 Aronson（1980）和 Newman（1980）独立发现，Newman（1983）研究了更一般形式并推进了工作。

对该模型寻找行波解 u(x, t) = U(z)、z = x − ct，考虑 (UU')' + cU' + U(1 − U) = 0，相平面系统为 U' = V、UV' = −cV − V² − U(1 − U)。我们对 U(−∞) = 1、U(∞) = 0 的波前解感兴趣，预期 U' < 0。第二个方程在 U = 0 处有奇点。引入新变量 ζ 使 U d/dz = d/dζ（即 dU/dζ = UV、dV/dζ = −cV − V² − U(1 − U)）以消除奇点。(U, V) 相平面的临界点为 (0,0)、(1,0)、(0, −c)。线性分析表明 (1,0) 与 (0, −c) 为鞍点，(0,0) 为稳定非线性结点。

由 11.2 节可知波可能在某特定点 z_c 出现切线不连续（即 U ≡ 0 对 z ≥ z_c）。相轨迹需从 (1,0) 出发在 U = 0 轴上某有限非零负 V 处终止。分析显示：若 0 < c < c_min，从 (1,0) 到 U = 0 的轨迹不存在（除非 V 趋于无穷的不现实情形）；当 c 增大到分岔值 c_min 时，存在从 (1,0) 到 (0, −c_min) 的唯一轨迹，意味着在波前 z_c（U = 0）处导数从 V = U' = −c_min 不连续地跳到 U' = 0；c 超过 c_min 后，从 (1,0) 到 (0,0) 的轨迹始终存在，且波解满足 U → 0、U' → 0（z → ∞）。图 13.4 给出各情形的相轨迹与典型波解示意。

对 (b) 中连接 (1,0) 到 (0, −c) 的轨迹，存在精确解：相轨迹是直线 V = −c_min (1 − U)，c_min 适当选取。代入相平面方程验证得 c_min = 1/√2。回到相平面方程 U' = V 并使用 V = −(1 − U)/√2，得 U' = −(1 − U)/√2，由 U(−∞) = 1 积分得 U(z) = 1 − exp((z − z_c)/√2)（z < z_c）、U = 0（z > z_c），其中 z_c 是波前，可任意选取（与 (13.44) 中 a 的作用相同）。该分析可推广到扩散系数为 uᵐ（m ≠ 1）或更一般 D(u)（需满足一定条件）的情形。

本节末尾，作者简要总结了 Fisher–Kolmogoroff 方程及其对一般 f(u)（规范化使 f(0) = 0 = f(1)、f'(0) > 0、f'(1) < 0）的推广结果：方程 (13.31) 在适当初值下可演化出 0 < u < 1 的行波前解，其波速 c ≥ c_min = 2√(f'(0))，通常计算所得解的波速为 c_min。对 Fisher–Kolmogoroff 方程 (13.4)，用非量纲化 (13.5) 回到量纲下波速 c* = 2√(kD)；取 D = 10⁻⁹ ~ 10⁻¹¹ cm²·sec⁻¹、k = O(1 sec⁻¹) 等典型生物值，传播速度为 O(2 × 10⁻⁴·⁵ ~ 10⁻⁵·⁵ cm·sec⁻¹)，覆盖 1 mm 量级距离的时间为 O(5 × 10²·⁵ ~ 10³·⁵ sec)，远小于纯扩散时间 O(10⁷ ~ 10⁹ sec)。正是反应与扩散的耦合大大增强了行波浓度变化所传递信息的效率。

本节还顺带提及：当对流项自然地作为守恒律的推广出现时，会得到形如 ∂u/∂t + ∂h(u)/∂x = f(u) + ∂²u/∂x² 的方程（"h(u)"由 Goldstein 与 Murray（1959）等研究，在离子交换柱和色谱中产生），其左端为标准守恒形式，对流"速度"为 h'(u)。非线性对流对反应扩散方程的解可有戏剧性影响。当扩散相对于对流可忽略时，解可呈现类激波行为（Murray 1968, 1970a,b, 1973）。

为分析带非线性对流的行波前存在条件，考虑 h'(u) = ku、f(u) 为 logistic 的简单但非平凡情形，方程化为 ∂u/∂t + ku ∂u/∂x = u(1 − u) + ∂²u/∂x²。k = 0 时回到 (13.6)。k ≠ 0 时寻找行波解 u(x, t) = U(z)、z = x − ct，代入得 U'' + (c − kU)U' + U(1 − U) = 0，边界条件 (13.58)。相轨迹 dV/dU = (−(c − kU)V − U(1 − U))/V，奇点为 (0,0) 和 (1,0)。要求存在 0 ≤ U ≤ 1、U'(z) ≤ 0 的单调解（即位于 U ≥ 0、V ≤ 0 象限连接两个奇点的相轨迹）。线性分析给出 c ≥ 2 即可使 (0,0) 为稳定结点、(1,0) 为鞍点。该方程及其行波解由作者与 R.J. Gibbs 研究（Murray 1977），可证对所有 c ≥ c(k) 都存在行波解，其中 c(k) = 2（2 > k > −∞）或 k/2 + 2/k（2 ≤ k < ∞）。可见此处 c = 2 仅在 k 的有限范围内是下界，更精确的下界由上式给出。

(13.59) 关于 (0,0) 的线性化得 dV/dU = (−cV − U)/V，特征值 e± = (−c ± √(c² − 4))/2。要求 U ≥ 0 需特征值为实数，即 c ≥ 2。此时 0 > e+ > |e−|，原点为稳定结点，对大 z 有 V/U 比例关系。当 a ≠ 0 时 dV/dU → e+/e−（z → ∞）；a = 0 时亦由对应特征值决定。(13.59) 有一个精确解 V = −(k/2)U(1 − U)，对应 c = k/2 + 2/k。结合 (c² − 4)^{1/2} 的分段表达和 e± 的分段结果，对 |k| ≥ 2 得 dV/dU|_{U=0} = −k/2 = e−/e+，从而由 (13.62) 得 dV/dU → e−（z → ∞），即 (13.64) 给出 c = k/2 + 2/k（|k| ≥ 2）。

对 k < 2 在 z → −∞ 附近线性化关于 (1,0) 得特征值 E± = (−(c − k) ± √((c − k)² + 4))/2；E+ > 0 > E−，(1,0) 为鞍点。z → −∞ 时 U → 1 − O(exp(E+z))，故 dV/dU → E+(c, k)（z → −∞）。由 (13.66) 知 dE+(k)/dk > 0（c ≥ 2 时），因此 U 充分接近 1 时 dV/dU 随 k 增大而增大。于是对 U 充分接近 1，相平面轨迹 V(U, c, k) 满足 V(U, c = 2, k) < V(U, c = 2, k = 2)（k < 2）。由此与 (13.59) 在 d 点的值比较可推出矛盾，进而证明 c ≥ 2（k < 2）。综合得 (13.60) 完整结果。

k ≠ 0 时可将 (13.55) 改写为突出非线性对流贡献的形式。k > 0 时令 ε = 1/k²、y = x/|k| = ε^{1/2} x，方程化为 u_t + u u_y = u(1 − u) + ε u_yy（(13.69)）；k < 0 时形式相同（(13.70)）。对 (13.55) 的解 u(x, t)，u(ky, t) 是 (13.69) 的解；U(x − ct) 是 (13.59) 满足边界条件的解时，U(ky − ct) 是 (13.69) 的解且波速 λ = c/k = c ε^{1/2}。用 (13.60) 的波速估计，(13.69) 的行波解在 λ ≥ λ(ε) = c(k)/|k| = c(ε^{−1/2}) ε^{1/2} 时存在，得 λ(ε) = 2 ε^{1/2}（ε > 1/4）或 1/2 + 2ε（1/4 ≥ ε > 0）。

ε → 0 时 (13.69) 化为 u_t + u u_y = u(1 − u)，其解可以是不连续的（弱解或激波解，Murray 1970a 详细讨论）。ε 小时波陡化为类激波解；而 (13.70) 在同样边界条件下不会出现不连续解。图 13.5 与 13.6 给出对两个方程分别在不同 ε 下的数值行波解。注意 (13.69) 中 ε 小时解出现不连续（图 13.5(b)），(13.70) 中 ε 小时波陡化但不出现不连续。

本节最后提到 Satsuma（1987）对标量密度依赖反应扩散方程精确解的研究，他所发展的方法新颖，可能有更广适用性。Hosono（1986）关于此类方程单调行波解的存在性与稳定性也直接相关。

13.5 多平稳态动力学模型中的波：昆虫种群的扩散与控制（Waves in Models with Multi-Steady State Kinetics: Spread and Control of an Insect Population）

酶反应系统中（参见第 6 章）的摄取函数或第 1 章中引入的种群增长-相互作用函数 f(u) 常常有不止两个平稳态，即 (13.31) 中的 f(u) 可能有三个或更多正零点。此类 f(u) 对应的波现象与前几节截然不同。实际的例子是云杉芽虫（spruce budworm）的增长函数，其空间齐次情形在第 1 章 1.2 节中详细讨论。该模型中具体的无量纲 f(u) 为 f(u) = ru (1 − u/q) − u²/(1 + u²)，其中 r、q 是与实际现场参数相关的无量纲参数（见 (1.17)）。对正参数 r、q 的某个范围，f(u) 形如图 1.5（重绘于图 13.7(a)）。回想平稳态的数目与大小对 r、q 的依赖；典型曲线再次示于图 13.7(b)。在无扩散（空间齐次）的情形下，可以有三个正平稳态：u₁ 和 u₃ 线性稳定，u₂ 不稳定；u = 0 也不稳定。

1.2 节中 u₁ 对应于芽虫的"避难所"状态，u₃ 对应于"爆发"状态。本节考虑两个问题：(i) 当我们纳入芽虫的空间扩散后，虫害或爆发如何传播；(ii) 是否能用分析结果提出控制策略防止爆发扩散。为回答这两个问题，假设芽虫通过线性扩散扩散，研究其行波可能性。实际问题显然是二维的，但这里讨论一维情形，因为它已能给出合理的回答，并至少能提出二维模型必须面对的问题；事实上二维模型并无本质上的新概念困难。所考虑的模型来自 (13.31)：∂u/∂t = f(u) + ∂²u/∂x²，f(u) 典型地如图 13.7(a)。

设行波解 u(x, t) = U(z)，z = x − ct，则 U 满足 U'' + cU' + f(U) = 0。相平面系统 dV/dU = (−cV + f(U))/V 有四个奇点 (0,0)、(u₁,0)、(u₂,0)、(u₃,0)。对 (U = 0 和 U = uᵢ, i = 1,2,3) 线性化 (13.74) 并作标准线性相平面分析得奇点类型：(0,0) 与 (u₂,0) 是 f'(0) > 0 与 f'(u₂) > 0 时稳定螺旋（c² < 4f'）/稳定结点（c² > 4f'）；(uᵢ, 0)（i = 1, 3）因 f'(uᵢ) < 0 而对所有 c 均为鞍点。c < 0 时 (0,0) 与 (u₂,0) 变为不稳定（类型相同）。c 变化时各种可能的相平面轨迹出现，整个存在性研究可由 Fife 与 McLeod（1977）严格处理；本节采取的直观方法虽未严格证明波的存在性，但强烈地表明其存在。

举例而言，若 c² > 4 max[f'(0), f'(u₂)]，则 (0,0) 与 (u₂,0) 为稳定结点。相平面图（13.8(a)）给出可能的奇点连接：若用 d₁ 包含原点和鞍点 (u₁, 0)，与图 13.1(b) 类似，可合理假设存在从 U(−∞) = u₁ 到 U(∞) = 0 的波解，且对所有 c ≥ 2√(f'(0)) 存在。其它域内类似地存在其他行波解。

特别地，关注 (u₁, 0) 与 (u₃, 0)（均为鞍点）。特征值 λ_{1,2} = (−c ± √(c² − 4f'(uᵢ)))/2（i = 1, 3，f'(uᵢ) < 0），对应特征向量 e_{i1}、e_{i2}。c 增大时特征向量趋于向 U 轴靠近。c 变化时相轨迹变化；尤其是图 13.8(a) 中标 T₁ 和 T₃ 的轨迹会变。由连续性论证可合理地预期，若 f'(u₁) 与 f'(u₃) 在适当范围，存在唯一的 c = c 使 T₁ 与 T₃ 轨迹相接，从而得到连接 (u₁, 0) 与 (u₃, 0) 的相路径（13.9(a)），对应波解（13.9(b)），波以唯一速度 c 传播，依赖于 f(u)。此情形下 U(−∞) = u₃、U(∞) = u₁，正是芽虫问题所要考虑的。

设初始 u = u₁（芽虫种群处于稳定的避难所态），但某有限域内局部增至 u₃（局部爆发）。为研究爆发扩散可能性，可更简单地提出代数上更简单的问题：是否存在从 u = u₁ 区域到 u = u₃ 区域的行波前解（如图 13.9(b)），若存在其波速与方向如何。由上面讨论预期存在这样的波。若 c > 0 波向 u₁ 区域推进，爆发扩散；c < 0 则不仅不扩散反而被压缩。

c 的符号可由 (13.73) 乘以 U' 并从 −∞ 到 ∞ 积分得到 c ∫ U'² dz = −∫ U' f(U) dz = −∫{u₁}^{u₃} f(U) dU。由于乘以 c 的项总为正，c 的符号由 ∫ f(u) du ⪌ 0，亦即面积 A₃ > A₁ 时 c > 0，A₃ < A₁ 时 c < 0）。该积分的符号等价于图 13.10 中比较面积 A₁ 与 A₃：若 A₃ > A₁ 则 c > 0，爆发扩散到避难所区域，此时称 u₃ 占优，t → ∞ 时 u → u₃ 处处成立；反之 A₃ < A₁ 时 c < 0，u₁ 占优，u → u₁，爆发被消除。}^{u₃} f(u) du 的符号决定（按源文记号 c ⪌ 0 对应 ∫_{u₃}^{u₁

从虫害控制的角度看，若虫害爆发正在扩散，我们希望知道如何改变局部条件使波被遏制或反向。需在局部将芽虫增长动力学改变，使新面积 A₁ > A₃。若 f(u) 的两个最大零点 u₂ 与 u₃ 彼此更接近即可实现此点。由图 13.7(b) 知这可通过减小 (13.71) 中的无量纲参数 q 来达成。芽虫模型的无量纲化（参见 1.2 节）将 q 与环境的基本承载力 K_B 相关联；故 q 的实际减小可通过例如喷洒带状区域以降低树叶承载力来实现，从而造成"虫害断裂带"，使 u₁ 占优，波速 c 不再为正。实际问题是"断裂带"需要多宽才能阻止爆发穿透；该问题需要仔细建模考虑，因为方程的抛物（扩散）特性使前缘很长。该思想与第一卷第 13 章中关于遏制狂犬病扩散的讨论密切相关，相关方法可直接应用于此处芽虫扩散的"断裂带"问题。

第 6 章曾简要描述过 (6.120) 形式的可作为生化开关的动力学。该机制下，从一个平稳态出发的足够大扰动可使系统跃迁到另一平稳态。一个重要的实验例子是钙刺激钙释放（calcium-stimulated-calcium-release）机制：当 Ca²⁺ 浓度被扰动到某阈值以上时，会引发被封存的钙进一步释放（即系统迁移到另一平稳态）。这发生于某些受精两栖类卵膜上的钙位点（下一节讨论一个真实例子）。除释放钙外，该膜还会重新封存钙。记 Ca²⁺ 浓度为 u，可建立动力学速率律 du/dt = A(u) − r(u) + L，其中 L 表示小量泄漏，A(u) 是钙的自催化释放，r(u) 是其重新封存。假设钙的重新封存服从一级动力学、自催化钙的产生在高 Ca²⁺ 时饱和，由此得反应动力学模型 (13.82) du/dt = L + k₁u²/(k₂ + u²) − k₃u = f(u)，其中各 k 与 L 为正参数。若各 k 满足特定关系（参见第 6 章末习题 3），f(u) 对 L 充分小可有三个正平稳态。该可激动力学情形下 f(u) 的形式示于图 13.11(a)。可激过程有密切相关的两种。

L = 0 时有三个平稳态（两个稳定，一个不稳定）。L 从零增大时先有三个正平稳态 uᵢ(L)（i = 1, 2, 3），u₁ 和 u₃ 线性稳定，u₂ 不稳定；L 增大到阈值 L_c 时 u₁ 与 u₂ 先合并再消失。故若初值 u = u₁，足够大的 L 脉冲可使平稳态跃迁到 u₃（两个稳定态中较大的）并保持。该定量分析并非简单，Kath 与 Murray（1986）在蝴蝶翅图式生成的模型机制中进行了处理（参见第一卷第 3 章）。

第二种可激性是 L 固定、f(u) 形如图 13.11(b) 曲线 du/dt (= f(u))。箭头方向表示给定浓度扰动下 u 如何变化：所有 0 < u < u₂ 时 u → u₁，所有 u > u₂ 时 u → u₃，u₂ 为阈值浓度。前一种情形以 L 为分岔参数，此处以扰动幅度（与 u₂ 相关）为分岔参数。

实际的钙刺激钙释放过程很复杂，(13.81) 及其定量形式 (13.82) 只能是合理的简化抽象。合理地，可对 f(u) 作进一步简化抽象（只要保持其定性动力学行为与所需零点数），即用具有三个正零点的三次式 f(u) = A(u − u₁)(u₂ − u)(u − u₃) 替代，A 为正常数且 u₁ < u₂ < u₃。该三次式定性如图 13.11(a) 中 0 < L < L_c 时的曲线。

对具此反应动力学的反应扩散方程 ∂u/∂t = A(u − u₁)(u₂ − u)(u − u₃) + D ∂²u/∂x²（未重整化以突出 A 与 D 的作用），可假设其有形如 u(x, t) = U(z)、z = x − ct、U(−∞) = u₃、U(∞) = u₁ 的波前解，代入得 L(U) = DU'' + cU' + A(U − u₁)(u₂ − U)(U − u₃) = 0。

为寻找精确解，尝试让 U 满足一个更简单的方程（其解呈指数行为）而同时又能满足 (13.85)。设 U' = a(U − u₁)(U − u₃)，其解在 z → ∞ 时指数地趋于 u₁ 和 u₃，符合所需行为。代入 (13.85) 得 L(U) = (U − u₁)(U − u₃) [(2Da² − A) U − (Da²(u₁ + u₃) − ca − Au₂)]。令 L(U) = 0 须有 2Da² − A = 0、Da²(u₁ + u₃) − ca − Au₂ = 0，由此得 a = (A/(2D))^{1/2}、c = (AD/2)^{1/2} (u₁ − 2u₂ + u₃)。由此方程 (13.86) 的解即为 (13.85) 的解。求解 (13.86) 得 U(z) = (u₃ + K u₁ exp[a(u₃ − u₁) z]) / (1 + K exp[a(u₃ − u₁) z])，K 为任意常数（确定 z 平面原点的位置）。U(−∞) = u₃、U(∞) = u₁。c 的符号由 uᵢ (i = 1, 2, 3) 相对大小决定：u₂ 大于 u₁、u₃ 平均值时 c < 0，否则为正。这与用 (13.80) 积分公式加 (13.83) 的三次 f(U) 所得结果相同。

(13.83) 及其某些推广由 McKean（1970）研究，产生于神经动作电位传播的简单模型背景（参见第一卷第 1 章）。(13.83) 有时称为约化 Nagumo 方程，与第 7.5 节中讨论的 FitzHugh–Nagumo 神经动作电位模型相关。

13.6 两栖类卵上的钙波：Medaka 卵上的激活波（Calcium Waves on Amphibian Eggs: Activation Waves on Medaka Eggs）

两栖类卵的皮层（cortex）是包裹卵的一种膜壳。紧接受精后、第一次卵裂前，Ca²⁺ 的若干化学波扫过皮层。波的起始点附近的卵顶部是动物极（animal pole），由精子入卵点决定；底部是植物极（vegetal pole）。波从精子入卵点发出。每一道波是发育中某重大事件的前兆，且每道波后伴随一个机械事件。此类 Ca²⁺ 波称为激活波。图 13.12(a) 展示了硬骨鱼 Medaka 卵上此类钙波的前进过程，数据来自 Gilkey 等（1978）的实验。本节构造的模型是化学波的简化机制，源自 Cheer 等（1987）和 Lane 等（1987）关于脊椎动物卵皮层波的论文。Lane 等（1987）还基于分段线性方法给出一些解析结果，与完整非线性系统的数值模拟吻合良好。图 13.12(d) 显示伴随钙波的机械表面波。

这里为 Ca²⁺ 构造的简单模型基于钙动力学是可激的这一事实，使用上节描述的钙刺激钙释放机制。假设 Ca²⁺ 在卵的皮层（表面）上扩散，于是得到一个反应扩散模型，其反应与扩散都发生在球面上。Ca²⁺ 波前实际是在表面传播的环，其数学描述只涉及一个独立变量 θ（从球顶部量起的极角），0 ≤ θ ≤ π。动力学涉及通过钙刺激钙释放机制从表面位点释放钙。这里的小泄漏是少量 Ca²⁺ 扩散到卵内部所致。存在一个钙浓度阈值，触发表面位点的钙"倾卸"。捕获可激动力学与一些已知事实的唯象模型由 (13.82) 给出。采用 (13.83) 中更简单的三次动力学 caricature，得到模型反应扩散系统 (13.89)：∂u/∂t = f(u) + (D/R²)(∂²u/∂θ² + cot θ ∂u/∂θ)，其中 f(u) = A(u − u₁)(u₂ − u)(u − u₃)，A 为正常数，R 为卵的半径（模型中的参数）。

若假设 Ca²⁺ 在卵表面均匀处于较低的平稳态 u₁，扰动超过阈值 u₂ 则 u 趋向高平稳态 u₃；扰动小于 u₂ 则返回 u₁。因此存在一个触发阈值，高于它则 u → u₃。

考虑 (13.89) 的可能波解。若 cot θ 项不在方程中，已知方程有 (13.84) 类型的波前解：u(θ, t) = U(z)，z = Rθ − ct，U(−∞) = u₃、U(∞) = u₁。当然球形卵问题中 t = 0 起时间时 z 不能趋于 −∞；并且 cot θ 项在方程中。然而为直观理解波在球面上传播时的行为，可作如下论证：对每个固定 θ 假设有形如 (13.91) 的波前解，代入 (13.89) 得 DU'' + (c + (D/R) cot θ) U' + A(U − u₁)(u₂ − U)(U − u₃) = 0。对固定 θ 此方程与 (13.85) 形式相同，只是 c 替换为 c + (D/R) cot θ。因此可合理地预期，卵表面上波前型解的传播速度为 (13.93) 形式 c = (AD/2)^{1/2} (u₁ − 2u₂ + u₃) − (D/R) cot θ。

(13.93) 的物理含义是：波从动物极（θ = 0）向植物极（θ = π）移动时波速变化。在 0 < θ < π/2 区域 cot θ > 0，波在上半球较慢；在 π/2 < θ < π 区域 cot θ < 0，波在（北半球的）下半球较快。这可由反应扩散方程 (13.89) 与一维版本 (13.83) 的扩散项比较看出：若波从 u = u₃ 区域向 u = u₁ 区域传播，∂u/∂θ < 0；在动物半球 cot θ > 0，cot θ ∂u/∂θ < 0 意味着扩散过程被有效减弱（而扩散对波的传播至关重要），故波在上半球减速；下半球 cot θ ∂u/∂θ > 0，波加速。这与直观上一致：上半球波前需不断扩张其周长，下半球则相反。

(13.93) 表明对球面上的表面波，不可能对所有 θ 都存在 c > 0 的行波解 (13.89)；这显然依赖于参数（参数需细致地依赖于空间）。

回到真实生物学：所表明的是，钙刺激钙释放机制的简化模型可给出卵表面上的行波钙波前型解。结合实验与各相关时间的量级，可估算相关参数。然而真实卵中前缘行为与模型卵之间存在一个重要的定性差异：真实卵中波在植物半球减速，而模型中加速。Cheer 等（1987）得出的一个重要推论是：皮层性质的非均匀性克服了球面波传播的自然加速趋势；由 (13.93) 知 AD 与 uᵢ (i = 1, 2, 3) 必随 θ 变化。该波速公式在参数在球面上缓慢变化时也成立。因此可以解析地确定模型性质中的定性行为以实现卵上正确的波传播性质，从而推断卵皮层性质可能的参数变化。图 13.12(b) 给出 Cheer 等（1987）使用上述非均匀参数模型的数值结果。读者可参考原文详细讨论生物、全模型与所得生物结论。第一卷第 6 章中引入并详细讨论生物图式生成的新力学化学方法，本节与 Cheer 等（1987）、Lane 等（1987）即是其例子。

13.7 具有扩散可变性的入侵波速（Invasion Wavespeeds with Dispersive Variability）

昆虫、种子、动物、疾病等对新领地的殖民化是生态学和流行病学中极为重要的问题。设计例如生物控制方案时至少需要对此过程有所了解。Kot 等（1996）的论文与该问题尤为相关。这里虽限于讨论连续模型，离散增长与扩散模型也同等重要。本章讨论的这类模型已被广泛用于估计入侵速度；参见 Shigesada 与 Kawasaki（1997）关于入侵问题的专著。

简单标量方程的连续模型在实际世界中有局限：种群的每个成员并不一定以相同方式扩散，总存在某种变异性。本节讨论 Cook（Julian Cook，个人通讯 1994）对这一问题的开创性贡献——他重新考察了经典 Fisher–Kolmogoroff 模型，研究个体扩散的变异性对入侵波速可能产生的影响。这一考察过去几十年间研究此标量方程的所有学者都完全忽视了。变异性对入侵速度的影响是相当出人意料的、直观上完全不可见的。

从基本的一维 Fisher–Kolmogoroff 方程开始：种群以 logistic 方式增长，在齐次环境中以恒定扩散系数 D、本征线性增长率 r、承载力 K 扩散。由 13.2 节分析，波速（即入侵速度）为 2√(rD)，即 (13.13) 中的最小速度。考虑种群分为扩散者与不扩散者，两个亚种群完全自由交配，所有新生个体以相同、固定概率为扩散者。该模型严格说并非单物种模型，但因其与经典 Fisher–Kolmogoroff 模型的密切联系仍归入本章。

将种群分为扩散者 A 与不扩散者 B。模型系统取为 ∂A/∂t = D ∂²A/∂x² + r₁(A + B)[1 − (A + B)/K]、∂B/∂t = r₂(A + B)[1 − (A + B)/K]。D 是扩散亚种群的扩散系数（与全种群的平均扩散率不同），K 是环境承载力，r 是每头总种群的本征增长率。新生个体为扩散者的概率为 p = r₁/(r₁ + r₂)。若 r₂ = 0，全种群都扩散，系统化为标准 Fisher–Kolmogoroff 方程。

Cook 指出该模型中个体要么是扩散者（恒定 D）、要么是不扩散者（D = 0），是离散变异性情形。尽管模型基于 logistic 增长，与修改后 Fisher–Kolmogoroff 方程类似，分析可在更一般增长函数下进行；入侵速度受到类似影响但不改变一般原则。

作为分析第一步，对系统无量纲化：u = A/K、v = B/K、T = Rt、X = √(R/D) x，其中 R = r₁ + r₂。p = r₁/(r₁ + r₂) 给出新生个体为扩散者的概率。系统化为 ∂u/∂T = ∂²u/∂X² + p(u + v)[1 − (u + v)]、∂v/∂T = (1 − p)(u + v)[1 − (u + v)]。

按常规寻找行波解：u = U(X − CT)、v = V(X − CT)、Z = X − CT，C 为波速（C > 0 时波向 X 增大方向移动）。代入得 Z 的一阶常微分方程组：−CU_Z = U_ZZ + p(U + V)[1 − (U + V)]、−CV_Z = (1 − p)(U + V)[1 − (U + V)]。寻找满足 U + V → 1（Z → −∞）与 U = V = 0（Z → ∞）的行波解。设 W = U_Z，(13.99) 与 (13.100) 化为 (U, V, W) 中的一阶方程组。行波解对应于连接两个平稳态（具体为 (0, 0, 0) 和 (U₀, 1 − U₀, 0)）的异宿轨道。

在零平稳态 (0, 0, 0) 附近，对线性化系统分析得使解不取负、且使 U、V 分量同号（对应方向为异宿轨道在趋于 (0, 0, 0) 时的特征向量方向）的条件。令 W = U_Z，(13.99) 与 (13.100) 关于 (0, 0, 0) 线性化（对应波前且此处拥挤对繁殖的影响可忽略）得 d/dZ [U, V, W]^T = M [U, V, W]^T，其中 M 为给定 3×3 矩阵。|M − λI| = 0 给出三次方程 λ[Cλ² + (C² + 1 − p)λ + C] = 0，化简为 λ = 0 与 Cλ² + (C² + 1 − p)λ + C = 0。

将 (13.102) 改写为关于 C 的二次方程 λC² + (1 + λ²) C + (1 − p) λ = 0，得 C(λ) = (−(1 + λ²) ± √((1 + λ²)² − 4(1 − p) λ²)) / (2λ)。这是 C 关于 λ 的色散关系。图 13.13 示意画出 C(λ) 的多支结构。

对 (13.105) 两个根关于 λ 的极值可由 (13.104) 求导得 dC/dλ = 0 推出 C² + 2λC + 1 − p = 0；与 (13.104) 联立得极值点出现在 λ = ±1。由此得到正波速 C 的两个可能范围：0 ≤ C ≤ 1 − √p = C₁ 和 C₂ = 1 + √p ≤ C ≤ ∞。与经典 Fisher–Kolmogoroff 分析相比，第一范围不存在。

为确定哪个范围对应实际的行波解，需检查 (0, 0, 0) 附近解的行为（保证 U、V 行为正确即保持非负）。先考虑低范围 (0, C₁) 中的小 C 情形，由 (13.102) 展开得 λ₁ = −C/(1 − p) + O(C³)、λ₂ = −(1 − p)/C + pC/(1 − p) + O(C³)。代入 (13.101) 对应的线性方程组求特征向量分量：由 (1 − p) e₁ + (pC²/(1 − p) + O(C⁴)) e₂ = 0 得 e₁ 与 e₂ 必须反号（p ≤ 1）。这意味着小 C 时 U、V 沿特征向量方向趋于平稳态时符号相反，不构成对 U、V 都有意义的非负解。随 C 在 (0, C₁) 中变化，特征值与特征向量连续变化，但 (13.109) 无非平凡解使 e₁ = 0 或 e₂ = 0，因此由连续性知 (0, C₁) 内不存在生态学上合理的行波解。

对高范围 (1 + √p, ∞) 中的 C 大情形，展开 (13.105) 得 λ₁ = −1/C + O(1/C³)、λ₂ = −C + p/C + O(1/C³)。类似分析表明对应特征向量的 U-、V-分量同号（除去 λ₀ = 0 对应的不可接受特征向量）。因此只有在 C 的高范围中才存在容许的解轨迹。这意味着 C₂ = 1 + √p 是现实解波速的下界，对应于标准 Fisher–Kolmogoroff 方程分析中的最小波速 2√(RD)（量纲形式）。

行波本质上是种群增长波，因此虽然扩散者负责波的传播，波上每一点扩散者与不扩散者都按 (13.94) 增长。沿行波解的扩散者与不扩散者的相对大小可由方程的特殊形式求出。将 (13.99) 与 (13.100)（令 W = U_Z）耦合系统通过引入 Q = U + V、U_Z = P 解耦为：U_Z = P、P_Z = −CP − pQ(1 − Q)、Q_Z = P − (1 − p) Q(1 − Q)/C。在 P − Q 平面中，特征向量论证以 Q 正、P 负于 (0, 0) 处为依据。

P − Q 平面上所需轨迹的形状由 dP/dQ = (−CP − pQ(1 − Q)) / (P − (1 − p) Q(1 − Q)/C) 给出。沿轨迹"向后"走 U 的总变化（即 Z = −∞ 处到 Z = ∞ 处）为 U₀ = ∫₀¹ (P dQ) / (P − (1 − p) Q(1 − Q)/C)。利用 dP/dQ 的表达式化为 P/(P − (1 − p) Q(1 − Q)/C) = p − (1 − p)/C · dP/dQ，代入积分得 U₀ = ∫₀¹ [p − (1 − p)/C · dP/dQ] dQ = p。

因此对任意波速 C，Z = −∞ 处 U 必须为 p，U(∞) = 0。即波后远处扩散者比例为 p，这与直觉一致。进一步由 U(Q) = ∫₀^Q (P dQ')/(P − (1 − p) Q'(1 − Q')/C) = pQ − (1 − p) P/C，得 U/Q = p − (1 − p) P/(C Q)。由于 P 为负，扩散者比例在所有点（除 P = 0 处）都高于 p，且越接近波前 P 越负，扩散者比例越高，符合直觉。

利用 l'Hôpital 法则可进一步考察轨迹趋于 (0, 0, 0) 时的斜率：由 (13.113) 得 (dP/dQ)² + (C − (1 − p)/C) dP/dQ + p = 0。由于 dP/dQ 必须为负，要求 C > √(1 − p)。此条件在 C 的高范围 (C₂, ∞) 中恒满足，低范围 (0, C₁) 中不成立，因此再次确认 (13.107) 给出的 C 范围在容许性上的结论。

Cook 的数值求解 (13.97) 显示解快速收敛到接近预测最小速度的行波解。例如 p = 1.0, 0.5, 0.1, 0.05, 0.01 时理论最小波速分别为 1 + √p = 2.00, 1.70, 1.33, 1.22, 1.10，对应数值波速为 2.01, 1.77, 1.34, 1.22, 1.10。

13.8 物种入侵与范围扩张（Species Invasion and Range Expansion）

物种的空间扩散在生态学上极为重要。Elton（1958）的经典著作列举了大量例子，此后又记录了许多较新的例子。从巴西向美国西南蔓延的杀人蜂是较近期的突出例子，Vancouver 岛南部的美洲牛蛙扩散是更新的例子。Skellam（1951）的开创性论文奠定了理论方法的基础，他实际上使用了 Fisher–Kolmogoroff 方程 (13.4) 的线性化形式（包含扩散与 Malthus 增长即指数增长）。他对麝鼠的范围扩张尤为感兴趣，发现入侵波速近似为 2√(rD)，并证明范围随时间线性扩张。Shigesada 与 Kawasaki（1997）的专著讨论了多种具体入侵（哺乳动物、植物、昆虫、流行病等），给出了若干主要模型，研究了齐次与异质空间环境、捕食-被捕食、竞争等多种种间相互作用下的行波型入侵。

基本上，当个体运动的尺度远小于观察尺度时，连续模型是合理的起点。一个模型与数据良好结合的范例是加州海岸的加州海獭（Enhydra lutris）再入侵。Lubina 与 Levin（1988）使用 Fisher–Kolmogoroff 方程 (13.4) 结合现有数据。海獭种群在 20 世纪初因过度捕猎被认为几近灭绝，1911 年通过国际条约保护时认为已经灭绝。1914 年在 Big Sur 附近发现约 50 头，此后种群数量及其领地均增加。Lubina 与 Levin（1988）充分记录这一再入侵的一个有趣方面：它本质上是沿海岸线的一维现象。他们能估计 (13.4) 中的参数并表明基本波速 2√(rD)（r 为线性增长率、D 为扩散系数）给出极佳结果。常速下范围随时间线性增长（从约 1938 年至 1984 年的再入侵数据基本支持）。北方入侵 D = 13.5 km²/yr、南方 D = 54.7 km²/yr、估计种群增长率 r = 0.056/yr，得到北方和南方波速分别为 1.74 km/yr 和 3.4 km/yr，与 1938–1972 年间的观测值 1.4 km/yr 和 3.1 km/yr 以及 1973–1980 年南方 3.8 km/yr 吻合。他们有说服力地论证北南入侵速度差异并非来自方程中的对流项，而是来自参数的栖息地变化。

回到上节对变扩散模型的结果并置于物种领地入侵背景下：已表明对 (13.97)，为使形如 (13.98) 的行波解生态上合理（即非负），波速 C 必须满足 C ≥ C₂ = 1 + √p，其中 p = r₁/(r₁ + r₂) 是新生个体为扩散者的概率。量纲下 c ≥ c₂ = √(RD) (1 + √p)，其中 c 是量纲下波速，D 是扩散亚种群的扩散率，R 是本征增长率。图 13.14 给出最小波速作为扩散者概率 p 的函数，并与经典 Fisher–Kolmogoroff 结果比较。

当然并未证明该行波解可由某些初值产生（如 Kolmogoroff 等 1937 对 Fisher–Kolmogoroff 方程所做），但作者认为合理初值会收敛到最小波速 (13.123) 的解，正如 Fisher–Kolmogoroff 方程一样。

考虑两个特例。p = 1 时所有个体以相同扩散系数扩散，系统化为 Fisher–Kolmogoroff 方程，波速下界 c₂ = 2√(DR)，与 (13.123) 在 p → 1 时一致。p ≪ 1 时扩散者极少；保持 D 固定，(13.123) 的下界为 c₂ ∼ √(RD) (1 + 0(p^{1/2}))，恰为所有个体以固定 D 扩散时下界的一半。这与直觉相反（图 13.14）：仅含极少量扩散者的种群，其波速与全种群以相同速率扩散的波速相差不超过两倍。自然环境因素可轻易造成此效应。

最初的直觉是若种群中扩散者极少，入侵范围应很小，极限下为零。当扩散者数量极低时，连续扩散假设不再成立，随机效应将占主导。但即使在到达此情形之前，波速仍不会接近零。

Cook 模型及其分析的主要点或许是：仅需少量扩散者就能以与全种群扩散时相近的速度驱动入侵。这显然有重要的生态学含义。如上节指出，入侵波本质上是繁殖波——种群一旦大于零，模型中的繁殖项就开始起作用，既产生扩散者也产生不扩散者。可将快速扩散者视为"播种"了不扩散者的繁殖：它们是繁殖波的驱动力。

Cook 还研究了其模型的其他方面与修改，例如扩散率方差对入侵波的影响、种群增长中的 Allee 效应（最小可生存种群；参见第 1 章讨论）以及两亚种群以不同速率扩散的情形。他进行了广泛的数值模拟以确认解析结果，并将基本思想应用于使用相关随机游动建模运动的其他方程，表明其波速的主要结果不仅限于经典扩散模型。

Lewis 与 Schmitz（1996）的工作与 Cook 的工作直接相关。他们也考虑具有可分流动与静止状态（含状态间切换）的生物入侵问题，对扩散和繁殖都考虑切换，证明即使转移率无穷小也能发生快速入侵。

Shigesada 等（1995）的论文（参见 Shigesada 与 Kawasaki 1997 专著）尤为相关。他们考察了多种物种（如北美的英格兰麻雀、欧洲椋鸟、日本稻水象甲）的范围扩张。为研究范围扩张，他们使用的模型之一是 Skellam（1951）经典工作中提出的两维 Fisher–Kolmogoroff 方程的线性化标量形式（即增长为 Malthus 指数增长）。他们从量纲下方程 ∂u/∂t = ∇²u + εu 出发，u 为局部种群密度、空间径向对称。在 δ 函数初值 u(r, 0) = N₀ δ(r)（代表局部引入物种）下解为 u(r, t) = (N₀/(4πDt)) exp(εt − r²/(4Dt))。

从物种空间扩散的实际角度看，扩张范围取决于少数个体的入侵。Shigesada 等（1995）指出可能存在一个最小可检测密度 u，低于此密度种群实际不可检测。这意味新引入物种在开始扩张其栖息地范围前存在一个事实上的等待期。由解 (13.125)，原点附近对小 t 的种群密度很快降到阈值 u 之下；但由于 (13.124) 中的指数增长项（使解含 εt 因子），u 开始增加并最终穿越阈值 u。这段滞后或建立期是引入种群到其规模穿过阈值水平之间的时间。令 u = u、r = r 即得到 r 关于 t 的关系 (13.126)：r = 2t[εD + (D/t) ln(4πDt u/N₀)]^{1/2}。

引入无量纲量 R = √(ε/D) r、T = εt、γ = εN₀/(Du)，得无量纲范围-时间关系 (13.128)：R = 2T[1 + (1/T) ln(γ/(4πT))]^{1/2}，仅依赖于无量纲参数 γ。当 (1/T) ln(γ/(4πT)) ≪ 1 时，范围随时间线性扩张，R* ≈ 2T。Shigesada 等（1995）进一步发展了具有分散殖民地（由长程扩散者启动）的物种入侵与范围扩张模型；这些模型实质上是具有可变扩散的入侵模型，重要的是他们将解析结果与真实数据关联，得到良好相关。

使用阈值和径向对称线性扩散反应产生入侵前沿的思想由 Murray（1981）应用于完全不同的生物学情形——蝴蝶翅上眼斑的发育。他也将该模型应用于其他非径向对称情形，相关应用详见第一卷第 3 章。

本章个人批注

本章是反应扩散理论中"行波分析"专题的一次系统化纵览，从最经典的 Fisher–Kolmogoroff 方程出发，把单物种行波解的存在性、最小波速、渐近解、稳定性、密度依赖扩散、精确解、多平稳态情形、可激动力学（钙波）以及变异性对入侵波速的影响连贯展开。我特别注意到几个对理解反应扩散行波机制至关重要的"反直觉"现象：(1) Fisher–Kolmogoroff 方程的最小波速由线性化前缘给出 c_min = 2√(kD)，恰好也是紧支集初值演化的渐近速度；紧凑的 2√(kD) 公式同时是色散关系的极值点和 Kolmogoroff–Petrovsky–Piscounov 紧支集定理的产物，这并非巧合——线性稳定性 + 紧支集 + 强非线性 Logistic 项的组合恰好把行波速度锁在最小值上。(2) 13.4 节中精确解（(13.46) 给出 q = 1 时的 c = 5/√6 ≈ 2.04）虽然数学上严格，但定量波形与渐近解差异显著；这警示精确解往往"找到的是某个解而非最相关的解"，必须与相平面分析结合。(3) 13.5 节中波速方向由 ∫ f(u) du 的符号决定，相当于在相平面上比较两个稳定态间"势垒"两侧的面积——这给出了通过调控参数（例如减小 q）改变 f(u) 形状从而遏制虫害扩散的几何解释。(4) 13.6 节中 cot θ 项对球面波速的影响是几何性的：上半球波前周长扩张减慢、下半球周长收窄加速；这解释了模型中"反方向"加速与真实卵中减速的差异需要皮层非均匀性来补偿——这是 Cheer 等 1987 论文的一个关键推论。(5) 13.7 节 Cook 模型的核心反直觉结论是：仅需少量扩散者即可驱动接近全种群扩散时的入侵速度，差距不超过两倍；这与"如果只有几个扩散者，入侵范围应接近零"的直觉相悖；其机制是"快扩散者播种了不扩散者的繁殖"，入侵波本质是繁殖波而非纯粹的扩散波。这是本章中我个人觉得最有生态学洞察力的结果——它把"入侵"从单一种群扩散的标量模型扩展到亚种群结构层面，并展示了看似无关的离散变异性可在波速上产生显著影响。我对其中几处"亚种群相对大小"的解析推导（U₀ = p 的积分恒等式）尤为欣赏，其技巧性虽强但结论简洁。整体上，本章的处理风格是"经典理论 + 精确解 + 真实生物学应用"三位一体——13.1–13.3 是 Fisher–Kolmogoroff 经典理论的精炼版（包含色散关系、渐近解与稳定性的标准三件套），13.4 给出精确解及其陷阱，13.5 把 f(u) 从双稳态扩展到多稳态并给出几何的面积法则，13.6 用钙波把"轴对称几何的 cot θ 效应"与真实发育生物学挂钩，13.7–13.8 跳到生态入侵的实际问题。行文逻辑上从"理论 → 模型 → 应用 → 反直觉推广"是层层递进的，但每节末尾往往不收束于证明而是把读者引向"该方法如何用于更现实的模型"，留出继续研究的余地。我对第 12 章与本章的区别需要再确认——第 12 章讲的是极限环振荡产生的波（Oscillator-Generated Wave Phenomena），本章讲的是单物种反应扩散方程的传播波（Travelling Wavefronts）；两者都属"波"，但驱动力和数学工具截然不同：第 12 章涉及耦合振子的相波（phase wave），本章则纯由反应-扩散耦合驱动。

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

第 12 章讨论了由极限环振荡器产生的波现象与中心模式生成器，研究重点是耦合振子的同步与相波传播——其核心数学对象是相方程与振荡器网络；本章则转入反应扩散理论中的"传播波"（travelling wavefronts）专题，把第 12 章的"波现象"从振子驱动切换到"反应 + 扩散"双驱动：驱动力的本质是稳定/不稳定齐次态之间的非线性切换加上线性（或密度依赖）扩散的协同；本章因此可以看作从第 12 章的"动力学波"走向第 14 章"分形"的过渡——前几章积累的连续介质建模与分析技术在本章的标量反应扩散方程中达到第一个高峰（最小波速公式、色散关系、渐近解、稳定性、精确解），并为第 14 章将要展开的"分形几何"应用提供了生物-物理接口（例如通过波浪式分布解释斑图）。作者也在 13.8 节末尾把目光引向第一卷第 3 章中蝴蝶翅图式发育与第 6 章中力学-化学方法的应用，提示本章建立的行波机制是后续"图式生成"与"生物波在异质介质/异质皮层中传播"两章的重要预修。位置上看，本章是第三部分"波与空间中的图式"的支柱性章节：第 12 章是行波分析的动力学/振子入口，本章是反应扩散行波的完整工具箱，第 14 章则把这些工具进一步推向分形几何与生物形态学。