第二章：反应扩散系统与空间斑图形成

学习笔记

2.1 生物学中斑图的角色

位置信息概念

发育生物学中关于斑图形成与分化的一个重要现象学概念是由 Wolpert (1969) 提出的位置信息 (positional information) 理论。该理论认为细胞被预先编程为对某种化学物质（形态发生素 morphogen）的浓度产生反应，并据此分化为不同类型的细胞，如软骨细胞等。

位置信息的核心观点： - 化学预构型 (chemical prepattern) 观点将胚胎发育过程分离为多个步骤 - 第一步是创造形态发生素浓度的空间预构型 - 一旦预构型建立，形态发生就是"从属过程" - 位置信息不依赖于建立空间预构型的具体机制

Turing 的反应扩散理论

Turing (1952) 在其经典论文中提出了形态建成的化学预构型理论，即反应扩散理论。该理论认为在特定条件下，化学物质可以通过反应和扩散产生稳态的非均匀空间斑图。反应扩散理论如今已拥有大量文献，是一个独立的研究领域。

关键洞察： 扩散通常被视为稳定过程，但 Turing 指出在某些条件下扩散可以导致不稳定——这是一个新颖且深刻的概念。

类比说明

为了直观理解扩散如何导致不稳定，Turing 提出了一个形象的比喻： - 将一片干草比作反应场 - 蚱蜢代表抑制剂 - 火代表激活剂 - 如果抑制剂（蚱蜢）的扩散系数远大于激活剂（火），火蔓延的区域将受到限制 - 如果两者扩散系数相同，则无法形成空间斑图

2.2 反应扩散 (Turing) 机制

基本方程

反应扩散机制的控制方程为：

\[\frac{\partial \mathbf{c}}{\partial t} = \mathbf{f}(\mathbf{c}) + D \nabla^2 \mathbf{c} \qquad(2.1)\]

其中： - $\mathbf{c}$ 是形态发生素浓度向量 - $\mathbf{f}$ 代表反应动力学 - $D$ 是正常数扩散系数对角矩阵

对于双化学物种系统 A(r,t) 和 B(r,t)：

\[\frac{\partial A}{\partial t} = F(A,B) + D_A \nabla^2 A \qquad(2.2a)$$ $$\frac{\partial B}{\partial t} = G(A,B) + D_B \nabla^2 B \qquad(2.2b)\]

Turing 的核心思想

Turing (1952) 的思想： 如果在没有扩散（$D_A = D_B = 0$）时，A 和 B 趋向于线性稳定的均匀稳态，那么在特定条件下（$D_A \neq D_B$），扩散驱动的失稳可以产生空间非均匀斑图。

典型反应扩散系统

1. Schnakenberg (1979) 系统：

\[F(A,B) = k_1 - k_2 A + k_3 A^2 B \qquad(2.3a)$$ $$G(A,B) = k_4 - k_3 A^2 B \qquad(2.3b)\]

2. Gierer-Meinhardt (1972) 激活子-抑制子系统：

\[F(A,B) = k_1 - k_2 A + \frac{k_3 A^2}{B} \qquad(2.4a)$$ $$G(A,B) = k_4 A^2 - k_5 B \qquad(2.4b)\]

3. Thomas (1975) 底物抑制系统：

\[F(A,B) = k_1 - k_2 A - H(A,B) \qquad(2.5a)$$ $$G(A,B) = k_3 - k_4 B - H(A,B) \qquad(2.5b)$$ $$H(A,B) = \frac{k_5 AB}{k_6 + k_7 A + k_8 A^2} \qquad(2.5c)\]

无量纲化

引入特征长度尺度 $L$，定义：

\[u = A\sqrt{\frac{k_3}{k_2}}, \quad v = B\sqrt{\frac{k_3}{k_2}}, \quad t^* = \frac{D_A t}{L^2}, \quad \mathbf{x}^* = \frac{\mathbf{x}}{L}\]

