第 9 章：受扰与耦合振子及黑洞（Perturbed and Coupled Oscillators and Black Holes）

9.1 振子中的相位重置（Phase Resetting in Oscillators）

随着已知生物振子数量的丰富及其普遍重要性的认可，一个自然的问题是外部扰动会对振子后续振荡产生什么影响。Winfree 在 1960 年代研究果蝇 Drosophila melanogaster 蛹周期性羽化时，在生物学语境下提出了这个看似简单却深刻的问题——这是开创性工作。自那时起，外部扰动对振子、空间耦合振子、与扩散过程耦合振子等的惊人性质被陆续发现（第 1 卷第 12 章与第 1 章）。Winfree 发展了一套新的生物时间几何概念理论，提出了许多富有挑战性的数学问题。Winfree (2000) 的开创性著作全面讨论了该领域并给出详尽文献，提供了大量生物实例说明理解这些效应对解释现象的必要性。

心脏的周期起搏点是一个重要的振子并被广泛研究，特别是外加扰动的影响。Jalife 与 Antzelevitch (1979) 处理心脏组织中的起搏点活动，其结果将在 9.4 节讨论；Krinsky (1978) 讨论心脏波型心律失常；Winfree (1983a,b) 讨论突发性心脏死亡的拓扑学方面。也有在跑步者与马的呼吸—跨步神经同步控制方面有趣的工作（参见 Hoppensteadt 1985）。

第 7.5 节表明神经细胞在某些条件下可呈现规则周期性发放。鉴于神经信号传递的关键重要性，研究外部刺激对这种振荡的影响显然极具意义。Best (1979) 的工作与本节及随后三节尤为相关：他将动作电位传播的一个被广泛接受的模型——FitzHugh–Nagumo 模型（第 7.5 节方程 (7.39)）——置于周期性脉冲下，展示了本章讨论的若干重要现象。

心脏起搏点问题及某些类型的心脏衰竭可能与受扰振子相关的波现象有关（例如 Winfree (1983b) 的综述文章）。这一领域与心脏衰竭相关的研究已持续多年，是本节内容的核心动机。然而本节的结果与结论是相当普适的、事实上与模型无关的——尽管为教学便利我们使用具体模型。

作为引言，简要描述关于果蝇蛹近似 24 小时节律性羽化的实验观察。在蛹期发育中经历变态，最终羽化为成虫。若让变态中的蛹在典型的昼夜光暗循环中自由放置，蝇在约 24 小时内成批羽化，每次持续约 6–8 小时。若将这种蛹置于完全黑暗环境中，蝇仍以几乎完全相同的方式羽化；图 9.1 汇总了多次实验的累积结果。

在暗环境中，蛹受到一束短暂光脉冲作用时，蝇的羽化时间或相位发生移动。换言之，潜在的生物钟出现了相位移动。相位移动取决于光脉冲的施予时间 T 及其持续时间即总光通量 D（单位 ergs/cm²）。我们关心光脉冲后羽化时间 T_E；T_E 取决于 T 与 D。若 T 时刻给药 D=0，相位变化显然为零；T_E = 24 − T 小时。Winfree (1975) 给出大量 T、D 变化下 T_E 测定的实验结果。此处需注意的一个重要事实是：存在一个临界剂量 D，若在特定时间 T 施予，将导致不再有周期性羽化，而是连续羽化。换言之周期性行为被破坏。同样令人惊讶的是导致这种结果所需的剂量非常小（参见 Winfree 1975）。这些实验结果基本表明存在一个临界相位和刺激使潜在周期性行为即生物钟被摧毁。这对一般振子具有重要意义。

本节及随后三节主要讨论生物振子、刺激与时间对周期性行为的影响以及实验证据与意义。在果蝇实验中振子的刺激-时间-响应空间存在奇点（或多个奇点），在奇点处振子简单地停止或行为变得不可预测。远离该奇点，行为则相对可预测。9.4 节将描述其他刺激实验——心脏组织上的实验——它们呈现相似的相位奇点行为。

在进行具体分析前——对所选教学例子分析非常简单——先以简单单摆说明相位重置与刺激-时间-相位奇点现象。设单摆以周期 ω 摆动，从摆锤在右侧最高点 S 时刻起度量零相位 t=0。则每当 t = nω（n 整数）摆锤再次到达 S。规则振荡中对摆锤施加冲量显然会扰动周期摆动。冲量或刺激后摆锤最终再次呈现简谐运动，但摆锤到达 S 的时刻不再是 t = nω 而是 t = t_s + nω，其中 t_s 为某常数。换言之相位被重置。若在摆锤正好摆到最低点时施予刺激，刺激量恰当时可使单摆完全停止。即在恰当的相位或时间施加恰好的刺激可完全停止振荡；这就是上面果蝇实验中刺激-相位-响应空间的奇点。

设某振子由向量状态变量 u 描述，满足微分方程组 du/dt = f(u, λ)，其中 f 为非线性速率函数，λ 表示振子参数。为视觉清晰与代数简洁，设 (9.1) 描述一个仅含两物种 x、y 的极限环振子。极限环轨迹通常是物种平面上一条简单闭合曲线 γ（如图 9.2(a)）。通过适当变量变换可将该极限环变换为闭合轨迹是圆的振子，振子状态本质上由角 θ 描述——相位——其原点在圆上某任意点。极限环以速度 v = dθ/dt 遍历，每完整遍历一次 θ 增加 2π。

这样一个极限环系统的简单例子是 dr/dt = R(r)、dθ/dt = Φ(r)，其中 R(r) > 0 对 0 < r < r₀，R(r) < 0 对 r > r₀，R(r₀) = 0，Φ(r₀) = 1。这些条件意味着 (9.2) 有唯一吸引极限环 r = r₀、dθ/dt = 1（第 3 章提到一个特别简单的情形：R(r) = r(1 − r)、Φ(r) = 1，其解可平凡给出）。关于 r₀ 归一化后可取极限环为 r = 1。图 9.2(b) 展示典型相平面极限环解。图 9.2(c) 展示点 x_P = x_S + cos θ(t) 随 t 变化的典型周期行为，其中 x_S 是图 9.2(a) 中的稳态，标出了 θ 的等价值点。

极限环解可被可视化为围绕圆的运动这一事实已被 Winfree (2000) 在环动力学的总题下以直观方式发展。其拓扑方面产生一些出人意料的结果和新概念。这里仅考虑主题的基本元素，但足以展示某些重要概念。

以生理振子建模为背景，设想某事件如心跳发生在某个特定相位值处，可归一化为 θ = 0。起搏点在周期中循环，于该特定相位（即时间）发放，然后在该周期的一部分处于不应期，随后再次发放，如此反复。用振子的环或圆概念，可将起搏点视为以常速沿环运动的点，每当通过相位 θ = 0 的位置即发放一次。

虽然从时间角度 t 线性增加，在特定时间（周期的整数倍）起搏点发放。为理解外部刺激对振子相位重置的基本概念，我们取最简单的非平凡极限环振子系统 dr/dt = r(1 − r)、dθ/dt = 1 作为教学例子，其相位 θ(t) = θ₀ + t（模 2π）；参见图 9.2(b)。基于此讨论两种基本类型的相位重置，即 Type 1 与 Type 0。

Type 1 相位重置曲线（Type 1 Phase Resetting Curves）

先扰动相位使控制方程变为 r = 1、dθ/dt = 1 + v(θ, I)，其中 v(θ, I) 表示对角速度 dθ/dt 施加的扰动即刺激。I 是表示施加在振子上的冲量大小的参数。再取一个 Winfree (1980) 用过的简单但非平凡的 v：dθ/dt = 1 + I cos 2θ，其中 I 可正可负。若 t = 0 时刻施予刺激 I 并持续时间 T，则对 (9.6) 积分给出新相位 φ 关于刺激起始时的旧相位 θ 的关系。由 (9.6)，从 θ 到 φ 积分 (1 + I cos 2s)⁻¹ ds = T，可得解析解：|I| < 1 时 tan φ = A tan[T B + tan⁻¹(A⁻¹ tan θ)]，|I| = 1 时 tan φ = 2T + tan θ，|I| > 1 时 tan φ = A(|K| + 1)/(|K| − 1) 若 |tan θ| > A，tan φ = A(|K| − 1)/(|K| + 1) 若 |tan θ| < A，其中 K = ((A + tan θ)/(A − tan θ)) exp(2TB)，A = ((1 + I)/(1 − I))^(1/2)，B = (|1 − I²|)^(1/2)。这些显式给出新相位 φ 关于旧相位 θ 及刺激强度 I 与持续时间 T 的函数。因此施予刺激会导致振子的相位移动；换言之重置了相位。t > T 时振子恢复 dθ/dt = 1 但现在有相位移动。这意味着振子将在不同时间发放但以与之前相同的相位发放；后续周期当然也与刺激前相同。即图 9.2(c) 那样的周期"波"形仅仅被沿时间轴平移一段。

