第八章：Belousov-Zhabotinskii反应

学习笔记

8.1 Belousov反应与Field-Körös-Noyes（FKN）模型

8.1.1 历史背景

Belousov-Zhabotinskii反应（以下简称BZ反应）是一个重要的振荡反应，由俄罗斯生物化学家Boris Belousov于1951年首次发现。Belousov在研究铈离子催化柠檬酸被溴酸氧化的过程中，观察到了铈离子浓度比例的周期性振荡现象。当使用铁离子（ferroin）作为催化剂时，振荡表现为鲜明的颜色变化：Fe²⁺状态时为砖红色，Fe³⁺状态时为亮蓝色。这一发现最初被学术期刊编辑拒绝发表，因为当时普遍认为振荡反应根本不可能存在。直到1959年，Belousov才在一份俄罗斯医学会议的论文集中发表了一篇简短注释。后来Zhabotinskii（1964）继续研究这一反应，使其为西方世界所知。Belousov的工作在他去世后于1980年被追授列宁奖。

8.1.2 BZ反应的生物学意义

BZ反应与真实生物振荡器具有重要的类比关系。虽然这种类比不是在分子层面上的，但在数学方程形式上，BZ反应的近似方程与细胞聚集过程中周期性信号传导的模型完全相同（例如盘基网柄菌Dictyostelium discoideum）。Winfree（1987）还将其巧妙地应用于心室三维活动的建模，Tyson（1991）将其用于细胞周期的分子生物学研究。因此，BZ反应不仅是一个化学反应，更是研究生物振荡器的原型系统。

8.1.3 反应的基本机制

原始Belousov反应的基本机制包括：在酸性介质中，以铈（Ce³⁺/Ce⁴⁺）为催化剂，用溴酸离子BrO₃⁻氧化丙二酸。反应可分为两个主要过程，由溴离子浓度[Br⁻]决定哪个过程占主导：

过程I（[Br⁻]高时）：Br⁻被消耗，铈离子主要处于Ce³⁺状态。
过程II（[Br⁻]低时）：Ce³⁺转变为Ce⁴⁺，同时Ce⁴⁺反应重新生成Br⁻，并恢复Ce³⁺状态。

整个序列不断重复，从而产生观察到的振荡。

8.1.4 FKN模型的关键化学反应

BZ反应涉及众多化学反应，但可以合理地简化为5个关键反应，这就是著名的Field-Körös-Noyes（FKN）模型。模型中的化学物质定义为：

\[ X = HBrO_2, \quad Y = Br^-, \quad Z = Ce^{4+}, \quad A = BrO_3^-, \quad P = HOBr \qquad(8.1) \]

五个关键反应为：

\[ \boxed{ \begin{aligned} \text{反应1:} \quad & A + Y \xrightarrow{k_1} X + P \\ \text{反应2:} \quad & X + Y \xrightarrow{k_2} 2P \\ \text{反应3:} \quad & A + X \xrightarrow{k_3} 2X + 2Z \\ \text{反应4:} \quad & 2X \xrightarrow{k_4} A + P \\ \text{反应5:} \quad & Z \xrightarrow{k_5} fY \end{aligned} } \qquad(8.2) \]

其中前两个反应近似对应过程I，后三个反应近似对应过程II。f是化学计量因子，通常取0.5。

8.1.5 三组分动力学方程

根据质量作用定律，从反应(8.2)可得三维动力学方程组：

\[ \boxed{ \begin{aligned} \frac{dx}{dt} &= k_1 ay - k_2 xy + k_3 ax - k_4 x^2 \\ \frac{dy}{dt} &= -k_1 ay - k_2 xy + fk_5 z \\ \frac{dz}{dt} &= 2k_3 ax - k_5 z \end{aligned} } \qquad(8.3) \]

这个振荡系统有时被称为"Oregonator"模型，因为Field等人的研究是在俄勒冈大学进行的。

8.1.6 无量纲化

对方程组(8.3)进行无量纲化。参照Tyson（1985）的做法，引入无量纲变量：

\[ x^* = \frac{x}{x_0}, \quad y^* = \frac{y}{y_0}, \quad z^* = \frac{z}{z_0}, \quad t^* = \frac{t}{t_0} \qquad(8.4) \]

