第 2 章:空间图样形成与反应扩散系统(Spatial Pattern Formation with Reaction Diffusion Systems)
2.1 图样在生物学中的作用(Role of Pattern in Biology)
胚胎学关心胚胎从受精到出生的形成与发育过程。胚胎发育是顺序性的,依循一份"发育总图"(ground plan)展开——在人类这份总图大致在第 5 周时即已铺就。形态发生(morphogenesis),即胚胎学中我们主要关心的部分,研究的是图样与形态的发育。发育总图究竟如何建立,以及产生空间图样(从而指定各种器官)的机制是什么,至今仍然未知。
本章与接下来若干章主要讨论能够在多种发育情境下作为可能图样形成过程而被提出的机制。波动现象当然也产生空间图样,但那是时空图样。本章关心的是稳态的、空间异质的图样。本章引入并分析反应扩散图样形成机制,主要关注发育生物学。§2.7 关注图样形成的生态学侧面,并由此提出一种可能的害虫控制策略;其数学分析虽与胚胎学情形不同,但直接适用于许多胚胎学情境。
作者想回答(或更现实地说,是想得到一些启示)的问题不胜枚举。例如,细菌能形成复杂图样、狼群标记领地,这些过程之间是否存在共同的图样形成原理?空间图样在生物医学中无处不在,理解它们如何形成无疑是最基本的科学挑战之一。在本书余下篇幅里,我们将研究多种多样的图样形成机制。
受精后细胞开始分裂。当分裂达到一定规模后,核心问题就是:均一的细胞团如何在空间上组织起来,从而使顺序性的发育过程得以推进。细胞按其在空间组织中的位置而分化(生物学意义上的分化)。细胞在胚胎中还会迁移。后一种现象是形态发生的重要元素,催生了一种在第 6 章详细讨论的图样与形态生成的新方法。
面对生物学中图样的丰富、多样与精美,很难不感到震撼与入迷。图 2.1 仅仅给出四例。这些图样以及无数其他图样究竟如何被"画"出来,至今仍无定论——虽然在多个方向上已取得显著进展,例如对果蝇胚胎早期图样化、黏菌空间图样化以及第 5 章讨论的细菌图样的研究。仅就图 2.1 的几幅图样所提出的图样问题,就已经相当多样。
作为图 2.1(c) 的脚注,需提一下蛾的触角。触角非常高效地收集一种被叫做 bombykol 的信息素化学分子——这是雌蛾释放以吸引雄蛾的化学气味物质。以蚕蛾为例:不能飞的雄蛾能够从远达一公里外探测到雌蛾释放的信息素,并沿浓度梯度向雌蛾方向移动。这种触角的过滤效率(它们事实上在"数"分子数)提出了一个与本书讨论的不同但有趣的数学生物学图样问题——如何设计这种过滤触角使之效率最高?这一具体问题——一个有趣的流体力学与扩散问题——Murray (1977) 已详细讨论过。
图样与形态在生物学中的根本重要性是不言自明的。我们的认识是:无论我们在动物界观察到什么图样,几乎可以肯定的是,产生该图样的过程至今仍是未知的。图样形成研究经常因模型中没有包含基因而受到批评。但同样的批评可以指向任何对一个复杂系统进行相对简单抽象建模的做法。应当记住,图样与形态的生成,特别是在发育过程中,通常与基因组水平相距甚远。当然基因起着关键作用,机制必须受基因控制;然而基因本身并不能"画"出图样,它们只是为图样生成提供了一份"蓝图"或"配方"。许多演化中的图样根本无法仅由基因信息预先决定。生物学的另一重大问题就是:基因信息如何被物理地转译为所需的图样与形态。发育生物学中大量的研究——无论实验的还是理论的——都致力于弄清早期发育中图样与形态形成的潜在机制。
Wolpert (1969,综述见 1971, 1981) 提出了一种被称为"位置信息"(positional information)的图样形成与分化现象学概念。他设想:细胞对某种化学物质(形态发生素,morphogen)的浓度具有预编程的反应,并据此分化成不同种类的细胞,比如软骨细胞。Wolpert (1977) 的总览性论文对动物中图样与形态的发育及其位置信息场景的概念与应用,给出了非常清晰且非技术性的说明。尽管它是一种现象学方法,没有涉及实际机制,但已催生了大量启发性的实验研究,其中很多与鸡胚中肢芽软骨图样化以及其他鸟胚(如鹌鹑和珍珠鸡)羽毛图样的发育相关(参看 Richardson et al. 1991 及其中的参考文献)。在发育中检索"位置信息"的文献,会得到数量惊人的参考文献。尽管这是一个简单而吸引人的概念,并已使我们对发育某些方面的认识取得重要进展,但它不是一个机制。
胚胎发生的"化学预图样"(chemical prepattern)观点把发育过程分成若干步骤:第一个关键步骤是建立形态发生素浓度的空间图样。"形态发生素"之名因此得名——它驱动形态发生。位置信息的概念依赖于化学预指定,使细胞能在化学浓度坐标系中读出自身位置,并据此分化、改变形状或迁移。因此,一旦预图样建立,形态发生即成为一种"从属"过程。位置信息不依赖于设置形态发生素空间预图样的具体机制。本章关心的是把反应扩散模型作为生成生物学图样的可能机制来研究。
形态发生的化学反应扩散理论基础——亦称形态发生的反应扩散理论——最早由 Turing (1952) 在那篇经典论文中提出。反应扩散理论如今已有海量文献,本身就是一个研究领域。
面对动物形态的复杂性,"位置信息"概念必然意味着细胞对"形态发生素图"有着极为复杂的解读。无论这是否真是位置信息、抑或只是细胞以某种方式对化学浓度细微差异作出反应,都无关紧要——当我们想到细胞有多么复杂即可。反应扩散所能形成的图样尺度可以非常小——图 2.11 给出的实验图样即为佐证。细胞复杂性的一个粗略量度是把细胞 DNA 分子(脱氧核糖核酸)每比特信息的重量约为 \(10^{-22}\),与之相比,电子束成像约为 \(10^{-10}\),磁带约为 \(10^{-5}\)。最复杂最紧凑的计算机芯片也完全不能与细胞相提并论。
从理论模型中可以引出一个重要观点:任何图样都自带其历史。试考虑以下简单的工程类比(Murray et al. 1998)——它类比的是我们在试图理解一个生物学过程中的角色。说一座桥需要一千吨钢、少了会太弱、多了会过于刚性,这是一回事;而指导工人如何最有效地把构件拼装起来,则是另一回事。在形态发生中,例如可以设想:参与组织形成与形变的细胞具备足够的能力——只要给它们正确的原料和初始指令,便可被说服去构建任何我们想要的元素。这是许多追寻完整可预测理解的人所抱的希望。然而很有可能:所有这些精密细胞活动的整体效应,会对发育过程中事件的序列极为敏感。作为科学家,我们应该关注的是:如何利用我们与"工人"交流的有限机会,从而把实验引导向可接受的结果。
针对任何具体的生物图样形成过程所提出的单个模型,乃至所有这些模型加在一起,都算不上一个完整模型。在一些被广泛研究的问题中(如发育肢芽中的图样化),每个模型都只照亮了过程的不同侧面,我们现在能说出一份完整模型所必须包含的重要概念元素。这些研究帮助我们看到自身知识的不足,并指出了富有成果的实验方向。这些理论构建对实验界的冲击是其关键检验。
作为本节小结,必须再次强调:对图样化场景的数学描述——包括现象学描述——并非解释。这一点虽被普遍接受,却常被遗忘。
2.2 反应扩散(Turing)机制(Reaction Diffusion (Turing) Mechanisms)
Turing (1952) 提出:在一定条件下,化学反应与扩散可以共同作用,产生化学或形态发生素浓度的稳态异质空间图样。在 Volume I 第 11 章中我们推导了反应扩散机制的控制方程——即式 (11.16)——这里以如下形式考虑:
其中 \(\mathbf{c}\) 是形态发生素浓度向量,\(\mathbf{f}\) 代表反应动力学,\(D\) 是正常数扩散系数的对角矩阵。本章主要关注两种化学物质 A(r, t) 和 B(r, t) 的模型。方程组形式为
其中 \(F\) 和 \(G\) 是反应动力学,通常为非线性。
Turing (1952) 的思想既简单又深刻。他说,如果在无扩散时(实际上 \(D_A = D_B = 0\))A 和 B 趋于一个线性稳定的均匀稳态,那么在 \(D_A \neq D_B\) 的某些条件下,空间非齐次图样可以通过扩散驱动的失稳而演化出来。扩散通常被视为稳定化过程——这正是此概念的新颖之处。要直观地理解扩散如何能够失稳,请考虑下面这个虽然不切实际但富有启发性的类比。
考虑一片干燥的草地,其中有许多蚱蜢,它们在感到温暖时会通过排汗产生大量水分。现在假设草在某个点被点燃,一个火焰前锋开始传播。我们可以把蚱蜢视为"抑制剂",把火视为"激活剂"。如果没有水分来扑灭火焰,火就会简单地蔓延过整片草地,造成一片均匀的焦黑区。设想:蚱蜢暖和到一定程度后能产生足够的水分浸湿草,于是当火焰蔓延到这样一块预先浸湿的区域时,草就不再燃烧。空间图样的场景如下:火开始蔓延——它是"反应物"之一(激活剂),具有"扩散"系数 \(D_F\)。当蚱蜢(抑制性"反应物")在火焰前锋之前感受到它的来临时,便迅速跳到火焰前方更远的地方——即它们的"扩散"系数 \(D_G\) 远大于 \(D_F\)。蚱蜢随后大量出汗,产生足够水分以阻止火蔓延至该浸湿区域。如此,焦黑区域被限制在某个有限空间内,其大小取决于"反应物"(火与蚱蜢)的"扩散"系数以及各种"反应"参数。如果不是只有一个初始火源,而是在整片草地上随机散布着许多火源,可以想象此过程最终将产生一个空间异质的稳态分布——焦黑区域与未焦黑区域交错,以及蚱蜢的空间分布——围绕每个火源都会发生上述过程。如果蚱蜢与火焰前锋以同样的速度"扩散",则这种空间图样无法演化。如何构造其他类比是显然的;§2.3 会给出另外的例子,Murray (1988) 在《科学美国人》的文章里也举了一个。
下面我们用反应与扩散的形态发生素语言来描述这一过程,推导反应动力学与扩散系数上需要满足的条件,并给出我们可能预期到的空间图样的类型。这里先简短记下两个特别简单的假想系统与一个实验实现的例子,它们都能满足 Turing 图样形成系统的条件。当然现在已有很多其他系统被用于空间图样化的研究,它们有不同程度的实验可实现性。鉴于 Volume I 第 8 章和上一章对 Belousov–Zhabotinskii 反应的广泛讨论,应当特别注意它。尽管已发现许多其他真实反应系统,它仍然是主要的实验系统。
最简单的系统是 Schnakenberg (1979) 反应,在 Volume I 第 7 章中已讨论。参考系统形式 (2.2),其动力学为
其中 \(k\) 为正常数速率常数。这里 A 通过 \(F(A,B)\) 中的 \(k_3 A^2 B\) 项自催化产生。这是原型反应扩散系统之一。另一个有影响力的系统是 Gierer 与 Meinhardt (1972) 提出的激活剂-抑制剂机制,自那以后被广泛研究和使用。该系统在 Volume I 第 6 章中已讨论:
这里 A 是激活剂,B 是抑制剂。\(k_3 A^2 / B\) 项同样是自催化的。Koch 与 Meinhardt (1994) 综述了 Gierer–Meinhardt 反应扩散系统在复杂结构的生物图样形成中的应用,并给出了该特定模型及其变体的广泛应用文献目录。
Thomas (1975) 在实验上研究过的真实底物抑制系统(在 Volume I 第 6 章中也详细描述)有
这里 A 与 B 分别是底物氧气与酶尿酸氧化酶的浓度。底物抑制通过 H 项中的 \(k_8 A^2\) 体现。因为 H 项为负,它有助于减少 A 和 B;减少的速率在 A 充分大时被抑制。基于 Field–Körös–Noyes (FKN) 模型动力学的反应扩散系统(参看 Volume I 第 8 章)由于具有实验验证理论的潜力,是一个特别重要的例子;具体文献会在下面相关位置给出。
在评论能够产生图样的反应动力学的类型之前,我们必须对 (2.2) 进行无量纲化,反应动力学取 (2.3) 至 (2.5) 的形式。作为示例,我们用 (2.2) 与 (2.3) 中的 \(F\)、\(G\) 详细推导——因其代数简洁且我们在 Volume I 第 7 章对其作过详细分析。引入 \(L\) 作为典型长度尺度,并设
无量纲化后的反应扩散系统(为代数方便去掉星号)为
其中 \(f\)、\(g\) 由这些方程定义。我们本可以通过设 \(\gamma^{1/2} r\) 和 \(\gamma t\) 替换 \(r\)、\(t\) 把 \(\gamma\) 吸收到新的长度与时间尺度中——这等价于定义长度尺度 \(L\) 使得 \(\gamma = 1\),即 \(L = (D_A / k_2)^{1/2}\)。我们保留 (2.7) 的特定形式,原因会在后面变得清楚——既为下一节的分析,也为后续章节中的应用。
对反应动力学 (2.4) 和 (2.5) 适当的无量纲化(见习题 1)给出
其中 \(a, b, \alpha, \rho, K\) 为正常数。如果在第一个系统的激活剂中加入自抑制,\(f\)、\(g\) 取为
其中 \(k\) 衡量自抑制的强度;可参看 Volume I 第 6 章 §6.7。
Murray (1982) 详细讨论了上述每个系统并就其作为图样生成器的相对优劣得出结论;他给出了一套研究任意两物种反应扩散系统的系统解析方法。对于大多数用作图样形成(解析)演示的反应扩散机制,最简单的——即 (2.7)——是已考虑过的系统中最为稳健的一个,并且幸运地也是最容易研究的。
所有这类反应扩散系统都可以无量纲化并缩放为如下通式:
其中 \(d\) 是扩散系数之比,\(\gamma\) 可有下列任一解释。
(i) \(\gamma^{1/2}\) 与一维空间域的线性大小成正比。在二维,\(\gamma\) 与面积成正比。这个含义非常重要,我们将在 §2.5 和第 3 章看到。
(ii) \(\gamma\) 代表反应项的相对强度。例如,\(\gamma\) 增大可能意味着反应序列中某个限速步骤的活性增强。
(iii) \(\gamma\) 增大也可以等价地被视为扩散系数比 \(d\) 的减小。
这种通式的特别优点在于:(a) 无量纲参数 \(\gamma\)、\(d\) 比有量纲参数有更宽的生物学解释;(b) 当我们考虑参数空间中特定空间图样出现的区域时,结果可以方便地在 \((\gamma, d)\) 空间里展示。这一方面由 Arcuri 与 Murray (1986) 利用。
系统 (2.2) 能否产生 Turing 型空间图样,关键取决于反应动力学 \(f\)、\(g\) 与 \(\gamma\)、\(d\) 的值。零等值线(null clines)的具体形式提供了基本的初步信息。图 2.2 展示了 \(f\)、\(g\) 由 (2.7)–(2.9) 定义时典型的零等值线。
尽管这些动力学有不同的化学动机与推导过程,但它们都等价于某种激活-抑制解释;只要配合反应物扩散的不均匀性,便能产生空间图样。空间激活-抑制的概念在 Volume I 第 11 章 §11.5 中已详细讨论,起源于一种积分方程公式:参看那里 (11.41) 式。在下一节我们可以看到,关于图样生成的动力学,关键之处包含在零等值线的形式以及它们在稳态邻域内如何相交。前一幅图展示了两种广义类型。Figures 2.2(b)、(c) 和 (d) 中的稳态邻域类似,代表一类;Figure 2.2(a) 是另一类。
这里应说明:还存在其他重要的零等值线类型我们不予考虑,例如具有多于一个正稳态的零等值线;Volume I 第 7 章中已讨论过这类动力学。具有这种动力学的反应扩散系统能生成更为复杂的空间图样:这里初始条件尤其重要。我们也不再讨论扩散系数为空间依赖或浓度/种群依赖的系统;这些在生态学情境中很重要(回顾第 1 章关于转基因生物扩散的讨论)。Volume I 第 11 章中简要考虑了密度依赖扩散的情形。后面我们会讨论一个扩散系数空间依赖的重要应用——在解剖学真实脑结构中模拟脑肿瘤扩散。
在模型构建中,把机制的图样动力学用示意图表示,并约定记号以表明自催化、激活、抑制、降解和不等扩散,常常是有用且直观的辅助手段。如果这样做,例如对 (2.10) 中第一个 (2.8) 式给出的激活剂-抑制剂动力学,我们可以采用图 2.3(a) 所示的约定。
不同的扩散系数(这里 \(d > 1\))的效果是把"局部激活、长程抑制"(local activation and lateral inhibition)的原型空间概念展示出来,如图 2.3(b) 与图 2.4(b) 所示。该一般性概念已在 Volume I 第 11 章中介绍。正是这种普适的空间行为对于空间图样化是必要的:蚱蜢与火的类比是一个明显的例子——火是局部激活,蚱蜢提供长程抑制。直观上很清楚:抑制剂的扩散系数必须大于激活剂的扩散系数。
局部激活与长程抑制的概念相当古老,至少可以追溯到 1885 年 Ernst Mach 提出的 Mach 带宽现象。这是一种当明暗带并置时产生的视觉错觉。图 2.4 给出这种现象的示意图,并展示了基于此的 Hermann 错觉的例子。
2.3 扩散驱动失稳的一般条件:线性稳定性分析与空间图样的演化(General Conditions for Diffusion-Driven Instability: Linear Stability Analysis and Evolution of Spatial Pattern)
反应扩散系统表现出"扩散驱动失稳"(有时也叫 Turing 失稳),如果其均匀稳态在无扩散时对小扰动稳定,而在有扩散时对小幅空间扰动不稳定。生物学中的"失稳"一词常出现在生态学语境下——均匀稳态对小扰动失稳,种群通常表现出某种时间振荡行为。我们这里关心的失稳是另一种完全不同的类型。驱动空间非齐次失稳的主要过程是扩散:机制决定了演化出来的空间图样。图样或模式如何被选择是分析的一个重要方面——这是本章(以及后续章节)讨论的主题。
下面推导对通式 (2.10) 来说,扩散驱动失稳以及空间图样启动的必要且充分条件。为把问题在数学上确定下来,需要边界条件与初始条件。这里取零通量边界条件与给定初始条件。数学问题由下式定义
其中 \(\partial B\) 是反应扩散区域 \(B\) 的封闭边界,\(\mathbf{n}\) 是 \(\partial B\) 的单位外法向。选择零通量边界条件有若干理由。最主要的一条是:我们关心的是图样的自组织;零通量条件意味着无外部输入。如果对 \(u\)、\(v\) 施加固定边界条件,空间图样化就可能是边界条件的直接后果——这将在 §2.7 的生态学问题中看到。§2.4 将对 (2.7) 给出的动力学在一维和二维具体情境下进行分析。
(2.11) 相关的均匀稳态 \((u_0, v_0)\) 是下列方程组的正解:
由于我们关心的是扩散驱动失稳,我们感兴趣的是该稳态关于纯空间依赖的小扰动的线性失稳。因此在没有空间变化时,均匀稳态必须是线性稳定的——我们先确定满足此条件的条件。这些条件在 Volume I 第 3 章中已推导;这里为方便起见并保持记号完整,我们再简要重导一遍。
无空间变化时 \(u\)、\(v\) 满足
按 Volume I 第 3 章同样方式在稳态 \((u_0, v_0)\) 附近线性化,令
(2.13) 对 \(|\mathbf{w}|\) 小的情形变为
其中 \(A\) 是稳定性矩阵。除非另作说明,此后 \(f\)、\(g\) 的偏导数都取在稳态处的值。我们寻找形如
的解,\(\lambda\) 是特征值。\(\mathbf{w} = 0\) 线性稳定当且仅当 \(\text{Re} \lambda < 0\),因为此时 \(\mathbf{w} \to 0\) 当 \(t \to \infty\)。把 (2.16) 代入 (2.15) 可知 \(\lambda\) 是
的根,从而
线性稳定(即 \(\text{Re} \lambda < 0\))在下列条件下得到保证:
由于 \((u_0, v_0)\) 是动力学参数的函数,这些不等式因而对参数施加了某些约束。注意对图 2.2 中的所有情形,在稳态邻域内有 \(f_u > 0\)、\(g_v < 0\);对图 2.2(a) 有 \(f_v > 0\)、\(g_u < 0\),而对图 2.2(b) 到 (d) 有 \(f_v < 0\)、\(g_u > 0\)。所以 \(\text{tr} A\) 与 \(\det A\) 可正可负:这里我们只关心满足 (2.19) 的条件与参数范围。
现在考虑完整的反应扩散系统 (2.11),同样在稳态附近线性化,结合 (2.14) 即 \(\mathbf{w} = 0\),得到
为在边界条件 (2.11) 下求解此方程组,我们首先定义 \(W(\mathbf{r})\) 为下列空间特征值问题的时间无关解:
其中 \(k\) 是特征值。例如,若区域是一维的,设 \(0 \leq x \leq a\),则 \(W \propto \cos(n\pi x / a)\)(\(n\) 为整数);它在 \(x = 0\) 与 \(x = a\) 处满足零通量条件。这种情况下特征值为 \(k = n\pi / a\)。于是 \(1/k = a / (n\pi)\) 度量波状图样:特征值 \(k\) 称为波数,\(1/k\) 与波长 \(\omega\) 成正比;此例中 \(\omega = 2\pi / k = 2a / n\)。此后我们在这一语境下把 \(k\) 称为波数。对于有限区域,由于 \(n\) 为整数,可能的波数是离散的。
设 \(W_k(\mathbf{r})\) 是波数 \(k\) 对应的特征函数。每个特征函数满足零通量边界条件。由于问题是线性的,我们寻找 (2.20) 的形如下列的解:
其中常数 \(c_k\) 由初始条件关于 \(W_k(\mathbf{r})\) 的 Fourier 展开确定。\(\lambda\) 是决定时间增长的特征值。把该形式代入 (2.20) 与 (2.21) 并消去 \(e^{\lambda t}\),对每个 \(k\) 有
我们要求 \(W_k\) 的非平凡解,因此 \(\lambda\) 由下列特征多项式的根决定:
把 \(A\)、\(D\) 用 (2.15)、(2.20) 的形式代入,可得 \(\lambda(k)\) 作为波数 \(k\) 的函数,是下列方程的根:
稳态 \((u_0, v_0)\) 线性稳定当且仅当 (2.23) 的两个解都满足 \(\text{Re} \lambda < 0\)。我们已经施加了稳态在没有任何空间效应时稳定的约束——即 \(\text{Re} \lambda(k^2 = 0) < 0\)。此情况下 (2.23) 化为 (2.17),而 \(\text{Re} \lambda < 0\) 的要求即 (2.19)。若稳态对空间扰动失稳,我们要求对某个 \(k \neq 0\) 有 \(\text{Re} \lambda(k) > 0\)。这可以在下列任一情况下发生:(2.23) 中 \(\lambda\) 的系数为负,或 \(h(k^2) < 0\) 对某个 \(k \neq 0\)。由 (2.19) 知 \(f_u + g_v < 0\),而 \(k^2 (1 + d) > 0\) 对所有 \(k \neq 0\),因此 \(\lambda\) 的系数
所以 \(\text{Re} \lambda(k^2)\) 能为正的仅有可能是 \(h(k^2) < 0\) 对某个 \(k\) 成立。从 (2.23) 的解可立即看出:
由于 (2.19) 要求 \(\det A > 0\),要使 (2.23) 中 \(h(k^2) < 0\) 唯一可能是 \((d f_u + g_v) > 0\)。而由 (2.19) 知 \((f_u + g_v) < 0\),这蕴含 \(d \neq 1\) 且 \(f_u\)、\(g_v\) 必须异号。因此除 (2.19) 外还需
由图 2.2 中的反应动力学可知 \(f_u > 0\)、\(g_v < 0\),于是 (2.19) 的第一个条件与 (2.24) 的最后一个不等式要求扩散系数比 \(d > 1\)。例如对激活剂-抑制剂机制 (2.8),这意味着抑制剂必须扩散得比激活剂更快——如上所述。
(2.24) 的不等式是 \(\text{Re} \lambda > 0\) 的必要而非充分条件。要使 \(h(k^2)\) 对某个非零 \(k\) 为负,最小值 \(h_{\min}\) 必须为负。由 (2.23) 对 \(k^2\) 求导可得
因此 \(h(k^2) < 0\) 对某个 \(k^2 \neq 0\) 成立的条件是
在分叉处 \(h_{\min} = 0\),要求 \(\det A = (d f_u + g_v)^2 / (4 d)\),因此对固定的动力学参数,临界扩散系数比 \(d_c (> 1)\) 是
的相应根。临界波数 \(k_c\) 由 (2.28) 给出
图 2.5(a) 显示 \(h(k^2)\) 随 \(k^2\) 变化的几种情形。
只要 \(h(k^2) < 0\),(2.23) 就有一个正的 \(\lambda\) 解,对应的波数范围就是使 \(h < 0\) 的相同范围。由 (2.23) 与 \(d > d_c\),不稳定波数 \(k_1^2 < k^2 < k_2^2\) 的范围由 \(h(k^2) = 0\) 的两个零点 \(k_1^2\) 与 \(k_2^2\) 给出:
图 2.5(b) 给出了典型的 \(\lambda(k^2)\) 对 \(k^2\) 的图。\(\lambda = \lambda(k^2)\) 这个表达式称为"色散关系"。我们将在下两节中更详细地讨论色散关系的重要性与用途。注意,在不稳定范围内,\(\text{Re} \lambda(k^2) > 0\) 在由 (2.25) 与 \(d > d_c\) 给出的波数 \(k_m\) 处取得极大。这意味着在 (2.22) 中对 \(\mathbf{w}\) 的求和中存在一个最快增长模式;这是我们接下来要利用的一个属性。
若考虑 (2.22) 给出的解 \(\mathbf{w}\),当 \(t\) 增大时主导贡献来自 \(\text{Re} \lambda(k^2) > 0\) 的模式,因为所有其他模式随时间指数衰减。从图 2.5 或 (2.29) 的解析形式可确定范围 \(k_1^2 < k^2 < k_2^2\)(此范围内 \(h(k^2) < 0\),从而 \(\text{Re} \lambda(k^2) > 0\)),于是由 (2.22) 对大 \(t\) 有
色散关系的分析与图谱因此极具信息量——它们立即表明哪些特征函数(即哪些空间图样)是线性不稳定并随时间指数增长的。必须牢记,对于有限域特征值问题,波数是离散的,因此 (2.29) 范围中只有特定的 \(k\) 是相关的;这一点稍后讨论其含义。
关键假设(也是实际发生的)是:(2.30) 中这些指数增长的线性不稳定特征函数最终会受反应扩散方程中非线性项的约束,从而涌现出一个最终的稳态空间非齐次解。这一假设的一个关键元素是动力学"约束集"(confined set)的存在性(见 Volume I 第 3 章)。直观上我们预期:如果动力学存在约束集,加入扩散后同一集合也包含解。这确实成立,且能被严格证明;见 Smoller (1983)。因此对特定机制的分析部分涉及证明正象限内约束集的存在。对于演化到有限振幅稳态空间图样的一般非线性分析仍然欠缺,但对 \(d\) 接近分叉值 \(d_c\) 的奇异摄动分析已经完成,确实得到非均匀的空间非齐次解(见例如 Lara-Ochoa 与 Murray 1983、Zhu 与 Murray 1995)。奇异摄动分析可在分叉附近任一临界参数附近进行。现在已经用各种特定反应扩散机制数值计算了大量空间非齐次解;数值方法已经相当标准。下一章给出的结果展示了反应扩散能够生成的图样的一些丰富性。
为回顾起见,下面汇总两物种反应扩散机制 (2.11) 生成空间图样的条件。请记得所有导数都是在稳态 \((u_0, v_0)\) 处取值,由 (2.19)、(2.24) 与 (2.26) 得
导数 \(f_u\)、\(g_v\) 必须异号:对图 2.2 的反应动力学,\(f_u > 0\)、\(g_v < 0\),所以 (2.31) 的第一与第三式蕴含扩散系数比 \(d > 1\)。
对交叉项 \(f_v\) 与 \(g_u\) 有两种可能,因为唯一的限制是 \(f_v g_u < 0\)。所以必有 \(f_v < 0\)、\(g_u > 0\) 或反过来。这对应于定性上不同的反应。两种情形如图 2.6 所示。回想促进其中一种增长的反应物为激活剂,另一种为抑制剂。图 2.6(a) 中,\(u\) 是激活剂,也自我激活;而抑制剂 \(v\) 不仅抑制 \(u\),也自我抑制。要使图样形成发生,抑制剂必须扩散得比激活剂快。图 2.6(b) 中,\(v\) 是激活剂,但仍自我抑制,并扩散得更快。两种情形间还有另一差别。图样沿与正特征值相关的不稳定流形增长。在图 2.6(a) 中,这意味着随图样增长两个物种在同处同时高或同时低密度(如图 2.6(c));在图 2.6(b) 中,\(u\) 在 \(v\) 低处密度高,反之亦然(如图 2.6(d))。相平面(仅反应项)在稳态附近的定性特征如图 2.6(e) 与 (f) 所示,分别对应这两种情形。图样处于同相或反相对生物学应用有根本性的影响。
为对这两种情形获得直观感受,考虑两个不同的生态学捕食者-被捕食者场景。在图 2.6(e) 中,设 \(u\)、\(v\) 分别代表被捕食者与捕食者。在高捕食者密度处被捕食者数量减少,而在低密度处其数量增加。稳态附近,被捕食者彼此受益——数量增加会被暂时放大。捕食者在捕食者-被捕食者比例高时数量减少,否则增加。另一个例子来自寄生虫学:设 \(v\) 是通过可移动宿主扩散的寄生虫,\(u\) 是受影响严重的较定居宿主。这些情形下稳态附近的相互作用如图 2.6(a) 与图 2.6(e) 中所示的局部零等值线与定性增长。
捕食者-被捕食者情形中扩散驱动失稳的一个必要条件是捕食者分散得比被捕食者快。这种情况下图样以图 2.6(c) 的形式形成。假设有一块被捕食者密度升高的区域。在没有扩散时这会被抑制——因为捕食者会暂时增加然后回到稳态。然而,当捕食者扩散时,本地的捕食者增加(源于被捕食者的增加)会部分地分散出去,因此不足以将被捕食者种群压回平衡。捕食者扩散会降低邻近区域被捕食者的密度。最终可能形成被捕食者与捕食者种群同时高的"团块"与两者都低的区域交错的图样。在寄生虫类比中,定居被捕食者(宿主)的成簇与寄生虫高密度区域相重合。宿主在某些区域也可维持高水平——因为寄生虫会持续扩散到附近的"死区",而该死区中这种宿主数量稀少。图样化的尺度依赖于扩散系数之比 \(d\)。
现在考虑图 2.6(b)、(d) 与 (f) 所示的第二类相互作用。同样取捕食者-被捕食者情形,但 \(u\) 为捕食者,\(v\) 为被捕食者。捕食者在此情形下是"自催化的"——因为当密度接近稳态时,捕食者密度的增加会被暂时放大;这是一种并不少见的情形。例如,捕食者密度增加可以提升捕食或繁殖效率。此情形与第一类的另一差别是:现在是被捕食者扩散得更快。
假设再次有一块被捕食者密度高的区域。在没有扩散时,捕食者数量会增加,并最终使两个种群都回到稳态。然而有可能发生:捕食者增长并将被捕食者种群压到稳态值以下(被捕食者的暂时增加足以触发捕食者自催化增长的启动)。这会造成被捕食者从相邻区域的净流入,相邻区域捕食者密度会下降(自催化反向起作用),从而使被捕食者种群得以超出稳态值增长。可以建立起一种图样——低捕食者/高被捕食者区域与被捕食者少而捕食者多的区域交错出现。实质上,自催化的捕食者既受益于本地高密度,也受益于附近捕食者稀少从而不断供应额外被捕食者的区域。被捕食者因扩散的随机性继续流向高捕食压力的区域。
若 (2.31) 的条件满足,则存在一个尺度(\(\gamma\))依赖的图样范围,其波数由 (2.29) 定义,是线性不稳定的。最初(指数)增长的空间图样是那些波数在 \(k_1\) 与 \(k_2\) 之间、由 (2.29) 确定的特征函数 \(W_k(\mathbf{r})\)——即 (2.30) 中所示的那些。注意尺度参数 \(\gamma\) 在这些表达式中起关键作用,我们将在下一节进一步讨论。通常我们预期动力学与扩散系数是固定的。在胚胎发生情形中,天然的变动参数就是 \(\gamma\),它反映胚胎(或更具体说胚胎区域,如发育中的肢芽)的大小。
无限区域中的扩散驱动失稳:连续特征值谱。 在有限区域,可能的波数 \(k\) 与对应空间波长是离散的,部分依赖于边界条件。在发育生物学中,胚胎在空间图样化阶段的大小相对于待形成的图样往往足够大,以至于"边界"在孤立特定波长中不能起主要作用——例如在后面章节中讨论的毛发、鳞片与羽毛原基的生成。因此在实用上,图样形成区域可视为无限大。这里我们说明如何在无限区域中确定不稳定特征值的谱——这比有限区域更容易。
我们从线性化系统 (2.20) 出发,寻找如下形式的解:
其中 \(\mathbf{k}\) 是波矢,大小 \(k = |\mathbf{k}|\)。代入 (2.20) 同样得
因此 \(\lambda\) 与波数 \(k\) 的色散关系仍由 (2.23) 给出。\(\text{Re} \lambda(k^2) > 0\) 的特征值范围仍由 (2.29) 给出。无限区域与有限区域情形的关键差别在于:只要 (2.29) 中 \(0 < k_1^2 < k_2^2\),就总存在一个空间图样(因为我们不受 (2.21) 特征值问题定义的离散 \(k^2\) 类限制)。因此在分叉时——即 \(k_c^2\)(由 (2.28) 给出)线性不稳定时——机制会演化出具有临界波长 \(\omega_c = 2\pi / k_c\) 的空间图样。一般而言,图 2.5(b) 中指数增长最大的波长将成为至少在一维中最终涌现的图样:但这并非总是如此,依赖于不稳定模式的数目与初始条件。下一章在生物学应用中我们会看到:有限区域与有效无限区域之间的区别具有重要的生物学含义——有限区域对允许的图样施加相当多的限制。
2.4 反应扩散机制中图样启动的详细分析(Detailed Analysis of Pattern Initiation in a Reaction Diffusion Mechanism)
这里以一个具体的两物种反应扩散系统为例进行详细分析。本节为后续在真实生物图样形成问题中的应用奠定基础。我们计算特征函数,得到启动空间图样所必需的具体参数条件,并确定指数增长的初始空间扰动的波数与波长。
我们研究最简单的反应扩散机制 (2.7),先在一维情形;即
动力学零等值线 \(f = 0\) 与 \(g = 0\) 在图 2.2(a) 中已示。均匀正解 \((u_0, v_0)\) 为
且在稳态处,
由于 \(f_u\) 与 \(g_v\) 必须异号,需要 \(b > a\)。由这些表达式,(2.31) 的条件要求
这些不等式在 \((a, b, d)\) 参数空间中定义一个区域,称为图样形成空间(Turing 空间),在此区域内机制对给定波数 \(k\) 的某些空间扰动不稳定——我们下面将确定这些 \(k\)。
考虑相关特征值问题 (2.21),并取区域为 \(x \in (0, p)\),\(p > 0\)。则
其解为
其中 \(A_n\) 为任意常数。特征值为离散波数 \(k = n \pi / p\)。每当 (2.34) 满足且存在波数 \(k = n \pi / p\) 落在 (2.29) 定义的范围中时,相应特征函数 \(W_n\) 线性不稳定。因此 (2.37) 中波长 \(\omega = 2\pi / k = 2p / n\) 的那些特征函数是按 \(\exp\{\lambda([n \pi / p]^2) t\}\) 随时间初始增长的。由 (2.29) 与 (2.34),不稳定波数的带为
用波长 \(\omega = 2\pi / k\) 表示,不稳定模式 \(W_n\) 的波长介于 \(\omega_1\) 与 \(\omega_2\) 之间,其中
注意 (2.38) 中尺度(以 \(\gamma\) 度量)的重要性。最小波数为 \(\pi / p\);即 \(n = 1\)。对固定的参数 \(a\)、\(b\)、\(d\),若 \(\gamma\) 充分小,(2.38) 表示范围内没有允许的 \(k\),于是 (2.37) 中没有任何模式 \(W_n\) 可被驱动失稳。这意味着 (2.30) 中 \(\mathbf{w}\) 的所有模式指数衰减,稳态稳定。尺度的这一重要作用我们将在下面更详细讨论。
由 (2.30),涌现出的空间非齐次解是不稳定模式之和,即
其中 \(\lambda\) 由 (2.23) 取以 (2.34) 导数的正解,\(n_1\) 是大于等于 \(p k_1 / \pi\) 的最小整数,\(n_2\) 是小于等于 \(p k_2 / \pi\) 的最大整数,\(C_n\) 由 \(\mathbf{w}\) 初始条件的 Fourier 级数分析决定。生物学语境中的初始条件带有一定的随机性,因而 Fourier 谱中必然包含全部 Fourier 模式——即 \(C_n\) 非零。我们因此可在此阶段假设 \(\gamma\) 充分大以确保允许的波数存在于 \(k\) 的不稳定范围内。在讨论可能涌现的图样之前,先给出对应的二维结果。
考虑二维区域 \(0 < x < p\)、\(0 < y < q\),矩形边界记为 \(\partial B\)。取代 (2.36) 的空间特征值问题现在是
其特征函数为
其中 \(n\)、\(m\) 为整数。线性不稳定的二维模式 \(W_k(x, y)\) 是那些波数 \(k\)(由最后那个方程定义)落在 (2.38) 以 \(a\)、\(b\)、\(d\) 定义的不稳定带内的那些。我们同样假设 \(\gamma\) 充分大使得不稳定波数范围至少包含一个可能的模式。现在不稳定空间图样解由 (2.30) 与 (2.42) 给出
求和遍及所有满足不等式的 \((n, m)\) 对,\(L\)、\(M\) 与前述 (2.38) 定义相同,\(\lambda(k^2)\) 是 (2.23) 取以 (2.34) 给出的 \(f\)、\(g\) 导数的正解。当 \(t\) 增大时,空间图样将演化,其最初由 (2.43) 中的模式构成。
现在考虑由 (2.40) 与 (2.43) 的不稳定解可能产生的空间图样的类型。假设区域大小(以 \(\gamma\) 度量)使得 (2.38) 的不稳定波数范围只允许波数 \(n = 1\):\(\lambda\) 关于波长 \(\omega = 2p / n\) 的色散关系如图 2.7(a) 所示。从 (2.37) 知此时唯一的不稳定模式是 \(\cos(\pi x / p)\),增长的不稳定由 (2.40) 给出为
其中 \(\lambda\) 是 (2.23) 二次方程的正解,导数 \(f_u, f_v, g_u, g_v\) 取 (2.34),\(k^2 = (\pi / p)^2\)。这里所有其他模式都随时间指数衰减。我们只能由初始条件确定 \(C_1\)。为对正在发生的事情获得直观理解,我们简单取 \(C_1\) 为 \((\varepsilon, \varepsilon)\)(\(\varepsilon\) 为某小正数)并考虑形态发生素 \(u\);即从上一方程与 \(\mathbf{w}\) 的定义 (2.14) 得
这个不稳定模式——它将随 \(t\) 增大而涌现为主导解——如图 2.7(b) 所示。换言之,这就是图 2.7(a) 的色散关系所预测的图样。
显然如果指数增长的解对所有时间都成立,则会蕴含 \(u \to \infty\) 当 \(t \to \infty\)。对机制 (2.32),动力学在正象限内有一个约束集,它界定了解的范围。所以上一方程的解必须有界且位于正象限内。我们假设此增长解最终会稳定为一个与图 2.7(b) 中所示单个余弦模式类似的空间图样。如前所述,对分叉附近任一参数(例如接近域大小 \(\gamma\) 的临界值,使某一单个波数刚好失稳,或接近临界扩散系数比 \(d_c\))的奇异摄动分析印证了这一点;完整非线性方程的大量数值模拟也印证了这一点。图 2.7(c) 是展示反应扩散机制空间图样结果的一种有用方式——阴影区代表浓度高于稳态值 \(u_0\),空白区代表浓度低于稳态值。正如我们将要看到的,这种简单的结果展示方式在把化学预图样理论应用于发育生物学的图样化问题中非常有用——按其设想,细胞依据所处形态发生素浓度是否高于(或低于)某一阈值而分化。
现在假设区域大小加倍。由于 \(\gamma\) 的定义代表尺度,在一维情形下 \(\gamma^{1/2}\) 与区域大小(这里即长度)成正比,因此加倍区域长度等价于把原 \(\gamma\) 乘以 4。这意味着色散关系与不稳定范围只是沿 \(k^2\) 轴(或 \(\omega^2\) 轴)平移。设原 \(\gamma = \gamma_1\)。(2.38) 定义的不稳定模式为波长 \(\omega\)(由 (2.39) 定义的那些)满足
让这对应于图 2.7(a) 所示情形——产生图 2.7(c) 的图样。现在令区域大小加倍。我们考虑与图 2.7 同样的区域,但把 \(\gamma\) 增大到 \(4 \gamma_1\)。这等价于保持 \(\gamma_1\) 不变但区域是图 2.7 中的 4 倍。我们选择前一种方式来表示尺度变化。等价的色散关系如图 2.8(a) 所示——它只是把图 2.7(a) 的原色散关系沿 \(\omega^2\) 轴平移,使被激发的不稳定模式现在有 \(\omega = p\),即 \(n = 2\)。等价的空间图样如图 2.8(b) 所示。在后续应用章节中我们将会看到,在呈现空间图样解时,仅通过改变 \(\gamma\) 体现尺度,是一个特别方便的方式。
由此例可看到图样形成过程在区域大小方面的运作方式。对给定的 \(\gamma = \gamma_1\),分析挑出一个基本波长——本例中为 \(n = 1\) 的波长。随着区域增长,它最终能容纳 \(n = 2\) 的图样,并依次容纳更高模式(图 2.8(c));同样地,如果区域足够小,显然存在一个 \(\gamma = \gamma_c\),使色散关系(在图 2.8(a) 中右移)甚至不能容纳 \(n = 1\) 的波长。这种情况下没有任何模式失稳,因此不能生成任何空间图样。空间图样存在的临界区域大小这一概念在发育生物学以及空间依赖的生态学模型中都很重要。
注意图 2.8(b) 中相同参数与零通量边界条件下的两个可能解。具体哪个解出现取决于初始条件的偏向性。它们的存在从发育生物学的位置信息语境下提出了一些概念性困难。若细胞在形态发生素浓度大于某阈值时分化,那么两种可能解对应的分化细胞图样显然不同。然而发育是顺序性的过程,自带其历史——前一阶段一般提示下一阶段。在反应扩散模型语境下,这暗示初始条件对某一种图样有偏向。
现在考虑二维问题,其色散关系使得不稳定模式由 (2.43) 给出。这里的情形不那么直接,因为对给定的 \(\gamma\)(代表尺度),实际不稳定的模式依赖于区域几何——以长度 \(p\) 与宽度 \(q\) 度量。参看 (2.43),首先注意:若宽度足够小(即 \(q\) 足够小),即使 \(m = 1\) 的第一模式也落在不稳定范围外,问题等价于一维情形。当宽度增大(即 \(q\) 增大)时,\(n \neq 0\)、\(m \neq 0\) 的真正二维模式变得不稳定,因为
落在不稳定波数范围内。图 2.9 示意了 (2.43) 描述的、在不同 \((n, m)\) 下随时间增长的空间图样。
规则平面镶嵌图样。 上图所示的线性图样源于 (2.41) 最简单的二维特征函数。较为复杂的区域需要求解
除非区域几何简单,否则分析很快变得相当复杂。即使是圆域,特征值也必须数值确定。然而对于某些对称区域,存在可平铺整个平面的初等解——正方形、六边形、菱形,以及细分后的三角形;它们由 Christopherson (1940) 找到。也就是说,我们可以用例如规则六边形瓦片覆盖整个平面(规则多边形的基本对称群为六边形、正方形与菱形,当然还有三角形——它们是前两者的子单元)。六边形图样,正如我们将要看到的,在许多真实发育情境中很常见——鸟皮上羽毛的分布就是一个例子(只需看看拔了毛的鸡皮即可)。另请参看下面的图 2.11,展示了多种实验获得的图样。因此我们需要的 \(\psi\) 解,其单元胞在边界上满足零通量条件,且是能覆盖平面的规则镶嵌之一。也就是说,我们需要的解具有胞周期性;这里的"胞"自然指镶嵌的单元。
方程 (2.46) 对六边形的解为
由 (2.46)(线性方程),\(\psi\) 在乘以任意常数意义下独立:这里选取的形式使 \(\psi = 1\) 在原点。当 \(k = n \pi\)(\(n = \pm 1, \pm 2, \ldots\))时,该解在六边形对称边界上满足零通量条件。图 2.10(a) 展示了该解可生成的图样类型。极坐标形式显示对六边形旋转的不变性——即按 \(\pi / 3\) 旋转不变——这是必须的。也就是说,
其中 \(H\) 是六边形旋转算子。
正方形的解为
其中 \(k = \pm 1, \pm 2, \ldots\) 且 \(\psi(0, 0) = 1\)。此解具有正方形旋转不变性,因为
其中 \(S\) 是正方形旋转算子。典型图样如图 2.10(b) 所示。
菱形的解为
其中 \(\phi\) 是菱形角,且 \(k = \pm 1, \pm 2, \ldots\)。此解在菱形旋转下不变;即
其中 \(R\) 是菱形旋转算子。示意性图样如图 2.10(c) 所示。
另一种胞周期解是正方形的一维版本:只在 \(x\) 方向有变化。其解形如
代表条纹状图样,如图 2.10(d) 所示。当然这些只是一维解 (2.37) 的简单形式。
当完整非线性方程以稳态附近的小随机扰动为初始条件数值求解时,线性理论对最终稳态的预测在一维情形下效果不错,尤其当不稳定模式具有较大波长(即较小波数)时。波数较大时,预测可靠性下降——二维结构更甚。由于我们研究的方程是线性的且在乘以常数时不变,我们可以有等价的解,即关于直线 \(u = u_0\) 的镜像;参看图 2.8(b)。因此演化出来的图样依赖于初始条件,最终图样倾向于最接近初始条件的那一个。从某种意义上说,空间图样关于初始条件有一个吸引域。再一次,在分叉附近,奇异摄动分析表明非线性图样与线性预测密切相关。但一般而言,远离分叉边界时线性预测要不可靠得多;参看下一章中展示的计算图样。除最简单的图样外,我们对二维与三维情形应仅把线性理论作为"反应扩散机制可以生成什么图样"的指引。然而线性理论确实决定了图样生成的参数范围。
图 2.10 展示了一些反应扩散方程在线性理论基础上能形成的规则图样。从数学上(当然也包括实验上)一个关键问题是:给定初始条件,哪一种会被形成?若形成了一种图样,改变哪些参数能使其变为另一种?为确定各种可能图样——六边形、菱形、正方形或条纹——中的哪个会稳定,我们必须超越线性理论进行弱非线性分析;即参数接近从均匀到非均匀的分叉边界。当做这样的非线性分析时,我们能确定使这些稳态空间非齐次解稳定的参数条件。这已由 Ermentrout (1991) 与 Nagorcka 与 Mooney (1992) 用多尺度奇异摄动分析对反应扩散方程完成。其他图样形成机制——细胞趋化与力学机制——由 Zhu 与 Murray (1995) 研究。后者比较了趋化系统与反应扩散系统的图样形成能力。Zhu 与 Murray (1995) 特别感兴趣的是确定能产生稳定条纹、点状、正方形与六边形的参数空间,以及它们的波长等空间特征。
对点状图样,他们也能确定哪种镶嵌点状排列会稳定。他们比较了不同模型的稳健性与灵敏度,并用大量方程的数值模拟加以确认。分析方法已经成熟,但细节相当复杂。Zhu 与 Murray (1995) 通过对方程的数值研究表明:从条纹到点状再到六边形图样的转变与反过来的路径——六边形失稳并最终变为条纹。六边形实际上被拉长而呈菱形特征,点状排成行并最终融合;这符合直觉。分叉附近的解析过程称为弱非线性稳定性分析,Wollkind et al. (1994) 给出了广泛综述。一般而言,相互作用动力学的形式在所获得的图样中起主要作用。三次相互作用倾向于支持条纹,而二次相互作用倾向于产生点状。当使用不同(非零通量)边界条件时,所获得的图样可能很不一样且更难预测。Barrio et al. (1999) 研究了这些条件与非线性在所获图样中的作用。凭藉大量数值模拟,他们推测这些反应扩散机制可能在鱼身上观察到的某些复杂图样中起作用。
迄今为止我们讨论的图样在某种意义上都是规则的——条纹、点状、六边形图样等等。反应扩散系统能生成大量不规则图样——我们将在下一章讨论一些实际例子时看到。Meinhardt (2000;亦见该文其他参考文献) 讨论了复杂图样,特别是反应扩散机制在基因激活图样上的应用——这是本书未涉及的主题。他还综述了其他本书未涉及的重要应用,例如植物形态学中的分枝结构。
在对反应扩散系统图样形成能力的分析中,这里假设反应物(即形态发生素)必须具有不同的扩散系数。在许多发育情境下,相同形态发生素在不同方向的扩散可能具有不同的值——即扩散是各向异性的。虽然我们这里不讨论,但如可预期,这对均匀态 Turing 失稳所形成的图样有显著影响(见习题 10)。
早在 1970 年代就已从许多数值研究中知道反应扩散系统能产生稳态有限振幅的空间图样。然而直到近十年,这种稳态图样(有时叫 Turing 结构或 Turing 图样)才在实验中发现。实验突破始于 1989 年;见 Ouyang et al. (1990, 1993)、Castets et al. (1990)、Ouyang 与 Swinney (1991)、Gunaratne et al. (1994)、De Kepper et al. (1994) 及其中的其他参考文献。最后两篇是获得整体图像的良好综述。后者也描述了 Turing 结构与行波相互作用时获得的复杂结构——它们可以非常复杂,伴随诸如时空间歇性与点分裂形成更复杂图样的现象。Ouyang 与 Swinney (1991) 在实验上演示了从均匀态到六边形最终再到条纹图样的转变;该转变与 Zhu 与 Murray (1995) 对反应扩散与细胞趋化图样形成机制的发现类似。自这些早期实验研究以来,Turing 图样已在几种相当不同的反应系统中发现;具体的化学与实验安排在原论文中详细给出。图 2.11 展示了 Gunaratne et al. (1994) 用氯-碘-丙二酸(chlorite-iodide-malonic acid)反应扩散系统在实验上获得的化学 Turing 图样。注意区域小、波长精确定义,介于 0.11 mm 到 0.18 mm;这确实在我们对许多形态发生情境所预期的范围内,并清楚地展示了反应扩散机制精细尺度图样化的潜力。Zhu 与 Murray (1994) 对其他图样生成器做了理论研究。在发育中鸡肢的情形下,与软骨形成相关图样化时期肢芽的宽度约为 2 毫米(见第 6 章关于肢芽图样化场景的讨论)。Wollkind 与 Stephenson (2000a, b) 对图样间的各种转变——包括图 2.11 中所示的黑眼图样——进行了详尽而全面的讨论。他们专门研究了在实验中所用的氯-碘-丙二酸反应系统,并将其结果与实验进行了重要比较。他们还把化学系统中对称破缺结构间的转变与相当不同科学语境中类似转变联系起来。
把反应扩散图样生成应用于具体的发育生物学问题,通常是在"预图样"理论的语境下:细胞依据所处形态发生素浓度水平进行分化。若空间图样相当分明(如上所述或具有较大的梯度),那么细胞执行指定功能所需的敏感性调谐较小;而图样变化或浓度梯度较小时则需要更精细的调谐。因此尝试获得一种对空间异质性的生物可解读的定量度量以比较不同机制,可能是有用的。下一节将讨论另一种生物上相关的方法。
Berding (1987) 引入了一种"异质性"函数,针对零通量边界条件下反应扩散系统生成的空间图样。假设通式 (2.10) 的机制在一维空间变量中扩散失稳,解在 \(t \to \infty\) 时演化为空间非齐次稳态解 \(U(x)\) 与 \(V(x)\)。由 (2.6) 中 \(\gamma\) 的定义(与域长的平方成正比),我们用 \(\gamma\) 度量区域大小,取 \(x \in (0, 1)\)。则 \((U, V)\) 满足无量纲方程
非负异质性函数定义为
它只依赖于系统的参数与区域尺度 \(\gamma\)。\(H\) 是一个"能量函数"。若用零通量边界条件 (2.51) 进行分部积分,
再利用 (2.51) 中关于 \(U''\)、\(V''\) 的表达式,可得
若不存在空间图样化,\(U\)、\(V\) 简单地是 \(f(U, V) = g(U, V) = 0\) 的均匀稳态解,于是 \(H = 0\),当然这也可从定义 (2.52) 得出。
由 (2.53) 我们看到尺度参数与扩散系数比如何出现在异质性的定义中。例如,假设域恰好支撑 \(\gamma = \gamma_1\) 下的单个波——在有量纲下,区域长度 \(L = L_1\)。如果域大小加倍到 \(2 L_1\),则可容纳两个波——直观上从 (2.52) 知 \(H\) 必须增大,因为有更多异质性。由于 \(\gamma \propto L^2\),由 (2.53) 知 \(H\) 简单地翻 4 倍。从胚胎学观点看,这意味着胚胎成长时我们预期会有越来越多的结构。图 2.18 后面将展示一个在增长区域中这种结构增加的例子。Berding (1987) 讨论了具体的应用并就异质性潜力比较了具体的反应扩散机制。
2.5 模式形成模型中的色散关系、Turing 空间、尺度与几何效应(Dispersion Relation, Turing Space, Scale and Geometry Effects in Pattern Formation Models)
我们首先就色散关系的一些一般性质作说明,然后用在上一节分析过的具体情形进一步挖掘。任何形态发生模型的空间图样形成主要是非线性现象。然而如我们所述,对一维中图样一个良好的指示可以由简单线性分析得到。要使空间图样形成,我们看到两个条件必须同时满足。第一,空间均匀态对小扰动稳定——即 (2.22) 中的所有 \(\lambda(k^2 = 0)\) 都满足 \(\text{Re} \lambda(k^2 = 0) < 0\);第二,只有特定空间范围的图样——即波长 \(k\) 在一定范围内——可以开始增长,\(\text{Re} \lambda(k^2 \neq 0) > 0\)。这些条件封装在 (2.5)(b) 与图 2.7(a) 那样的 \((\lambda, k^2)\) 或 \((\lambda, \omega^2)\) 形式的色散关系中。后者——例如——也表明:若空间扰动图样 \(k^2\) 很大(即波长非常小的扰动),则稳态仍线性稳定。因此色散关系立即给出各种大小图样的初始增长或衰减率。色散关系由图样形成机制的一般演化方程得到。Oster 与 Murray (1989) 给出了图样形成模型一般且非技术的、面向生物学的讨论。
由于线性特征函数方程 (2.36) 的解只是正弦与余弦函数,因此各种空间图样的"大小"由三角函数的波长度量;例如 \(\cos(n \pi x / p)\) 的波长 \(\omega = 2p / n\)。因此寻找增长的空间图样就变成看多少个正弦或余弦波能"装进"给定大小的区域中。二维情形类似,但波如何组合有更大的灵活性。
色散关系的一个非常重要用途在于:它立即表明图样能否增长,若能增长,其大小是多少。图 2.5(b) 与图 2.7(a) 中的曲线是生成空间图样的典型"纯形式"色散关系。后面我们将展示其他形式也是可能的,并意味着不同的图样形成场景。然而这些较少见,且对哪些图样会从中演化仍知之甚少。第 6 章中详细讨论的力学化学模型事实上能生成令人惊讶的丰富色散关系谱(见 Murray 与 Oster 1984),其中大多数无法由两物种或三物种反应扩散模型生成。
原型色散关系有上述两个基本特征:(i) 空间无特征态(\(k = 0\),\(\omega = \infty\))稳定——即非常大波长波的增长率为负;(ii) 存在一个小的波长带或"窗口"可以增长(即在有限域情形下,对有限整数 \(n\) 的不稳定"模式" \(\cos(n \pi x / L)\) 的有限带)。这些增长模式中有一个增长最快——最接近色散曲线峰值的那一个——其波数 \(k_m\) 由下式给出
严格来说 \(k_m\) 在有限域情形下可能不是允许的模式。这种情况下它是与分析所得 \(k_m\) 最接近的允许模式。
因此色散曲线表明:尽管空间均匀态稳定,但系统会放大特定空间范围的图样,只要它们被随机波动(在生物系统中总是存在)或发育早期图样的提示所激发。一般而言,模型参数中有一个被"调谐",直到色散曲线获得所示定性形状。例如,在图 2.5(b) 中,若扩散比 \(d\) 小于临界 \(d_c\),则对所有 \(k^2\) 都有 \(\text{Re} \lambda < 0\)。随着 \(d\) 增大,曲线上升,直到 \(d = d_c\) 之后在某波数 \(k_c\)(即波长 \(\omega_c = 2\pi / k_c\))处把头部推过轴线——此时该波长的余弦波可以开始增长,假设它是允许的特征函数。该临界波数由 (2.28) 给出,由 (2.27) 取 \(d = d_c\),我们得到临界波长
对示例 (2.32),有 4 个无量纲参数:动力学参数 \(a\)、\(b\),扩散系数比 \(d\),以及尺度参数 \(\gamma\)。我们集中讨论色散关系随 \(d\) 的变化,并展示了均匀稳态失稳时分叉值 \(d_c\) 如何存在,图样"大小"由 \(k_c\) 或 \(\omega_c\) 由上式确定。了解涉及所有参数的图样形成参数空间以及如何通过改变所选参数(或同时改变多个参数)进入该图样形成域,是非常有用的。显然参数越多,相应参数空间或 Turing 空间就越复杂。我们现在用 Volume I 第 3 章中描述的参数方法(用于确定振荡解可能存在的空间)的扩展,对模型 (2.32) 解析地确定该参数空间。这一方法由 Murray (1982) 发展并应用于若干反应扩散模型;它是一个适用于其他图样形成机制的一般过程。也可以用 Zhu 与 Murray (1995) 发展的数值方法。
若区域足够大,机制 (2.32) 生成空间图样的参数 \(a\)、\(b\)、\(d\) 条件由 (2.35) 给出,\(\gamma\) 通过 (2.38) 决定的可能不稳定模式进入考量。尽管 (2.35) 的不等式可能是任何反应扩散机制中我们能写出的最简单现实的一组,处理起来代数上仍然相当繁琐。对于比最简单动力学更复杂的情况,不可能解析地进行类似分析。因此我们从 Volume I §7.4 中 (7.24) 给出的稳态表示出发,以 \(u_0\) 为非负参数变量;即 \(v_0\) 与 \(b\) 用 \(a\)、\(u_0\) 由 (2.55) 表示
利用这些表达式,(2.35) 中的不等式涉及
现在把 (2.31) 给出的扩散驱动失稳条件表示为关于参数 \(u_0\) 的不等式;它们在参数空间中定义边界曲线。由第一个,
作为边界,其中(由于我们对边界曲线感兴趣)\(b = u_0 - a\) 来自稳态定义 (2.55),并且用 \(a\) 来自仅涉及 \(u_0\) 与 \(a\) 的不等式的表达式来替换。这些在 \((a, b)\) 空间中以参数形式定义一个域——让 \(u_0\) 取所有正的值即可;若不等式取等号,(2.57) 即以参数形式定义边界曲线。我们现在对 (2.31) 的每个条件做同样处理。
(2.31) 的第二个条件,用 (2.56),要求
自动满足。第三个条件要求
作为边界曲线。
(2.31) 的第四个条件更复杂一些。这里
左侧因式分解后蕴含
因此此不等式产生两条边界曲线,即
由 (2.57)–(2.60) 定义的曲线以及它们所围的区域共同定义了参数空间(或 Turing 空间)(见 Murray 1982),其中稳态可被扩散驱动失稳从而产生空间图样。正如 §2.4 中所指出的,(2.35) 的第一与第三条件要求 \(f_u\) 与 \(g_v\) 异号,这又要求 \(b > a\) 从而 \(d > 1\)。
现在剩下的是简单的绘图练习——让 \(u_0\) 取一系列正值并对给定的 \(d\) 计算对应的 \(a\)、\(b\)——即可得到由 (2.57)–(2.60) 定义的曲线。一般而言,结合 (2.57)–(2.60) 的不等式,五条曲线参与定义边界。这里——如常见的情况——其中若干条是冗余的,被另一条所覆盖。例如,对 (2.60) 中第一个,
因为我们考虑 \(u_0 > 0\),所以如果满足 (2.60) 中的第一个条件,(2.59) 自动满足。又因为 \(d > 1\),
所以 (2.57) 定义的曲线位于 (2.59) 定义的曲线之下;前者是下方的极限边界曲线——所以若使用 (2.60) 的第一式即可定义一个合适的域。此外,由于
不可能存在同时满足 (2.59) 与 (2.60) 中第二条曲线的域。
因此,对本机制我们只需要两条参数曲线——即 (2.57) 与 (2.60) 的第一式所定义者——Turing 空间由下式确定
我们知道当 \(d = 1\) 时没有 Turing 空间——即不存在可以产生空间图样的域。(2.61) 在 \(d = 1\) 时定义的曲线相互矛盾,因此不存在 Turing 空间。现在让 \(d\) 取大于 1 的值。对一个临界 \(d\)(记为 \(d_c\)),Turing 空间在 \(d > d_c\) 时开始出现。具体地,\(d = d_c = 3 + 2 \sqrt{2}\),由 (2.61) 通过确定使两条曲线都在 \(b = 1\) 处给出 \(a = 0\) 的 \(d\) 计算——此时两个不等式不再矛盾。该空间实际上由 \((a, b, d)\) 空间中的两个曲面定义。图 2.12 展示了 (2.32) 的 Turing 空间——即参数空间中反应扩散机制 (2.32) 可以生成空间图样的区域。例如,若 \(d = 25\),则位于由标为 C(\(d = 1\))的曲线与 \(d = 25\) 的曲线所围域内的任何 \(a\)、\(b\) 值都将导致扩散驱动失稳。若区域(以 \(\gamma\) 度量)足够大以容纳 (2.38) 与 (2.43) 定义的允许 \(k^2\),则空间图样将演化。
即使给定 \(d > 1\) 时 \(a\)、\(b\) 位于 Turing 空间内,这也并不保证机制会产生空间图样,因为尺度与几何起重要作用。取决于 \(\gamma\) 的大小以及机制实际作用的空间区域,不稳定特征函数(或模式)可能不是允许的解。这里色散关系的具体形式再次重要。具体而言,考虑由 (2.36) 定义的一维有限域问题。特征值——即波数——\(k = n \pi / p\)(\(n = \pm 1, \pm 2, \ldots\))是离散的。因此,参看图 2.5(b),除非色散关系在其不稳定模式范围中至少包含这些离散值中的一个,否则不能形成任何结构。因此我们必须在图 2.12 的 Turing 空间上叠加另一条表示尺度参数 \(\gamma\) 的轴。若把 \(\gamma\) 包含在 Turing 空间的参数中,则该空间不一定简单连通——因为若色散关系随 \(\gamma\) 变化在其不稳定模式中不包含任何允许的特征函数,则没有图样演化。我们来考虑这一方面并更详细地审视色散关系。
Turing 空间只涉及无量纲参数,它们是有量纲模型的合适参数组合。图 2.12 中的参数 \(a\)、\(b\)、\(d\) 由 (2.6) 给出
假设例如 \(d = 100\) 且 \(a\)、\(b\) 取与图 2.12 中 P 点相关的值——即机制不在图样形成模式中。进入图样形成域并无唯一方式;我们可以减小 \(a\) 或 \(b\),从而分别到达 Q 或 R。在有量纲意义上,例如通过适当改变 \(k_1\)、\(k_2\) 或 \(k_3\)——或同时改变——来减小 \(a\)。若改变的不是 \(k_1\),则也会影响 \(b\),因此我们必须同步关注 \(b\)。若仅改变 \(k_2\),则 Turing 空间中的路径在定性上类似从 P 到 S 的路径。若 \(d\) 也可变——这意味着 \(D_A\) 或 \(D_B\) 可变——我们可以通过简单增大 \(d\) 把 P 包络入图样形成区域。
从生物学观点解释这些结果,我们因此看到:产生图样是多种效应共同协调的结果,而非仅仅一种效应——因为我们可以通过改变若干参数中任一个来进入图样形成域。显然我们可以由多条路径到达空间中的特定点。等效效应(通过参数变化)产生相同图样的概念在解释和设计相关模型实验时是重要的。这在生物学中尚未被广泛认识。我们将在后续章节中讨论无量纲组合的实际用途在重要生物学应用中的体现。
简要回顾一下:对通式反应扩散系统 (2.10),色散关系由 (2.23) 实部最大的根 \(\lambda(k^2)\) 给出。空间不稳定模式存在的关键是函数
对某个 \(k^2 \neq 0\) 是否为负;见图 2.5(a)。记住 \(f\)、\(g\) 的导数都在稳态 \((u_0, v_0)\) 处取值,且 \(f(u_0, v_0) = g(u_0, v_0) = 0\),所以 \(h(k^2)\) 是关于 \(k^2\) 的二次式,其系数仅为动力学参数、扩散系数比 \(d\) 与尺度参数 \(\gamma\) 的函数。\(h\) 在 \(k = k_m\) 处的最小值 \(h_{\min}\) 对应 \(\lambda\) 中实部最大者——也就是具有最大增长因子 \(\exp[\lambda(k_m^2) t]\) 的模式。由 (2.25),或直接由上式,\(h_{\min}\) 给出为
空间稳定与不稳定模式之间的分叉是 \(h_{\min} = 0\)。成立时,存在一个临界波数 \(k_c\),由 (2.28) 或再次直接由 (2.62) 推导,当参数满足
随着参数在 Turing 空间中移动,我们可以通过让某一参数通过其分叉值(其他参数保持固定)来达到所需的等式——即 (2.64) 的第一式。在上一节以及图 2.5(b) 中,我们选择 \(d\) 作为变化的参数,对给定的 \(a\)、\(b\) 计算分叉值 \(d_c\)。在这种情况下,恰好在分叉处(即 \(h_{\min}(k_c^2) = 0\))时,一个具有波数 \(k_c\) 的单一空间图样被驱动失稳,或被激发,对 \(d = d_c + \varepsilon\)(\(0 < \varepsilon \ll 1\))。由 (2.64) 知该临界波数正比于 \(\sqrt{\gamma}\),因此我们可以通过改变 \(\gamma\) 来改变所启动的空间图样。这被称为模式选择,在应用中是重要的。
在有限域情形下,我们可以通过选择不稳定波数带宽窄且以期望模式为中心来孤立某个特定模式被激发或驱动失稳。设动力学中参数固定,\(d = d_c + \varepsilon\)(\(0 < \varepsilon \ll 1\))。由 (2.64) 我们得到对应特定 \(k\) 的合适 \(\gamma\) 近似为
其中在分叉处(有时称为边际动力学状态)的动力学参数满足 (2.64) 的第一式。因此通过变化 \(\gamma\),我们可以孤立我们希望激发的任何模式。图 2.13(a) 给出了典型情形。Arcuri 与 Murray (1986) 对更复杂的 Thomas (1975) 机制在这种情况下进行了广泛的 Turing 空间分析。注意图 2.13(a) 中随着 \(\gamma\) 增大 \(h_{\min}\) 变得更负,正如 (2.63) 所示。
现在保持 \(\gamma\) 与动力学参数固定,让 \(d\) 从其分叉值 \(d_c\) 增大。由 (2.63),对 \(d\) 充分大 \(h_{\min} \sim - (d / 4) (\gamma f_u)^2\),于是 \(\lambda \to \infty\) 随 \(d \to \infty\)。不稳定模式波数带由 (2.62) 中 \(h(k^2)\) 的零点 \(k_1\)、\(k_2\) 所界定。这由 (2.66) 给出(或直接由 (2.62)),由此得
因此对固定尺度,不稳定模式波数有上界,从而空间图样波长有下界。图 2.13(b) 示意了对 Thomas (1975) 系统 (2.8) 的典型情形。
所有动力学参数都固定时,每对参数 \((d, \gamma)\) 定义 (2.62) 中 \(h(k^2)\) 的唯一抛物线,它进而指定一组不稳定模式。因此我们可以把 \((d, \gamma)\) 平面视为被分成若干区域,其中特定模式或一组模式是扩散性不稳定的。当存在多个不稳定模式时,由 (2.5)(b) 形式的色散关系形式,显然存在一个具有最大增长率的模式,因为对某 \(k_m^2\) 而言 \(\text{Re} \lambda\) 有最大。由 (2.23),正特征值 \(\lambda_+(k^2)\) 给出为
对波数 \(k_m\) 取最大值
正如我们已经指出的,预测是最快增长的 \(k_m\) 模式将占主导,因此将演化为稳态非线性图样。这对低阶模式是合理的预测。原因是:高阶模式下,非线性引起的相互作用比仅最简单模式线性不稳定时更为复杂。因此用 (2.68) 我们可以映射出 \((d, \gamma)\) 空间中特定模式从而特定图样将演化的区域;见 Arcuri 与 Murray (1986)。图 2.14(a) 与 (b) 分别显示了由线性理论与完整非线性 Thomas 系统的 Thomas (1975) 系统在一维中的映射结果;图 2.14(c) 显示了对应于图 2.14(b) 的空间形态发生素图样。
此类参数空间的一个重要用途是衡量机制对随机参数变化的稳健性。例如,用图 2.14(b),假设生物条件导致 \((d, \gamma)\) 参数对给出某点 P(位于演化为 4-模式的区域)。任何模型的一个关键性质是它对真实世界中不可避免的随机扰动的敏感程度。由图 2.14(b) 我们可以看到若在发育序列中需要 4-模式图样时有多少余地。当然这个 \((d, \gamma)\) 空间只是需要考虑的相关空间之一,因为任何机制还涉及其他参数。因此在评估稳健性或模型灵敏度时,我们还必须考虑涉及所有动力学参数的 Turing 空间的大小与形状。不同的反应扩散系统,\((d, \gamma)\) 空间在定性上可能不会太不一样。但肯定不同的是 Turing 空间的大小与形状——正是这个空间提供了比较模型相关稳健性的另一个有用准则。Murray (1982) 研究了这个具体问题并就相关机制进行了比较,得出关于较稳健机制的若干结论:Thomas (1975) 与 Schnakenberg (1979) 系统(分别由 (2.7) 与 (2.8) 给出)都具有相对较大的 Turing 空间,而 Gierer 与 Meinhardt (1972) 的激活剂-抑制剂模型(由 (2.9) 给出)的 Turing 空间相当小,意味着图样对小的参数变化具有相当大的灵敏度。在下一章关于具体生物学图样形成问题的讨论中,我们会涉及由色散关系形式与所采用的无量纲化所暗示的模型相关性的其他重要方面。
指定特定图样区域的参数空间都是在以稳态附近随机扰动为初始条件时得到的。即使在低阶模式,极性也能被偏向的初始条件明确影响。例如,我们可以在区域中心创造一个具有单个最大值的单峰图样,或在中心具有单个最小值的单峰图样;见图 2.8(b)。因此即使能孤立特定模式,初始条件仍能强烈影响极性。当多个模式可激发、且有一个由色散关系自然主导时,我们仍能通过适当的初始条件影响最终图样。若初始条件包含一个在不稳定带内且振幅充分大的模式,则该模式可以持续经过非线性阶段而主导其他不稳定模式,最终图样通常粗略地具有该波长。我们将在下一节更详细讨论这些。这些事实在生物学应用中也有高度相关的意义。
2.6 模式选择与色散关系(Mode Selection and the Dispersion Relation)
考虑一个典型的"纯形式"或最简单的色散关系——它给出增长因子 \(\lambda\) 关于波数或波长 \(\omega\) 的关系,例如图 2.5(b)——其中一波数带是线性不稳定的。再设区域是有限的,因此特征值谱是离散的。在上一节中我们看到几何与尺度在确定由线性理论预测的特定图样中起关键作用,这也由非线性系统的数值模拟所证实;亦见下一章中给出的结果。我们指出初始条件能起作用——例如确定图样的极性或特定图样是否涌现。若初始条件由均匀稳态附近的小随机扰动组成,那么很可能演化的图样是线性增长最大的那个。然而在许多发育问题中,图样启动的触发是尺度;下面几章中——特别是第 4 章——有几个例子。在其他发育情境中,从均匀稳态的扰动在空间域的一端启动,空间图样从那里发展,最终扩散到整个区域。因此对给定机制演化的具体图样可能临界地依赖于失稳如何启动。本节我们进一步研究这个问题,以建议在关键参数通过分叉值时,从哪些初始条件会演化出哪些图样。哪个图样会演化的问题,即模式选择,是一个不断重现的问题。下面的讨论虽然以反应扩散图样生成器为动机,但实际上相当一般,适用于任何产生类似类型色散关系的图样形成模型。
考虑基本色散关系 \(\lambda(\omega^2)\)(其中波长 \(\omega = 2\pi / k\),\(k\) 是波数),如图 2.15(a) 所示。现在取一维区域并依次考虑图 2.15(b)、(c) 与 (d) 所示的三种可能的图样启动方式。
首先考虑图 2.15(b) 的情形。这里的初始扰动在关于特征函数的展开中包含所有模式,因此图 2.15(a) 中不稳定波长带内的所有模式都被激发。具有最大 \(\lambda\)、\(\omega_2\) 的模式是增长最快的,它最终占主导。持续存在的稳态非齐次图样就是波长为 \(\omega_2\) 的那个。
在图 2.15(c) 中,我们设想区域以与生成空间图样所需时间相比缓慢的速率增长。本节后面将描述一个生长不微小的漫画系统;第 4 章会更详细地讨论生长对图样的影响——那里的相互作用至关重要。对小 \(L(t)\),区域不能容纳不稳定带内任何波长的波。当它达到 \(L_c\)(图样的临界域大小时),它能维持最小波长图样——波长 \(\omega_1\) 的图样。在 \(L(t)\) 增长到足以容纳下一个波数(即波长 \(\omega_1\))所需的时间内,\(\omega_1\) 的图样已经充分建立以主导非线性阶段。因此最终涌现的图样是基本波长 \(\omega_1\) 的图样。
行波启动的图样。 考虑图 2.15(d) 的情形——图样在区域一端启动;这里发生的事更微妙。我们预期最终图样将具有波长——位于色散关系预测的不稳定带内的某个波长。为看清如何一般地计算波数,让我们从无限一维区域与一般线性系统出发:
这里 \(J\) 是线性算子,例如与反应扩散方程线性形式相关的算子。对 (2.20) 而言,\(J = (\partial / \partial t) - \gamma A - D \nabla^2\) 而色散关系 \(\lambda(k)\) 如图 2.5(b) 或图 2.16(a) 那样的形式——其中 \(\omega\) 由 \(k\) 替换;即经典形式。线性系统 (2.69) 的一般解为
其中 \(A(k)\) 由初始条件 \(\mathbf{w}(x, 0)\) 的 Fourier 变换决定。由于我们关心最终结构而非瞬态,这里无需计算 \(A(k)\)。
设初始条件 \(\mathbf{w}(x, 0)\) 局限于原点 \(x = 0\) 附近的一个小有限区域,且图样从该区域向外传播。我们关心图 2.15(d) 第二幅所示的波动式图样生成。这意味着我们应考察远离原点的解的形态。换言之,我们应该关注 \(t \to \infty\) 时解的渐近形式,但 \(x\)、\(t\) 都大且 \(x/t = O(1)\)——即以速度 \(c = x / t\) 移动,处于"前沿"附近的位置。把 (2.70) 写成
此积分对 \(t \to \infty\) 的渐近估计通过把被积函数解析延拓到复 \(k\) 平面并用最陡下降法(见 Murray 1984 年专著第 3 章)给出
其中 \(J\) 是常数,\(k_0\)(现在为复数)由下式给出
解的渐近形式因此是
其中 \(K\) 是常数。
对大 \(t\),波的"前沿"大致是图样形成尾部与起始扰动的引导边缘之间的点——即 \(\mathbf{w}\) 既不增长也不衰减的点。因此该点满足
在"前沿"处波数为 \(\text{Re} k_0\),解的振荡频率 \(\omega\) 为
把"前沿"后所铺设的图样的波数记为 \(k_*\)。我们假设跨过"前沿"节点守恒,这意味着
(2.72)、(2.74) 与 (2.75) 三个方程现在确定 \(k_0\) 与我们感兴趣的量 \(c\) 与 \(k_*\)——分别是图样铺设速度与稳态图样波数。由于复变量的原因,这并不像看起来那样简单。Myerscough 与 Murray (1992) 针对细胞趋化系统的情形(也见第 4 章)使用了该技术,并发展了漫画色散关系以解析地求解这三个方程。他们把解析结果与精确的数值模拟做了比较;比较是良好的,并且由于得到了解析结果,定性上有用。该技术也由 Dee 与 Langer (1983) 数值地用于反应扩散机制。
增长区域中图样形成的动力学。 增长区域中图样随时间的演化可以相当复杂,特别是当区域增长与空间图样生成时间相当时,以及在二维或多维情形下。当尺度 \(\gamma\) 增大时,色散关系的形式可能具有高度相关的生物学含义——我们将在第 6 章中看到,例如考虑发育肢中软骨的形成。这里我们引入这一现象并讨论 \(\gamma\) 随时间增大时两类具体色散关系行为的一些含义。
增长区域中反应扩散方程的形式必须仔细推导,这将在第 4 章完成——第 4 章主要讨论图样化问题及增长对所形成图样的影响。这里我们仅考虑一个漫画形式以演示增长区域的某些时间依赖效应。对简化的漫画我们不应期望捕捉到所有可能的顺序空间图样——这确实如此。Crampin et al. (1999) 在对反应扩散系统在增长一维区域中的全面解析与数值研究中,把演化出的顺序图样分类。他们考虑不同的增长形式。他们用自相似性论证预测指数增长情形下频率加倍的发生,并展示增长可以作为一种增强图样稳健性的机制。Kulesa et al. (1996a, b)(也见第 4 章)证明牙齿原基(前体的牙齿)的顺序定位与颌的生长(由实验确定)密切相关。原基出现的正确顺序关键地依赖于图样形成过程与区域增长动力学的相互作用。Murray 与 Myerscough (1992) 在其关于蛇图样的研究中(我们在第 4 章详细讨论)采用了稍不同的处理增长区域的方法:他们考察了稳态方程(此情形为细胞趋化方程)解的分叉。
在图 2.8(a) 中我们看到随着尺度 \(\gamma\) 增大,色散曲线沿轴平移,连续激发波长更小的模式。图 2.16(a) 是这种行为的重示例。图 2.16(b) 是当尺度 \(\gamma\) 增大时色散关系可能具有的另一种行为。它们意味着增长区域中不同的图样生成场景。
首先考虑图 2.16(a) 的情形。这里对 \(\gamma = \gamma_1\),波长 \(\omega_1\) 的模式被激发并开始增长。随着区域增大,我们看到对 \(\gamma = \gamma_2\) 没有模式落在不稳定带内,因此图样衰减到空间均匀稳态。随着尺度进一步增大到 \(\gamma = \gamma_3\),我们看到具有波长 \(\omega_2\) 的图样被产生。因此图样形成实际上是一个离散过程——随着 \(\gamma\) 增大,越来越多结构被生成,但每次结构增加之间穿插着空间均匀的过渡阶段。图 2.17(a) 示意了 \(\gamma\) 增大时事件序列。
现在考虑图 2.16(b) 中色散关系关于尺度的依赖所暗示的行为。这里尺度的效果仅是扩大不稳定模式的带。主导模式随 \(\gamma\) 变化,因此从一种模式(在 \(\gamma = \gamma_1\) 时主导)连续演化到另一种模式(在 \(\gamma = \gamma_3\) 时主导)。图样这种动态发展如图 2.17(b) 所示。我们将在 §6.6 中看到图 2.16 与 2.17 的含义对发育肢中软骨图样如何形成有直接影响。
在将不同模型与实验比较时,不一定能在图样生成的给定时间选择来进行实验,因为我们通常不知道图样形成在胚胎发生中究竟在何时发生。当能够选定时,图样相似性是与理论比较的必要第一步。当不能选定时,图样的动态形式可能是重要的,并可能是决定哪种机制更合适的关键一步。我们将在第 4 章与第 6 章回顾这些评论。
图 2.18 展示了尺度 \(\gamma\) 增大时图样动力学形成的一个计算示例。
在这些模拟中,机制的图样生成时间小于代表性的增长时间,因为图样序列显然在破碎启动下一个图样之前就已经形成。这是图 2.16(b) 那样的色散关系行为的例子;即不存在空间均匀性的阶段。图 2.18(c) 所示的趋向周期加倍的现象是有趣的,迄今仍未完全解释。Arcuri 与 Murray (1986) 考虑了这一点以及增长区域中图样形成的其他方面。必须记住,后面这项研究是增长区域中反应扩散系统的漫画;参看第 4 章关于指数区域增长的精确公式化以及 Crampin et al. (1999) 的综合讨论。
2.7 单物种模型中图样的生成:云杉芽虫模型中的空间异质性(Pattern Generation with Single-Species Models: Spatial Heterogeneity with the Spruce Budworm Model)
我们上面看到,如果区域大小不够大——即 \(\gamma\) 太小——零通量边界条件下的反应扩散模型不能生成空间图样。零通量条件意味着反应扩散区域与外部环境隔离。现在我们考虑不同的边界条件,以考虑反应扩散区域外部区域的影响。具体而言,考虑单反应扩散方程的形式
并把模型放在生态学设定下;即 \(u\) 表示某物种的种群密度。这里 \(f(u)\) 是该物种的动力学,因此假设 \(f(0) = 0\)、\(f'(0) \neq 0\);若只有一个(正)稳态则 \(f(u_i) = 0\)(\(i = 1\)),若有三个则 \(i = 1, 2, 3\)。稍后我们将考虑种群动力学 \(f(u)\) 取云杉芽虫的形式——我们在 Volume I 第 1 章 §1.2 中详细研究过,它有三个稳态(如图 1.5(b) (Volume I))。扩散系数 \(D\) 是相关物种分散效率的度量。
我们首先考虑一维问题,区域 \(x \in (0, L)\),其外部对该物种完全不利。这意味着在区域边界上 \(u = 0\)。我们考虑的数学问题是
其中 \(u_0\) 是初始种群分布。我们要回答的问题是:这样的模型是否能支撑空间图样。
在空间均匀情形下 \(u = 0\) 与 \(u = u_2\) 不稳定,\(u_1\) 与 \(u_3\) 是稳定稳态。在没有扩散时,动力学意味着 \(u\) 趋于两个稳定稳态之一,究竟哪一个依赖于初始条件。在空间情形下,因此我们预期 \(u(x, t)\) 试图从 \(u = 0\) 增长,除了边界外。因为在边界 \(u_x \neq 0\),扩散的效应意味着有 \(u\) 从区域 \((0, L)\) 流出的通量。所以对小 \(u\) 而言存在两个竞争效应:动力学的增长与边界的损失。作为第一步,我们考察通过对 \(u = 0\) 线性化得到的线性问题。相关形式由 (2.77) 得
我们寻找如下形式的解:
由观察知它满足 \(x = 0, L\) 处的边界条件。代入 (2.78) 并比较 \(\sin(n \pi x / L)\) 的系数得 \(\lambda = [f'(0) - D (n \pi / L)^2]\),于是解为
其中 \(a_n\) 由初始条件 \(u_0(x)\) 的 Fourier 级数展开确定。本分析中我们不需要 \(a_n\)。由 (2.79) 我们看到 \(u\) 表达式中的主导模式是具有最大 \(\lambda\) 的那个,即 \(n = 1\),因为
因此若主导模式随 \(t \to \infty\) 趋于零,则其他所有模式也是如此。我们于是得到 \(u = 0\) 线性稳定的条件
以有量纲方式 \(D\) 的单位是 \(\text{cm}^2 \text{s}^{-1}\)(或 \(\text{km}^2 \text{yr}^{-1}\) 或其他我们感兴趣的尺度),\(f'(0)\) 的单位是 \(\text{s}^{-1}\)(因为对 \(u\) 小,\(f(u) \approx f'(0) u\),即线性出生率),合起来使 \(L_c\) 的单位为厘米。因此如果域大小 \(L\) 小于临界大小 \(L_c\),则 \(u \to 0\) 当 \(t \to \infty\) 且无空间结构演化。扩散系数越大,临界域大小越大;这与以下观察一致:\(D\) 增大时,从该区域流出的通量也增大。
增长区域中空间结构的情景是:随着区域增长且 \(L\) 刚过 \(L_c\) 时,\(u = 0\) 变得不稳定,第一模式
开始随时间增长。最终非线性效应开始起作用,\(u(x, t)\) 趋于一个稳态空间非齐次解 \(U(x)\),由 (2.77) 可知它由下式确定
其中撇号表示对 \(x\) 的导数。因为 \(f(U)\) 是非线性的,我们一般不能得到 \(U\) 的显式解。
由 (2.77) 与 (2.81) 中的空间对称性(\(x \to -x\) 保持方程不变),我们预期解关于中点 \(x = L / 2\) 对称。由于 \(u = 0\) 在边界,我们假设中点是最大值,记为 \(u_m\),此处 \(U' = 0\);现在可参看图 2.19(a)。将 (2.81) 乘以 \(U'\) 并对 \(x\) 从 0 到 \(L\) 积分,我们得到
因为 \(U = u_m\) 时 \(U' = 0\)。把原点变到 \(L / 2\) 比较方便,使 \(U'(0) = 0\)、\(U(0) = u_m\);即令 \(x \to x - L / 2\)。则
积分得
它隐式地给出解 \(U(x)\);典型解如图 2.19(b) 示意。\(u = 0\) 在 \(x = \pm L / 2\) 处的边界条件与上式给出
我们于是隐式地获得 \(u_m\) 作为 \(L\) 的函数。\(u_m\) 对 \(L\) 依赖关系的实际确定必须数值地完成。注意被积函数在 \(w = u_m\) 处有奇点,但由于平方根,它是可积的。一般 \(u_m\) 随 \(L\) 增大而增大,如图 2.19(b) 所示。
云杉芽虫的空间图样化。 现在考虑云杉芽虫的模型——其动力学我们在 Volume I 第 1 章中推导过。这里用 (1.8) 表示 \(f(u)\),(2.77) 变为
其中正常数 \(r\)、\(q\) 与该模型中有量纲参数定义的无量纲量相关;\(q\) 与环境容纳量成正比,\(r\) 与线性出生率成正比、与捕食强度成反比。种群动力学 \(f(u)\) 草图见图 2.20(a),其中参数处于给出三个正稳态 \(u_1\)、\(u_2\)、\(u_3\) 的参数域内,第一与第三线性稳定,第二个不稳定。以 (2.82) 定义的 \(F(u)\) 代入 (2.84),我们把 \(u_m\) 作为域大小 \(L\) 的函数。Ludwig et al. (1979) 数值上计算了它,形式如图 2.20(b) 所示;存在另一个临界长度 \(L_0\) 使得对 \(L > L_0\) 有多于一个解。我们下面分析此现象。
从生态学观点看,我们想知道的是最大种群能处于"爆发"(outbreak)状态的临界域大小 \(L_0\);即图 2.20(a) 中 \(u_m > u_2\)。这由 (2.84) 的数值积分确定,如图 2.20(b) 所示。当 \(L > L_0\) 时由图 2.20(b) 看到有三个可能的解,对应不同的 \(u_m\)。种群在庇护与爆发区的那些是稳定的,另一个——中间的——不稳定。哪个解出现依赖于初始条件。稍后我们将考虑该模型在芽虫控制中可能的生态学用途。在此之前,我们描述一种解析地确定临界域大小的有用技术。
确定临界域大小与最大种群的解析方法。 当对给定 \(L\) 有三个可能的 \(u_m\) 时,\(u_m(L)\) 的数值计算并不完全平凡。由于临界域大小 \(L_0\)(维持爆发的临界值)是我们在实际应用中所需要的重要且有用的量之一,我们现在利用 Lions (1982) 描述的一个思想,导出一个临时性的解析方法来获得它。
稳态问题由 (2.81) 定义。让我们重标度问题使域为 \(x \in (0, 1)\),设 \(x \to x / L\),则等价 \(U(x)\) 由下式确定
从图 2.19 知解定性上像正弦曲线。重标度后 \(x \in (0, 1)\),解因此定性上像 \(\sin(\pi x)\)。这意味着 \(U'' \approx -\pi^2 U\),于是上式蕴含
我们对使上式有 \(U\) 的三根的 \(L\) 值感兴趣;这对应于图 2.20(b) 中 \(L > L_0\) 的情形。因此要确定近似 \(L_0\),我们只需画上式如图 2.21,并确定使三解存在的 \(L\)。
对固定扩散系数 \(D\) 我们看到 \(U\) 的解如何随 \(L\) 变化。当 \(L\) 从 \(L \approx 0\) 增大时,第一个临界 \(L\)(\(L_c\))是当直线 \(D \pi^2 U / L^2\) 与 \(f(U)\) 相交时给出——即 \(D \pi^2 / L^2 = f'(0)\),与 (2.80) 一致。当 \(L\) 进一步增大,我们可以确定临界 \(L_0\),当 \(D \pi^2 U_0 / L_0^2\) 在图 2.21 中点 P 处与曲线 \(f(U)\) 相切时。简单地说就是确定使
有双重正根的 \(L\)。作为习题(习题 7),确定 \(L_0\) 作为 \(r\)、\(q\)、\(D\) 的函数,其中 \(f(U)\) 由 (2.85) 给出。对任何给定的 \(L\),该过程也近似地确定最大 \(U\)。从图 2.21 我们显然可以通过此过程得到与图 2.20(b) 类似的图形。这个简单过程对确定临界域大小是相当一般的——既适用于从零稳态分叉的结构,也适用于种群动力学有多个正稳态时能维持较大种群的域。
2.8 标量种群相互作用扩散方程中的空间图样与对流:生态学控制策略(Spatial Patterns in Scalar Population Interaction Diffusion Equations with Convection: Ecological Control Strategies)
此类模型的实际应用中感兴趣的域通常是二维的,因此我们必须考虑 (2.76)。同时,考虑到害虫昆虫,外部区域一般不是完全不利的,因此 \(u = 0\) 在边界上是过于严格的限制。这里我们简要考虑一维与二维问题,其中外部区域并非完全不利且存在恒定风向。这在许多昆虫扩散情形中很常见,并能以主要方式改变种群的空间分布。
为代数简单起见,假设二维域是一个矩形区域 \(B\),定义为 \(0 \leq x \leq a\)、\(0 \leq y \leq b\),面积为 \(A\)。完全不利的问题由下式给出
按上节同样步骤,对小 \(u\) 我们得到线性化问题的解为
因此临界域大小涉及 \(a\) 与 \(b\) 的任意组合,使
由于
我们得到空间图样化存在的一个不等式估计
一般二维域的估计由 Murray 与 Sperb (1983) 得到。显然这里的数学问题是:对所考虑的空间域求最小特征值。
在上述所有标量模型中获得的空间图样都只有单个最大值。在完全不利的边界条件下,这些是唯一能生成的图样类型。然而对两物种反应扩散系统,我们已看到能生成更多样的图样。自然要问:在单物种模型中,以一维为背景,能否得到类似的多峰图样?现在我们展示这样的图样如何能发生。
设有一个恒定风向 \(\mathbf{w}\),它对种群 \(u(\mathbf{r}, t)\) 的守恒方程贡献对流通量 \((\mathbf{w} \cdot \nabla) u\)。再设外部环境并非完全不利,于是合适的边界条件是
其中 \(\mathbf{n}\) 是域边界 \(\partial B\) 的单位法向。参数 \(h\) 衡量不利程度:\(h = \infty\) 意味着外部完全不利,而 \(h = 0\) 意味着封闭环境——即零通量边界。我们后面简要考虑后一情形。数学问题因此为
边界条件 (2.91) 与给定初始分布 \(u(\mathbf{r}, 0)\)。这里我们考虑一维问题,遵循 Murray 与 Sperb (1983) 的分析,他们也处理了二维类似以及此类问题更一般化的方面。
我们简要考虑的问题是一维系统——它定义稳态空间非齐次解 \(U(x)\)。由 (2.91) 与 (2.92),因为 \((\mathbf{w} \cdot \nabla) u = w_1 u_x\),
其中 \(w_1\) 是 \(\mathbf{w}\) 的 \(x\) 分量,\(U(x)\) 的数学问题是
我们用相平面分析来研究该问题,设
我们寻找相平面轨迹——由 (2.93) 中边界条件——从下面两条线之一上的任一点连接到另一条上的任一点
相平面情形如图 2.22(a) 与 (b) 所示,下面我们来展示。参看图 2.22(a)。由 (2.94) 可知任一点 \((U, V)\) 处 \(dV / dU\) 的符号。在曲线 \(V = f(U) / w_1\) 上,\(dV / dU = 0\),当 \((U, V)\) 分别在曲线之上(\(V > 0\))与之下时 \(dV / dU\) 为正与为负。因此若从边界线 \(V = h U\) 上某点 P 出发,定性上轨迹如 \(T_1\) 所示(其上 \(dV / dU < 0\) 处处成立)。若从 S 出发——尽管轨迹起始时 \(dV / dU < 0\)——它与 \(dV / dU = 0\) 的线相交并穿过到 \(dV / dU > 0\) 的区域,因此轨迹向上转折。\(T_2\)、\(T_3\)、\(T_4\) 都是可能的情形,依赖于参数与解轨迹的起点。\(T_3\) 与 \(T_4\) 不是满足 (2.94) 的解轨迹,因为它们没有终止于边界曲线 \(V = -h U\) 上。\(T_1\) 与 \(T_2\) 是允许的解路径,每个都具有单一最大值 \(U_m\)(轨迹穿过 \(V = 0\) 轴处)。
现在我们必须把相应的域长 \(L\) 与这些解轨迹联系起来。具体地,让我们聚焦于轨迹 \(T_2\)。记解的 \(V > 0\) 部分为 \(V_+(U)\),\(V < 0\) 部分为 \(V_-(U)\)。若我们现在对 (2.94) 的第一个方程从 \(U_Q\) 到 \(U_Q'\)(即 \(T_2\) 轨迹两端的 \(U\) 值)积分,则对应 \(T_2\) 代表的解的域长为
因此对每个允许的解轨迹我们可以得到相应解域的大小。解 \(U(x)\) 关于 \(x\) 的定性形式可由相轨迹推出——因为我们在其上每处知道 \(U\) 与 \(U'\)——并由上式计算域大小。对应于图 2.22(a) 的情形,\(U(x)\) 中只能有单个最大值。然而由于风对流项,对称性不再像上节那样——解关于 \(x\) 不再对称。
现在假设外部不利程度降低——即 (2.95) 中 \(h\) 减小——使边界线如图 2.22(b) 所示。按与图 2.22(a) 中解轨迹相同的方法进行,我们看到对应于轨迹 \(T\) 的解可能存在。画出相应的解 \(U(x)\),我们看到在域中现在有两个最大值:见图 2.22(c)。但在这种情况下,我们实际上是在把若干可能的解拼在一起。参看图 2.22(b),我们看到一条可能的解由轨迹 \(T\) 从 A 到 \(B_1\) 段代表。它具有单一最大值,对应域长 \(L_1\) 由 (2.96) 的等价式给出。因此若我们把域大小限制为 \(L_1\),这就是相应的解。然而若允许更大的 \(L\),则从 \(B_1\) 到 \(B_2\) 的延续是可能的——于是轨迹 \(A B_1 B_2\) 对应 (2.93) 的一个解。进一步增大 \(L\),我们可以包括轨迹的其余部分到 \(B_3\)。因此若域足够大,多峰解是可能的。对应于解路径 \(T\) 的长度 \(L\) 严格地按上述方法由 (2.96) 的等价式得到。
因此对足够小的 \(h\),当轨迹绕 \((U, V)\) 相平面中点 \(u_2\) 盘旋时,越来越复杂的结构是可能的。当然对这样的解存在 \(w_1 \neq 0\) 是必要的。若 \(w_1 = 0\),则解关于 \(U\) 轴对称,因此不可能有螺旋解。因此恒定风对复杂图样的形成是必要的。它也影响图样存在的临界域大小。一般结果与进一步分析由 Murray 与 Sperb (1983) 给出。
害虫控制策略。 现在考虑害虫控制问题。森林芽虫问题非常是一个二维空间问题。正如我们在 Volume I 第 1 章 §1.2 中所指出的,一个好的控制策略是把种群维持在庇护水平。正如我们在 Volume I §1.2 中所示的,战略上有利的是若能改变 (2.85) 中的动力学参数 \(r\)、\(q\) 使得只存在一个正稳态。但这不是真正生态上可行的。然而对更现实的空间问题,我们有进一步的、可能的方法把害虫水平保持在庇护范围内:即确保它们的空间域大小不允许处于爆发区。该论证对二维域同样成立,但为了示例我们首先考虑一维情形。
参看图 2.20(b)。如果空间区域被分成大小为 \(L < L_0\) 的若干区域——即最大 \(u_m\) 始终小于庇护种群水平 \(u_1\)——我们将达到目标。因此一个可能的策略是:把该区域用条带状喷洒,使未喷洒区域施加一个有效 \(L < L_0\),如图 2.23(a) 所示:分隔喷洒区域的实线竖线是完全不利外部的边界。
当然实际上不可能消灭所有偏离未处理区域的害虫——所以更现实的模型是采用 (2.91) 形式的边界条件,其中一些昆虫可以在处理区域外生存。需要解的关键数学问题是确定昆虫"间断"宽度 \(L_b\)。这必须使得相邻未处理区域的贡献不能通过间断扩散来传递足够数量的昆虫,以启动相邻斑块的爆发——即使 \(L < L_0\)(独立情形下的临界大小)。典型的种群分布定性上如图 2.23(a) 中虚线所示。
二维类似是显然的,但优化问题的求解更为复杂。首先必须对边界条件 (2.91) 确定能维持害虫爆发的临界域 \(A_0\)。然后必须确定喷洒条带的宽度。这不是平凡问题,但当然是一个可能的问题。这些问题的初步研究由 Ben-Yu et al. (1986) 完成。
虽然我们集中于芽虫问题,但这些技术与控制策略同样适用于其他害虫。昆虫扩散领域提出了若干非常重要的生态学问题,例如控制正在席卷美国西部的杀人蜂(见例如 Taylor 1977)以及非洲的蝗灾。Levin(见例如 1981a, b)已经对这些问题以及与空间异质生态学模型相关的其他问题进行了真实而实际的研究。"间断"控制策略以防止疾病流行的空间扩散将在第 13 章关于狂犬病空间扩散的讨论中详细讨论——该策略现已在使用。
在扩散驱动失稳的一个有趣的生态学应用中,Hastings et al. (1997) 调查了西方杂色蛾(western tussock moth)的爆发——曾被假设为蛾与寄生蜂之间捕食者-被捕食者相互作用的结果。他们分析的模型在定性上是一个两物种系统,其中被捕食者不能移动。我们曾在 Volume I 第 13 章 §13.7 中看到,让一定比例的被捕食者定居会引起反直觉的结果。Hastings et al. (1997) 也通过考虑一个相当一般的、典型捕食者-被捕食者相互作用(其中被捕食者不扩散)的系统而得到反直觉的结果。他们的新分析非常简单但极富启发性。他们的方法让人想起 §1.6 中关于可激发波的部分——在该部分中通过确定相平面中波轨迹是否闭合来确定波速等——并围绕每个空间位置存在三个可能稳态的可能性(其中两个稳定),可能存在一个跳跃或陡峭奇异区段将它们连接起来。然后他们反直觉地证明:被捕食者的空间分布在爆发域的边缘处具有最高密度。这种现象在西方杂色蛾爆发中已被观察到。理论在预测反直觉行为以及随后实验或观察确认中的作用是特别重要的——因为通常情况下(事实上经常如此)大量观察事实令人困惑而非有启发。Kareiva (1990) 的文章与理论-数据关系特别相关。
2.9 反应扩散系统中空间图样的不存在性:一般与特殊结果(Nonexistence of Spatial Patterns in Reaction Diffusion Systems: General and Particular Results)
零通量的标量一维反应扩散系统通常形式为
直观上唯一稳定的解是空间齐次解 \(u = u_0\)(即 \(f(u) = 0\) 的稳态解):若 \(f(u) = 0\) 有多于一个的稳定解,则哪个会得到取决于初始条件。可以证明任何空间非齐次稳态解都是不稳定的(分析在本书第一版中给出:它涉及估计特征值)。这个结论不能完全推广到多于一维空间变量的标量方程——正如 Matano (1979) 在 \(f(u)\) 有两个线性稳定稳态的情形中所表明的。然而能获得的空间图样依赖于具体域边界(具体来说是非凸边界)。例如,一个有足够窄颈部的哑铃形域是一个例子。图样依赖于通过颈部扩散性输送足够物质通量的困难程度——这种通量是把一个稳态改变为另一个稳态以实现齐次性所必需的。
我们在 §2.3 与 §2.4 中看到,零通量边界条件下的反应扩散系统能在参数与动力学满足适当条件时生成丰富的空间图样谱:关键点是反应物的扩散系数必须不同。这里我们证明:对一般多物种系统,若扩散充分大,则图样化会被破坏。这正是我们可能直观预期的——但若扩散系数不相等则并不显然。我们现在证明这一点。分析(如我们所示)给出了另一个涉及机制动力学弛豫时间的条件——这当然不是显然的。
在讨论多物种多维理论之前,从教学法上考虑一般一维两物种反应扩散系统是有益的:
零通量边界条件与初始条件
其中 \(u_0'(x)\) 与 \(v_0'(x)\) 在 \(x = 0, 1\) 上为零。定义能量积分为
这——除 1/2 外——是 (2.52) 中引入的异质性函数。对 \(E\) 关于 \(t\) 求导
并从 (2.98) 代入(对 \(x\) 求导后),分部积分得
由于零通量条件,被积项为零。
现在定义量 \(d\) 与 \(m\) 为
其中 \(\max_{u, v}\) 意味着 \(u\)、\(v\) 取所有可能解值时的最大。若我们愿意,可以用涉及 \(f\)、\(g\) 导数的某种范数定义 \(m\);对我们的结果而言这并不关键。由 \(dE/dt\) 的方程与这些定义,可得
其中我们用了如下结果(附录 A 中推导)
对 \(v\) 有类似不等式。
由不等式 (2.102) 我们现在看到:若 (2.101) 中最小扩散系数 \(d\) 充分大使得 \((4 m - 2 \pi^2 d) < 0\),则 \(dE/dt < 0\),这意味着 \(E \to 0\) 当 \(t \to \infty\)(因为 \(E(t) \geq 0\))。由 (2.100) 中 \(E\) 的定义,这意味着 \(u_x \to 0\)、\(v_x \to 0\),从而 \(t \to \infty\) 时 \(u\)、\(v\) 中具有空间齐次性。结果并不精确,因为有许多合适的 \(m\) 选择;(2.101) 仅是一个例子。结果的目的只是表明扩散有可能衰减所有空间异质性。我们下面简要评论此结果的生物学含义。
我们现在证明对一般反应扩散系统的类似结果。考虑
其中 \(u\)(分量 \(u_i\),\(i = 1, 2, \ldots, n\))是浓度或种群向量,\(D\) 是正扩散系数 \(D_i\)(\(i = 1, 2, \ldots, n\))的对角矩阵,\(f\) 是非线性动力学。我们证明的结果对带某些交叉扩散项的扩散矩阵同样成立;为简单起见这里只处理 (2.104)。\(u\) 的零通量边界条件与初始条件为
其中 \(\mathbf{n}\) 是域 \(B\) 边界的单位外法向。与前面一样,我们假设所有解 \(u\) 对所有 \(t \geq 0\) 有界。实用上,如果反应动力学存在约束集,则这实际上得到保证。
我们现在推广前面的分析;参看上面相应的等价步骤是有帮助的。定义能量 \(E(t)\) 为
其中范数
设 \(d\) 为矩阵 \(D\) 的最小特征值,在对角矩阵情形下即所有物种的最小扩散系数。现在定义
其中 \(u\) 取所有可能解值,\(\nabla_u\) 是关于 \(u\) 的梯度算子。
对 (2.106) 中 \(E(t)\) 求导,用分部积分、边界条件 (2.105) 与原系统 (2.104),用 \(\langle a, b \rangle\) 表示 \(a\)、\(b\) 的内积,我们得到
在附录 A 中我们证明当 \((\mathbf{n} \cdot \nabla) u = 0\) 在 \(\partial B\) 上时
其中 \(\mu\) 是
的最小正特征值,\(\phi\) 是标量。在 (2.108) 中用 (2.109) 得
因此再一次——若最小扩散系数充分大——这意味着 \(\nabla u \to 0\),从而所有空间图样随 \(t \to \infty\) 趋于零。
Othmer (1977) 指出:(2.101) 与 (2.107) 中定义的参数 \(m\) 是反应速率对 \(u\) 变化敏感性的度量——因为 \(1/m\) 是该机制的最短动力学弛豫时间。另一方面 \(1 / (2 \mu d)\) 是最长扩散时间的度量。于是结果 (2.110)(即 \(1/m > 1/(2 \mu d)\))意味着:若动力学的最短弛豫时间大于最长扩散时间,则所有空间图样化将随 \(t \to \infty\) 消失。该机制将由纯动力学支配。记住后者的解可以包括极限环振荡。
考虑一维情形,取一个典型的胚胎学感兴趣域,比如 \(L = O(1 \text{ mm})\)。取 \(d = O(10^{-6} \text{ cm}^2 \text{ s}^{-1})\),结果 (2.110) 蕴含若动力学的最短弛豫时间 \(1/m > L^2 / (2 \pi^2 d)\)——即 \(O(500 \text{ s})\) 的时间——则会出现齐次性。
考虑一般系统 (2.104),重标度使长度尺度为 1,扩散系数相对于某个 \(D_1\) 标度。现在回到 (2.10) 中那种形式——其中尺度 \(\gamma\) 出现在形式 \(\gamma f\) 的动力学前面。这对 (2.110) 条件的影响是现在产生 \(\gamma m - 2 \mu < 0\) 作为稳定性要求。我们立即从此形式看到:存在一个临界 \(\gamma\)(与域面积成正比——在一维即 (长度)\(^2\)),低于它则不能存在任何结构。当然这与我们在 §2.3 与 §2.4 中得到的结果类似。
我们应当重申,这里的结果给出的是定性界限,而不是针对模型机制相关各种参数的估计。合适的 \(m\) 的求值并不容易。在 §2.3 与 §2.4 中我们导出了当动力学为特定类别时,参数之间的具体定量关系以给出空间结构解。然而本节中的一般结果适用于所有类型的动力学——振荡型或其他——只要解有界。
本章中我们主要处理了反应或种群相互作用动力学——它们在我们所考虑的参数空间限制区域内,在没有扩散时不表现出振荡行为。我们可以问:当振荡动力学与扩散耦合时,能得到什么类型的空间结构?第 1 章中我们看到当区域无限时这种组合可以产生行波波列。若区域有限,我们可以预期某种"晃动"型的规则性——这是空间与时间上不稳定模式存在性的反映。这其实能发生,但并非总是如此。值得指出的一种情形是:经典 Lotka–Volterra 系统具有相等的物种扩散系数。Murray (1975) 证明在有限域中所有空间异质性必须消失(见习题 11)。
现在已有若干其他图样形成机制,超越反应扩散-趋化系统。Wittenberg (1993) 是关于发育中自组织模型最透彻的综述之一。他详细描述了这些模型并批判性地综述了若干种不同的机制,包括反应扩散-趋化系统、力学化学机制以及细胞自动机模型。
在下一章中我们将讨论若干具体的实际生物图样形成问题。在后续章节中我们将描述其他能生成空间图样的机制。一个被广泛研究的重要系统是用于生成细菌以及黏菌 amoebae 中聚集图样的反应扩散-趋化机制——其模型我们在 Volume I 第 11 章 §11.4 中推导过。用与上面我们对扩散驱动失稳所讨论的完全相同类型的分析,我们能展示这些模型方程中空间图样如何产生以及在什么参数条件下会发生(见习题 9)。如上所述,这些趋化系统随着对细菌图样兴趣的高涨而日益重要——这正是第 5 章被纳入的原因。我们将在第 4 章考虑增长区域对图样化影响时讨论细胞趋化机制的其他相当不同的应用——例如蛇身上观察到的复杂图样。
本章个人批注
本章是整本 Murray II 的技术核心:把所有后续章节关于具体生物图样(毛发、鳞片、羽毛、牙原基、软骨、血管网络、脑肿瘤、伤口愈合、流行病、狼领地等)的数学工具——反应扩散机制的 Turing 失稳——在数学上彻底展开。9 节分四个层次展开:(1) §2.1 在概念层面交代为何"图样"是发育生物学的核心问题,并介绍"位置信息"与"化学预图样"的现象学概念;(2) §2.2-§2.4 给出最简 Schnakenberg 模型的数学分析:方程无量纲化、零通量边界下的色散关系、Turing 空间的几何刻画、规则镶嵌(正方形/六边形/菱形/条纹)特征函数的显式构造;(3) §2.5-§2.6 把分析扩展到增长域、模式选择与不稳定带的几何解释;(4) §2.7-§2.9 转向生态学应用(云杉芽虫、害虫"间断"控制)与一般性"图样不存在性"判据(Othmer 弛豫时间比 \(1/m\) vs \(1/(2\mu d)\))。
几个让我特别留意的技术点:(a) 临界波数 \(k_c^2 = \gamma (d_c f_u + g_v) / (2 d_c)\)——尺度参数 \(\gamma\) 显式出现,因此"波数选择"本质上是"尺度选择"——这是后面第 4 章增长域中频率加倍、§6.6 软骨图样等应用的根基;(b) §2.3 的色散关系最大增长率 \(\text{Re} \lambda(k_m^2)\) 在 \(h(k^2) = 0\) 范围内由 \(k_m^2\) 取极值给出——这条曲线是判断"哪种图样会涌现"的统一判据;(c) §2.4 末 Berding 异质性函数 \(H\) 给出"在胚胎增长过程中图样复杂度单调上升"的数学表达——这是发育生物学的关键定性观察;(d) §2.9 的反定理告诉我们:当 \(d\) 充分大(扩散时间远短于弛豫时间)时,任何反应扩散机制的空间图样都会消失——这是反应扩散机制的"对立面",提醒我们图样化是有条件而非普遍的。
特别值得比较的概念是 §2.3 末尾对两种"激活剂-抑制剂"拓扑(\(f_v < 0, g_u > 0\) vs \(f_v > 0, g_u < 0\))的区分:前者导致"同相"图样(高低与高低共区域),后者导致"反相"图样(高 \(u\) 与低 \(v\) 共区域)。这两种拓扑对生物学的意义重大——例如在捕食者-被捕食者模型中,前者对应"捕食者扩散更快"产生同相团块;后者对应"被捕食者扩散更快"产生反相斑块。这种拓扑分类在后续章节的应用中会反复出现。
本章习题部分(§2.10 起)已经在原文末尾给出,包括:(1) 推导 (2.4) 与 (2.5) 的无量纲化;(2) 完整处理 \(f(u, v) = u^2 / v - b u\) 激活剂-抑制剂系统的 Turing 空间;(3) 含激活剂自抑制的 Turing 空间参数刻画;(4) 圆域的临界半径 \(R_c\);(5) 极限 \(d \to \infty\) 的不稳定模式上界;(6) 渔业区最优化 \(H > (\pi / 2) [D / (r - E)]^{1/2}\);(7) 解析地确定云杉芽虫临界长度 \(L_0\);(8) Lotka–Volterra 系统的 Lyapunov 论证(与 Murray 1975 等价);(9) 黏菌系统 \(n_t = D_n n_{xx} - \chi (n a_x)_x\), \(a_t = h n - k a + D_a a_{xx}\) 的无量纲化与色散关系推导;(10) 各向异性扩散 \(d_1, d_2, d_3, d_4\) 的临界比 \(d_c\);(11) 异质性能量分解 \(H = H_G^2 + H_S\);(12) 反应扩散线性化的特征多项式与必要条件 \(d_i\) 不全相等;(13) Jacobian 对称时 \(w \cdot A w < -\delta w \cdot w\) 与 \(E(t) = \int w \cdot w\) 衰减证明。
按 skill 约束我不涵盖习题——这些是源文为读者准备的练习题而非章节内容。
与上下章的衔接(一段话)
第 1 章把多物种反应扩散的时空图样(追逐-逃避行波、目标图案、螺旋波、可激发波、BZ 反应、灰松与红松竞争流行病学)作为引子,展示了"多一个物种"后时空行为的丰富性——但其重点在时间维度上的波与振荡。第 2 章是这种多物种数学在空间维度上的彻底展开:把"反应扩散"机制——特别是 Turing 失稳——从方程本身到色散关系到 Turing 空间几何到模式选择,全部系统化。本章为此后所有具体应用章节(§3 章动物毛皮图样、§4 章增长域、§5 章细菌趋化、§6 章力学图样形成、§8 章血管网络、§11 章脑肿瘤、§12 章神经图样等)铺设了共同的分析框架。可以说,本章是"工具箱",第 3 章起是"工具用出来"。
第 3 章直接是本章数学工具在动物毛皮图样(豹斑、斑马纹、鲑鱼斑点、蛇鳞、贝壳花纹、长颈鹿斑块等)上的应用:作者会展示 §2.3-§2.5 的色散关系与 Turing 空间如何与具体生物图样对接,以及观察中"斑点排列的规则性"为何能被预测。第 2 章与第 3 章的衔接因此是"理论工具"与"应用案例"的衔接——后者把前者的方程系数、尺度参数 \(\gamma\)、扩散系数比 \(d\) 与具体生物体的尺寸、形态发生素的扩散能力联系起来,并把 §2.4 的规则镶嵌特征函数(正方形、六边形、菱形)与动物毛皮中的点状/条纹图样对应起来。