我们关注不同刺激强度（即 I 及其持续时间 T 依赖）下 φ 关于 θ 的相位重置曲线。关于形如 (9.5) 的刺激，一个重要点是 dφ/dθ > 0 对所有 I、T、θ 成立。这表明在较晚的旧相位施加冲量导致较晚的新相位。这可由对 (9.8) 关于 θ 求导立即得到 dφ/dθ 严格正表达式。画出新相位 φ 关于旧相位 θ 的曲线，即得相位重置曲线，典型形态如图 9.3。无论刺激 I 如何，新相位 φ 的值覆盖完整相位周期，此处为 0 到 2π。换言之任何 0 < φ ≤ 2π 的新相位可通过选取适当的 0 < θ ≤ 2π 旧相位和刺激 I 来获得。对给定的 I 和 T，新相位 φ 由旧相位 θ 唯一确定。这称为 Type 1 相位重置曲线，其特征是 dφ/dθ > 0 对所有 0 < θ ≤ 2π：整个周期的平均梯度为 1——故得此名。虽图 9.3 中 I 在 0 < θ < π 内是相位提前，另一个振子可能呈现相位延迟。重点是 Type 1 重置中 dφ/dθ > 0。

若 I 足够强（即 |I| > 1），相速度 dθ/dt = 1 + v(θ, I) 在某些相位可能为负：对于 (9.6)，这发生在 1 + |I| cos 2πθ < 0 的 θ。这意味着刺激期间存在相位吸引子与相位排斥子，即 dθ/dt = 0 且 d[dθ/dt]/dθ 分别负和正处（参见第 1 章单种群模型的稳定性分析）。刺激不会持续无穷长时间，因此刺激移除后振子恢复周期循环——当然相位由 (9.8) 确定而不同。

Type 0 相位重置曲线（Type 0 Phase Resetting Curves）

考虑同一极限环 (9.4) 但现在施加使解离开极限环 r = 1 的刺激 I。具体地，取 I 为沿 y 轴方向的冲量（如图 9.4）。分析对任何扰动均成立但代数更复杂且易于掩盖主要观点。规定 I > 0 为图 9.4 所示情形：即 0 < θ < π/2 时新相位 φ 小于旧相位 θ 且新位置一般有 r = ρ ≠ 1。我们想求新相位 φ 关于旧相位 θ 和刺激 I。由图，ρ cos φ = cos θ、ρ sin φ + I = sin θ，消去 ρ 即得 φ = φ(θ, I) 的隐式表示；这是 (φ, θ, I) 空间中的三维曲面。如后所见，将该曲面投影到 (I, θ) 平面特别值得关注。在此之前先构造等价于图 9.3 的相位重置曲线，即各种刺激 I 下新相位 φ 关于旧相位 θ：这些是曲面 φ = φ(θ, I) 在固定 I 下向 (φ, θ) 平面的投影。

由 (9.10)，tan φ = tan θ − I/cos θ，它给出给定 I 下 φ 关于 θ。先设 0 < I < 1。由图 9.4 定性可知对 0 < θ < π/2 与 3π/2 < θ < 2π，φ < θ；对 π/2 < θ < 3π/2，φ > θ。因此定性相位重置曲线如图 9.5(a)：它与零刺激对角线在 θ = π/2、3π/2 相交。定量细节此处不重要。由 (9.11) 对 θ 求导得 (1 + tan²φ) dφ/dθ = 1 + tan²θ − I sin θ/cos²θ = 1 − I sin θ/cos²θ。容易验证对 0 < I < 1，dφ/dθ > 0 对所有 θ。因此在相位重置曲线上，0 < I < 1 时 dφ/dθ > 0 对所有 θ 如图 9.5(a)。与图 9.3 的曲线比较，它们拓扑等价，故图 9.5(a) 是 Type 1 相位重置曲线。−1 < I < 0 时同样成立。

现考虑 I > 1。由 (9.12)，存在 dφ/dθ < 0 的 θ 区间。参见图 9.4(b) 并让 P 沿圆移动。S 永不进入上半平面。即当 θ 在整个 2π 周期内变化时，φ 至少不会取 (0, π) 范围内的任何相位；事实上确切范围可由 (9.11) 或 (9.12) 容易算出。这种情况的相位重置曲线定性如图 9.5(b)。该曲线与图 9.5(a) 的曲线拓扑不等价。所有 I > 1 的相位重置曲线都与 Type 1 重置曲线拓扑不同。形如图 9.5(b) 的相位重置曲线——即旧相位 θ 取遍 (0, 2π) 内所有相位值时，新相位 φ 仅取完整周期范围的一个子集——称为 Type 0 重置曲线。这种曲线上 dφ/dθ < 0 对某 θ 范围：曲线上平均梯度为 0，这正解释了该类重置曲线之得名。另注意 Type 0 重置曲线不能仅由相位刺激得到。刺激 I < −1 时也得到同样的 Type 0 重置曲线。

另一种清晰展示 Type 0 重置曲线到 Type 1 重置曲线在 I 通过 I = 1 时分岔的方式是画出 (9.11) 给出的新相位 φ 关于旧相位 θ。图 9.6 给出代表性 0 < I < 1、I = 1、I > 1 的例子。

尽快重置生物钟是患时差综合征时每个人想做的事。Winfree 不容置疑地表明了昼夜节律对光的极高敏感性。他引入了研究内时钟的全新方法。关于生物钟重置，人与果蝇本质上无别，恰当时间给予恰当光刺激可重置人的生物钟。从 Winfree (1975) 的文章可确定何时施加强阳光以重置跨多个时区飞行后的生物钟：例如从西雅图到伦敦有 8 小时时差，约下午 1 点用 15 分钟强烈阳光照射眼睛应可重置生物钟；问题是在英国不易得到那阳光。反向则在西雅图下午 5 点需要强烈阳光，显然也有同样问题。Winfree (1982, 2000) 讨论了人体生物钟及睡眠时间安排，提示理解这些节律可能有实用医学和精神病学意义。

9.3 黑洞（Black Holes）

上一节分析表明，刺激 I 由 0 增大时，在 I 通过 I = 1 时相位重置类型出现显著分岔。即对 I = 1 相位重置存在奇点。为清晰看出物理上发生什么，必须将 (9.11) 给出的 φ = φ(θ, I) 曲面投影到 (I, θ) 平面，对各种 φ。即构造曲线 I = (sin θ − cos θ tan φ)，对 0 ≤ φ ≤ 2π 内各种 φ。这虽是初等曲线绘制的练习，需要相当仔细的微积分运算，结果示意如图 9.7。先考虑旧相位范围 0 ≤ θ ≤ π，设此时 φ ≠ π/2、3π/2。无论 φ 取何值，所有曲线均过点 I = 1、θ = π/2，因该处 cos θ tan φ = 0 且 I = sin(π/2) = 1 对所有 φ。所有 π/2 > φ > 0、2π > φ > 3π/2 的曲线在 θ = 0 轴上相交于 I = −tan φ。对 π/2 < φ < 3π/2，由 (9.13) 的微积分分析得所示曲线。特殊值 φ = π/2、3π/2 给出经过 θ = π/2 的竖直奇点线（取极限 φ → π/2 或观察 φ 趋近 π/2 时的常 φ 相曲线）。处理完 θ 范围 (0, π) 后类似处理 (π, 2π)，整体图像如图 9.7。需注意的重要事实是：存在两个奇点 S1 与 S2，进入每个奇点的曲线包含 (0, 2π) 内每个相位的常相位曲线；一种情况下 φ 递增的曲线按顺时针排列，另一种按逆时针。

现考虑图 9.7 的重要含义。设拥有这样一个振子并在给定相位 θ 施予刺激 I。只要 |I| < 1，可由 I 和旧相位 θ 读出新相位，且结果唯一。对所有 |I| > 1，给定旧相位 θ，新相位也唯一确定。但在这种情况下，给定 I 可对应两个不同旧相位 θ 的同一新相位 φ。前者（参见图 9.5）是 Type 1 相位重置，后者是 Type 0 相位重置。

现取特定刺激 I = 1 在相位 θ = π/2 施加于振子；图 9.7 中结果点为奇点 S1，它没有与某一特定相位 φ 对应，而是对应整个范围 0 ≤ φ ≤ 2π。换言之该特定刺激在特定相位下的效果不确定。这些奇点 S1 与 S2 即刺激-相位空间中的黑洞，是刺激结果未知的点。若 I 非恰为 1 而接近 1，则所有相位 φ 通过该奇点，结果显然微妙。从实用观点看这种刺激对生物振子的效果不可预测。然而从数学上看，若精确刺激 I = 1 严格在 θ = π/2 施加，则没有结果新相位 φ。这正是在简单单摆情形下，当摆锤正好经过竖直位置时施以恰好算出的冲量会发生的。在实践中让真实单摆完全停止显然很困难，即使能接近数学上算出的条件，所得相位结果也远非显然。