其中特征尺度为：

\[ \boxed{ \begin{aligned} x_0 &= \frac{k_3 a}{k_4} \approx 1.2 \times 10^{-7} \, M, \quad y_0 = \frac{k_3 a}{k_2} \approx 6 \times 10^{-7} \, M \\ z_0 &= \frac{2(k_3 a)^2}{k_4 k_5} \approx 5 \times 10^{-3} \, M, \quad t_0 = \frac{1}{k_5} \approx 50\, s \end{aligned} } \]

无量纲参数为：

\[ \boxed{ \begin{aligned} \varepsilon &= \frac{k_5}{k_3 a} \approx 5 \times 10^{-5}, \quad \delta = \frac{k_4 k_5}{k_2 k_3 a} \approx 2 \times 10^{-4} \\ q &= \frac{k_1 k_4}{k_2 k_3} \approx 8 \times 10^{-4}, \quad f \approx 0.5 \end{aligned} } \qquad(8.4) \]

代入后得到简洁的无量纲系统（省略星号）：

\[ \boxed{ \begin{aligned} \varepsilon \frac{dx}{dt} &= qy - xy + x(1 - x) \\ \delta \frac{dy}{dt} &= -qy - xy + 2fz \\ \frac{dz}{dt} &= x - z \end{aligned} } \qquad(8.5) \]

8.2 FKN模型的线性稳定性分析与极限环解的存在性

8.2.1 定常态求解

令方程组(8.5)左边为零，可得非负定常态：

\[ \boxed{(0, 0, 0) \quad \text{或} \quad z_s = x_s, \quad y_s = \frac{2fx_s}{q + x_s}} \qquad(8.7) \]

其中

\[ x_s = \frac{1}{2}\left[(1 - 2f - q) + \sqrt{(1 - 2f - q)^2 + 4q(1 + 2f)}\right] \qquad(8.7) \]

另一个非零定常态为负值，故只考虑上述正解。

8.2.2 原点的稳定性

在原点(0, 0, 0)处线性化，稳定性矩阵A的特征值λ满足：

\[ |A - \lambda I| = 0 \Rightarrow \lambda^3 + \lambda^2\left(1 + \frac{q}{\delta} - \frac{1}{\varepsilon}\right) - \lambda\left[\frac{1}{\varepsilon}\left(1 + \frac{q}{\delta}\right) - \frac{q}{\delta}\right] - \frac{q(1 + 2f)}{\varepsilon\delta} = 0 \qquad(8.8) \]

特征根的乘积为 \(\frac{q(1 + 2f)}{\varepsilon\delta} > 0\)，因此至少有一个正实部根，原点始终线性不稳定。

8.2.3 正定常态的稳定性与Hopf分岔

在正定常态\((x_s, y_s, z_s)\)处线性化，特征方程为：

\[ \lambda^3 + A\lambda^2 + B\lambda + C = 0 \qquad(8.8) \]

其中

\[ \boxed{ \begin{aligned} A &= 1 + \frac{q + x_s}{\delta} + \frac{E}{\varepsilon} \\ E &= 2x_s + y_s - 1 = \frac{x_s^2 + q(x_s + 2f)}{q + x_s} > 0 \\ B &= \frac{q + x_s}{\delta} + \frac{E}{\varepsilon} + \frac{(q + x_s)E + y_s(q - x_s)}{\varepsilon\delta} \\ C &= \frac{(q + x_s)E - 2f(q - x_s) + y_s(q - x_s)}{\varepsilon\delta} = \frac{x_s^2 + q(2f + 1)}{\varepsilon\delta} > 0 \end{aligned} } \qquad(8.9) \]

根据Routh-Hurwitz条件，所有特征根具有负实部的必要充分条件是：

\[ AB - C > 0 \qquad(8.10) \]

定义 \(\phi(\delta, f, \varepsilon) = AB - C = \frac{N\delta^2 + M\delta + L}{\varepsilon\delta^2}\)，其中N > 0。稳定性边界由 \(\phi = 0\) 给出。

8.2.4 Hopf分岔条件