可得无量纲形式：

\[u_t = \gamma (a - u + u^2 v) + \nabla^2 u \qquad(2.7a)$$ $$v_t = \gamma (b - u^2 v) + d \nabla^2 v \qquad(2.7b)\]

其中： - $d = D_B/D_A$（扩散系数比） - $\gamma = L^2 k_2/D_A$（尺度参数） - 参数 γ 的含义： - (i) $\gamma^{1/2}$ 与一维空间域的线性大小成正比 - (ii) γ 表示反应项的相对强度 - (iii) 增加 γ 等价于减小扩散系数比 d

局部激活与侧向抑制

反应扩散系统的核心空间概念是局部激活与侧向抑制。对于激活子-抑制子系统，抑制子的扩散速度必须快于激活子：

\[u_t = a - bu + \frac{u^2}{v} + \nabla^2 u \qquad(2.8)$$ $$v_t = u^2 - v + d\nabla^2 v\]

这个概念可以追溯到 1885 年 Ernst Mach 的 Mach 带研究。

2.3 扩散驱动失稳的通用条件

数学 formulation

考虑一般形式：

\[u_t = \gamma f(u,v) + \nabla^2 u, \quad v_t = \gamma g(u,v) + d\nabla^2 v \qquad(2.11)\]

边界条件为零通量： $$\frac{\partial u}{\partial n} = \frac{\partial v}{\partial n} = 0, \quad \mathbf{r} \in \partial B$$

线性稳定性分析

步骤 1：无空间变化的均匀稳态稳定性

令 $w = \begin{pmatrix} u - u_0 \\ v - v_0 \end{pmatrix}$，线性化得：

\[w_t = \gamma A w \qquad(2.15)\]

其中 $A = \begin{pmatrix} f_u & f_v \\ g_u & g_v \end{pmatrix}$。

特征方程： $$\lambda^2 - \gamma(\text{tr} A)\lambda + \gamma^2 |A| = 0 \qquad(2.17)$$

线性稳定性条件（$\text{Re}\lambda < 0$）：

\[\text{tr} A = f_u + g_v < 0, \quad |A| = f_u g_v - f_v g_u > 0 \qquad(2.19)\]

步骤 2：包含扩散的线性稳定性

\[w_t = \gamma A w + D \nabla^2 w \qquad(2.20)\]

其中 $D = \begin{pmatrix} 1 & 0 \\ 0 & d \end{pmatrix}$

定义空间特征值问题：

\[\nabla^2 W + k^2 W = 0, \quad \frac{\partial W}{\partial n} = 0 \qquad(2.21)\]

寻找形式解 $w(\mathbf{r},t) = \sum_k c_k e^{\lambda t} W_k(\mathbf{r})$，代入得特征方程：

\[\lambda^2 + \lambda[k^2(1+d) - \gamma(f_u+g_v)] + h(k^2) = 0 \qquad(2.23)\]

其中： $$h(k^2) = dk^4 - \gamma(d f_u + g_v)k^2 + \gamma^2 |A| \qquad(2.24)$$

扩散驱动失稳条件

对于扩散驱动失稳，需要： 1. 无扩散时稳定：$\text{Re}\lambda(k^2=0) < 0$ ✓ 2. 有扩散时不稳定：$\text{Re}\lambda(k^2 \neq 0) > 0$ 对某些 $k$

由于 $f_u + g_v < 0$ 且 $k^2(1+d) > 0$，唯一可能是 $h(k^2) < 0$。

关键条件：

\[d f_u + g_v > 0 \Rightarrow d \neq 1 \qquad(2.25)\]

完整条件集：

\[f_u + g_v < 0, \quad f_u g_v - f_v g_u > 0 \qquad(2.31a)$$ $$d f_u + g_v > 0, \quad (d f_u + g_v)^2 - 4d(f_u g_v - f_v g_u) > 0 \qquad(2.31b)\]

色散关系