上述概念（Winfree 1970；亦见 2000）显然适用于任何内生振子，故结果与含义相当普适。可呈现 Type 1 与 Type 0 相位重置的生物振子的一个关键特征是：在其旧相位-刺激空间中存在某冲量与相位对应黑洞。或许最重要的应用是：对这样的振子存在一种刺激，若在特定相位施加，将彻底湮灭振荡。黑洞存在性的连续性论证是：若刺激连续增大时在某特定值出现 Type 1 到 Type 0 重置的转换，则在转换的相位与刺激值处存在黑洞。

现考虑实际振子中黑洞与湮灭的实验证据。

9.4 实际生物振子中的黑洞（Black Holes in Real Biological Oscillators）

现已有多份记录详实的实验案例支持 Type 0 相位重置与适当刺激在正确相位下对基本振荡的湮灭——均如上述理论预测。除本节讨论的案例外，还有 Taddei-Ferretti 与 Cordella (1976) 在 Hydra attenuata 测得的 Type 0 相位响应曲线；Pinsker (1977) 关于 Aplysia 突发神经元在突触输入扰动下的工作——同样是 Type 0 情形；Guttman 等 (1980) 显示了乌贼轴突膜神经元振子中的湮灭。在详细描述一个实验案例前，先给出 Best (1979) 对第 7 章 7.5 节 Hodgkin–Huxley 模型中黑洞存在的直接验证，该模型建模了乌贼巨轴突空间钳位膜的振荡。

Hodgkin 与 Huxley (1952) 对乌贼轴突空间钳位神经元发放的模型（方程 (7.37)、(7.38)）呈现极限环振荡。第 7.5 节已说明该模型的 FitzHugh–Nagumo 简化模型有极限环周期行为。Best (1979) 对完整 Hodgkin–Huxley 模型 (7.37)、(7.38) 数值研究在极限环振荡出现的电流范围下进行。然后他对振子施加电压变化（作为刺激），旨在实验实现其结果。他发现了如所预期的 Type 1 与 Type 0 相位重置曲线；图 9.8 给出每种的一个例子。

Best (1979) 模拟中黑洞或零空间的存在由 Type 1 到 Type 0 重置曲线的转变（随电压刺激增大）所指示。这可从图 9.8 预期：9.8(a) 与 (b) 拓扑不同因而由某分岔状态分隔。由于任何数值模拟固有的近似性，不能像图 9.7 那样确定单个奇点。取而代之，存在奇点周围一个区域——黑洞或零空间——其中在适当旧相位范围内施加适当扰动后新相位不确定。图 9.9 示意 Best (1979) 的结果。除阴影区域外，对给定的旧相位 θ 与刺激 I 存在唯一重置相位 φ；注意存在一个 (I, θ) 子空间，可对单个 I 与两个 θ 得同一 φ。还应注意新相位值在一个完整周期中以顺时针方式围绕 Hole 1 变化，以逆时针方式围绕 Hole 2 变化。关于刺激-旧相位等高线图（如 9.9）需记住的一个关键特征是：等高线收敛到黑洞，一个为正刺激、一个为负刺激。

黑洞的另一个关键性质是：若内生振子在适当相位受到临界刺激，振荡直接消失。Best (1979) 在 Hodgkin–Huxley 模型上演示了这一点；结果如图 9.10。注意内生振荡的湮灭。Guttman 等 (1980) 在实验上表明：浸于弱钙溶液的空间钳位轴突中的重复发放可被在周期中特定时间施加的恰当强度刺激停止。

Jalife 与 Antzelevitch (1979) 在心脏起搏点细胞的规则周期跳动上做了类似工作——这与心脏起搏点直接相关。他们使用狗、猫、小牛的心脏组织并对基本振荡施加电刺激。实验得到呈现 Type 1 与 Type 0 重置曲线的相位重置曲线；图 9.11 给出部分结果。

由图 9.11 的重置曲线可预期存在刺激持续时间的转换值或多个值，能摧毁振荡——即内生心脏振子的零空间或黑洞。事实确是如此，结果如图 9.12(a)。图 9.12(b) 给出刺激接近转换值（介于 9.11(a) 与 (b) 值之间）的重置曲线。图 9.12(c) 展示狗心脏纤维中刺激对规则振荡的摧毁。

A. T. Winfree 多年来研究突发性心脏死亡的可能原因及其与起搏点振子拓扑（时间与空间）的关系。虽然心脏收缩涉及起搏点，但偏离正常常涉及循环收缩波的出现而非发放机制的中断。在纤颤情况下，当心律失常使心脏看起来像一把蠕动的蠕虫时，可能对奇点或黑洞出现的透彻理解有助于阐明该问题。Winfree (1983b) 的《科学美国人》文章专门关注突发性心脏衰竭的拓扑学。第 1 卷第 1 章讨论与此类心脏问题直接相关的螺旋旋转波。空间分布的振子在空间非齐次外加刺激下相关的数学问题显然富有挑战性且引人入胜，并具有重大生物学意义。

9.5 耦合振子：动机与模型系统（Coupled Oscillators: Motivation and Model System）

生态学、流行病学、发育生物学等领域中生物振子与周期过程的出现是公认的事实。在大量情形中振子必然以某种方式耦合以获得所需输出。刚刚看到理解振子扰动影响的重要性。这里我们考虑振子耦合的若干效应并描述用于研究此类问题的关键分析技术之一。耦合极限环振子被数学家广泛研究多年，分析问题远非平凡。毫不奇怪它们共同表现的现象范围远大于任何单个振子（参见 Winfree 2000）。该课题是日益增长的研究方向，不仅在生物学中，也在非线性动力学总题下。已观察到的许多过程至今仅被部分理解。本章其余部分主要关注同步过程及其破坏。这些同步现象可能是相位锁定、频率协调等，均源自极限环振子的交互耦合。这里我们限于两个振子的耦合并仅考虑弱耦合；基本上沿用 Neu (1979) 的分析。第 12 章将考虑当模型化某些游泳脊椎动物的神经排列时一链耦合振子相关的重要现象。

考虑数学问题前，先简述我们研究的具体模型系统的一个实验动机。Marek 与 Stuchl (1975) 研究了用不同参数（即不同周期振荡）耦合两个 Belousov–Zhabotinskii 反应系统的效应。他们将每个反应置于独立搅拌反应器中，通过公共多孔壁间的物质交换耦合。他们观察到：若自主振子频率几乎相同，则相位差随时间趋于常值——这称为相位锁定。但若自主频率差过大，则相位锁定不持续而代之以耦合系统有长时间缓慢变化的相位差区间，被非常短时间内的快速涨落分隔。下面给出的分析将解释这些现象。

耦合极限环振子的分析研究中，并不需要详细知道所建模的具体系统。鉴于上述实验，我们以 Belousov 反应系统为例。设极限环振子相同且各自单独由方程组 dx_i/dt = F(x_i, y_i)、dy_i/dt = G(x_i, y_i)（i = 1, 2）控制，其中非线性函数 F、G 表示振子动力学（它们可以例如为 (8.27) 右端函数，第 8.4、8.5 节讨论的 Belousov 反应两反应物模型之一，或第 3 章 3.3 节 (3.18) 的捕食者-猎物交互动力学）。假设 (9.14) 的解呈现周期 T 的稳定极限环行为：x_i = X(t + ψ_i)、y_i = Y(t + ψ_i)（i = 1, 2），其中 ψ_i 为任意常数。故 X(t + ψ_i + T) = X(t + ψ_i)、Y(t + ψ_i + T) = Y(t + ψ_i)（i = 1, 2）。

为以方便（如后所见）但仍普适的方式形式化弱耦合，考虑无量纲模型系统： dx₁/dt = F(x₁, y₁) + ε{k(x₂ − x₁) + λf(x₁, y₁)}, dy₁/dt = G(x₁, y₁) + ε{k(y₂ − y₁) + λg(x₁, y₁)}, dx₂/dt = F(x₂, y₂) + εk(x₁ − x₂), dy₂/dt = G(x₂, y₂) + εk(y₁ − y₂), 其中 0 < ε ≪ 1，k > 0 为耦合常数。ε = 0 时振子解耦，方程简化为 (9.14)。形式 (9.16) 的普适性来自 λ 项。若 ε ≠ 0 且 λ = 0，两振子相同且有如 (9.15) 的解耦解。若 ε ≠ 0 且 λ ≠ 0，两不同振子被耦合，(9.16) 前两方程中的 ελ 项是独立极限环振子的一部分。我们所选的具体耦合（由 k 项代表）正比于差 x₁ − x₂ 与 y₁ − y₂。对 Marek 与 Stuchl (1975) 的实验，这反映了存在质量传递。对相互作用种群，可视为物种质量传递、扩散通量近似。事实上当考虑栖息地间影响对种群动力学影响时，常以这种方式纳入：它将总体空间效应纳入而不需要扩散项，否则将模型变为偏微分方程系统；这些后面考虑。

