第 6 章：反应动力学（Reaction Kinetics）

6.1 酶动力学：基本酶反应

生化反应在所有活体生物中不断发生，其中多数涉及被称为酶的蛋白质——它们以极高的效率起催化作用。酶只对特定的化合物（即底物）起选择性反应；例如，红细胞中的血红蛋白就是一种酶，而与它结合的氧就是它的底物。酶在调节生物过程中有重要作用，例如作为激活剂或抑制剂。要理解酶的作用必须研究酶动力学，即反应速率的研究、各反应物的时间行为，以及影响它们的条件。本章将讨论若干模型反应机制——它们能映照大量真实反应——以及若干一般性的反应现象与相应的数学实现。

最基础的酶促反应之一最早由 Michaelis 和 Menten 在 1913 年提出：底物 S 与酶 E 反应形成复合物 SE，SE 进一步被转化生成产物 P 并释放酶。其反应机制可写为

S + E ⇌ SE → P + E

其中 k1、k−1、k2 是与各反应速率相关的常数参数。双箭头表示该反应可逆，单箭头表示反应只能单向进行。整个机制表达了 S 通过酶 E 的催化被转化为 P：1 分子 S 与 1 分子 E 结合形成 1 分子 SE，SE 最终产生 1 分子 P 和 1 分子 E。

质量作用定律（Law of Mass Action）指出，反应的速率与反应物浓度的乘积成正比。用小写字母表示浓度（s = [S], e = [E], c = [SE], p = [P]），对上述机制应用质量作用定律，便得到一组非线性反应方程：

ds/dt = −k1 es + k−1 c de/dt = −k1 es + (k−1 + k2) c dc/dt = k1 es − (k−1 + k2) c dp/dt = k2 c

其中 k 为速率常数，是应用质量作用定律时的比例常数。完整的数学表述还需要初始条件——取反应开始时的值：s(0) = s0, e(0) = e0, c(0) = 0, p(0) = 0。这样便给出了浓度作为时间函数、从而反应速率作为时间函数的解。任何反应动力学问题中我们只关心非负浓度。

dp/dt 方程与前三个方程解耦，一旦 c(t) 求出便可通过积分得到 p(t)，因此只需分析前三个方程。在该机制中，酶 E 是催化剂，只起促进作用，其总浓度（游离的与结合的）保持恒定。这一酶守恒律也可由前三个方程中第 2、3 个相加得到：de/dt + dc/dt = 0，即 e(t) + c(t) = e0。由此方程组简化为只关于 s 和 c 的两个方程：

ds/dt = −k1 e0 s + (k1 s + k−1) c dc/dt = k1 e0 s − (k1 s + k−1 + k2) c

通常的处理方法是假设复合物 c 的形成有一个非常快的初始阶段，之后便基本处于平衡，即 dc/dt ≈ 0；由上式的第二个方程得到 c 关于 s 的代数关系

c(t) = e0 s(t) / [s(t) + Km], Km = (k−1 + k2) / k1

代入第一个方程便得到

ds/dt = −k2 e0 s / (s + Km)

Km 称为 Michaelis 常数。由于传统上认为酶的浓度远小于底物浓度，因此在这一初始瞬态阶段底物浓度被认为基本不变。这种近似下动力学可由上式加初始条件 s = s0 描述，称为拟稳态（pseudo-steady state）近似或准稳态（quasi-steady state）近似。解得关于 s(t) 的隐式解

s(t) + Km ln s(t) = s0 + Km ln s0

代入 c 关于 s 的关系可得到 c(t)。该 c(t) 表达式并不满足 c(0) = 0，但或许大多数时候是合理的近似——这是通常应用该方法时的信念。事实上在许多实验情形下是够用的，但并不总是够用。

实际上该系统涉及两个时间尺度：一个是 t = 0 附近的初始瞬态尺度，另一个是较长的、底物显著变化的尺度——在该尺度上酶复合物可由 c = e0 s / (s + Km) 合理近似。这种基本推理引出几个重要问题：(i) 初始瞬态有多快？(ii) 在哪些参数范围内上述近似足够好？(iii) 若酶浓度不显著小于底物浓度，应如何处理？其他问题在后续章节讨论。

作为第一步，必须对该系统无量纲化。无量纲化的方法有多种。一个关键的无量纲量是时间，因为上述基本假设依赖于初始瞬态有多短。拟稳态分析的标准做法是引入以下无量纲量：

τ = k1 e0 t, u(τ) = s(t) / s0, v(τ) = c(t) / e0 λ = k2 / (k1 s0), K = (k−1 + k2) / (k1 s0) = Km / s0, ε = e0 / s0

若 ε ≪ 1 这是合理的无量纲化。代入原方程与初始条件，便得到传统拟稳态近似的无量纲方程组

du/dτ = −u + (u + K − λ) v ε dv/dτ = u − (u + K) v u(0) = 1, v(0) = 0

由 (6.12) 可知 K − λ > 0。由反应机制可知最终稳态为 u = 0 且 v = 0——即底物与底物-酶复合物浓度均为零。我们关心反应的时间演化，因此需求该非线性系统的解（不能简单闭式求解）。但可定性观察解的形态：τ = 0 附近，du/dτ < 0，故 u 从 1 开始下降；dv/dτ > 0，v 从 0 上升并继续上升直至 v = u/(u + K)，此时 dv/dτ = 0，但 u 仍在下降。v 达到最大值后便随 u 一同下降趋于零，u 则单调趋于零。无量纲酶浓度 e/e0 先从 1 下降再回升至 1。典型解如图 6.1 所示。

6.2 瞬态时间估计与无量纲化

在生物学中普遍认为，酶卓越的催化效力反映在反应中所需酶的浓度远小于底物浓度。在 Michaelis–Menten 模型的无量纲形式中即 ε = e0/s0 ≪ 1。然而如上所述，e0/s0 ≪ 1 并非普遍成立。Segel（1988）以及 Segel 和 Slemrod（1989）通过引入新的无量纲方式扩展了传统分析，使其既包含此情形也覆盖 e0/s0 = O(1) 的情形。本节描述他们的分析。

首先需要对两个时间尺度——快瞬态尺度 tc 与慢时间尺度 ts（底物显著变化的时间）——做估计。在初始瞬态中，复合物 c(t) 迅速增加而 s(t) 变化不大，因此由 dc/dt 方程（令 s(t) = s0）可估出该快时间尺度 tc = 1 / [k1 (s0 + Km)]。长尺度 ts 的估计：取底物的最大可能变化 s0 除以底物变化率的最大值（令 s = s0），得 ts ≈ (s0 + Km) / (k2 s0)。

拟稳态近似有效的一个前提是初始快瞬态时间远小于底物显著变化的长时间尺度，即 tc ≪ ts。结合 tc、ts 的表达式，这意味着参数须满足 k2 e0 / [k1 (s0 + Km)^2] ≪ 1。拟稳态近似的另一个要求是 s(t) 的初始条件可取 s(0) = s0，即在快瞬态中底物的消耗量 Δs(t) 仅为 s0 的小部分：|Δs/s0| ≪ 1。Δs(t) 的过估由最大消耗速率 k1 e0 s0 乘以 tc 得到，除以 s0 后给出 ε = e0 / (s0 + Km) ≪ 1。这比前述条件更强，因此是保证拟稳态近似的条件。可见即使 e0/s0 = O(1)，只要 Km 足够大该条件仍可满足，正如许多反应中的情形。

由于无量纲化取决于所关注的时间尺度，从两个时间尺度 tc 与 ts 中选哪个取决于想看哪个区段：用 tc 关注的是 t = 0 附近，用 ts 关注的是 s(t) 显著变化的长时间段。涉及两个这样的时间尺度的问题一般是奇异摄动问题，对此有标准方法（Murray 1984 的小书便是介绍之一）。本节随后对这种问题作详细的奇异摄动分析。

使用快瞬态尺度 tc 时，引入以下无量纲变量和参数：

u(τ) = s(t) / s0, v(τ) = (s0 + Km) c(t) / (e0 s0) τ = t / tc = k1 (s0 + Km) t, Km = (k−1 + k2) / k1 ρ = k−1 / k2, σ = s0 / Km, ε = e0 / (s0 + Km)

代入原方程与初始条件得

du/dτ = ε [−u + (σ / (1 + σ)) u v + (ρ / ((1 + ρ)(1 + σ))) v] dv/dτ = u − (σ / (1 + σ)) u v − v / (1 + σ) u(0) = 1, v(0) = 0

用长时间尺度 ts 时，可将时间写为 T = (1 + ρ) t / ts = (1 + ρ) k2 e0 t / (s0 + Km) = ε (1 + ρ) k2 t，引入系数 (1 + ρ) 仅是为代数方便。结合上述无量纲形式但用 T 替代 τ，原方程变为

du/dT = −(1 + σ) u + σ u v + (ρ / (1 + ρ)) v ε dv/dT = (1 + σ) u − σ u v − v

需要记住我们研究的是原方程 (6.7)；三个方程组 (6.13)、(6.21) 与 (6.23) 实际上是同一个系统，差别只在于无量纲化的方式。小参数 ε 出现的位置不同决定了所用的解析方法。本节中我们以一种具体方式做奇异摄动分析，并展示如何得到对所有时间统一有效的近似解。