特征值 $\lambda(k^2)$ 称为色散关系。当 $h(k^2) < 0$ 时，$\lambda$ 为正，对应的模式线性不稳定。不稳定波数的范围为 $k_1^2 < k^2 < k_2^2$：

\[k_{1,2}^2 = \frac{\gamma}{4d}\left[(df_u+g_v) \mp \sqrt{(df_u+g_v)^2 - 4d|A|}\right] \qquad(2.29)\]

最快增长模式在 $k = k_m$ 处：

\[k_m^2 = \frac{\gamma(df_u+g_v)}{2d} \qquad(2.26)\]

2.4 斑图启动的详细分析

Schnakenberg 系统详解

以 Schnakenberg 系统为例：

\[u_t = \gamma(a - u + u^2 v) + u_{xx} \qquad(2.32a)$$ $$v_t = \gamma(b - u^2 v) + d v_{xx} \qquad(2.32b)\]

均匀稳态：

\[u_0 = a + b, \quad v_0 = \frac{b}{(a+b)^2} \qquad(2.33)\]

在稳态处的偏导数：

\[f_u = \frac{b-a}{a+b}, \quad f_v = (a+b)^2 > 0 \qquad(2.34a)$$ $$g_u = \frac{-2b}{a+b} < 0, \quad g_v = -(a+b)^2 < 0 \qquad(2.34b)$$ $$f_u g_v - f_v g_u = (a+b)^2 > 0 \qquad(2.34c)\]

Turing 空间

由条件 (2.31) 可得：

\[0 < b - a < (a+b)^3 \qquad(2.35a)$$ $$d(b-a) > (a+b)^3 \qquad(2.35b)$$ $$[d(b-a) - (a+b)^3]^2 > 4d(a+b)^4 \qquad(2.35c)\]

一维模式

对于域 $x \in (0,p)$，边界条件为零通量：

\[W_n(x) = A_n \cos\left(\frac{n\pi x}{p}\right), \quad k = \frac{n\pi}{p} \qquad(2.37)\]

不稳定模式数由参数范围决定。

二维模式

对于矩形域 $0 < x < p, 0 < y < q$：

\[W_{n,m}(x,y) = C_{n,m} \cos\frac{n\pi x}{p} \cos\frac{m\pi y}{q} \qquad(2.42)\]

波数 $k^2 = \pi^2(n^2/p^2 + m^2/q^2)$

规则平面镶嵌斑图

Christopherson (1940) 发现的解可铺满平面的规则镶嵌包括：

六边形： $$\psi(x,y) = \cos kx + \cos\left(k\frac{\sqrt{3}y + x}{2}\right) + \cos\left(k\frac{\sqrt{3}y - x}{2}\right) \qquad(2.47)$$

正方形： $$\psi(x,y) = \cos kx + \cos ky \qquad(2.48)$$

菱形： $$\psi(x,y) = \cos kx + \cos[k(x\cos\phi + y\sin\phi)] \qquad(2.49)$$

一维条纹： $$\psi(x,y) = \cos kx \qquad(2.50)$$

异质性函数

Berding (1987) 引入了异质性函数来量化空间非均匀程度：

\[H = \int_0^1 (U^2 + V^2)dx \geq 0 \qquad(2.52)\]

物理意义：H 表征了空间斑图的结构程度。

2.5 色散关系、Turing 空间、尺度与几何效应

色散关系的基本性质

原型色散关系具有两个基本特征： 1. 空间均匀态（$k=0$）是稳定的 2. 存在有限的波长带（不稳定模式窗口）

临界波长：

\[\omega_c = \frac{2\pi}{k_c} = 2\pi\left[\frac{d_c}{\gamma^2(f_u g_v - f_v g_u)}\right]^{1/4} \qquad(2.54)\]

Turing 空间分析

对于 Schnakenberg 系统，Turing 空间由参数曲线定义：

曲线 1： $$a > \frac{1}{2}u_0(1-u_0^2), \quad b = \frac{1}{2}u_0(1+u_0^2) \qquad(2.57)$$