9.6 振荡相位锁定：萤火虫同步（Phase Locking of Oscillations: Synchronisation in Fireflies）

生物振子耦合时可产生极为丰富的现象：节律分裂、相位锁定、牵入等。耦合振子的数学具有挑战性且可能非常复杂。相关文献众多，从非常抽象到非常实用（如跑步者同步与人体睡眠-觉醒周期 Strogatz 1986）。Winfree (1987) 的生物钟优美著作（科学和视觉上）讨论了生物钟在众多领域中的应用，重点在昼夜节律与相位重置；亦见 Winfree (1987) 的 When Time Breaks Down 主要讨论心脏节律。Glass 与 Mackey (1988) 的入门书给出大量与生物钟相关的节律现象例子；应用主要在生理学。Strogatz 与 Stewart (1993) 的非技术性《科学美国人》文章描述多种有趣的同步现象，包括萤火虫同步闪烁的惊人照片；这是下面讨论的话题。

Glass 与同事（Guevara 等 1981、Guevara 与 Glass 1982、Keener 与 Glass 1984）研究了一个数学与生物学紧密相关的有趣实用例子。该模型与下面讨论的密切相关，但产生第 2 章讨论的非线性差分方程。设 φ_i 为 delta 函数刺激施加于系统前的相位——即旧相位。新相位（这里在离散模型中记为 φ_{i+1}）由 φ_{i+1} = g(φ_i) + 2πτ (mod 2π) 给出，其中 g 是 φ_i 的函数（后来由特定系统的实验确定），τ 是相对于周期长度的归一化刺激周期。该方程的稳定稳态对应相位锁定；初步分析类似第 2 章。Guevera 与 Glass (1982) 用该方程研究受周期性 delta 函数刺激的基本方程 (9.4) 的牵入。在参数范围内他们也发现了混沌行为。即使如此简单的方程也能产生大量复杂与意外的解。该工作特别有趣的是：在分析之后 Guevara 等 (1981) 实际在自发跳动胚胎心脏细胞制备中测量了重置曲线。在宽范围的振幅与频率下，他们能预测观察到的节律，包括相位锁定与混沌节律。这是数学建模与分析如何促进我们理解重要复杂生物现象的优秀例子。

生物医学科学中耦合振子的建模挑战似乎无穷无尽。这里我们仅触及主题。9.7–9.9 节将讨论两个振子的弱耦合并用奇异摄动技术分析。另一个涉及许多振子的例子将在第 12 章讨论。

一个突出且有名的耦合生物振子视觉例子是大量静止萤火虫（Pteroptyx malaccae）的周期闪光相位锁定同步。雄萤发光以吸引雌萤，雌萤飞行寻找闪光特别吸引的雄萤。许多实验研究量化了单只萤火虫（及其他昆虫）变化光闪周期性的能力。Buck 与 Buck (1976) 的《科学美国人》文章是萤火虫同步主题的优秀入门。Buck (1988) 综述这些萤火虫同步节律闪光的生物学文献。Hanson (1978) 的早期工作表明：单只萤火虫能改变其光发射振子的相位（由内生神经起搏点控制），并能牵入或同步于发光光——只要其周期在萤火虫自然周期约 0.9 秒附近。若人造光刺激周期离自然周期过远则不能牵入。一些萤火虫在牵入上更擅长；Pteroptyx malaccae 似乎是大师，能变化频率近 15%。已有多个萤火虫同步的数学模型：Rinzel 与 Ermentrout (1983)、Ermentrout (1991)、Mirollo 与 Strogatz (1990)。这里我们仅讨论 Rinzel 与 Ermentrout (1983) 的简单但有效的模型。

如 9.1 节，设萤火虫振子在时刻 t 的相位为 θ(t)，其自然频率为 ω。即无任何外部刺激时相位满足 dθ/dt = ω。具体地，假设振子在 θ = 0 发放即发光。设外部相位为 θ_e(t)，频率为 ω_e，故满足 dθ_e/dt = ω_e。萤火虫试图将其频率与外部刺激同步，过慢则加速，过快则减速。实现此功能的简单模型是 dθ/dt = ω + I sin(θ_e − θ)，其中参数 I > 0。刺激大小 I 是萤火虫改变频率效率的度量。若 θ_e 领先 θ（0 < ω_e − ω < π），则 θ' > ω，萤火虫加速其相位。若 θ' < ω，萤火虫减速。形式 (9.19) 是相位重置方程 (9.5) 的特殊情形。一种基于 θ − θ_e(t) 函数但更复杂情形的相似假设用于第 12 章 12.3 节。

当关心何时发生同步时（如第 12 章 12.3 节所示），考虑相位差 φ = θ − θ_e 的方程。由 (9.18) 与 (9.19) 得 dφ/dt = dθ_e/dt − dθ/dt = ω_e − ω − I sin φ，φ(t) = θ_e(t) − θ(t)。引入新变量 τ = I t、δ = (ω_e − ω)/I，方程变为 φ' = dφ/dτ = δ − sin φ。无量纲参数 δ 有明确物理意义：是外部频率与自然频率差相对于刺激强度 I 的度量。如 9.2、9.3 节所见，I 的大小至关重要。

我们关心 (9.22) 的稳态解及其稳定性。若有稳定稳态解 φ_s > 0，则由 (9.20) 外部刺激相位 θ_e 总领先萤火虫相位 θ 常量。萤火虫振子因此与刺激相位锁定但在刺激后即刻闪光。若 δ = 0，φ = 0 是 (9.22) 的解，此时若 φ = 0 稳定，两振子零相位差闪烁，故为齐声。牵入问题归结为 (9.22) 的稳态及其稳定性。

第 1 章已看到所需做的只是画出 (9.22) 右端图、读出稳态并注意梯度正负。图 9.13 示意 δ ≥ 0 的主要解可能性。稳态稳定性由稳态处梯度决定——梯度负则稳定，正则不稳定。对该模型，萤火虫与刺激相位锁定时仅有一个稳定稳态，条件 −1 < δ < 1。−1 < δ < 0 的情形与图 9.13 相似但 δ < 0。

图 9.13(b) 中稳定稳态为 0 < φ₁ < π，故萤火虫必须增加频率以相位锁定。若 δ > δ_c = 1 如图 9.13(c)，则不能跟上，相位差 φ 持续增加直至达到 2π 周期重新开始。后一种情形为相位漂移。9.13(c) 中 φ' > 0 且非常数，意味着相位漂移持续增大但速率不均匀。这与 Hanson (1978) 的实验结果一致。

用此模型可作多个预测：δ > 1 时的非均匀相位漂移仅为一个。关键预测是：刺激的相位锁定可能当且仅当外部频率 ω_e 满足 ω − I ≤ ω_e ≤ ω + I，该不等式给出牵入的刺激频率范围。同样，刺激强度 I 重要。若由实验知道刺激频率范围，可计算 I 并由 (9.20) 与 (9.22) 预测相位锁定相位差 φ_s = θ_e − θ_s = sin⁻¹((ω_e − ω)/I)，−π/2 ≤ φ_s ≤ π/2。

当 −1 ≤ δ ≤ 1 时，相位锁定萤火虫振子的有量纲周期 T 由 (9.22) 中 φ 变化 2π 所需时间得到：T = (1/I) ∫₀^{2π} dφ/(δ − sin φ)，化简得牵入周期 T = 2π / [I (δ² − 1)^(1/2)] = 2π / [(ω_e − ω)² − I²]^(1/2)。δ → ±1 时周期变为无穷大；即无牵入，正如图 9.13(c) 所示。现在可见萤火虫如何能同步其周期发光脉冲：例如一只有较强 I 的萤火虫牵入另一只时，群刺激以单频增长直至全部被牵入。较强 I 意味着较小 δ，故"起搏点"与附近频率者频率差小。若整群现以同频闪烁，对绕圈雌萤要选定领头者必较困难！

萤火虫 Pteroptyx malaccae 不是唯一的萤火虫种类，但它似乎在能与外部刺激相位锁定上最为灵活。看到简单模型如何能捕获一些实验结果后，我们应考察能反映更多生物学的模型。适应由正弦函数（如 (9.19)）支配的假设过于简单。第 12 章在处理 N 个振子时再次用此假设。一个更合适的耦合方程代替 (9.19) 应为 dφ/dt = ω + h(φ)，其中 h(φ) 是其参数的周期函数但不必对称。Ermentrout (1991) 给出了结合更多 Pteroptyx malaccae 适应特征并用更一般耦合函数的更精细模型。他还数值模拟了一群萤火虫在萤火虫树中如何趋近同步。上述模型展示了同步但除 δ = 0 外有永久相位滞后。这是其作为萤火虫树模型的缺点——实际中萤火虫树有几乎无相位滞后的同步。Ermentrout (1991) 也讨论了此方面并提出一种可能机制。