在结束无量纲化的话题前，应提一下酶过量（ε 不小）的情形。这出现在多种酶反应中，也出现在一个完全不同的情境中——抗原引起的 T 细胞增殖。De Boer 和 Perelson（1994）研究过该系统：这里的"底物"是复制中的细胞，"酶"是抗原递呈细胞上的位点，"复合物"是结合的 T 细胞与抗原递呈细胞。其动力学表示为

E + S ⇌ C → E + 2S

Borghans 等（1996）研究了酶过量的该反应系统，将上述分析扩展得到了一个统一有效的渐近解。他们将 s 方程替换为底物总量 s̄(t) = s(t) + c(t) 的方程（即将 s 与 c 方程相加），并与 c 方程（用 c 与 s̄ 表示）一起，结合边界条件得到

d s̄/dt = −k2 c dc/dt = k1 [(e0 − c)(s̄ − c) − Km c] s̄(0) = s0, c(0) = 0

分析稍复杂但概念相似。他们导出了等效拟稳态近似成立的条件并讨论若干例子，包括一般性捕食者-猎物模型。

6.3 Michaelis–Menten 拟稳态分析

本节对上述某一无量纲方程组做奇异摄动分析。原则上该方法可应用于任何一种，为具体起见，我们对 (6.13) 做详细教学性分析，解释该方法的背景与用法，从而得到对 0 < ε ≪ 1 的极好近似解（渐近解）。事先应重申：特定的无量纲化 (6.12) 只是若干种选择中的一种。下一节中我们将分析一个稍复杂的实际酶反应系统，使用更一般的公式——该系统中 e0/s0 并不小，但 Km 足够大使得 (6.20) 中的 ε 很小。Frenzen 和 Maini（1988）研究的实际反应就是这种情形，他们用的分析方法与本节类似。

考虑方程组 (6.13)。设简单地用 Taylor 展开形式寻找 u 与 v 的正则展开

u(τ; ε) = Σ εn un(τ), v(τ; ε) = Σ εn vn(τ)

代入 (6.13) 并令 ε 同次幂相等，便得到关于 un(τ)、vn(τ) 的逐阶微分方程序列——这等价于假设 u、v 是 ε → 0 时的解析函数。O(1) 方程为

du0/dτ = −u0 + (u0 + K − λ) v0 0 = u0 − (u0 + K) v0 u0(0) = 1, v0(0) = 0

这种做法已有难处：第二个方程是代数方程，且不满足初始条件（若 u0 = 1，v0 = 1/(1 + K) ≠ 0）。若解此 O(1) 方程组

v0 = u0 / (u0 + K), du0/dτ = −λ u0 / (u0 + K)

积分得 u0(τ) + K ln u0(τ) = A − λ τ；要求 u0(0) = 1 得 A = 1，因此

u0(τ) + K ln u0(τ) = 1 − λ τ v0(τ) = u0(τ) / (u0(τ) + K)

这与 (6.11)、(6.9) 的解相同。然而该解并非对所有 τ ≥ 0 一致有效的近似解，因为 v0(0) ≠ 0。这并不奇怪：方程组 (6.25) 是在 (6.13) 中令 ε = 0 得到的；它只有一个积分常数（来自 u 方程），故无法同时满足 u0 与 v0 的初始条件。

(6.13) 中小参数 ε 乘在导数项上，提示这是一个奇异摄动问题。一类典型的奇异摄动问题可立即识别：在令 ε = 0 时方程组阶数降低；降阶系统在一般情形下无法满足所有初始条件。奇异摄动技巧是确定此类方程组在 ε 小时渐近解的重要且强大的方法。渐近解通常作为精确解的极好近似。关于若干关键技巧的实用且初等的讨论可参见 Murray（1984）的渐近分析书。下面我们详细描述奇异摄动方法的理念与实际技术，并推导 (6.13) 在 0 < ε ≪ 1 时的渐近解。主要目的是指出在实际情形中何时可以忽略 ε 项。

由于解 (6.26)（特别是 v0(τ)）不满足初始条件（且高阶项也不能补救），必须断定 u(τ; ε) 与 v(τ; ε) 中至少有一个在 ε → 0 时不是 ε 的解析函数。在导出 (6.25) 时，我们暗中假设 v 是解析的。由于忽略了 ε dv/dτ 而使 v(0) = 0 不可满足，因此在分析中——至少在 τ = 0 附近——必须保留该项。因此在 τ = 0 附近更合适的时间尺度是 σ = τ/ε 而非 τ：此时 ε dv/dτ = dv/dσ。从 σ = τ/ε 变换看，τ = 0 附近的小邻域被放大了，固定 0 < τ ≪ 1 在 ε → 0 时有 σ → ∞，即 τ = 0 的极小邻域对应 σ 的极大区域。用此变换分析 (6.13) 在 τ = 0 附近的行为，再分析远离 τ = 0 的解，最后将两部分合并得到对所有 τ ≥ 0 一致有效的解。

在变换 σ = τ/ε, u(τ; ε) = U(σ; ε), v(τ; ε) = V(σ; ε) 下，方程 (6.13) 变为

dU/dσ = −ε U + ε (U + K − λ) V dV/dσ = U − (U + K) V U(0) = 1, V(0) = 0

若令 ε = 0 并做正则摄动展开

U(σ; ε) = Σ εn Un(σ), V(σ; ε) = Σ εn Vn(σ)

则得

dU0/dσ = 0 dV0/dσ = U0 − (U0 + K) V0 U0(0) = 1, V0(0) = 0

这与原方程 (6.28) 阶数相同。解为

U0(σ) = 1 V0(σ) = (1 / (1 + K)) (1 − e^{−(1 + K) σ})

该解不能期望对所有 τ ≥ 0 成立；若成立则意味着 dv/dσ = ε dv/dτ 对所有 τ 都是 O(1)。它仅在 0 ≤ τ ≪ 1 范围内是有效的奇异（内层）解；外层解 (6.26) 在所有不紧邻 τ = 0 的 τ 范围内有效。令 ε → 0，对任何固定 0 < τ ≪ 1，σ → ∞。因而在 ε → 0 极限下，应有 τ → 0 时的解 (6.26) 等于 σ → ∞ 时的解 (6.31)：即奇异解的 σ → ∞ 极限匹配非奇异解的 τ → 0 极限。这正是奇异摄动理论中匹配（matching）的实质。事实上从 (6.31) 与 (6.26) 可看到

lim_{σ→∞} [U0(σ), V0(σ)] = [1, 1/(1+K)] = lim_{τ→0} [u0(τ), v0(τ)]

图 6.1 示意了 u(τ)、v(τ) 及由 (6.6) 的无量纲形式 e/e0 = 1 − v(τ) 给出的无量纲酶浓度。τ = 0 附近 O(ε) 厚的薄层有时被称为边界层，是解发生极快变化的 τ 区域。从 (6.31) 可得