曲线 2： $$a < \frac{1}{2}u_0\left(1-\sqrt{2}-\frac{u_0^2}{\sqrt{d}}\right) \qquad(2.60)$$

关键发现

临界扩散比： $d_c = 3 + 2\sqrt{2}$，当 $d > d_c$ 时 Turing 空间开始存在
模式选择： 通过改变 γ（域大小）可以选择激发特定模式
鲁棒性： Thomas 系统和 Schnakenberg 系统具有较大的 Turing 空间，而 Gierer-Meinhardt 系统较小

尺度与几何效应

在有限域中，可能的波数是离散的。只有当色散关系的不稳定模式带包含允许的离散波数时，斑图才能形成。

参数空间 $(d, \gamma)$： 该平面上可以划分出不同模式对应的区域。

2.6 模式选择与色散关系

三种斑图启动方式

1. 随机扰动启动 - 初始扰动包含所有模式 - 最快增长模式（色散曲线峰值）最终主导

2. 生境生长启动 - 域生长缓慢 - 第一个被激发的模式保持主导地位 - 存在"间隙"区域，斑图暂时消失

3. 行波启动 - 斑图从域的一端开始形成 - 最终斑图的波长由行波前端的运动决定

行波启动的数学分析

考虑线性系统 $Jw = 0$，寻找 $w \propto e^{ikx + \lambda t}$ 形式解。

通解： $$w(x,t) = \int A(k) \exp[ikx + \lambda(k)t] dk \qquad(2.70)$$

在 $x/t$ 为 O(1) 的渐近区域，积分由最速下降法给出。

生境生长动力学

情况 A： 色散关系随 γ 变化存在"间隙" - 斑图形成是离散过程 - 中间存在空间均匀态

情况 B： 色散关系随 γ 连续变化 - 不稳定模式数量连续增加 - 模式连续演化

2.7 单物种模型的斑图生成：云杉芽虫模型

单物种反应扩散方程

\[u_t = f(u) + D\nabla^2 u \qquad(2.76)\]

云杉芽虫动力学

\[f(u) = ru\left(1 - \frac{u}{q}\right) - \frac{u^2}{1+u^2} \qquad(2.85)\]

该动力学具有三个正稳态：$u_1$（ refuge，稳定）、$u_2$（不稳定）、$u_3$（ outbreak，稳定）。

临界域大小

从零稳态失稳的临界大小：