9.7 奇异摄动分析：预备变换（Singular Perturbation Analysis: Preliminary Transformation）

(9.16) 一般而言难以分析。即使数值上也不易看出解行为对各参数的依赖，特别是在非相同自主振子情形。由于许多感兴趣情形中耦合是弱的，且在许多生物学应用中预期如此，我们利用 0 < ε ≪ 1 这一事实并使用奇异摄动理论（参见 Murray 1984 关于基本技术的简要教学性讨论）。

每个振子有其自己的极限环解，可由 x − y 相平面中的闭合轨迹 γ 表示。我们可以新坐标系用此曲线作为局部坐标系基础。可以周期 T 的相位 θ 表征周期极限环（从 0 变到 T 一周），用垂直于 γ 的距离 A 度量对 γ 的扰动；在 γ 上 A = 0。代数上证明用此特征在我们的耦合振子分析中尤为方便。故代替 (9.15) 作为自主极限环解有 x_i = X(θ_i)、y_i = Y(θ_i)（i = 1, 2），其中 X(θ_i) 与 Y(θ_i) 是 θ_i 的 T 周期函数。注意 θ_i 与 t 关系为 dθ_i/dt = 1。

用相位与垂直于极限环的扰动表示接受周期极限环解的相平面系统解的思路可由以下（虽人造但仍具教益）例子说明。考虑微分方程组 dx₁/dt = x₁(1 − r) − ωy₁、dy₁/dt = y₁(1 − r) + ωx₁，r = (x₁² + y₁²)^(1/2)，其中 ω 为正常数。相平面分析（见附录 A）显示 (0, 0) 是唯一奇点且为不稳定螺旋，逆时针旋转。可找到一个约束集（取 r 大并注意在这样的大圆上轨迹向量 (dx₁/dt, dy₁/dt) 指向内），故由 Poincaré–Bendixson 定理存在极限环周期解，表示为 (x₁, y₁) 平面中的闭合轨道 γ。改用极坐标 (r, θ)，x₁ = r cos θ、y₁ = r sin θ，系统 (9.28) 变为 dr/dt = r(1 − r)、dθ/dt = ω。极限环轨迹 γ 即 r = 1。解如图 9.14(a)。极限环渐近稳定：由 (9.30)，对 r = 1 的任何扰动都将衰减，r 简单地卷回 r = 1，因 dθ/dt > 0 故为逆时针。在该例中若对极限环的扰动为 r < 1，由 (9.30) r 增大；若扰动为 r > 1，r 减小并趋于轨道 r = 1。此处 r = 1 等价于轨道 γ，A（即垂直距离）即 r − 1。由 (9.30) 关于 A (= r − 1) 与 θ 的微分方程组为 dA/dt = −A(1 + A)、dθ/dt = ω。当然可精确积分 (9.30) 得 r(t) = r₀eᵗ / [(1 − r₀) + r₀eᵗ]、θ(t) = ωt + θ₀，其中 r(0) = r₀、θ(0) = θ₀，并由此 x₁(t) = r(t) cos θ(t)、y₁(t) = r(t) sin θ(t)。t → ∞ 时 r(t) → 1（故 A(t) → 0）且 x₁ → cos θ、y₁ → sin θ，等价于 (9.27) 中的 X(θ)、Y(θ)：它们是 θ 的 2π 周期函数。此处遍历 γ（即 r = 1）的速率 dθ/dt = ω 由 (9.31) 给出。图 9.14(b) 示意一般情形。θ = 0 取在某点 P，相位在沿 γ 逆时针方向遍历一周时增加 2π。

现考虑两个振子，每个有其自主闭合极限环轨道 γ_i（i = 1, 2）。耦合的效应将改变每个的轨道与相位。我们可以以局部坐标项来表征每个的效应：用参数化 γ_i 上点的相位 θ_i 和垂直于原极限环轨道的扰动 A_i。回想耦合振子系统 (9.16) 中我们关心弱耦合故 0 < ε ≪ 1。ε = 0 时每个振子有其极限环解，按相位可写为 (9.27)：x_i = X(θ_i)、y_i = Y(θ_i)（i = 1, 2）。θ_i = t + ψ_i ⇒ dθ_i/dt = 1。

我们预期 O(ε) 耦合的效应是使 (9.27) 给出的轨道 γ_i 产生 O(ε) 位移。因此可见合适的变量变换是从 (x_i, y_i) 到局部变量 A_i 与新相位 θ_i（i = 1, 2），其中 A_i 是垂直于轨道 γ_i 的距离。为激励我们将用的具体变量变换，参见图 9.15。无耦合时，轨迹 γ 以平行于 γ 的速度 (dx/dt, dy/dt) 遍历。以相位 θ 表示（相位随轨道遍历单调增加），由 (9.34) 与 (9.35) 该速度等于 (X'(θ), Y'(θ))，其中撇号表示关于 θ 的导数。该速度向量因耦合受扰，轨道 γ 将被位移。该位移在每点可由它被位移的垂直距离描述，记为图中向量 A。由于 A 是速度 (X'(θ), Y'(θ), 0) 与垂直于 (x, y) 平面单位向量 (0, 0, 1) 的向量积，故得 A = A(X'(θ), Y'(θ), 0) × (0, 0, 1) = (A Y'(θ), −A X'(θ), 0)。

现考虑 0 < ε ≪ 1 的 (9.16) 系统，假设自主轨道被 O(ε) 扰动。那么从 (x_i, y_i) 到 (A_i, θ_i) 的合适变量变换用 (9.34) 与 (9.36) 为 x_i = X(θ_i) + εA_i Y'(θ_i)、y_i = Y(θ_i) − εA_i X'(θ_i)（i = 1, 2）。此处用 εA 代替 A 以强调 ε 在我们分析中是小量，故轨道扰动也是小量。

9.8 奇异摄动分析：变换后系统（Singular Perturbation Analysis: Transformed System）

现将变量变换 (9.37) 用于耦合系统 (9.16)（0 < ε ≪ 1）。即将 (9.37) 用于右端并对 ε 作 Taylor 展开：代数复杂繁琐，但简洁有趣且重要的最终结果值得，不仅为本章展示的结果，还为第 12、13 章将讨论的两个显著现象。我们进行足够代数以展示如何将 (9.16) 表示为关于变量 θ_i 与 A_i 的方程；不过要跟随代数细节建议单独用纸笔（希望跳过代数的读者可跳到 (9.45) 但后续会引用此处的某些定义）。

下面各函数（主要是 X、Y）的参数除另行说明或强调外均为 θ₁。(9.16) 第一式用 (9.37) 变为： dx₁/dt = X' dθ₁/dt + εY' dA₁/dt + εA₁ Y'' dθ₁/dt = F(X, Y) + εA₁[Y' F_X(X, Y) − X' F_Y(X, Y)] + εk[X(θ₂) − X(θ₁)] + ελf(X, Y) + ε²k[A₂ Y'(θ₂) − A₁ Y'(θ₁)] + ε²λA₁[Y' f_X(X, Y) − X' f_Y(X, Y)] + O(ε³),

第二式变为： dy₁/dt = Y' dθ₁/dt − εX' dA₁/dt − εA₁ X'' dθ₁/dt = G(X, Y) + εA₁[Y' G_X(X, Y) − X' G_Y(X, Y)] + εk[Y(θ₂) − Y(θ₁)] + ελg(X, Y) + ε²k[A₁ X'(θ₁) − A₂ X'(θ₂)] + ε²λA₁[Y' g_X(X, Y) − X' g_Y(X, Y)] + O(ε³)。

ε = 0 时由 (9.14) 与 (9.34) 有 X'(θ₁) = F(X, Y)、Y'(θ₁) = G(X, Y)。现将 (9.38) 乘以 X'(θ₁) 并加 Y'(θ₁) 乘 (9.39)，得 (X'² + Y'²) dθ₁/dt + εA₁(X'Y'' − Y'X'') dθ₁/dt = [X' F(X, Y) + Y' G(X, Y)] + εA₁{X'Y'[F_X(X, Y) − G_Y(X, Y)] − X'² F_Y(X, Y) + Y'² G_X(X, Y)} + εk{X'[X(θ₂) − X(θ₁)] + Y'[Y(θ₂) − Y(θ₁)]} + ε²kA₂[Y'(θ₂)X'(θ₁) − X'(θ₂)Y'(θ₁)] + ελ[X' f(X, Y) + Y' g(X, Y)] + ε²λA₁{X'Y'[f_X(X, Y) − g_Y(X, Y)] − X'² f_Y(X, Y) + Y'² g_X(X, Y)} + O(ε³)。