当参数满足 \(B = C/A\) 时，发生Hopf分岔。临界情况下特征根为：

\[ \lambda = -A, \quad \pm i\sqrt{B} \quad \text{当} \quad B = \frac{C}{A} \qquad(8.13) \]

分岔附近的小振幅极限环周期为：

\[ \boxed{T = \frac{2\pi}{\sqrt{C/A}} = 2\pi\sqrt{\frac{A}{C}}} \qquad(8.14) \]

8.2.5 稳定性边界与参数域

由 \(\phi(\delta, f, \varepsilon) = 0\) 得二次方程：

\[ N\delta^2 + M\delta + L = 0 \qquad(8.10) \]

对于正定常态的不稳定（即产生振荡），参数必须满足：

\[ 0 < \delta < \frac{1}{2N}\left[-M + \sqrt{M^2 - 4LN}\right] \qquad(8.15) \]

这要求L < 0，由此可得f的临界值。当 \(q = 8 \times 10^{-4}\) 时，不稳定参数域为：

\[ \boxed{\frac{1}{4} \approx \,_1f_c < f < \,_2f_c \approx \frac{1 + \sqrt{2}}{2}} \qquad(8.17) \]

8.3 FKN模型的非局部稳定性

8.3.1 confined set（有界集）的概念

线性稳定性分析只给出局部信息。要证明极限环的存在，需要找到满足以下条件的有界闭曲面S（confined set）：

\[ \boxed{\mathbf{n} \cdot \frac{d\mathbf{r}}{dt} < 0, \quad \mathbf{r} \in S} \qquad(8.19) \]

其中n是S的外法向量。这意味着一旦解进入S，就永远不会离开。

8.3.2 有界矩体的构造

Hastings和Murray（1975）给出了严格证明。取一个包围定常态\((x_s, y_s, z_s)\)的矩形箱体，其六个面为：

\[ x = x_1, \, x = x_2; \quad y = y_1, \, y = y_2; \quad z = z_1, \, z = z_2 \qquad(8.20) \]

通过分析每个面上的向量场，可确定有界集边界：

\[ \boxed{ \begin{aligned} x &= q, \quad x = 1 \\ y &= \frac{2fq}{q + 1}, \quad y = \frac{f}{q} \\ z &= q, \quad z = 1 \end{aligned} } \qquad(8.20) \]

在典型参数值下，该有界集确实存在，且定常态\((x_s, y_s, z_s)\)位于其内部。结合局部不稳定性和有界集的存在，运用Poincaré-Bendixson定理的推广，可以证明极限环解的存在。

8.4 松弛振荡器：BZ反应的近似分析

8.4.1 松弛振荡器的特征

观察BZ反应的浓度-时间曲线，某些部分变化非常迅速（特别是[Br⁻]的突升和突降），这是松弛振荡器的典型特征。松弛振荡器是指极限环的某些部分被快速跨越的振荡器，其数学特征是微分方程组中存在小参数在关键位置。

8.4.2 简单松弛振荡器模型

考虑简单的松弛振荡器系统：

\[ \boxed{ \begin{aligned} \varepsilon \frac{dx}{dt} &= y - f(x) \\ \frac{dy}{dt} &= -x, \quad 0 < \varepsilon \ll 1 \end{aligned} } \qquad(8.21) \]

当 \(f(x) = \frac{1}{3}x^3 - x\) 时，这就是著名的Van der Pol振荡器。

8.4.3 周期计算

沿y = f(x)的nullcline（分支AB和CD），慢变过程满足：