\[L_c = \pi\sqrt{\frac{D}{f'(0)}} \qquad(2.80)\]

维持 outbreak 的临界大小 $L_0$： 可通过数值方法确定。

解析近似方法

将问题重新标度到 $x \in (0,1)$：

\[DU'' + L^2 f(U) = 0, \quad U(0) = U(1) = 0 \qquad(2.86)\]

通过几何分析可确定 $L_c$ 和 $L_0$。

公式汇总表

编号	公式名称	方程	条件/说明
(2.1)	反应扩散基本方程	$\frac{\partial \mathbf{c}}{\partial t} = \mathbf{f}(\mathbf{c}) + D \nabla^2 \mathbf{c}$	Turing (1952)
(2.2)	双物种系统	$\frac{\partial A}{\partial t} = F + D_A \nabla^2 A, \quad \frac{\partial B}{\partial t} = G + D_B \nabla^2 B$	$D_A \neq D_B$
(2.7)	无量纲形式	$u_t = \gamma(a-u+u^2v)+\nabla^2 u, \quad v_t = \gamma(b-u^2v)+d\nabla^2 v$	Schnakenberg 系统
(2.19)	无扩散稳定性	$\text{tr}A = f_u+g_v < 0, \quad	A
(2.23)	色散关系	$\lambda^2 + \lambda[k^2(1+d)-\gamma(f_u+g_v)] + h(k^2) = 0$	特征方程
(2.24)	h(k²) 定义	$h(k^2) = dk^4 - \gamma(df_u+g_v)k^2 + \gamma^2	A
(2.25)	关键条件	$df_u + g_v > 0 \Rightarrow d \neq 1$	扩散比必须不等
(2.26)	最快增长波数	$k_m^2 = \frac{\gamma(df_u+g_v)}{2d}$	h取最小值处
(2.29)	不稳定波数范围	$k_{1,2}^2 = \frac{\gamma}{4d}\left[(df_u+g_v) \mp \sqrt{(df_u+g_v)^2-4d	A
(2.31)	Turing 失稳完整条件	$f_u+g_v<0,\quad f_ug_v-f_vg_u>0,\quad df_u+g_v>0,\quad (df_u+g_v)^2>4d(f_ug_v-f_vg_u)$	必要充分条件
(2.33)	Schnakenberg 稳态	$u_0 = a+b, \quad v_0 = \frac{b}{(a+b)^2}$
(2.34)	稳态偏导数	$f_u=\frac{b-a}{a+b},\quad f_v=(a+b)^2,\quad g_u=\frac{-2b}{a+b},\quad g_v=-(a+b)^2$
(2.35)	Schnakenberg Turing 条件	$0<b-a<(a+b)^3,\quad d(b-a)>(a+b)^3$	参数约束
(2.37)	一维本征函数	$W_n = A_n\cos(n\pi x/p)$	零通量边界
(2.42)	二维本征函数	$W_{n,m} = C_{n,m}\cos(n\pi x/p)\cos(m\pi y/q)$	矩形域
(2.47)	六边形解	$\psi = \cos kx + \cos(k\frac{\sqrt{3}y+x}{2}) + \cos(k\frac{\sqrt{3}y-x}{2})$	平面镶嵌
(2.52)	异质性函数	$H = \int_0^1 (U^2+V^2)dx$	Berding (1987)
(2.54)	临界波长	$\omega_c = 2\pi\left[\frac{d_c}{\gamma^2(f_ug_v-f_vg_u)}\right]^{1/4}$
(2.76)	单物种 RD 方程	$u_t = f(u) + D\nabla^2 u$
(2.80)	临界域大小	$L_c = \pi\sqrt{\frac{D}{f'(0)}}$	零稳态失稳
(2.85)	芽虫动力学	$u_t = ru(1-u/q) - \frac{u^2}{1+u^2} + Du_{xx}$

关键概念总结

Turing 失稳机制

均匀稳态在无扩散时线性稳定
扩散驱动使系统失稳，产生空间斑图
必要条件：扩散系数不等 ($d \neq 1$)，且 $df_u + g_v > 0$
充分条件：完整的 Turing 条件 (2.31)

模式选择机制

色散关系确定不稳定模式范围
域大小 (γ) 影响可包含的模式
初始条件影响最终斑图极性
行波启动产生不同的波长选择

生物应用

动物毛皮斑图：豹斑、斑马条纹、长颈鹿斑块
发育生物学：软骨形成、羽毛原基
生态学：害虫控制、生物入侵

重要发现

非常小的动物应该颜色均匀（域太小无法形成斑图）
非常大的动物也应该颜色均匀（斑图太细密）
中等大小的动物呈现最丰富的斑图
斑图的时间进程可以被"冻结"并被生长扭曲

参考文献

Turing, A.M. (1952). The chemical basis of morphogenesis. Phil. Trans. R. Soc. B, 237, 37-72.
Schnakenberg, J. (1979). Simple chemical reaction systems with limit cycle behaviour. J. Theor. Biol., 81, 389-400.
Gierer, A. & Meinhardt, H. (1972). A theory of biological pattern formation. Kybernetik, 12, 30-39.
Thomas, D. (1975). Artificial enzyme membranes. J. Coll. Int. Sci., 51, 228-239.
Murray, J.D. (1982). Parameter space in reaction-diffusion theory and morphogenesis. J. Math. Biol., 15, 89-96.
Arcuri, P. & Murray, J.D. (1986). Pattern sensitivity to boundary conditions in reaction-diffusion models. J. Math. Biol., 25, 89-109.
Zhu, J. & Murray, J.D. (1995). Analysis of pattern formation in reaction-diffusion models. Phys. Rev. E, 52, 177-195.

本笔记基于 J.D. Murray 的《Mathematical Biology II》第二章内容整理

编号	公式名称	方程	条件/说明
(2.1)	反应扩散基本方程	\(\frac{\partial \mathbf{c}}{\partial t} = \mathbf{f}(\mathbf{c}) + D \nabla^2 \mathbf{c}\)	Turing (1952)
(2.2)	双物种系统	\(\frac{\partial A}{\partial t} = F + D_A \nabla^2 A, \quad \frac{\partial B}{\partial t} = G + D_B \nabla^2 B\)	\(D_A \neq D_B\)
(2.7)	无量纲形式	\(u_t = \gamma(a-u+u^2v)+\nabla^2 u, \quad v_t = \gamma(b-u^2v)+d\nabla^2 v\)	Schnakenberg 系统
(2.19)	无扩散稳定性	$\text{tr}A = f_u+g_v < 0, \quad	A
(2.23)	色散关系	\(\lambda^2 + \lambda[k^2(1+d)-\gamma(f_u+g_v)] + h(k^2) = 0\)	特征方程
(2.24)	h(k²) 定义	$h(k^2) = dk^4 - \gamma(df_u+g_v)k^2 + \gamma^2	A
(2.25)	关键条件	\(df_u + g_v > 0 \Rightarrow d \neq 1\)	扩散比必须不等
(2.26)	最快增长波数	\(k_m^2 = \frac{\gamma(df_u+g_v)}{2d}\)	h取最小值处
(2.29)	不稳定波数范围	$k_{1,2}^2 = \frac{\gamma}{4d}\left[(df_u+g_v) \mp \sqrt{(df_u+g_v)^2-4d	A
(2.31)	Turing 失稳完整条件	\(f_u+g_v<0,\quad f_ug_v-f_vg_u>0,\quad df_u+g_v>0,\quad (df_u+g_v)^2>4d(f_ug_v-f_vg_u)\)	必要充分条件
(2.33)	Schnakenberg 稳态	\(u_0 = a+b, \quad v_0 = \frac{b}{(a+b)^2}\)
(2.34)	稳态偏导数	\(f_u=\frac{b-a}{a+b},\quad f_v=(a+b)^2,\quad g_u=\frac{-2b}{a+b},\quad g_v=-(a+b)^2\)
(2.35)	Schnakenberg Turing 条件	\(0<b-a<(a+b)^3,\quad d(b-a)>(a+b)^3\)	参数约束
(2.37)	一维本征函数	\(W_n = A_n\cos(n\pi x/p)\)	零通量边界
(2.42)	二维本征函数	\(W_{n,m} = C_{n,m}\cos(n\pi x/p)\cos(m\pi y/q)\)	矩形域
(2.47)	六边形解	\(\psi = \cos kx + \cos(k\frac{\sqrt{3}y+x}{2}) + \cos(k\frac{\sqrt{3}y-x}{2})\)	平面镶嵌
(2.52)	异质性函数	\(H = \int_0^1 (U^2+V^2)dx\)	Berding (1987)
(2.54)	临界波长	\(\omega_c = 2\pi\left[\frac{d_c}{\gamma^2(f_ug_v-f_vg_u)}\right]^{1/4}\)
(2.76)	单物种 RD 方程	\(u_t = f(u) + D\nabla^2 u\)
(2.80)	临界域大小	\(L_c = \pi\sqrt{\frac{D}{f'(0)}}\)	零稳态失稳
(2.85)	芽虫动力学	\(u_t = ru(1-u/q) - \frac{u^2}{1+u^2} + Du_{xx}\)