由 (9.40) X' F(X, Y) = X'²、Y' G(X, Y) = Y'²，故上式变为 R²(1 + εΓ A₁) dθ₁/dt = R² + ε[R²Ω A₁ + R²k r + R²k V + R²γ λ] + O(ε²)，其中 R² = X'² + Y'² ≠ 0，R²Γ = X'Y'' − Y'X''，R²γ = X' f(X, Y) + Y' g(X, Y)，R²r = −XX' − YY'，R²Ω = X'Y'[F_X(X, Y) − G_Y(X, Y)] − X'² F_Y(X, Y) + Y'² G_X(X, Y)，R²V = X'(θ₁)X(θ₂) + Y'(θ₁)Y(θ₂)。

两边除以 R²(1 + εΓ A₁) 并对 0 < ε ≪ 1 作级数展开得 dθ₁/dt = 1 + ε[{Ω(θ₁) − Γ(θ₁)} A₁ + λγ(θ₁) + kr(θ₁) + kV(θ₁, θ₂)] + O(ε²)。

类似地，对 (9.38) 乘以 Y'(θ₁) 减去 (9.39) 乘以 X'(θ₁)，得 A₁ 方程。dA₁/dt = Φ(θ₁) A₁ + kU(θ₁, θ₂) + λφ(θ₁) + εΨ(A, θ) + O(ε²)，其中 R²U(θ₁, θ₂) = X(θ₂)Y'(θ₁) − Y(θ₂)X'(θ₁)，Φ、φ、Ψ 均确定：A 与 θ 是向量 (A₁, A₂) 与 (θ₁, θ₂)。我们唯一需要确切形式的函数是 (9.44) 与 (9.41) 给出的 U(θ₁, θ₂) 与 V(θ₁, θ₂)。

现对 (9.16) 第三、第四方程做同样处理，发现变换到 (A_i, θ_i) 依赖变量的效果是用以下系统替代模型耦合振子系统 (9.16)： dA₁/dt = Φ(θ₁) A₁ + kU(θ₁, θ₂) + λφ(θ₁) + εΨ₁(A, θ) + O(ε²), dθ₁/dt = 1 + ε[{Ω(θ₁) − Γ(θ₁)} A₁ + λγ(θ₁) + kr(θ₁) + kV(θ₁, θ₂)] + O(ε²), dA₂/dt = Φ(θ₂) A₂ + kU(θ₂, θ₁) + εΨ₂(A, θ) + O(ε²), dθ₂/dt = 1 + ε[{Ω(θ₂) − Γ(θ₂)} A₂ + kr(θ₂) + kV(θ₂, θ₁)] + O(ε²)。

函数 V 与 U 由 (9.41) 与 (9.44) 给出，后续将引用。函数 Φ、φ、Γ、γ、Ω、Ψ₁、Ψ₂、r 的确切形式对后续分析非必要，重要的是它们都是 θ₁、θ₂ 的 T 周期函数且 Φ 满足关系 ∫₀^T Φ(σ) dσ < 0。最后这个关系来自原始极限环解（解耦振子的）是稳定的这一事实；这里稍作题外证明。

解耦振子的极限环稳定性条件：当 k、λ 为零时振子解耦。我们要保持 ε ≠ 0 因为我们准备用变换 (9.37) 研究受扰极限环振子。关于变量 A 与 θ，由 (9.45) 与 (9.46)，k = λ = 0 时的控制系统为 dA/dt = Φ(θ) A + O(ε)、dθ/dt = 1 + O(ε)，其中 Φ 是 θ 的 T 周期函数。对 O(1) 时间第二方程给出 θ ≈ t，第一方程变为 dA/dt = Φ(t) A + O(ε)，从 t 到 t + T 积分得 A(t + T)/A(t) = [1 + O(ε)] exp(∫_t^{t+T} Φ(σ) dσ) = [1 + O(ε)] exp(∫₀^T Φ(σ) dσ)，其中积分限已因 Φ 为 T 周期函数而改变。未受扰极限环为 A ≡ 0、θ = t + ψ。故极限环稳定当且仅当 (9.48) 的所有解有 A(t) → 0 当 t → ∞。由 (9.49) 可见若 ∫₀^T Φ(σ) dσ < 0 则 A(t + T) < A(t) 对所有 t 成立，故 A(t) → 0 当 t → ∞。

9.9 奇异摄动分析：双时间展开（Singular Perturbation Analysis: Two-Time Expansion）

(9.45) 与 (9.46) 中的 O(ε) 项在长时间（实际上 O(1/ε)）后才有影响。这提示将 (9.45) 与 (9.46) 中 A_i 与 θ_i 寻找如下形式的渐近解（ε → 0）：A_i ∼ ⁰A_i + ε ¹A_i、θ_i ∼ ⁰θ_i + ε ¹θ_i，其中 A 与 θ 是快时间 t 与慢时间 τ = εt 的函数。换言之仅当 τ = O(1) 时（即 t = O(1/ε)）ε 效应才显现。（参见 Murray 1984 关于双时间展开过程的初等阐述。）现在所有时间导数 d/dt = ∂/∂t + (dτ/dt) ∂/∂τ = ∂/∂t + ε ∂/∂τ，(9.45) 与 (9.46) 系统变为偏微分方程系统。

本节其余代数也相当复杂。最终结果即 (9.69) 在 9.10 节是必需的，那里导出耦合振子的一个重要结果。现将 (9.50) 与 (9.51) 代入 (9.45) 与 (9.46) 并令 ε 各幂次相等得方程层级： O(1)： ∂⁰A₁/∂t − Φ(⁰θ₁) ⁰A₁ = kU(⁰θ₁, ⁰θ₂) + λφ(⁰θ₁), ∂⁰θ₁/∂t = 1, ∂⁰A₂/∂t − Φ(⁰θ₂) ⁰A₂ = kU(⁰θ₂, ⁰θ₁), ∂⁰θ₂/∂t = 1。