(dV/dτ)|{τ=0} ∼ ε^{-1} (dV0/dσ)| ≫ 1} = ε^{-1

原方程 (6.13) 中从第二个方程与边界条件也可直接看到这一点。

为以系统奇异摄动方式推进，先寻求全系统 (6.13) 的外层解为正则级数展开 (6.24)。逐阶方程为

O(1): du0/dτ = −u0 + (u0 + K − λ) v0, 0 = u0 − (u0 + K) v0 O(ε): du1/dτ = u1 (v0 − 1) + (u0 + K − λ) v1, dv0/dτ = u1 (1 − v0) − (u0 + K) v1

这些方程在 τ > 0 时有效。其解含未定的积分常数（每一阶一个），须由这些解在 τ → 0 的极限与奇异解在 σ → ∞ 的极限相匹配来确定。

对解的奇异部分（在 0 ≤ τ ≪ 1 有效），将 (6.29) 代入 (6.28) 并令 ε 同次幂相等得

O(1): dU0/dσ = 0, dV0/dσ = U0 − (U0 + K) V0 O(ε): dU1/dσ = −U0 + (V0 + K − λ) V0, dV1/dσ = (1 − V0) U1 − (V0 + K) V1

等等。这些解须满足 σ = 0（即 τ = 0）处的初始条件

U0(0) = 1, Un≥1(0) = 0, Vn≥0(0) = 0

此时奇异解是唯一确定的。在一般的奇异摄动问题中并非如此（Murray 1984）。匹配要求选取外层方程 (6.32) 解中未定的积分常数，使得对 ε 的各阶

lim_{σ→∞} [U(σ; ε), V(σ; ε)] = lim_{τ→0} [u(τ; ε), v(τ; ε)]

由 (6.32) 与之前的讨论，外层解为

u0(τ) + K ln u0(τ) = A − λ τ, v0(τ) = u0(τ) / (u0(τ) + K)

其中 A 是待匹配的积分常数。(6.33) 第一个方程加 (6.34) 的解前面已给出 (6.31)。由匹配过程

lim_{σ→∞} V0(σ) = 1/(1+K) = lim_{τ→0} v0(τ) ⇒ v0(0) = 1/(1+K) = u0(0) / (u0(0) + K) ⇒ u0(0) = 1 ⇒ A = 1

于是对 0 < ε ≪ 1 的 O(1) 阶一致有效的渐近解为：τ > 0 用 (6.26)，0 < τ ≪ 1 用 (6.31)；奇异部分更自然地用 0 ≤ τ/ε < ∞ 表示。可继续由 (6.33) 求 U1(σ)、V1(σ) 及由 (6.32) 求 u1(τ)、v1(τ)，依次求到任意 ε 阶；虽然方程都是线性的，解却越来越复杂。这样便得到对所有 τ ≥ 0 在 0 < ε ≪ 1 时一致有效的非线性动力学 (6.13) 的渐近解。总结 O(1) 阶结果为

u(τ; ε) = u0(τ) + O(ε), u0(τ) + K ln u0(τ) = 1 − λ τ v(τ; ε) = V0(σ) + O(ε), V0(σ) = (1/(1+K)) (1 − e^{−(1+K)τ/ε}), 0 < τ ≪ 1 v(τ; ε) = v0(τ) + O(ε), v0(τ) = u0(τ) / (u0(τ) + K), 0 < ε ≪ τ

在大多数生物应用中 0 < ε ≪ 1，因此只需计算 O(1) 项；O(ε) 项的贡献可忽略。

为完成原动力学问题 (6.3)、(6.4) 的分析，写无量纲产物与游离酶浓度

z(τ) = p(t) / s0, w(τ) = e(t) / e0

用 (6.36) 关于 u、v 的结果，结合 (6.5) 与 (6.6) 得

z(τ) = λ ∫_0^τ v(τ') dτ', w(τ) = 1 − v(τ)

底物-酶复合物 v(τ; ε) 的快速变化发生在 τ = O(ε) 这段无量纲时间内，对应的量纲时间 t 也极短，量级为 O(1/(k1 s0))，在许多实验情形下不可测量。因此在许多实验中奇异解 v(τ) 永远观测不到。相关的解是 O(1) 外层解 u0(τ)、v0(τ)（即 (6.26)），由 (6.13) 令 ε = 0 并只满足 u(τ) 的初始条件得到。换言之，我们说复合物 v(τ) 的反应实质上处于稳态，或数学上 ε dv/dτ ≈ 0：v 反应极快，几乎始终处于平衡。这就是通常的 Michaelis–Menten 拟稳态假设。

(6.13) 的形式一般是

du/dτ = f(u, v), dv/dτ = ε^{−1} g(u, v), 0 < ε ≪ 1

这立即表明若 g(u, v) 不接近零，则 dv/dτ ≫ 1。所以 v 反应比 u 反应快得多。v 反应极快达到拟稳态，意味着在 τ = O(1) 时它实质上处于平衡，模型机制可由

du/dτ = f(u, v), g(u, v) = 0, u(0) = 1

近似。若解代数方程 g(u, v) = 0 得 v = h(u)，则

du/dτ = f(u, h(u))

这就是底物浓度的速率或摄取方程。生物过程的许多建模都依赖于对摄取函数 f(u, h(u)) 的定性假设。

生物学上有兴趣的是反应速率，即底物浓度 u(τ) 求出后 du/dτ 的值。实验中通常在多个时间点测量量纲底物浓度 s(t)，外推至 t = 0，求初始速率 [ds/dt]_{t=0} 的大小 r。由于时间测量几乎都在 τ ≫ ε 即 t ≫ 1/(k1 s0)（量级为秒），相应的解析速率由非奇异（外层）解给出。由 (6.36) 第一式，O(1) 阶 0 < ε ≪ 1 的速率 r0 为

r0 = (du0/dτ)|_{τ=0} = λ u0(0) / (u0(0) + K) = λ / (1 + K)

量纲上，利用 (6.12)，O(1) 阶反应速率 R0 为

R0 = k2 e0 s0 / (s0 + Km) = Q s0 / (s0 + Km), Km = (k−1 + k2) / k1, Q = [R0]_{max} = k2 e0

其中 Q 是反应的最大速率，Km 是 (6.9) 中的 Michaelis 常数。这便是从生物学角度通常所求的基于拟稳态假设的速率。由 (6.13) 底物的精确初始速率为 [du/dτ]{τ=0} = −1，复合物的为 [dv/dτ] = 1/ε。

当一种底物的摄取被描述为 Michaelis–Menten 摄取时，指的就是形如 (6.41) 的反应速率，图 6.2 给出其示意。实际随时间变化的反应速率是量纲形式的外层解 du0/dτ。即 Michaelis–Menten 摄取由方程

ds/dt = − Q s / (Km + s)

控制。这只是 (6.39) 的量纲形式（也与 (6.10) 相同）。对 s ≪ Km 摄取对 s 是线性的，右端近似为 −Q s / Km。最大速率 Q = k2 e0 取决于产物反应 SE → P + E 的速率常数 k2，k2 称为反应机制 (6.1) 的限速步骤。

拟稳态假设尽管有用且重要，但因忽略了 (6.13) 中的 ε dv/dt 项并把实验结果套到一个不能同时满足所有初始条件的理论上，仍有所失。利用 Michaelis–Menten 理论与实验结果所能确定的是如图 6.2 的曲线，从而给出最大速率 Q 和 Michaelis 常数 Km 的值。这并不确定全部三个速率常数 k1、k−1、k2，仅确定 k2 与三者间的一个关系。要确定全部三者需要对 τ = O(ε) 时的实验测量。然而通常拟稳态假设下的摄取速率（即 Michaelis–Menten 理论）就够用了。

6.4 自杀底物动力学

自杀底物系统（mechanism-based inhibitor）是一类具有相当实验兴趣的酶系统（如 Seiler 等 1978, Walsh 1984），由 Walsh 等（1978）所表示：

E + S ⇌ X → Y → E + P, k4 ↓ Ei

其中 E、S、P 分别表示酶、底物与产物，X 与 Y 是酶-底物中间物，Ei 是失活的酶，k 为正常速率常数。在该系统中 Y 可走两条路径之一：以 k3 速率到 E + P 或以 k4 速率到 Ei。这两条路径的速率之比 k3/k4 称为分配比（partition ratio），记为 r。这两条路径在该反应的时间尺度上都视为不可逆（Waley 1980）。S 被称为自杀底物（suicide substrate），因为它像底物一样与酶的活性位点结合，但酶却将其转化为抑制剂并不可逆地失活酶——故酶"自杀"。

自杀底物之所以重要是因为它提供了一种针对特定酶失活的方法。它们在药物施用中尤为有用，因为它们的常见形态无害，只有被指定的酶才能将其转化为抑制剂形式。例：自杀底物曾被研究用于治疗抑郁（单胺氧化酶抑制剂，Seiler 等 1978）、癫痫（脑 GABA 转氨酶抑制剂，Walsh 1984）以及某些肿瘤（鸟氨酸脱羧酶抑制剂，Seiler 等 1978）。

自杀底物动力学曾由 Waley（1980）和 Tatsunami 等（1981）考虑，他们关心的是决定底物是否在所有酶失活之前耗尽。Waley 建议关键因子是 rµ，其中 µ 是酶初浓度与底物初浓度之比即 e0/s0（即 (6.12) 中的 ε）。Tatsunami 等（1981）则认为决定因子是 (1+r)µ。当 (1+r)µ > 1 底物耗尽；当 (1+r)µ < 1 所有酶失活；二者相等时两者同时发生。Burke 等（1990）用奇异摄动分析做了深入研究，本节的分析即沿用之。兴趣在于 e0/s0 不小的情形——这在 Waley 和 Tatsunami 等使用拟稳态近似时被实际地假定了。从上节经验可知，随 e0/s0 增大近似的有效性下降。Duggleby（1986）指出事实上 e0/s0 并不小。因此必须用奇异摄动技巧，但此时须用与 (6.20) 等价的无量纲化，而非 (6.12)，因为此情形下小量是 e0/(s0 + Km)，而 e0/s0 = O(1)。

由 (6.43) 用质量作用定律得速率方程

d[S]/dt = −k1 [E][S] + k−1 [X] d[E]/dt = −k1 [E][S] + k−1 [X] + k3 [Y] d[X]/dt = k1 [E][S] − k−1 [X] − k2 [X] d[Y]/dt = k2 [X] − k3 [Y] − k4 [Y] d[Ei]/dt = k4 [Y] d[P]/dt = k3 [Y]

其中 [ ] 表示浓度，t 为时间。典型实验初始条件为

E(0) = e0, S(0) = s0, X(0) = Y(0) = Ei(0) = P(0) = 0

(6.49) 仍解耦；求出 [Y] 后 [P] 由积分得出。

利用酶的守恒律可进一步降阶：将 (6.45)–(6.48) 相加得

d/dt {[E] + [X] + [Y] + [Ei]} = 0 ⇒ [E] + [X] + [Y] + [Ei] = e0

用此式消去 [E]，得到降阶系统

d[S]/dt = −k1 (e0 − [X] − [Y] − [Ei]) [S] + k−1 [X] d[X]/dt = k1 (e0 − [X] − [Y] − [Ei]) [S] − (k−1 + k2) [X] d[Y]/dt = k2 [X] − (k3 + k4) [Y] d[Ei]/dt = k4 [Y]

无量纲形式

无量纲化有多种方式。由于 e0/s0 = O(1)，采用与 6.2 节中 (6.20) 等价的过程。引入变量

[S] = s0 s, [X] = (e0 s0 / (s0 + Km)) x, [Y] = e0 y, [Ei] = e0 ei

其中 Km = (k−1 + k2) / k1。快瞬态时间尺度（参照 (6.20)）取为

τ = t / tc = t k1 (s0 + Km)

准稳态时间尺度为

T = (1 + ρ) t / ts = t ε (k−1 + k2) (1 + ρ)

其中 ρ 如 (6.66)，且 ε = e0 / (e0 + Km)。在尺度 τ 下，快瞬态阶段的方程 (6.53)–(6.56) 化为

ds/dτ = ε [−s + (σ/(1+σ)) s x + s y + s ei + (ρ / ((1+ρ)(1+σ))) x] dx/dτ = s − (σ/(1+σ)) s x − s y − s ei − x/(1+σ) dy/dτ = (σ / ((1+σ)^2 (1+ρ))) x − (ψ / (1+σ)) y dei/dτ = (φ / (1+σ)) y

其中 σ = s0/Km, ρ = k−1/k2, ψ = (k3 + k4) / (k−1 + k2), φ = k4 / (k−1 + k2)。

初始条件 (6.50) 由 (6.57) 变为 s(0) = 1, x(0) = 0, y(0) = 0, ei(0) = 0。方程 (6.62)–(6.65) 是 (6.21) 的等价形式，给出奇异（内层）解。

在 T 尺度下，非奇异（外层）拟稳态阶段的速率方程为

ds/dT = −s [(σ + 1) − σ x − (σ + 1) y − (σ + 1) ei] + (ρ / (1+ρ)) x ε dx/dT = s [(σ + 1) − σ x − (σ + 1) y − (σ + 1) ei] − x ε dy/dT = (σ / ((1+σ)(1+ρ))) x − ψ y ε dei/dT = φ y

参数 ε、σ、ρ、ψ、φ 由 (6.66) 给出。这些是 (6.23) 在此处的等价形式。

渐近技术与解

利用 0 < ε ≪ 1（(6.61) 中的 ε），用上节详细讨论的奇异摄动技巧求解。有一些分析上重要的差异（不仅是代数上更复杂）。

内层或奇异解：从快瞬态阶段方程 (6.62)–(6.65) 加初始条件 (6.67) 出发；因 ε 小，对每个变量 s、x、y、ei 取 Taylor 展开

s(τ) = s^(0)(τ) + ε s^(1)(τ) + ε^2 s^(2)(τ) + ...

代入并按 ε 同次幂相等得

ds^(0)/dτ = 0, dy^(0)/dτ = −(ψ/(1+σ)) y^(0)

加上 (6.67) 给出唯一解 s^(0)(τ) ≡ 1, y^(0)(τ) ≡ 0。同理 (6.65) 给出 O(1) 阶 dei^(0)/dτ = −(φ/(1+σ)) y^(0) = 0，故 ei^(0)(τ) ≡ 0（因 ei(0) = 0）。最后代入 (6.65) 得

dx^(0)/dτ = s^(0) − s^(0) y^(0) − s^(0) ei^(0) − x^(0)/(1+σ) − σ s^(0) x^(0)/(1+σ)

代入上面的解后化为 dx^(0)/dτ = 1 − x^(0)，加 x^(0)(0) = 0 得 x^(0)(τ) = 1 − e^{−τ}。

为得到 y 与 ei 的非零解，需要至少确定 O(ε) 项 y^(1)(τ) 与 ei^(1)(τ)。这需要匹配 O(ε) 系数。注意由 (6.61) 与 (6.66)

ε = e0 / [s0 (1 + Km/s0)] = (e0/s0) (σ/(1+σ))

故 σ = (s0/e0) ε + O(ε^2)。因 s0/e0 = O(1)，有 σ = O(ε)。引入 σ 的相似变量 σ = ε p，p 是 O(1) 的常数。把 (6.79) 代入 (6.64) 并令 O(ε) 项相等得

dy^(1)/dτ = (p / (1+ρ)) x^(0) − ψ y^(1)

已知 x^(0)(τ) = 1 − e^{−τ}，可解此线性方程

y^(1)(τ) = (p / (ψ(1+ρ))) [(1 − e^{−ψτ})/ψ + (e^{−ψτ} − e^{−τ})/(ψ − 1)]

由 (6.65) 中 O(ε) 系数匹配可得 dei/dτ 关于 y^(1) 的方程；稍作代数得

ei^(1)(τ) = (φ p / (1+ρ)) [τ/ψ + (e^{−τ} − 1)/(ψ − 1) + (1 − e^{−ψτ})/(ψ^2 (ψ − 1))]

在求 ei^(1)(τ) 时假设 φ = O(1)。若 φ = O(ε)，则须用另一相似变量 q = ε φ，结果发现 ei^(1)(τ) 与 ei^(2)(τ) 与上式结果相同。可类似地求高阶项。例如 (6.62) 的 O(ε) 项给出

s^(1)(τ) = −τ/(1+ρ) + (ρ/(1+ρ)) (e^{−τ} − 1)

所有这些解都满足初始条件 (6.67)。

外层或拟稳态解：现在求长时间尺度的解——给出拟稳态近似。然后匹配两个时段的解。长时间尺度解在一般情况下不满足初始条件；未定的积分常数通过 6.3 节的匹配法求定。

求 (6.68)–(6.71) 形式的解

s(T) = s^(0)(T) + ε s^(1)(T) + ε^2 s^(2)(T) + ...

对每个变量 s、x、y、ei 都做这种展开。代入并按 ε 同次幂相等。这里须用匹配渐近展开法解未定的积分常数——即内层解 τ → ∞ 极限匹配外层解 T → 0 极限。

取 O(1) 项，记住 σ = ε p = O(ε)（由 (6.79)）：(6.69) 给出

0 = s^(0) − s^(0) y^(0) − x^(0) − s^(0) ei^(0)

假设 ψ = O(1) 时 (6.71) 给出 y^(0) = 0；合起来得

x^(0) = s^(0) (1 − ei^(0))

类似地得 y^(1) = (p / (ψ(1+ρ))) x^(0)。

要进一步匹配各系数，需要确定每项的量级。实验上有两种基本不同的结果：底物全部耗尽或酶全部失活。这对应 φ = O(1) 且 ψ = O(1)，以及 ψ = O(1) 且 φ = O(ε)（参数关系见 (6.66)）。须对每种约束情形分别求解。

情形 1：ρ = O(1), ψ = O(1), φ = O(1)。这是所有速率常数同量级的情形。假设 φ = O(1) 时，(6.68) 联合 (6.79)、(6.84)、(6.86) 给出

ds^(0)/dT = −(1/(1+ρ)) s^(0) (1 − ei^(0))

由 (6.69) 联合 (6.79)、(6.84)、(6.86)、(6.87) 得

dei^(0)/dT = (φ ρ / (ψ(1+ρ))) s^(0) (1 − ei^(0))

将两式相除并积分得

ei^(0)(T) = (1/β) (B − s^(0)(T))

其中 B 为积分常数，β = ψ/(φ ρ)。B 用 6.3 节详细讨论过的匹配条件确定（即 s^(0)(T)、x^(0)(T)、y^(0)(T)、ei^(0)(T) 在 T → 0 时的值必须与 s^(0)(τ)、x^(0)(τ)、y^(0)(τ)、ei^(0)(τ) 在 τ → ∞ 时的值匹配）。已知 s^(0)(τ) ≡ 1, x^(0)(τ) = 1 − e^{−τ}, y^(0)(τ) ≡ 0, ei^(0)(τ) ≡ 0，故 O(1) 外层解的匹配条件为 s^(0)(T) → 1, x^(0)(T) → 1, y^(0)(T) → 0, ei^(0)(T) → 0 (T → 0)。这给出 B = 1；代入 (6.88) 后积分并用 (6.92) 的条件得

s^(0)(T) = (1 − β) / (1 − β e^{T (1 − 1/β) / (1+ρ)}) ei^(0)(T) = 1 − s^(0)(T) / β

情形 2：ρ = O(1), ψ = O(1), φ = O(ε)。假设 φ = O(ε) 时

s^(0)(T) = e^{−T/(1+ρ)}, ei^(0) = 0, ε ei^(1)(T) = (1 − e^{−T/(1+ρ)}) / β

仍用与内层解匹配得到。

在内外层解中都可以继续求 (6.72) 与 (6.84) 中 ε 高阶项。解会越来越复杂，但每一步方程都是线性的。对大多数实际用途，第一非零项已足够准确。

对所有时间的一致有效解：求得快瞬态与拟稳态阶段的解后，可用 Kevorkian 和 Cole（1996）所述的简单方法得到对所有 t ≥ 0 有效的复合解——即把内层解的第一项与外层解的对应项相加，再减去它们的公共部分（内层解 t → ∞ 的极限，即外层解 T → 0 的极限）。例如 s 的内层解为 s^(0)(τ) = 1，情形 2 的外层解为 s^(0)(T) = e^{−T/(1+ρ)}；上述两个极限都是 1，故复合解为

s0_comp = 1 + e^{−T/(1+ρ)} − 1 = e^{−t/ts} = e^{−ε (k−1+k2) t}

对其他解作同样处理，可得两组复合解（情形 1 与情形 2），对所有时间有效：

情形 1：

s0_comp(t) = (1 − β) / (1 − β e^{t (1−1/β) / ts}) ei0_comp(t) = (1 − s0_comp) / β x0_comp(t) = s0_comp (1 − ei0_comp) − e^{−t/tc} y0_comp(t) = 0 ε y1_comp(t) = (σ / (ψ(1+ρ))) [(e^{−ψ t/tc} − ψ e^{−t/tc}) / (ψ − 1) + s0_comp (1 − ei0_comp)]

情形 2：

s0_comp(t) = e^{−t/ts}, ei0_comp(t) = 0 ε ei1_comp(t) = (1 − s0_comp) / β x0_comp(t) = s0_comp − e^{−t/tc}, y0_comp(t) = 0 ε y1_comp(t) = (σ / (ψ(1+ρ))) [(e^{−ψ t/tc} − ψ e^{−t/tc}) / (ψ − 1) + s0_comp]

其中 β = ψ/(φ p)，σ、ρ、ψ 如 (6.66)。

注意区分情形 1 与情形 2 的重要参数是 β。β < 1 时情形 1 成立，β > 1 时情形 2 成立。这个 β 实际上就是 Tatsunami 等（1981）所称的 (1 + r)µ。上述表达式表明 β < 1 时 ε 一阶 ei → 1 (T → ∞)，而 β > 1 时 ε 一阶 s → 0 (T → ∞)。这直接关系到失活酶的量——即本节开头讨论的。

数值解与解析解的比较：求出无量纲系统的渐近近似解后，可与 Burke 等（1990）所得数值解比较以显示其准确性。Burke 等在量纲系统 (6.53)–(6.56) 上做数值解；为便于比较将无量纲浓度乘以相应的尺度因子后再绘图。复合解的前两项与数值解的比较见图 6.3，复合解在内层区段远比以前解更准确。图 6.4 给出中间浓度 X 与 Y 的数值解与复合解比较；每种 X、Y 都使用情形 2 复合解的第一项。这些中间结果比任何拟稳态方法都更准确，因为该方法纳入了到达拟稳态之前中间时间导数的变化。

上述结果表明，解析解是 (6.43) 所代表的自杀底物系统动力学的很好近似。它看起来并不比基本酶反应 (6.1) 复杂多少，但如我们所见，分析要复杂得多。这里发展的方法在估计中间浓度（X 与 Y）上特别有用——以前的分析无法做到。该方法最重要的结果或许是：以动力学参数解析地得到解；这些解可用 Waley（1985）和 Duggleby（1986）所描述的方法来估计参数。当方程是刚性的（小的参数乘以微分方程组中的导数项）时，解析解特别重要——因为数值解此时难以精确计算。

6.5 协同现象

在模型机制 (6.1) 中，1 分子酶与 1 分子底物结合——即酶有 1 个结合位点。许多酶有不止一个底物结合位点。例如血红蛋白（Hb）——红细胞中的携氧蛋白——对氧（O2）分子有 4 个结合位点。当一个酶分子在一处结合一个底物分子之后，还能在另一处结合另一个底物分子，这种酶与底物之间的反应被描述为协同（cooperative）。这种协同现象很常见。

另一种重要的协同行为是：一个具有多个结合位点的酶，一处结合一个底物分子会影响另一处结合其他底物分子的活性。这种不同且特异结合位点之间的间接相互作用称为变构（allostery）或变构效应（allosteric effect），表现出该效应的酶称为变构酶。若一个底物在一处结合增加了另一处的结合活性，则该底物为激活剂；若减少了活性则为抑制剂。变构反应动力学的详细数学分析可在 Murray（1977）一书中简略找到，在 Rubinow（1975）一书中更为详细；后者还给出了酶动力学的图论方法。

作为协同现象的实例，考虑一个酶有 2 个结合位点的情况，并计算等效的拟稳态近似和底物摄取函数。模型由酶分子 E 与底物分子 S 结合形成单结合的底物-酶复合物 C1 构成；C1 不仅能分解为产物 P 和酶 E，还能与另一个底物分子结合形成双结合的底物-酶复合物 C2；C2 分解为产物 P 和单结合复合物 C1。该模型的反应机制为

S + E ⇌ C1 → E + P S + C1 ⇌ C2 → C1 + P

其中 k 为相应的速率常数。

用小写字母表示浓度，对上述机制应用质量作用定律得

ds/dt = −k1 s e + (k−1 − k3 s) c1 + k−3 c2 dc1/dt = k1 s e − (k−1 + k2 + k3 s) c1 + (k−3 + k4) c2 dc2/dt = k3 s c1 − (k−3 + k4) c2 de/dt = −k1 s e + (k−1 + k2) c1 dp/dt = k2 c1 + k4 c2

初始条件 s(0) = s0, e(0) = e0, c1(0) = c2(0) = p(0) = 0。酶的守恒律由 (6.99) 中第 2、3、4 个方程相加并用初始条件得到：dc1/dt + dc2/dt + de/dt = 0 ⇒ e + c1 + c2 = e0。p(t) 方程仍解耦；求出 c1、c2 后可积分得到。因此利用 (6.101)，最终系统为

ds/dt = −k1 e0 s + (k−1 + k1 s − k3 s) c1 + (k1 s + k−3) c2 dc1/dt = k1 e0 s − (k−1 + k2 + k1 s + k3 s) c1 + (k−3 + k4 − k1 s) c2 dc2/dt = k3 s c1 − (k−3 + k4) c2

初始条件同 (6.100)。

与往常一样将系统无量纲化。若 e0/s0 ≪ 1，写

τ = k1 e0 t, u = s / s0, v1 = c1 / e0, v2 = c2 / e0 a1 = k−1 / (k1 s0), a2 = k2 / (k1 s0), a3 = k3 / k1 a4 = k−3 / (k1 s0), a5 = k4 / (k1 s0), e = e0 / s0

则 (6.102) 变为

du/dτ = −u + (u − a3 u + a1) v1 + (a4 + u) v2 = f(u, v1, v2) ε dv1/dτ = u − (u + a3 u + a1 + a2) v1 + (a4 + a5 − u) v2 = g1(u, v1, v2) ε dv2/dτ = a3 u v1 − (a4 + a5) v2 = g2(u, v1, v2)

加上初始条件 u(0) = 1, v1(0) = v2(0) = 0，构成适定的数学问题。

该问题与 6.5 节中分析的 Michaelis–Menten 问题 (6.13) 一样，在 0 < ε ≪ 1 时是奇异摄动问题。完整的内外层解可按 6.3 节所述方法类似求得（作者留作练习）。但这里关心的是：在实验可测量的时间区段（τ ≫ ε）中底物浓度 u 的摄取函数的形式。所以只需 O(1) 阶外层解——即在 (6.104)–(6.107) 中令 ε 项为零：

du/dτ = f(u, v1, v2) g1(u, v1, v2) = 0 g2(u, v1, v2) = 0

后两个方程是代数方程，解出 v1、v2：

v2 = a3 u v1 / (a4 + a5) v1 = u / (a1 + a2 + u + a3 u^2 (a4 + a5)^{−1})

代入 f(u, v1(u), v2(u)) 即得 u 的摄取方程（速率方程）：

du/dτ = f(u, v1(u), v2(u)) = −u [a2 + a3 a5 u (a4 + a5)^{−1}] / [a1 + a2 + u + a3 u^2 (a4 + a5)^{−1}] = −r(u) < 0

无量纲反应速率即为 r(u)。量纲上，利用 (6.103)，对 0 < e0/s0 ≪ 1 的 Michaelis–Menten 速率 R0(s0) 由 (6.108) 得

R0(s0) = |ds/dt|_{t=0} = e0 s0 / (k2 K'_m + k4 s0) · Km K'_m / (Km K'_m + K'_m s0 + s0^2)

其中 Km = (k2 + k−1) / k1, K'_m = (k4 + k−3) / k3。Km 和 K'_m 是机制 (6.98) 的 Michaelis 常数，相当于 (6.41) 中的 Michaelis 常数。

反应速率 R0(s0) 见图 6.5。若部分参数为零会有一个拐点：例如 k2 = 0 时从 (6.109) 可清楚看到，因为此时对小的 s0 有 R0 ∝ s0^2。这种协同行为的一个好例子是血红蛋白与氧的结合；实验测得的摄取曲线非常像图 6.5 中的下曲线。肌红蛋白（Mb）——红肌纤维中大量存在的一种蛋白——只有一个氧结合位点，其摄取是 Michaelis–Menten 形式（图 6.5 中也示出以作比较）。

当怀疑酶反应中存在协同现象时，常作 Hill 图。其基本假设是反应速率或摄取函数形式为

R0(S0) = Q s0^n / (Km + s0^n)

其中 n > 0 通常不是整数；这常称为 Hill 方程。解出 s0^n：

s0^n = R0 Km / (Q − R0) ⇒ n ln s0 = ln Km + ln (R0 / (Q − R0))

Hill 图是 ln (R0 / (Q − R0)) 对 ln s0 的图；其斜率给出 n，且在 Hill 方程为摄取动力学的有效描述时该斜率为常数。若 n < 1, n = 1 或 n > 1 分别称为负协同、零协同或正协同。尽管 Hill 方程作为 Michaelis–Menten 意义上的反应速率定量形式是合理的，但产生它的具体反应机制并不太现实：本质上它就是 (6.1)，但需要 E + nS 一步结合形成复合物。若 n 不是整数这不太可能（虽然可以作为一种化学计量形式）。若 n 是整数且 n ≥ 2，该反应是三分子或更高阶的——除非中间反应极快相当于把若干步压缩在一起，否则这种反应不会发生。

尽管在反应机制的含义上有这些缺点，Hill 方程等经验速率形式在建模中极其有用。毕竟我们希望从模型中得到对现象背后的动力学与机制的一些理解。第一步是找到一个在生物上合理且能定性描述行为的模型。细节的修正和完善留待以后。

6.6 自催化、激活与抑制

许多生物系统内建有反馈控制。这些控制非常重要，我们必须知道如何建模。下一章（生物振子）中会描述一个反馈控制必不可少的领域。代谢反馈控制系统的理论模型与动力学综述见 Tyson 和 Othmer（1978）。本节描述若干较重要的反馈控制类型。基本上，反馈指反应序列中某一步的产物对序列中其他反应步骤有影响——该影响通常是非线性的，可能是激活或抑制。下一章会给出具体反应机制的实际例子。

自催化（autocatalysis）是一种化学物质参与自身生产的过程。一个简单的教学例子是

A + X ⇌ 2X

其中 1 分子 X 与 1 分子 A 结合形成 2 分子 X。若 A 保持恒定浓度 a，对该反应应用质量作用定律得反应速率

dx/dt = k1 a x − k−1 x^2 ⇒ x(t) → xS = k1 a / k−1

其中 x = [X]，xS 是 t → ∞ 时的非零终态稳态。零稳态不稳定（从形式上即可看出）。该自催化反应表现出强反馈，其中"产物"抑制反应速率。显然须有某种逆反应（k−1 ≠ 0）。这正是第 1 章所讨论的 logistic 增长的化学等价。

若反应系统改为

A + X ⇌ 2X, B + X → C

即 X 在生成 C 时被消耗。该机制呈现一个简单分岔——下面说明。设 B 与 A 同为恒定浓度 a 与 b，则

dx/dt = k1 a x − k−1 x^2 − k2 b x = (k1 a − k2 b) x − k−1 x^2

这里 k1 a 是 x 的单位生产率，k2 b 是单位损失率。由 (6.114) 可看出：若 k1 a > k2 b 则稳态 x = 0 不稳定，x(t) → xS = (k1 a − k2 b) / k−1 > 0 (t → ∞)，且稳定；若 k1 a < k2 b 则 x = 0 稳定（不等式意味着损失率大于生产率故不奇怪）。数学上此时仍有另一个稳态但为负（不现实）且不稳定。该反应的简单分岔见图 6.6，图中以参数 k1 a − k2 b 给出稳态 xS；分岔发生在 k1 a − k2 b = 0 处，稳定性由一个稳态变为另一个。

为下一章（生物振子）作准备，经典 Lotka（1920）反应机制——他提出该机制作为假想模型振子——是自催化的另一例子：

A + X → 2X X + Y → 2Y Y → B

其中 A 维持恒定浓度 a。前两个反应是自催化的。质量作用定律给出

dx/dt = k1 a x − k2 x y dy/dt = k2 x y − k3 y

加上无量纲变量

u = k2 x / k3, v = k2 y / (k1 a), τ = k1 a t, α = k3 / (k1 a)

方程变为

du/dτ = u (1 − v) dv/dτ = α v (u − 1)

这是 3.1 节中详细讨论的 Lotka–Volterra 方程 (3.4)。解 u、v 关于时间是周期的；但如我们在那里所见，它们是结构不稳定的。

在几乎所有生物过程中我们都不知道发生的具体生化反应。然而我们常常知道改变某种已知反应物或改变运行条件的定性影响。因此在建模时，将这些已知的定性行为直接纳入模型机制通常更有效率也更有启发性。正是这种模型机制在解释、揭示基本过程中如此有用，并在宽广的生物医学问题中作出有用的预测。既然我们知道如何将反应序列表示为微分方程系统，便可构造直接将这些定性行为结合进浓度微分方程的模型。此时微分方程系统本身就构成了模型。

设有一个微分方程系统——其模型可经渐近方法（如上节所讨论）简化为由以下无量纲机制控制的关键两个要素

du/dτ = a / (b + v) − c u = f(u, v) dv/dτ = d u − e v = g(u, v)

其中 a、b、c、d、e 为正常数。该模型的生物解释是：u 激活 v（通过 d u 项），u 与 v 都按其浓度一阶动力学降解——这是 −c u 与 −e v 项。a/(b + v) 项表示 v 对 u 的生产的负反馈——v 增加时 u 的生产减少，进而间接减少 v 自身。v 越大，u 的生产越小。这是反馈抑制的一个例子。

可容易证明该机制存在一个稳定的正稳态。相关稳态 (u0, v0) 是

f(u0, v0) = g(u0, v0) = 0 ⇒ v0 = d u0 / e, u0^2 + (e b / d) u0 − (a e) / (c d) = 0

的正解。方程 (6.116) 与第 3 章详细分析过的形式完全相同。线性稳定性由线性化 Jacobian（或反应矩阵或稳定性矩阵，相当于第 3 章的群落矩阵）的特征值 λ 决定，特征方程为

det[∂f/∂u − λ, ∂f/∂v; ∂g/∂u, ∂g/∂v − λ] |_{u0, v0} = 0

即 λ^2 + (c + e) λ + [c e + (c d u0) / (b + v0)] = 0 ⇒ Re λ < 0

故 (u0, v0) 线性稳定。它也是全局吸引的稳态——在 (u, v) 平面上可构造一个矩形约束集使其边界上向量 (du/dt, dv/dt) 指向内部。

若干具体模型系统被提议作为某些基本生物现象（如振子行为、胚胎发育中的图式形成、哺乳动物的毛皮图样等）的机制。我们在后续章节详细研究其中一些。本节简要地看两个。

Thomas（1975）机制基于涉及底物氧和尿酸在酶 uricase 存在下反应的具体反应。氧（v）与尿酸（u）的经验速率方程的无量纲形式可写为

du/dt = a − u − ρ R(u, v) = f(u, v) dv/dt = α (b − v) − ρ R(u, v) = g(u, v) R(u, v) = u v / (1 + u + K u^2)

其中 a、b、α、ρ、K 为正常数。基本上 u 与 v 以恒定速率 a 与 αb 供应，按其浓度一阶降解，且都以 ρ R(u, v) 的速率在该反应中被消耗。R(u, v) 的形式表现为底物抑制：对给定的 v，u 小时 R(u, v) = O(u v)，对 u 线性；当 u 大时 R(u, v) = O(v/(K u))。因此 u 小时 R 随 u 增加而增加，u 大时 R 随 u 增加而减少。这就是底物抑制的含义。参数 K 衡量抑制的强度。由图 6.7（给出 R(u, v) 关于 u 的函数）可见摄取速率对小 u 呈 Michaelis–Menten 形式，在 u = 1/√K 处达最大值，然后随 u 增加而下降。R(u, v) 取最大值时的浓度以及实际最大速率都随抑制增强（即 K 增大）而减小。

如第 3 章中相互作用的种群模型那样，在 (u, v) 相平面上画出反应动力学的零线总是有信息量的。这里 (6.117) 的零线为

f(u, v) = 0 ⇒ v = (a − u) (1 + u + K u^2) / (ρ u) g(u, v) = 0 ⇒ v = α b (1 + u + K u^2) / [ρ u + α (1 + u + K u^2)]

其示意图见图 6.8。视参数不同，可有一个或三个正稳态。虽然这些零线是针对具体的底物抑制机制，但它们对一般底物抑制模型是相当典型的，特别是 f = 0 零线（也参见图 6.9）。稳态稳定性的问题将在下章详细且较为一般地讨论。在现阶段，从 (u, v) 相平面上看零线即可对稳定性得到直观印象。考虑图 6.8 中存在三个稳态 P1、P2、P3 的情形，具体先看 P1(u1, v1)。沿通过 P1 的线 v = v1 移动，跨过 f = 0 零线时记下 f(u, v1) 的符号：f = 0 零线左侧 f > 0，右侧 f < 0。所以在 P1 处 v = v1（常数）时 ∂f/∂u < 0。若这是 u 的非耦合标量方程 du/dt = f(u, v1)，则按 1.1 节讨论，P1 局部线性稳定。但 (6.117) 中 u 方程并未解耦——耦合可能有不稳定的影响。

继续考虑 P1；用同样论据沿通过 P1 的线 u = u1 跨过 g = 0 零线。可见 ∂g/∂v < 0，故在 P1 处 du/dt = g(u1, v) 局部 ∂g/∂v < 0——同样理由支持 P1 线性稳定的直观判断。这些论据在下章中严格发展，证实直观判断是对的。类似地可直观推断 P3 也是稳定的。对 P2 用同样符号论据：v = v2 时由图 6.8 可见 ∂f(u, v2)/∂u > 0，故 P2 不稳定。当存在单一稳态 S 时需仔细分析（见第 7 章）。

不作任何分析也可以清楚：必然存在某些参数范围使得稳态唯一，另一些范围使得有三个稳态。因此做有信息量的分析需要确定这两种情形各自的参数域。虽然原则上简单——联立代数方程 f(u, v) = g(u, v) = 0 求正稳态——但代数上通常困难，需要数值完成。这样的分析会产生一些有趣结果，我们将在 6.7 节详细讨论。

另一模型机制——代数上比 Thomas 系统 (6.117) 简单——是如下假想但生物上合理的反应方案

du/dt = a − b u + u^2 / (v (1 + K u^2)) = f(u, v) dv/dt = u^2 − v = g(u, v)

其中 a、b、K 为常数。这是激活剂（u）—抑制剂（v）系统，是 Gierer 和 Meinhardt（1972）所提模型的无量纲形式。它已被用于多种建模情境，我们在后续章节中指出。本处 u 的自催化生产由 u^2 / (v (1 + K u^2)) 项表达，但当 u 大时该项饱和到 1/(K v)。抑制剂 v 由 u 激活（第二个方程），但它抑制其激活剂的生产——因为 u^2 / (v (1 + K u^2)) 随 v 增加而下降。(6.118) 中的零线 f = 0 和 g = 0 见图 6.9。注意图 6.8 与图 6.9 零线之间的定性相似性，尤其在稳态附近以及大 u 时——我们稍后考虑其含义。下一章中我们引入其他反应系统，第 8 章中详细讨论一个具有相当实验重要性和生物相关性的具体系统。

对一般系统

du/dt = f(u, v) dv/dt = g(u, v)

若 ∂g/∂u > 0 则 u 是 v 的激活剂；若 ∂f/∂v < 0 则 v 是 u 的抑制剂。依赖于具体动力学，反应物可以仅在某一浓度或参数范围内表现为激活剂。因此存在许多分岔现象的可能性，这些现象在生物上有重要意义（我们将在本书后面看到）。

由于反应动力学与相互作用种群之间在数学上的相似，我们也预期观察到如 3.8 节中讨论的阈值现象。事实确实如此——(6.117) 模型系统在参数使稳态位于 S 或 S′（图 6.10）时表现出类似的阈值行为。3.8 节的分析在此直接适用。

现在我们可以开始构建模型反应以纳入多种反应动力学行为——自催化、激活与抑制等——因为定性上知道需要什么。举例：设有细胞以对化学物质 S 的局部浓度水平做出反应，激活某基因使细胞产生产物 G；设产物以饱和方式自催化生产并按一阶动力学降解——即按其浓度线性降解。用小写字母表示浓度，一个定性体现所有这些要求的关于 g 的速率方程为

dg/dt = k1 s + k2 g^2 / (k3 + g^2) − k4 g = f(g)

其中 k 为正常数。该模型具有一些有用的生物开关性质，在第 3 卷（Volume II）第 3 章中讨论产生生物空间图式的模型时我们将考虑并使用它。

至此已清楚：对 n 个反应的反应动力学研究产生 n 阶一阶微分方程组

dui/dt = fi(u1, …, un), i = 1, …, n

这与相互作用种群模型中出现的同一类一般系统（具体为第 3 章的 (3.43)）形式上相同。在那里我们只关心非负解，这里也一样——因为 u(t) 是浓度向量。所有分析稳态（即 fi(u1, …, un) = 0 的解）稳定性的方法都可适用。极限环、阈值现象等的条件在此也都成立。

反应物之间的相互作用细节与相互作用种群之间的相互作用细节在形式与动机上当然不同。在生物系统中由于所需的微分方程模型阶数一般更高因而更复杂。但如我们所见，酶催化剂的存在常常为此提供生物上合理的依据，使得阶数大大降低。例如产生如下无量纲方程的系统

dui/dt = fi(u1, …, un), i = 1, 2 εi dui/dt = fi(u1, …, un), i = 3, …, n 0 < εi ≪ 1, i = 3, …, n

对实际几乎所有目的，在 ε 足够小时化简为二阶系统

dui/dt = fi(u1, u2, u3(u1, u2), …, un(u1, u2)), i = 1, 2

其中 i = 3, …, n 的 fi(u1, …, un) = 0 是代数方程，解出 un≥3 作为 u1 与 u2 的函数。这正是拟稳态近似向高阶系统推广的形式，也正是我们广泛研究二反应物动力学模型的理由。数学上最后一个方程是 εi → 0（对所有 i）时 (6.122) 非奇异解的 O(1) 渐近系统。生物上这一般就是我们所需要的，因为通常是机制相对长时间的行为主导生物发展。

6.7 多稳态、蘑菇与孤立圈

图 6.8 表明可能存在多个正稳态。从单一稳态到三个稳态的转变发生在模型某参数通过分岔值时。图 6.11 示意了发生这种转变的典型情形。例如参考图 6.9 与 (6.118) 的动力学，稳态将定性地如图 6.11(a) 那样随参数变化——抑制参数 K 起 p 的作用。

设现在随参数 k 变化，us 关于 p 的曲线以这样的方式变化：在 k 的某个范围内曲线的定性形式如图 6.11(b)。固定 k，对 p1 < p < p2 有三个稳态，分别在 BC、CD、DE 三条分支上。这等价于图 6.8 的三稳态情形。由上节讨论可知位于 CD 分支上的稳态线性不稳定；这一点将在下章证明。

图 6.11(b) 中 (us, p) 图的形式提示了滞后（hysteresis）的可能性（回顾 1.1 节）。如所述，设位于 ABC 与 DEF 分支上的稳态稳定。现将参数 p 从 p < p1 缓慢增加到 p > p2。p 达到 p2 前，us 仅沿 ABC 分支相应值增加。p 通过 p2 时，us 突然变化跳到 EF 分支上；p 继续增加时 us 由该分支相应值给出。现将 p 缓慢减小。此时 us 保持在下方 FED 分支上直到 p 达到 p1——该分支上解稳定。在 p1 处 us 突然跳到上方 BA 分支。这是典型的滞后环。p 增加时路径为 ABCEF，p 减小时路径为 FEDBA。

蘑菇：若 (us, p) 的变化不是图 6.11(a) 那种，而是一种 us 随 p 增大而增大的简单形式（图 6.12(a)），则到三稳态的转变如图所示。可以想象：图 6.12(a) 中的简单曲线演化为图 6.12(b) 那种蘑菇状——在 p 空间中有两个多稳态区——是容易理解的。图 6.12(b) 的蘑菇状 (us, p) 关系在两个不同的 p 范围内有三条稳态。位于 CD 与 GH 分支上的稳态不稳定。有两个等价的滞后环（如图 6.11(b)），即 BCED 与 IHFG。

孤立圈：图 6.12(c) 所示情形——一个分离的"断片"区域——是图 6.12(b) 情形的自然推广。这种解行为称为孤立圈（isola）。在 CDG 分支上的解不稳定。该情形所代表的物理情境与蘑菇情形相当不同。首先没有通常的滞后：us 在参数 p 从 p < p1 增加到 p > p2 时只是保持在 ABIJ 分支上；在多稳态区域 p1 < p < p2 中返回扫掠也仍保持在该分支上。孤立圈是解分支的孤立闭曲线，只能由非线性方程产生。

仍参考图 6.12(c)，若 us 位于 BI 分支上，则只有给 us 一个有限扰动使之移入 DFG 分支上稳态的吸引域，才能移到另一个稳定分支 DFG 上。各种可能情形现已清楚。

通过简单地操作曲线（如我们刚才做的）就可以预测相当复杂的解行为。多稳态的出现并不难想象——只要有合适的动力学。Dellwo 等（1982）给出了一般理论，以解析方式描述一类孤立圈的结构——即那些在某参数趋于临界值时缩为一点的孤立圈。立即出现的问题是：例如孤立圈是否能在真实世界存在。孤立圈已在多种真实实际情形中发现，包括化学反应；早期综述见 Uppal 等（1976），其他参考文献见 Gray 和 Scott（1986）的论文。

Gray 和 Scott（1983, 1986）提出了一个简单模型动力学系统，能展现——除其他之外——带蘑菇和孤立圈的多稳态：它涉及连续搅拌釜式反应器（CSTR）中的自催化。该系统包括两个反应物 X、Y（浓度 x、y）的假想反应，具体机制见图 6.13。过程包括三分子自催化步骤 X + 2Y → 3Y，描述该过程的特定方程组为

dx/dt = k0 (x0 − x) − k1 x y^2 − k3 x dy/dt = k0 (y0 − y) + k1 x y^2 + k3 x − k2 y

其中 k 为正常速率常数。合适的无量纲化为

u = x / x0, v = y / x0, t* = t k1 x0^2 c = y0 / x0, a = k0 / (k1 x0^2), b = k3 / (k1 x0^2), d = k2 / (k1 x0^2)

则方程变为（省略星号以简化记号）无量纲系统

du/dt = a (1 − u) − u v^2 − b u = f(u, v) dv/dt = a (c − v) + u v^2 + b u − d v = g(u, v)

现含四个无量纲参数 a、b、c、d。

此处只关心稳态 us、vs——即 f(u, v) = g(u, v) = 0 的解。稍作代数得

us (1 + c − us)^2 = a [(1 + d/a)^2 (1 − us) − b us] (a + d)^2 / a^2

这是一个三次方程

us^3 − 2 (1 + c) us^2 + [(1 + c)^2 + (a + d)^2 / a + b (a + d)^2 / a^2] us − (a + d)^2 / a = 0

由于该三次方程有三次符号变化，由 Descartes 符号规则（见附录 B）可能有三个正解。对参数的大值或小值某些解析解可渐近地求出。但完整的图形只能数值地得到（Gray 和 Scott 1986）。典型结果如图 6.14 所示。对该反应及其复杂行为以及解析与数值结果的良好综述见 Gray（1988）。

当然数学上可以构造越来越复杂的解行为并假想展现它们的反应。所以此时关键的问题是：是否有真实反应过程展现这些有趣现象（蘑菇和孤立圈）。已经证明无机碘酸盐-亚砷酸反应在适当条件下具有所需的动力学。Ganapathisubramanian 和 Showalter（1984）已令人信服地证实了这一点；他们的模型与实验结果将在下面描述。尽管这不是酶或生物反应，但它表明具有蘑菇和孤立圈解行为的真实反应机制确实存在。鉴于生物过程的丰富性与复杂性，若这种反应系统不存在于生物医学中将是难以置信的。带着这种信念，下面我们描述该无机反应的元素并给出相关实验结果。

碘酸盐-亚砷酸反应：双稳定性、蘑菇、孤立圈

连续流搅拌釜式反应器中的碘酸盐-亚砷酸反应可由两个复合反应描述：

IO3^- + 5 I^- + 6 H^+ → 3 I2 + 3 H2O I2 + H3AsO3 + H2O → 2 I^- + H3AsO4 + 2 H^+

总反应 (6.128) + 3 × (6.129) 为

IO3^- + 3 H3AsO3 → I^- + 3 H3AsO4

反应 (6.128) 的速率比 (6.129) 慢，因此它是整个过程 (6.129) 的限速步骤。设 (6.128) 的速率为 R，实验上已确定其经验形式为

R = −d[IO3^-]/dt = (k1 + k2 [I^-]) [I^-] [H^+]^2 [IO3^-]

定量描述连续流搅拌釜式反应器中碘酸盐-亚砷酸反应的简单模型机制，由关于碘化物 I^- 与碘酸盐 IO3^- 的速率方程（加上适当的流项与衰减项）构成：

d[I^-]/dt = R + k0 [I^-]_0 − (k0 + k3) [I^-] d[IO3^-]/dt = −R + k0 [IO3^-]_0 − (k0 + k3) [IO3^-]

其中 k0、k3 为正常数，[I^-]_0 与 [IO3^-]_0 分别为流入浓度，R 由 (6.131) 给出。

若写

X = [I^-], Y = [IO3^-] X0 = [I^-]_0, Y0 = [IO3^-]_0 k1 = k1 [H^+]^2, k2 = k2 [H^+]^2

则稳态 Xs、Ys 由以下方程的解给出

0 = R + k0 X0 − (k0 + k3) X 0 = −R + k0 Y0 − (k0 + k3) Y R = (k1 + k2 X) X Y

由此可得关于 Xs 的三次多项式

k2 (k0 + k3) Xs^3 + [k1 (k0 + k3) − k2 k0 (X0 + Y0)] Xs^2 + [(k0 + k3)^2 − k1 k0 (X0 + Y0)] Xs − k0 (k0 + k3) X0 = 0

k1 与 k2 的值已由实验确定，X0、Y0 与 [H^+] 可给定，故由 (6.134) 可确定 k1 与 k2。图 6.15 给出由三次方程 (6.135) 数值求得的正稳态碘化物浓度 Xs 关于 k0 + k3 的函数（对应不同 k3 值）。

将上述碘酸盐-亚砷酸反应模型与完整反应系统比较可得良好的定量结果。图 6.15 表明蘑菇与孤立圈多稳态行为是可能的。证明这种行为存在还需要实验确认。Ganapathisubramanian 和 Showalter（1984）也已做到这一点；他们的实验结果见图 6.16。注意这些实验结果与该碘酸盐-亚砷酸反应模型机制所得结果的比较。图 6.16 的结果清楚地显示图 6.11 与 6.12 提示的各种滞后行为。

本章个人批注

本章是 Murray 第一卷中最技术化的一章，几乎没有任何生物学应用的具体细节——纯粹是数学建模技巧的展示。但恰恰因为这一点，它的方法论价值对本书后续章节至关重要。第 7 章（生物振子）会大量复用这里的奇异摄动与拟稳态分析，第 11、12 章（反应扩散）会复用这些无量纲化思想。

6.1–6.3 的 Michaelis–Menten 拟稳态分析是教学性的：作者先给出传统拟稳态假设，再严格用奇异摄动法说明该假设何时成立、何时失效。关键的小参数是 ε = e0/(s0 + Km)，不是 e0/s0——这是 Segel 和 Slemrod 修正传统分析的核心洞见。许多教科书错误地以 e0/s0 ≪ 1 作为拟稳态条件，但实际只需 e0/(s0 + Km) ≪ 1；当 Km 足够大时即使酶与底物同量级，条件仍能成立。这一点对建模的实用意义是：拟稳态假设比一般想象的应用面更广。

6.4 的自杀底物分析比 6.3 复杂得多——因为有两条平行的反应路径（到产物与到失活酶）以及一个关键参数 β = ψ/(φ p) ≈ (1 + r)µ 决定是底物先耗尽还是酶先失活。值得记住的是情形 1（β < 1）对应"酶全失活"、情形 2（β > 1）对应"底物全耗尽"。β 实际上就是 Tatsunami 等（1981）的 (1 + r)µ——这是一个跨研究的同一参数的不同命名。

6.5 的协同现象给出 Hill 图的来源。Hill 方程 R0 = Q s0^n / (Km + s0^n) 是个方便的唯象工具，但机制上不现实（n 步同时结合）。这种"机制不现实但拟合数据好"的形式在生物建模中反复出现——作者自己也承认"经验速率形式极其有用"。Hemoglobin（n ≈ 2.8）vs Myoglobin（n = 1）的对比是经典例子。

6.6 给出两种激活剂-抑制剂模型的雏形：Thomas（底物抑制，R(u, v) = uv / (1 + u + K u^2)）和 Gierer–Meinhardt（u^2 / (v(1 + K u^2))）。零线的定性图（图 6.8、6.9）在一定参数范围内可形成 S 形曲线，从而产生 3 个正稳态——这是后续讨论分岔、滞后、多稳态的几何基础。Gierer–Meinhardt 模型是 Murray 第二卷中形态发生理论的核心。

6.7 给出蘑菇和孤立圈——多稳态的更复杂拓扑。Gray–Scott 模型（CSTR 中 X + 2Y → 3Y）是经典范例。重点是：数学上构造多稳态不难，难的是在真实反应中验证。碘酸盐-亚砷酸反应是少数几个实验证实蘑菇和孤立圈都存在的反应——这是一个无机化学反应，与生物无关，但作者以此暗示这种行为可能在生物系统中也存在。

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

第 5 章把差分方程模型用于心理学（婚姻互动），是一次跨域应用的展示。第 6 章回归到本书主线——生化反应的数学建模。本章是第 3 章（相互作用种群）、第 4 章（温度依赖性别决定）之后的"工具章"，建立了拟稳态近似、奇异摄动分析、零线分析等核心技术。

承接第 5 章结尾的"建模生成理论"方法论，第 6 章提供了具体的数学技巧库：6.1–6.4 的 Michaelis–Menten 与自杀底物分析展示奇异摄动法；6.5 的协同现象给出 Hill 图；6.6 的激活剂-抑制剂系统把反应动力学还原为第 3 章已分析过的二维 ODE 系统形式。这些技巧并非孤立——它们为第 7 章（生物振子与开关）做完整准备：振子分析中需要的就是拟稳态降阶 + 线性稳定性分析 + 极限环判定。第 7 章的 Goodwin 振子、Tyson 振子都基于本章建立的二维可激活-可抑制反应机制框架。

向后看，第 6 章末尾的二维系统形式 f(u, v)、g(u, v) 实际上是第 7–9 章反复出现的"标准形式"。6.7 的多稳态、蘑菇、孤立圈为第 8 章 BZ 反应、第 9 章耦合振子铺垫。第 11 章（反应扩散）将本章的"反应项"扩展为"反应 + 扩散"项，零线方法会被保留但加上扩散维度。