\[ f'(x)\frac{dx}{dt} \approx -x \qquad(8.22) \]

可从此积分得到x(t)。由于对称性，周期T（到O(1)）为：

\[ \boxed{T = 3 - 2\ln 2 \approx 1.614} \qquad(8.23) \]

对于Van der Pol振荡器。

8.5 BZ反应极限环振荡的松弛模型分析

8.5.1 降阶分析

利用 \(\varepsilon \ll \delta\) 的事实，令 \(\varepsilon \frac{dx}{dt} \approx 0\)，由(8.5)第一式得：

\[ \boxed{x = x(y) = \frac{1}{2}\left[(1 - y) + \sqrt{(1 - y)^2 + 4qy}\right]} \qquad(8.25) \]

代入原系统，得到关于y和z的二维系统：

\[ \boxed{ \begin{aligned} \delta \frac{dy}{dt} &= 2fz - y[x(y) + q] \\ \frac{dz}{dt} &= x(y) - z \end{aligned} } \qquad(8.26) \]

8.5.2 Nullcline的渐近近似

利用 \(0 < q \ll 1\)，z的nullcline（来自dz/dt = 0）为：

\[ \boxed{ z = x(y) \approx \begin{cases} \frac{qy}{1-y}, & 1 - y \gg q \\ y - 1, & y - 1 \gg q \end{cases} } \qquad(8.27) \]

y的nullcline（来自δdy/dt = 0）为：

\[ \boxed{ z \approx \begin{cases} \frac{y(1-y)}{2f}, & q \ll 1-y \ll 1 \\ \frac{qy}{f}, & y-1 \gg q \end{cases} } \qquad(8.28) \]

8.5.3 关键点坐标

nullcline的极值点为：

\[ \boxed{ \begin{aligned} &z_{\text{max}} = \frac{1}{8f} \quad \text{在} \quad y_{\text{max}} = \frac{1}{2} \\ &z_{\text{min}} = \frac{q(3 + 2\sqrt{2})}{2f} \quad \text{在} \quad y_{\text{min}} = \frac{2 + \sqrt{2}}{2} \end{aligned} } \qquad(8.29, 8.30) \]

周期轨上的关键点A、B、C、D坐标为（当 \(_1f_c < f < \,_2f_c\) 时）：

\[ \boxed{ \begin{aligned} &y_A \approx \frac{1}{8q}, \quad z_A = \frac{1}{8f}; \quad y_B \approx \frac{2+\sqrt{2}}{2}, \quad z_B \approx \frac{q(3+2\sqrt{2})}{2f} \\ &y_C \approx q(3+2\sqrt{2}), \quad z_C \approx \frac{q(3+2\sqrt{2})}{2f}; \quad y_D \approx \frac{1}{2}, \quad z_D = \frac{1}{8f} \end{aligned} } \qquad(8.31) \]

8.5.4 参数域与nullcline构型

根据f的值，nullcline有三种构型：

\[ \boxed{ \begin{aligned} &\text{当} \frac{1}{4} < f < \frac{1+\sqrt{2}}{2}: \quad \text{存在极限环（如图8.4a）} \\ &\text{当} f < \frac{1}{4}: \quad z_{\text{max}}位于定常态右侧（如图8.4b） \\ &\text{当} f > \frac{1+\sqrt{2}}{2}: \quad z_{\text{min}}位于定常态上方（如图8.4c） \end{aligned} } \qquad(8.32, 8.33, 8.34) \]

8.5.5 周期近似计算

周期由两部分组成：

\[ \boxed{T \approx T_{AB} + T_{CD}} \qquad(8.35) \]

其中

\[ T_{AB} = \int_{z_A}^{z_B} [x(y) - z]^{-1} dz \approx -\ln[4(3 - 2f + 2\sqrt{2})q] \qquad(8.37) \]

\[ \boxed{ T_{CD} \approx \frac{4f-1}{2f-1}\ln\frac{4f-1}{2f} - 2 \quad \text{（当} f \to 0.5\text{时，} T_{CD} \to 2\ln 2 - 1 \approx 0.386\text{）} } \qquad(8.39) \]

8.5.6 理论值与实验值的比较

表8.1总结了理论计算值与实验观察值的对比：

量	计算值	实验值
周期	183-228 s	110 s
**[Br⁻]₍ᵦ₎ = [Br⁻]临界	1.7×10⁻⁵[BrO₃⁻]	2×10⁻⁵[BrO₃⁻]
**[Br⁻]₍ᴄ₎ = [Br⁻]跳跃值	0.3[Br⁻]临界	0.3[Br⁻]临界
**[Br⁻]₍ₐ₎ = [Br⁻]最大值	1.6×10⁻³[BrO₃⁻] = 90[Br⁻]临界	3[Br⁻]临界

考虑到反应的复杂性和简化模型的诸多近似，理论与实验的一致性相当好。

公式汇总表

基本定义

符号	含义	典型值
\(X = HBrO_2\)	亚溴酸	-
\(Y = Br^-\)	溴离子	-
\(Z = Ce^{4+}\)	铈(IV)离子	-
\(A = BrO_3^-\)	溴酸离子	-
\(f\)	化学计量因子	0.5

无量纲参数

参数	定义	典型值
\(\varepsilon\)	\(k_5/(k_3 a)\)	\(5 \times 10^{-5}\)
\(\delta\)	\(k_4 k_5/(k_2 k_3 a)\)	\(2 \times 10^{-4}\)
\(q\)	\(k_1 k_4/(k_2 k_3)\)	\(8 \times 10^{-4}\)

核心方程

方程	描述
\(\varepsilon \frac{dx}{dt} = qy - xy + x(1-x)\)	x的演化方程
\(\delta \frac{dy}{dt} = -qy - xy + 2fz\)	y的演化方程
\(\frac{dz}{dt} = x - z\)	z的演化方程

稳定性判据

条件	物理意义
\(\phi(\delta,f,\varepsilon) = N\delta^2 + M\delta + L < 0\)	定常态不稳定（产生振荡）
\(AB - C > 0\)	Routh-Hurwitz稳定性条件
\(\frac{1}{4} < f < \frac{1+\sqrt{2}}{2}\)	f的参数域（振荡存在）

关键公式

公式	描述
\(T = 2\pi\sqrt{A/C}\)	Hopf分岔周期近似
\(x(y) = \frac{1}{2}[(1-y) + \sqrt{(1-y)^2 + 4qy}]\)	准稳态关系
\(T \approx -\ln[4(3-2f+2\sqrt{2})q] + \frac{4f-1}{2f-1}\ln\frac{4f-1}{2f} - 2\)	松弛振荡器周期

重要结论

BZ反应是化学振荡的原型：虽然它是化学反应而非生物反应，但它的方程与生物振荡器（如细胞周期、盘基网柄菌信号传导）的数学模型具有相同的形式。
FKN模型成功捕捉了振荡本质：五步简化机制结合质量作用定律，导出的三维常微分方程组能够很好地描述BZ反应的振荡行为。
Hopf分岔机制：当参数穿越 bifurcation surface \(\phi = 0\) 时，定常态失稳产生极限环振荡。
松弛振荡器特性：由于 \(\varepsilon \ll 1\)，系统表现出快速-慢速交替的松弛振荡特征，可以用分段近似方法解析求解。
理论与实验的定量一致：尽管模型经过多重简化，周期、临界浓度等关键量级的理论预测与实验测量吻合良好，验证了数学模型的有效性。

8.6 实验观察与数值模拟验证

8.6.1 极限环振荡的数值解

图8.2展示了FKN模型方程组(8.5)的数值解与实验观测到的BZ反应振荡之间的对比。使用参数值 \(f = 0.3\)，\(\delta = 1/3\)，\(q = 5 \times 10^{-3}\)，\(\varepsilon = 0.01\) 进行数值模拟，结果表明：

数值解产生的极限环轨迹与实验测量的浓度-时间曲线在形态上高度一致
振荡的振幅和周期与实验值可比
模型成功捕捉了弛豫振荡的典型特征：快速跳转与慢速演化交替

这种理论与实验的符合不是偶然的，而是建立在对反应机理的深入理解和合理的模型简化基础之上。Field等人（1972）的详细化学动力学分析为模型参数的确定提供了可靠依据。

8.6.2 阈值行为与双稳态

当参数f处于特定范围时，FKN模型不仅表现出极限环振荡，还可能出现阈值行为和双稳态现象。在可逆反应条件下，模型可以存在三个正定常态，其中两个是稳定的。这对应于生物学中常见的"开关"行为：一个小的扰动可以使系统从一个稳态跳转到另一个稳态。

此外，数值研究还发现了周期性爆裂（periodic bursting）和混沌（chaos）行为。Györgyi和Field（1991）提出的BZ反应模型展示了确定性混沌的存在。Barkley等人（1987）则演示了包括周期性爆裂、滞后现象和周期-混沌序列在内的复杂行为。

8.7 反应-扩散系统与空间模式

8.7.1 从均匀系统到空间非均匀系统

本章中我们只考虑了均匀或充分搅拌的系统，即假设反应器中各处的化学物质浓度处处相同。然而，当化学物质能够扩散时，系统会变得完全不同。扩散与反应的耦合产生了一系列令人惊叹的空间-时间模式现象，包括：

螺旋波（spiral waves）
靶心波（target patterns）
化学湍流（chemical turbulence）
相位波（phase waves）

这些空间模式的存在是BZ反应持续吸引生物学和物理学研究人员的重要原因之一。它们为生物形态建成的空间自组织研究提供了重要的实验模型和理论框架。

8.7.2 BZ反应对生物学的启示

BZ反应的研究对生物学产生了深远的影响。它不仅是一个"原型振荡器"，更重要的是，它证明了一个开放的、远离平衡的化学系统可以表现出自组织的时空模式。这一发现对理解以下生物学问题具有重要启发意义：

胚胎发育中的周期性现象：细胞分化和形态建成的时空组织
心脏节律：心室肌细胞的动作电位传播与心律失常
神经网络活动：脑皮层中的行波和同步振荡
细胞周期：细胞分裂的调控网络

Winfree（1987）将BZ反应的数学模型应用于心室三维电活动的建模，Tyson（1991）将其应用于细胞周期分子机制的研究，这些都表明BZ反应模型已经超越了化学本身，成为理解生命系统中振荡现象的重要工具。

公式汇总表（续）

阈值与临界值

公式	描述
\(y_{\text{max}} = \frac{1}{2}\)	y最大值位置
\(z_{\text{max}} = \frac{1}{8f}\)	z最大值
\(y_{\text{min}} = \frac{2+\sqrt{2}}{2}\)	y最小值位置
\(z_{\text{min}} = \frac{q(3+2\sqrt{2})}{2f}\)	z最小值

分岔边界

参数范围	系统行为
\(f < \frac{1}{4}\)	nullcline如图8.4b，无振荡
\(\frac{1}{4} < f < \frac{1+\sqrt{2}}{2}\)	nullcline如图8.4a，存在极限环
\(f > \frac{1+\sqrt{2}}{2}\)	nullcline如图8.4c，无振荡

时间尺度

量	表达式	典型值
\(t_0\)	\(1/k_5\)	50 s
\(T_{AB}\)	\(-\ln[4(3-2f+2\sqrt{2})q]\)	O(ln(1/q))
\(T_{CD}\)	\(\frac{4f-1}{2f-1}\ln\frac{4f-1}{2f} - 2\)	O(1)
周期T	\(T_{AB} + T_{CD}\)	~183-228 s（计算值）

重要结论

BZ反应是化学振荡的原型：虽然它是化学反应而非生物反应，但它的方程与生物振荡器（如细胞周期、盘基网柄菌信号传导）的数学模型具有相同的形式。
FKN模型成功捕捉了振荡本质：五步简化机制结合质量作用定律，导出的三维常微分方程组能够很好地描述BZ反应的振荡行为。
Hopf分岔机制：当参数穿越 bifurcation surface \(\phi = 0\) 时，定常态失稳产生极限环振荡。
松弛振荡器特性：由于 \(\varepsilon \ll 1\)，系统表现出快速-慢速交替的松弛振荡特征，可以用分段近似方法解析求解。
理论与实验的定量一致：尽管模型经过多重简化，周期、临界浓度等关键量级的理论预测与实验测量吻合良好，验证了数学模型的有效性。
更广泛的意义：当考虑扩散效应时，BZ反应展现出螺旋波、靶心波等丰富的空间模式，为理解生物形态建成中的自组织现象提供了重要启示。

笔记整理基于 Murray - Mathematical Biology I，第8章，J.D. Murray著 参考文献：Field et al. (1972), Field & Noyes (1974), Tyson (1976, 1977, 1985, 1991, 1994), Winfree (1984, 1987, 2000), Hastings & Murray (1975), Györgyi & Field (1991), Barkley et al. (1987)