O(ε)： ∂¹A₁/∂t − Φ(⁰θ₁) ¹A₁ = {Φ'(⁰θ₁) ⁰A₁ + k ∂U/∂θ₁(⁰θ₁, ⁰θ₂) ¹θ₁ + k ∂U/∂θ₂(⁰θ₁, ⁰θ₂) ¹θ₂ + λφ'(⁰θ₁) ¹θ₁} + Ψ₁(⁰A, ⁰θ) − ∂⁰A₁/∂τ, ∂¹θ₁/∂t = [Ω(⁰θ₁) − Γ(⁰θ₁)] ⁰A₁ + λγ(⁰θ₁) + kr(⁰θ₁) + kV(⁰θ₁, ⁰θ₂) − ∂⁰θ₁/∂τ, ∂¹A₂/∂t − Φ(⁰θ₂) ¹A₂ = {Φ'(⁰θ₂) ⁰A₂ + k ∂U/∂θ₂(⁰θ₂, ⁰θ₁) ¹θ₂ + k ∂U/∂θ₁(⁰θ₂, ⁰θ₁) ¹θ₁} + Ψ₂(⁰A, ⁰θ) − ∂⁰A₂/∂τ, ∂¹θ₂/∂t = [Ω(⁰θ₂) − Γ(⁰θ₂)] ⁰A₂ + kr(⁰θ₂) + kV(⁰θ₂, ⁰θ₁) − ∂⁰θ₂/∂τ。

积分 (9.52) 第 2、第 4 式得 ⁰θ_i = t + ψ_i(τ)（i = 1, 2），其中 ψ_i(τ) 此处是 τ 的任意函数。代入 (9.52) 第 1、第 3 式得 ∂⁰A₁/∂t − Φ(t + ψ₁) ⁰A₁ = kU(t + ψ₁, t + ψ₂) + λφ(t + ψ₁), ∂⁰A₂/∂t − Φ(t + ψ₂) ⁰A₂ = kU(t + ψ₂, t + ψ₁)。

为得所需解，再考虑较简洁的方程：dx/ds − Φ(s) x = f(s)、dy/ds − Φ(s) y = U(s, s + χ)，其中 χ ≡ ψ₂ − ψ₁。记住 Φ、φ、U 都是 T 周期函数。各方程的补函数或齐次解为 exp[v(s)]，v(s) = ∫₀^s Φ(σ) dσ。由 (9.47) 知 v(s) < 0，这对解耦极限环振子稳定性必要。故补函数最终衰减到零。

可证 (9.56) 各方程有唯一 T 周期解。考虑如 (9.56) 第一式。其精确解为 x(s) = x(0) exp(∫₀^s Φ(σ) dσ) + ∫₀^s exp(∫_α^s Φ(σ) dσ) φ(α) dα。x 方程在 s 替换为 s + T 时不变。由于上述 x(s) 解是周期的，x(0) = x(T) = x(0) exp(∫₀^T Φ(σ) dσ) + ∫₀^T exp(∫_α^T Φ(σ) dσ) φ(α) dα，这是初值 x(0) 的方程，代入 (9.57) 得 (9.56) 第一式的唯一 T 周期解。类似地 (9.56) 第二式有唯一 T 周期解。记这些周期解为 x = p(s)、y = ρ(s, χ)。

以这些解 (9.58) 表示，(9.55) 的一般解为 ⁰A₁ = kρ(t + ψ₁, χ) + λp(t + ψ₁) + h₁(τ) exp[v(t + ψ₁)]，⁰A₂ = kρ(t + ψ₂, −χ) + h₂(τ) exp[v(t + ψ₂)]，其中 h₁、h₂ 是 τ 的任意函数。代入 (9.53) 的 θ 方程得 ∂¹θ₁/∂t = {kρ(t + ψ₁, χ) + λp(t + ψ₁)}{Ω(t + ψ₁) − Γ(t + ψ₁)} + λγ(t + ψ₁) + kr(t + ψ₁) + kV(t + ψ₁, t + ψ₂) − dψ₁/dτ + h₁(τ) exp[v(t + ψ₁)]{Ω(t + ψ₁) − Γ(t + ψ₁)}, ∂¹θ₂/∂t = kρ(t + ψ₂, −χ){Ω(t + ψ₂) − Γ(t + ψ₂)} + kr(t + ψ₂) + kV(t + ψ₂, t + ψ₁) − dψ₂/dτ + h₂(τ) exp[v(t + ψ₂)]{Ω(t + ψ₂) − Γ(t + ψ₂)}。

若 f(t) 是 T 周期函数可写 f(t) = μ + ω(t)，μ = (1/T) ∫₀^T f(s) ds，其中 ω(t) T 周期且均值为零。用此于 (9.60) 得 ∂¹θ₁/∂t = μ₁(χ) + ω₁(t, τ)、∂¹θ₂/∂t = μ₂(χ) + ω₂(t, τ)，其中 μ₁ = H(χ) + λβ − dψ₁/dτ、μ₂ = H(−χ) − dψ₂/dτ，且 β = (1/T) ∫₀^T {p(s)[Ω(s) − Γ(s)] + γ(s)} ds、H(χ) = (1/T) ∫₀^T {kρ(s, χ)[Ω(s) − Γ(s)] + kr(s) + kV(s, s + χ)} ds。

函数 ω₁(t, τ)、ω₂(t, τ) 由指数衰减项与均值为零的周期项组成，故积分 (9.62) 得 ¹θ₁ = μ₁(χ) t + W₁(t, τ)、¹θ₂ = μ₂(χ) t + W₂(t, τ)，其中 W₁、W₂ 有界。代入 (9.53) 第三式并用 (9.54) 与 (9.57)，得 ¹A₂ 方程 ∂¹A₂/∂t − Φ(t + ψ₂) ¹A₂ = {S₁(t, τ) μ₁ + S₂(t, τ) μ₂} t + B(t, τ)，S₁ = k ∂U/∂θ₁(t + ψ₂, t + ψ₁)，S₂ = k ∂U/∂θ₂(t + ψ₂, t + ψ₁) + kρ(t + ψ₂, −χ) Φ'(t + ψ₂)，其中 B(t, τ) 是另一函数由指数衰减项与 T 周期部分组成。

按通常的渐近奇异摄动方法，我们现在要求振幅 O(ε) 部分（即 ¹A₂）及其时间导数 ∂¹A₂/∂t 对所有时间有界——这保证级数解 (9.50) 的一致有效性。由 (9.66) 这要求 {S₁(t, τ) μ₁ + S₂(t, τ) μ₂} t 对所有时间有界。然而由于 S₁ 与 S₂ 是 T 周期函数（故不能 t → ∞ 时趋于零），得到有界性的唯一方式是 S₁(t, τ) μ₁ + S₂(t, τ) μ₂ ≡ 0。一般 S₁ 与 S₂ 是两个不同周期函数，故确保 (9.67) 的唯一方式是 μ₁ 与 μ₂ 都为零。即由 (9.63) 要求 μ₁ = H(χ) + λβ − dψ₁/dτ = 0、μ₂ = H(−χ) − dψ₂/dτ = 0，故 dψ₁/dτ = H(χ) + λβ、dψ₂/dτ = H(−χ)。记住 χ = ψ₂ − ψ₁，两式相减得 χ 的单个方程 dχ/dτ = P(χ) − λβ，其中 P(χ) = H(−χ) − H(χ)、χ = ψ₂ − ψ₁。

因变量 χ 即耦合引起的相移：常微分方程 (9.69) 控制该相移的时间演化。推导该方程是 9.8、9.9 节奇异摄动分析的主要目的。它也是本章及第 12 章都要利用的方程。

9.10 相移方程分析及应用：耦合 Belousov–Zhabotinskii 反应（Analysis of the Phase Shift Equation and Application to Coupled Belousov–Zhabotinskii Reactions）

由 (9.64) 的定义，函数 H(χ)、H(−χ) 是 T 周期函数。故相移方程 (9.69) 中的 P(χ) 也是 T 周期函数。若 χ = 0，P(0) = H(0) − H(0) = 0。由 (9.64) 中 H(χ) 的形式，其在 χ = 0 的导数为 H'(0) = (1/T) ∫₀^T {k ∂ρ/∂χ(s, 0)[Ω(s) − Γ(s)] + ∂V/∂θ₂(s, s)} ds。

函数 ρ(s, χ) 是 (9.56) 第二式的 T 周期解：∂ρ/∂s − Φ(s) ρ = U(s, s + χ)。对 χ 求导并在 χ = 0 处取值，得 ∂/∂s (∂ρ/∂χ)(s, 0) − Φ(s) (∂ρ/∂χ)(s, 0) = (∂U/∂θ₂)(s, s)。由 (9.44) 中 U 的定义可见 ∂U/∂θ₂ = 0，故该方程的唯一周期解为 ∂ρ/∂χ = 0。用 (9.41) 中 V 的定义得 ∂V/∂θ₂ = 1。取这些值 (9.70) 给出 H'(0) = (1/T) ∫₀^T k ds = k。故 P'(0) = −H'(0) − H'(0) = −2k。图 9.16(a) 示意典型 P(χ)。

现考虑相移 χ = ψ₂ − ψ₁ 的时间演化方程 (9.69)。相位锁定是相位差对所有时间为常数，即 χ = χ₀ 常数。相位锁定解由 (9.69) 中令 dχ/dτ = 0 给出，即 χ₀ 的方程 P(χ) − λβ = 0。

χ₀ 的线性稳定性由对 (9.69) 关于 χ₀ 线性化给出 d(χ − χ₀)/dτ ≈ P'(χ₀)(χ − χ₀)，故 χ₀ 稳定若 P'(χ₀) < 0，不稳定若 P'(χ₀) > 0。例如若耦合振子相同则 λ = 0（由 (9.16)），由图 9.16(a) χ = 0 是 P(χ) = 0 的解且其导数 P'(0) < 0。故 χ = 0 稳定。这意味着耦合同步相同振子。

现设耦合振子不相同，即 (9.16) 中 λ ≠ 0、ε ≠ 0。这种情况下 (9.71) 的稳态 χ₀ 取决于水平线 λβ 是否与图 9.16(a) 中 P(χ) 曲线相交。故若 min P(χ) = P(χ_m) < λβ < P(χ_M) = max P(χ)，0 ≤ χ₀ ≤ T 内至少有两个稳态解。参见图 9.16(b)，两个典型解 χ₁、χ₂ 由 (9.72) 分别不稳定和稳定，因由观察 P'(χ₁) > 0、P'(χ₂) < 0。这种情形下耦合振子系统经长时间（τ 大即 εt 大）后将演化到稳定极限环振荡，两振子间有常相移 χ₂。至于相移 χ 是否有两个以上稳态解取决于 P(χ) 形式；图 9.16 例子是说明最简单情形。

只要 λβ 落在 P(χ) 最大与最小之间（即 (9.73) 满足），(9.71) 的稳态解 χ₀ 连续依赖于 λβ。这是显然的：若在上下界 P(χ_m)、P(χ_M) 间连续移动 P = λβ 线，由图 9.16(b) 可见。例如 λβ 增大时振子间稳定稳态相移减小。

现设 λβ 使两解合并（于 P(χ_m) 或 P(χ_M)）并问当 λβ 使这些解不再存在时会发生什么。为研究此情形具体考虑 λβ 略小于临界值 (λβ)_c = P(χ_m) 的情形。λβ 略小于 P(χ_m) 时 P(χ) 与直线 λβ 不相交。设 λβ = (λβ)_c − δ² = P(χ_m) − δ²，0 < δ² ≪ 1。现将相位差方程 (9.69) 写为 dχ/dτ = {P(χ) − P(χ_m)} + {P(χ_m) − λβ} = {P(χ) − P(χ_m)} + δ²。

图 9.17(a) 示意 dχ/dτ 关于 χ 的函数（δ² > 0 时）。这与图 9.16(a) 的曲线相同，仅上移距离 −P(χ_m) + δ²，在我们情形下足以使 (9.75) 对 χ 无稳态解；即 dχ/dτ 关于 χ 的曲线不穿过 dχ/dτ = 0 轴。

(9.75) 提出的求解问题在 0 < δ ≪ 1 时是另一奇异摄动问题，可用标准奇异摄动技术处理（参见 Murray 1984）。然而不需要进行渐近分析即可看出解行为。先注意 δ = 0 时（即 λβ = P(χ_m)）有解 χ = χ_m + nT（n = 0, 1, 2, ...）因 P(χ) 是 T 周期函数。具体地，设从图 9.17(a) 中 χ ≈ χ_m 开始。由 (9.75) 及图 9.17，dχ/dτ = O(δ²) > 0，这意味着 χ ≈ χ_m + δ² τ，故对 τ = O(1) 时间 χ 相对 χ_m 变化不大。然而对所有 τ > 0，dχ/dτ > 0，故 χ 缓慢增大。当 τ 足够大，使 χ 显著偏离 χ_m 时，P(χ) − P(χ_m) 不再近似为零，此时 dχ/dτ = O(1) 且 χ 在 τ = O(1) 时间内可测变化，χ → χ_m + T。当 χ 接近 χ_m + T 时，P(χ) − P(χ_m) ≈ 0 又成立，dχ/dτ = O(δ²)。

定性图像现已清楚。解在解 χ = χ_m + nT 附近停留长时间 τ = O(1/δ²)，然后在 τ = O(1) 时间变化到下一解（δ = 0 时），再次停留长时间。过程以准周期方式重复，与我们之前得到的 T 周期极限环行为显著不同。解如图 9.17(b)。快速变化区域是奇异区，而大致常值区域是解的非奇异部分。该行为称为节律分裂。故当 λβ → (λβ)_c = P(χ_m) 时 χ 的解从相位同步分岔为节律分裂。

现已知 χ = ψ₂ − ψ₁ 的定性行为，可由 (9.68) 求 ψ₁、ψ₂：dψ₁/dτ = H(χ) + λβ、dψ₂/dτ = H(−χ)。x_i、y_i 的解由 x_i = X(t + ψ_i(τ))、y_i = Y(t + ψ_i(τ))（i = 1, 2）给出。振荡频率由 O(ε) 给出 dθ_i/dt = 1 + ε dψ_i/dτ = {1 + ε[H(χ) + λβ], i = 1; 1 + εH(−χ), i = 2}。

图 9.18(a) 示意节律分裂解。小尖峰是 χ 在 O(1) 时间的快速变化，长平台区对应缓慢变化解（相位差 ψ₂ − ψ₁ (= χ) 近似为常）。图 9.18(b) 示意该节律分裂下典型解 x_i 关于 t 的函数。注意当相位差 χ 从一个解 χ_m 经快速变化到下一 χ_m + T 时发生频率的突然改变。

耦合振子所展示的相位锁定到节律分裂的突然转换分岔现象正是 Marek 与 Stuchl (1975) 在实验中演示的（如 9.5 节所述）。即他们先观察到相位锁定。然后当参数改变使自主极限环频率足够不同时，他们观察到相位差长时间缓慢变化但被短时间的快速涨落所分隔。这正是上面讨论的节律分裂现象。

本章个人批注

9.1 节里最吸引我的是 Winfree 实验的核心地位：果蝇蛹羽化节律的相位-剂量响应曲面存在奇点——在某个临界剂量 D 与临界时间 T 下，原本稳定的 24 小时周期羽化被"摧毁"，转为连续羽化。这种"黑点"现象被作者用单摆类比讲得很清楚：摆锤过最低点时施加恰当冲量可完全停止摆动。这是一个物理学直觉的完美示范——把抽象的"奇点"落到简单机械装置上。对我来说这是阅读时最有画面感的一节。

9.2 节的 Type 1 与 Type 0 重置曲线分类对我而言是本章最重要的概念性发现。关键是 dφ/dθ 的正负：Type 1 全周期 dφ/dθ > 0（平均梯度 1），Type 0 在某些 θ 范围 dφ/dθ < 0（平均梯度 0）。两类曲线拓扑不等价——这是从 9.2 节 → 9.3 节 → 9.4 节的逻辑链条。9.3 节把 (9.11) 投影到 (I, θ) 平面是这一章的"几何高峰"：所有 φ 值的常相位曲线汇聚到两个奇点 S1、S2，进入 S1 的曲线按顺时针排列、进入 S2 的按逆时针。这种"所有相位曲线汇聚于一点"是真正的奇异性。

9.4 节把上述理论落到 Hodgkin–Huxley 与心脏组织实验上。Best (1979) 的乌贼巨轴突数值工作确认了 Type 1→Type 0 的转变与"零空间"区域。Jalife 与 Antzelevitch (1979) 在狗、猫、小牛心脏起搏细胞上做了类似实验，得到与理论一致的实验观察。Winfree 关于突发性心脏死亡的拓扑学观点把抽象的黑洞与临床重要现象联系了起来——这是从纯数学到生物医学应用的精彩过渡。

9.5–9.6 节将焦点从单振子扰动转到耦合振子。Marek 与 Stuchl (1975) 的 Belousov 双反应器实验是核心动机，萤火虫 Pteroptyx malaccae 同步闪烁是令人印象深刻的视觉例子。萤火虫模型 dθ/dt = ω + I sin(θ_e − θ) 简洁漂亮：相位差方程 φ' = δ − sin φ 的稳态分析给出相位锁定的明确条件 −1 < δ < 1 与失稳（相位漂移）的临界 δ_c = 1。这一简单模型却能定量预测实际萤火虫群的同步行为——这是数学建模力量的精彩展示。

9.7–9.10 节是本章的数学重头戏。变量变换 (9.37) 的动机是几何的：在每个振子极限环 γ_i 上引入局部坐标系（相位 θ_i + 垂直扰动 A_i）。奇异摄动的合理性来自 0 < ε ≪ 1：耦合是弱耦合，对极限环的扰动为 O(ε)。两时间展开（快时间 t + 慢时间 τ = εt）的标准技术是 1970 年代以来奇异摄动分析的核心套路；推导 (9.69) 这样的相移方程是该技术的高级应用。

9.10 节把 (9.69) 的相移方程转化为实际物理：λβ 控制振子失配程度。λβ 处于 P(χ) 最大与最小之间时存在稳定相锁定解；λβ 越界则相锁定被破坏、出现节律分裂。这正好解释了 Marek–Stuchl 实验中观察到的"长时间慢变 + 短时间快变"模式。

Murray 整章的写法有一种"先简单再复杂再回到具体生物" 的结构感：从 9.1 节的果蝇实验 → 9.2 节的简化极限环模型 → 9.3 节的黑洞几何 → 9.4 节的真实生物振子实验验证 → 9.5–9.6 节的耦合振子实例（Belousov 双反应器、萤火虫）→ 9.7–9.10 节的耦合振子一般分析。这种"教学-几何-实验-理论"的螺旋递进是教科书最好的样子。

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

第 8 章是"具体化学振子的详细分析"——以 BZ 反应和 FKN 模型为载体，把第 7 章的极限环/约束集工具搬到一个具体三阶化学反应上。8.1 节引入历史与定性图像，8.2–8.3 节做线性稳定性与全局有界性证明，8.4–8.5 节做弛豫振子的奇异摄动分析。第 9 章紧随其后，把视线从单一振子拓展到"受扰动"与"耦合"的振子系统：9.1–9.2 节建立相位重置的 Type 1/Type 0 分类，9.3–9.4 节把刺激-相位空间中的"黑洞"奇点落到 H-H 模型与心脏起搏细胞实验上，9.5–9.6 节引入耦合振子的具体例子（Belousov 双反应器、萤火虫同步），9.7–9.10 节用奇异摄动技术推导出耦合振子相位差的控制方程并应用到 Belousov 双反应器的"节律分裂"分岔。第 9 章在数学上比第 8 章更"几何"（黑洞、奇点、节律分裂都是几何语言），在生物学上比第 8 章更"耦合"（从单振子到多振子网络）。第 12 章将沿用本章 9.7–9.10 节的奇异摄动技术处理一链 N 个耦合振子（神经排列中的同步），形成跨章的工具复用。