第 5 章:细菌图样与化学趋性(Bacterial Patterns and Chemotaxis)
5.1 背景与实验结果(Background and Experimental Results)
研究细菌的理由十分直观:细菌是大量疾病的病原,也是绝大部分物质循环的执行者;对细菌复杂生物学的深入理解也将扩展其在其他领域的应用。本节关注的中心问题是:在局部相互作用下,细菌种群如何形成全局图样。多种实验条件下,许多菌株会在培养基上聚集形成(暂时)稳定的、规则度极高的宏观图样。要单凭实验解释这些图样并不容易,但基于已知生物学建立数学模型大体上能够解释。Berg 及其同事(Budsrene 和 Berg 1991, 1995 以及更早的文献)系统地研究了 Escherichia coli(E. coli)和 Salmonella typhimurium(S. typhimurium)这两类常见细菌形成的多种图样。E. coli 在人体肠道内大量存在,S. typhimurium 则常见于未完全烹饪的家禽与肉类;两者均能借长鞭毛(flagella)自推动运动(Berg 1983)。Berg 同时绘制了单个细菌随时间的运动轨迹,发现其近似于随机游走,因此可用常规 Fick 扩散描述其随机运动,其扩散系数已被实验测出。
许多细菌的关键特征是:在某些化学物质存在时,会朝高浓度的化学吸引剂(chemoattractant)或低浓度的趋避剂(repellent)作偏好运动;对该梯度的敏感性常随浓度水平而变。本章模型化工作主要关注化学趋性与这种敏感性问题。化学趋性的基本概念已在卷 I 第 11 章 4 节以及上一章讨论过;图样能否形成取决于细菌种群与化学动力学之间的恰当耦合,以及扩散性分散与趋化聚集之间的竞争。
Budrene 与 Berg(1991, 1995)针对 S. typhimurium 和 E. coli 完成了一系列实验:让菌落在三羧酸循环(TCA 循环)中间产物(尤其琥珀酸和延胡索酸)上生长或暴露其中,发现细菌可形成既有趣又高度规则的图样。两种实验方法对应三种图样形成机制:一种在液体培养基上、两种在半固体培养基(0.24% 水琼脂)上。所有实验中细菌均分泌天冬氨酸(aspartate)作为强力化学吸引剂。
液体实验(E. coli 与 S. typhimurium)。该组实验产生的图样相对简单、出现于分钟量级、持续约半小时后永久消失。两种图样由初始条件选择。最简情形:液体培养基中细菌均匀分布、TCA 中间产物少量添加;细菌在整个液面上形成大小相近的聚集群,图样常从某一区域起始并向外扩散(图 5.1(a))。第二种情形:细菌密度均匀、TCA 中间产物局部加于某点("原点"),随后细菌形成围绕原点的环状聚集群,环内随机分布(图 5.1(b))。液体实验中图样生成的时间尺度小于细菌繁殖时间,因此增殖对图样形成没有贡献;细菌对初始培养基中的任何化学物质(包括该刺激物)都不表现趋化反应;实验者 Berg(H.C. Berg, 个人通讯 1994)也确认流体动力学效应不引起观察到的图样。
半固体实验(E. coli 与 S. typhimurium)。最引人注目的图样见于半固体实验,特别是 E. coli。高密度接种的细菌放在含均匀分布刺激物的 petri 皿中(0.24% 水琼脂),刺激物(琥珀酸或延胡索酸)同时作为细菌的主要食物来源,浓度远高于液体实验。两三天后细菌经历 25–40 代繁殖,从接种点向外铺开,最终覆盖整个皿面,形成静止的高密度聚集群图样,聚集群之间为接近零密度的区域(图 5.2、图 5.3)。S. typhimurium 图样为同心环,可连续或点状断裂;E. coli 图样更复杂,相邻聚集群之间有明显的位相关系。常见的有向日葵型螺旋、径向条纹、径向点状、人字纹(chevrons)等。延时摄像揭示 S. typhimurium 与 E. coli 形成图样的运动学过程差异显著。
S. typhimurium 图样(图 5.2)。过程起始是低密度的菌苔从接种点向外铺展;稍后某一半径上(小于菌苔半径)出现一道高密度细菌环;再隔一段时间当菌苔进一步扩张时,在更大半径上出现第二道高密度环。环一旦形成即静止。环可能保持连续(图 5.2(a))或断裂成点状环(图 5.2(b))。相邻环之间的高密度聚集群没有明显的位相关系。
E. coli 图样(图 5.3)。这些更显著的图样在径向与角向相邻聚集群之间有明确的位相关系;该关系似乎由已存在的聚集群诱导后续聚集群的形成(Budsrene 和 Berg 1995)。取代初始菌苔的是一束高活跃度、高密度的可动细菌——称为"swarm ring"——从接种点向外扩张;环内细菌密度持续升高直至环失稳,部分细菌被"留下"成为聚集群。这些聚集群短时间内仍保持明亮且充满活跃细菌,但随后溶解、细菌重新加入 swarm ring。在聚集群原位置留下一团不知何故失去运动能力的细菌:这些不可动细菌即是图样的标记物。
图样机制关键取决于 swarm ring 扩张速度与聚集群溶解时间的相对关系。若溶解发生得快,聚集群被 swarm ring 拖着前进,留下径向条纹状不可动细菌(图 5.3(c)(d))。若溶解稍慢,溶解细菌重新加入 swarm ring,并在其重新加入处诱导新聚集群形成,导致径向点状图样。若溶解更慢,swarm ring 在溶解细菌重新加入前已失稳;环倾向于在已有聚集群之间形成新聚集群,造成向日葵螺旋型图样。
需要记住的是,与液体实验一样,petri 皿中放置的任何化学物质都不是化学吸引剂。但半固体实验图样的时间尺度长,足以容纳数代 E. coli 的繁殖,因此增殖在这里起重要作用。刺激物的消耗亦不可忽略,尤其在 swarm ring 图样中。由于所用底物中没有化学吸引剂,图 5.3 中 E. coli 图样不能用外部化学吸引剂来解释。然而化学吸引剂确实扮演了主要角色,因为细菌自身产生强力化学吸引剂天冬氨酸(Budsrene 和 Berg 1991)。此前一直假设 E. coli 和 S. typhimurium 中化学趋性现象只是为了引导细菌找到食物来源;正是 Budsrene 与 Berg(1991, 1995)的实验首次给出证据:细菌能产生并分泌一种化学吸引剂作为信号机制。这让人联想到粘菌 Dictyostelium discoideum,其细胞产生环腺苷酸(cAMP)作为化学趋性的聚集信号。本章模型主要聚焦扩散与朝向内源性产生的化学吸引剂的化学趋性,以及它们如何相互作用产生细菌图样。
所有图样都形成于二维 petri 皿上;显然若在三维培养基中图样复杂度会更高,但实验上拍摄三维图样极有挑战。本章只对二维图样建模,并展示能够反映生物学的模型确实能再现实验观察到的图样。化学趋性现象能产生复杂且多变的几何图样,单凭实验直觉难以理解个体细菌间的相互作用如何形成这些复杂几何。在这类情况下,当生物学直觉难以提供充分解释时,数学建模能起到重要作用——甚至是关键作用。要理解图样需回答许多问题,常涉及生物学假设与参数估计的细节:例如,仅靠朝向内源性化学吸引剂的扩散与化学趋性是否足以解释图样形成?化学吸引剂刺激物的定量作用是什么、其产生速率要多快?还有哪些图样可能?要获得这些图样应改变实验中的哪些关键要素?Tyson(1996)给出了关于生物学、模型与数值方案的综述。
本章详细讨论的内容基于一系列针对这些特定图样的理论研究(Woodward 等 1995, Tyson 1996, Murray 等 1998, Tyson 等 1999)。这些作者提出了紧贴已知生物学的数学模型。Woodward 等(1995)——与实验者 Berg 和 Budsrene 的合作——研究了 S. typhimurium 形成的较简单图样并给出观察到的自组织机制解释。Ben-Jacob 及其同事(Ben-Jacob 等 1995, 2000 等)也理论与实验地研究了多种细菌:他们得到的许多图样复杂而显著;图样显然依赖于参数值与实验条件。自然界中细菌面对多种条件(敌对与友善),为应对环境因素发展了相应策略。这些策略涉及协同通信并影响它们形成的图样类型。Ben-Jacob(1997 等)研究了化学趋性反馈等可能的通信过程的影响;意料之中地,模型化学趋性系统中纳入协同性后图样复杂度谱更广。Mimura 等(Mimura 和 Tsujikawa 1996, Matsushita 等 1998, 1999, Mimura 等 2000 等)系统的解析与广泛数值研究特别凸显了反应-扩散-化学趋性系统能展现的复杂解行为。例如 Mimura 等(2000)将图样分为五类,并提出除一类外,形态多样性均可由反应-扩散模型生成。Mimura 和 Tsujikawa(1996)考虑了带种群增长的扩散-化学趋性模型,并在扩散与化学趋性较小时推导出图样随时间演化的方程。本章讨论 S. typhimurium 和 E. coli 形成的特定细菌图样,并非常简短地介绍 Bacillus subtilis 所展现的、相当不同的图样。
这些细菌图样远比化学趋性菌株在含吸引剂的营养培养基上形成的图样复杂得多(如 Agladze 等 1993)。它们也不同于粘菌 D. discoideum 的聚集细胞行波——E. coli 和 S. typhimurium 等形成的结构仅在时间上稳定。化学趋性的空间图样潜力已在多种生物学背景下被利用。涉及化学趋性的数学模型(与反应-扩散模型(第 2、3 章)和机械化学模型(第 6、7 章)一起)只是空间图样发展积分-微分方程模型这一大领域的一部分。Keller-Segel(1970)最早提出针对粘菌 D. discoideum 的连续统图样形成机制,本书卷 I 第 11 章 4 节已讨论。Dallon 和 Othmer(1997)给出了更具生物学基础、适合细胞信号的离散聚集模型。Othmer 和 Schaap(1998)给出了粘菌发育中 cAMP 振荡信号传导的详尽综述。自 Keller 和 Segel(1970, 1971)的开创性工作以来,针对这些图样的建模工作已经很多——例如 Ben-Jacob 等(1995)在模型中加入了天冬氨酸产生的阈值行为与细胞分泌的废物场,得到了与实验观察的 E. coli 图样相似的空间图样。Brenner 等(1998)对半固体培养基中 E. coli swarm ring 形成的模型机制作了一维分析,从图样形成角度研究了方程各项的相对重要性并获得一些解析结果,例如导出了给定域中聚集群数关于模型参数的表达式。
化学趋性在多种实际现象中起重要作用,例如伤口愈合(第 10 章)、癌症生长(第 11 章)以及白细胞响应细菌性炎症的运动(如 Lauffenburger 和 Kennedy 1983, Alt 和 Lauffenburger 1987)。以往研究细胞种群非恒定的工作相对较少;例外之一是 Kennedy 和 Aris(1980)的行波模型——其中细菌既繁殖又死亡同时迁移。化学趋性(或机械化学理论中的趋触性 haptotaxis,第 6 章讨论)的存在会导致比单纯反应-扩散更广泛的图样种类。当图样过程中发生显著增长时,图样谱自然更广。许多细菌所形成的特定图样敏感地依赖于参数以及实验中的实际条件——包括初始条件。因此细菌图样可能的一项实际应用是作为污染的定量测量。
5.2 半固体实验中 E. coli 的模型机制(Model Mechanism for E. coli in the Semi-Solid Experiments)
本节构造 Budsrene 和 Berg(1991, 1995)实验中控制细菌图样形成过程的生物学机制。先考虑 E. coli 半固体实验的模型。实验中的关键角色是细菌(当然)、化学吸引剂(天冬氨酸)和刺激物(琥珀酸或延胡索酸),因此引入三个变量:细菌密度 n、化学吸引剂浓度 c、刺激物浓度 s。细菌作扩散运动、朝着化学吸引剂梯度作趋化运动、增殖、并变成不可动;不可动细胞可视为"死亡",便于建模。化学吸引剂扩散、被细菌产生和摄取;刺激物扩散并被细菌消耗。建模的第一步通常是写下"词方程"描述你认为正在发生的事情。基于 5.1 节描述的生物过程,建议下述由三个守恒方程构成的模型:
n 的变化率 = n 的扩散 + n 向 c 的趋化 + n 的增殖(生长与死亡) (5.1)
c 的变化率 = c 的扩散 + n 产生 c − n 摄取 c (5.2)
s 的变化率 = s 的扩散 − n 摄取 s (5.3)
建模的关键是如何量化各项。详细讨论见 Tyson(1996),她还通过广泛的文献综述估计各项参数。本章主要讨论她的工作以及 Tyson 等(1999)和 Murray 等(1998)的工作。
扩散。方程 (5.2) 和 (5.3) 中化学吸引剂与刺激物的扩散项直接明了。这些化学物质按简单 Fick 扩散,扩散系数分别为 D_c 和 D_s,实验估计值(Berg 个人通讯 1993)为 D_c ≈ D_s ≈ 9 × 10⁻⁶ cm² s⁻¹。所有参数估计列于表 5.1。假设细菌种群也按 Fick 方式扩散,但其扩散系数 D_n 的估计不太直接。Berg(1983)在书中给出由细胞个体运动推出的扩散常数表达式,得到 D_n ≈ 2 × 10⁻⁶ cm² s⁻¹。Phillips 等(1994)汇编了过去 10 年文献中报告的扩散系数:上限约 D_n ≈ 1.9 (±0.9) × 10⁻⁴ cm² s⁻¹,下限约 D_n ≈ 1–10 × 10⁻⁷ cm² s⁻¹。最近的测量则落在 D_n ≈ 1–3 × 10⁻⁶ cm² s⁻¹,与 Berg(1983)的理论值一致。
选择模型合适的扩散系数时需注意测量时细菌可用的空间维度。本章所建模的实验中 petri 皿内液体或琼脂混合物的深度约 1.8 mm。在 10 μm 毛细管阵列中受限运动的 E. coli 扩散系数为 5.2 × 10⁻⁶ cm² s⁻¹;在 50 μm 毛细管中则为 2.6 × 10⁻⁶ cm² s⁻¹(Berg 和 Turner 1990)。正如所预期,大毛细管中的扩散较慢。估计细胞或细菌(实际上除化学物质外的任何东西)的扩散系数始终是个难题(参见第 9、10、11、13、14 章)。Sherratt 等(1993b)在其真核细胞运动研究中给出了一种理论方法;Ford 和 Lauffenburger(1991)以及 Sherratt(1994)从受体层面动力学出发发展了估计细菌扩散系数的模型。
对本章建模而言,无需纳入受体层面的细节,因为实验者相信绝对化学吸引剂浓度并不影响细菌扩散系数。因此假设细胞以常数扩散系数 D_n = 2–4 × 10⁻⁶ cm² s⁻¹ 扩散。
化学趋性。化学趋性(卷 I 第 11 章)涉及生物朝浓度梯度方向的定向运动,类似于负扩散。但细胞扩散仅依赖于其密度梯度,而化学趋性依赖于细胞、化学吸引剂与化学吸引剂梯度之间的相互作用。守恒方程中化学趋性项的一般形式为化学趋性通量的散度:
其中 J_c 为化学趋性通量,χ(n, c) 为化学趋性响应函数(尚待确定),n 和 c 分别为细胞密度与化学吸引剂浓度。大量研究致力于找到 χ(n, c) 的生物学准确表达式;Ford 和 Lauffenburger(1991)综述了尝试过的主要函数形式。与扩散项类似,χ(n, c) 的形式要么从细胞行为的微观描述向上构造,要么通过宏观种群实验曲线拟合。各种方法的综合表明,Lapidus 和 Schiller(1976)的宏观形式
是一个不错的选择——其中 k 为参数。这与实验数据吻合最好,特别是 Dahlquist 等(1972)专门针对化学趋性响应函数形式设计的实验。有趣的是,其他研究者纳入全部受体层面复杂性并未比 Lapidus 和 Schiller(1976)的结果带来显著改进。受体模型的主要优势是参数可直接对应于细菌可实验观测到的物理化学性质。但种群研究不一定要纳入这些细节。本章建模与分析关注 E. coli 和 S. typhimurium 种群的整体行为,因此最适用宏观导出的化学趋性系数。基于上述形式,Woodward 等(1995)选取:
参数 k_1 和 k_2 由 Dahlquist 等(1972)的实验结果确定,给出 k_2 = 5 × 10⁻⁶ M,k_1 = 3.9 × 10⁻⁹ M·cm²·s⁻¹。
细胞增殖。增殖项涉及生长与死亡。由 Budsrene 和 Berg(1995),细胞以受琥珀酸供应影响的常速率生长。在半固体实验中刺激物(琥珀酸或延胡索酸)是细菌的主要碳源(营养物),而液体实验中营养由其他形式提供且不受限制。因此假设增殖项为:
其中 k 均为参数。直观上这一形式合理:表现为 logistic 增长,其承载能力依赖于营养 s 的供应。当细菌密度低于承载能力时 (5.6) 为正,种群增加;n 高于承载能力时该项为负,种群净减少。该形式隐含假设每细胞的死亡率正比于 n;另一种稍不合理的可能是恒定的每细胞死亡率。
化学吸引剂与刺激物的产生和消耗。模型含一个产生项(化学吸引剂产生)和两个消耗项(化学吸引剂摄取与刺激物摄取)。由于数据不足,必须依赖直觉判断化学吸引剂的产生与摄取的合理形式。对营养消耗,预期营养以与细胞出现速率成比例的速率从培养基中消失。由于细胞线性出生率取为 k_3 k_4 s²/(k_9 + s²),提示下述营养消耗形式:
其中 k 均为参数。消耗形式具有 sigmoid 特征。细胞对化学吸引剂的消耗可能也有类似 sigmoid 特征。但化学物质并非生长所必需,实验中产生的量可能很少。因此简单假设细胞遇到天冬氨酸分子即摄取,建议化学吸引剂消耗形式为:
参数 k_7 标定。化学吸引剂产生项也没有被详细测量过。仅知(个人通讯 Berg 1993)化学吸引剂的产生量随营养浓度增大并可能随时间饱和——这提示一个饱和函数。具体选择:
其中 k_5、k_6 为其他参数。另一种合理的不饱和可能形式是
事实上这两种形式都能产生所需的图样,需要进一步实验来区分或找到另一种函数。关键特征是 n 小时的行为——此处产生函数的导数必须为正。
半固体培养基中细菌图样形成的数学模型。将上述函数形式代入词方程系统 (5.1)–(5.3) 得到:
其中 n、c、s 分别为细胞密度、化学吸引剂浓度和刺激物浓度。模型有三个扩散系数、三个初始值(t = 0 时的 n、c、s)和九个 k 参数。其中部分参数已有估计,部分可合理估计,另一些以现有生物学知识尚无法获知。下面讨论参数估计。
液体培养基中细菌图样形成的数学模型。合理假设化学吸引剂的产生与细胞的化学趋性响应在液体与半固体实验中由相同函数控制;两组实验的差异更多在于时间尺度与刺激物的作用。如前所述,细胞在液体实验的时间过程中来不及增殖,因此该模型不含生长项;此外刺激物不是细胞的主要食物来源(外部提供),其消耗可忽略。这显示液体实验模型只是半固体实验模型的特例。消去细胞生长、化学吸引剂降解和刺激物消耗,留下较简单的三方程模型:
比半固体情形少了未知参数。最后一个方程与其他方程不耦合;在最简液体实验中(刺激物在培养基中均匀分布),第三个方程也可省去。
参数估计。如前所述,我们有部分参数值,可从文献推导其他估计,也有一些参数组合的估计。乘积 k_3 k_4 为最大瞬时生长率,常由世代时间 t_gen 确定:
E. coli 实验的世代时间约为 2 小时,相应瞬时生长率为 0.35/小时。组合 k_3 k_4/k_8 称为产量系数,实验者按细菌生成重量/底物消耗重量计算。组合 k_3 k_4/k_7 为细菌作为化学吸引剂(氮源)函数的产量系数。参数 k_1 和 k_2 由 Dahlquist 等(1972)的细胞漂移速度和化学趋性梯度测量计算。
表 5.1 列出了从文献得到的维度参数估计值(E. coli 与 S. typhimurium 模型 (5.11)–(5.16))。其余 k 参数尚无估计。但因为后续将对方程作无量纲分析,只需对参数的特定组合有估计即可。这些估计将随同数值解的图例一起给出。
直觉解释图样形成机制。在分析任何生物学问题的模型前,直观地思考特定情形下会发生什么是很有教益的。记住感兴趣的域是有限的——实验 petri 皿的域。考虑完整模型 (5.11)–(5.13)。在 (5.13) 中摄取项是扩散方程的汇,因此随时间趋于无穷,营养浓度 s 趋于零。这又由 (5.11) 和 (5.12) 意味着最终细胞生长和化学吸引剂产生都趋于零,而化学吸引剂的消耗与细胞死亡继续。因此细胞密度与化学吸引剂浓度也随时间趋于零。该模型中唯一稳态是 (n, c, s) = (0, 0, 0)。但这并非我们关心的情况。其含义是不能在非零均匀稳态附近作通常的线性分析;必须考虑方程的动力学解。
对液体实验的模型系统 (5.14)–(5.16) 而言。最后一个方程意味着最终刺激物空间均匀,因为只是经典扩散方程,会随时间抹平所有空间异质性。容易看出存在均匀稳态 (n, c, s) = (n_0, 0, 0),其中 n_0 是细胞的初始浓度,是另一个可由实验变化的参数。由 (5.15) 源项始终为正,c 将无限增长。此时浓度高到一定程度会显著降低 (5.14) 中的化学趋性响应,最终简单扩散主导、解变成时间无关且空间均匀。因此对液体实验也必须考虑解的动力学演化。在均匀初始条件的扰动下,化学吸引剂浓度将持续增长。要从模型分析中得到任何有意义的结果,必须在扰动与化学趋性响应饱和之间的时间窗口内寻找图样。
直接观察液体模型 (5.14)–(5.16) 的物理扩散-化学趋性系统如何产生高密度细胞聚集群的出现与消失是直观的。t = 0 时细胞开始分泌化学吸引剂,因细胞随机分布,某些区域化学吸引剂浓度高于其他区域。由于化学趋性,这些细胞浓度较高的群吸引邻近细胞,增加局部细胞密度并降低周围密度。新加入聚集群的细胞也产生化学吸引剂,使局部浓度的增长速率高于周围低密度细胞的贡献。于是细胞密度与化学吸引剂浓度的峰与谷被放大。但这并非全部,因为细胞与化学物的扩散也参与,产生分散效应,抑制聚集的化学趋性过程或抹平峰谷——或首先阻止它们形成。这就是经典的局部激活与侧向抑制情形,而哪种过程(聚集或分散)占优,依赖于通过 n_0 的各种参数与初始条件之间的精细关系。
5.3 液相模型的直观分析(Liquid Phase Model: Intuitive Analysis)
5.2 节关于液体实验及其模型系统的讨论表明:随机排列的点状图样可能短期内出现在液体培养基实验中,但最终聚集群消退、均匀性再现。原因大概在于化学趋性响应的饱和。由于细胞的化学吸引剂产生没有任何形式的降解或抑制来对消,皿中化学吸引剂量持续增加,结果化学趋性响应最终饱和、扩散占优。
通常的均匀稳态附近线性分析在此不可用,必须发展不同的分析来研究图样形成动力学。Tyson 等(1999)发展的方法本质上是直观的而非精确的,但仍有信息量且具定性预测能力,能解释化学趋性模型与实验中随机或环形排列点状瞬态图样如何出现。同时给出一些数值解与解析预测作比较。所有分析与模拟均采用零通量边界条件,反映实验情形。
从液体实验的最简模型开始——半固体相模型中细胞增殖、化学吸引剂降解、刺激物均匀分布(既不被消耗也不降解,因此只是另一个参数)置零。则 (5.14)–(5.16) 化为:
为简化,所作分析针对一维情形,∇² = ∂²/∂x²。虽然针对一维域分析,结果可只作小修改扩展到二维(如同第 2 章反应扩散图样)。无量纲化方程时设:
去掉星号作代数简化后,无量纲方程为:
量 n_0 和 s_0(本质上可由实验变化的参数)分别为平均初始细胞密度与刺激物浓度。由于刺激物 s 既不消耗也不降解,w = 1。另外,液体实验中细胞种群既无生长也无死亡,存在守恒:
其中 u_0 为平均初始无量纲细胞密度,当 n_0 为初始均匀密度时 u_0 = 1。无量纲参数值见表 5.2,由表 5.1 的维度参数值计算;参数 μ 未知。实验的初始条件为非零均匀细胞密度和零化学吸引剂浓度。希望找到 (5.20) 和 (5.21) 的空间异质解,初期增长、随时间衰减。
(5.20) 和 (5.21) 的非平凡(即 u_0 ≠ 0)空间独立解为
初始条件 u(x, 0) = 1, v(x, 0) = 0。假设 (5.20) 和 (5.21) 的初始条件为关于初始细胞密度的 O(ε) 小随机扰动,寻找如下形式的解:
其中 0 < ε ≪ 1,k 为对应于随机初始条件 Fourier 级数的波数。为近似实际实验情形(初始化学吸引剂浓度恰为零)设 g(0) = 0;为说明取 f(0) = 1。寻找叠加在时间增长解上的空间变化解。
由于在有限域上寻找零通量边界条件的解,只有正弦(余弦)解涉及整数模式 m(参见第 2 章 2.4 节),它们与波数 k 通过
关联,其中 l 为无量纲域长。将 (5.23) 代入 (5.20) 和 (5.21) 并按通常方法在 ε 上线性化,对每个 k 得到 O(ε) 方程:
其中 τ = μ + 1 + t(注 τ_0 = τ_{t=0} = μ + 1 > 0),F(τ) ≡ f(t),G(τ) ≡ g(t)。(5.25) 右端第二项的系数是唯一依赖于化学趋性参数 α 的。
解析问题在于如何确定 F(τ) 和 G(τ) 的解行为。由 (5.25) 可见 τ → ∞ 时 G(τ) 的系数趋于零,F(τ) 的解简化为衰减指数。一旦发生这种情况,(5.26) 的解也给出衰减指数。因此解形式 (5.23) 解释了随时间最终图样的消失。现在考虑从初始扰动产生的空间图样增长。
在 τ 接近 τ_0(即 t 小)时,可结合 (5.25) 和 (5.26) 给出关于细胞密度图样振幅 F(τ) 的单一二阶常微分方程:
该方程有基于合流超几何函数的精确解,但对看清解行为实质上没什么用——而这正是我们想要的。改用启发性与定性推理。否则若一开始就直接数值求解系统,也难以看清实际发生的事。先假设 (5.27) 的二阶常微分方程系数变化速率远慢于函数本身及其导数。这让我们在 τ 的小区间内把 (5.27) 与常系数二阶方程比较。记 (5.27) 的系数为 D(τ) 和 N(τ):
如前述,N(τ) 中最后一项是唯一包含无量纲参数 α(含化学趋性参数)的项。参数 μ 也仅显含于该项,但其影响也包含在 τ 的表达式中,因此其效应不易分离。注意 D(τ) 对所有 τ > 0 为正,N(τ) 在 τ 接近 τ_0 = 1 + μ(即维度时间 t = 0 处)时可正可负可零。对足够大的 τ,N(τ) > 0。暂且视 N(τ) 和 D(τ) 为常数,则 (5.28) 的解记为 F̃ 为形式上的:
其中 L 为积分常数。在 τ 的小区间内可认为 N(τ) 和 D(τ) 近似常数。关于解 F̃,由于 D(τ) > 0 对所有 τ 成立,有 Re(λ_−) < 0 对所有 τ。但 Re(λ_+) 的符号可随 N(τ) 的符号而变。
此时主要关注化学趋性系数 α 和波数 k(或模式 m)如何改变解。考虑增大 α 的效应。若 α 足够大,N(τ) 对 τ 的小值(含 τ_0)为负。τ 增大时,N(τ) 经过零变为正。对 λ_+ 的效应是使本征值实部对 τ 足够小为正、对 τ 较大为负。λ_+ 经过零的 τ = τ_crit 与 N(τ) 变为零的点相同。因此对 τ 小,F̃ 的一项为增长指数;对 τ 较大两项均为衰减指数。因此预测 α 有去稳作用:随着 α 增大图样增长的可能性增加,因为它使 N(τ) 在 τ < τ_crit 时更负。当然从 (5.20) 可以预测 α 的去稳效应,这里只是其去稳影响的量化。
回顾模式 m² = k²l²/π²。对足够大的 k²,(5.29) 中 N(τ) 对所有 τ 为正,得到严格衰减的解。这引出预测:最低频模式最不稳定,且不会看到频率大于以下值的模式:
通过解 N(τ_0) = 0 给出。随时间增加不稳定的模式越来越少;τ → ∞ 时仅 0 附近缩小区间内的模式仍不稳定。任意时刻最快增长的波数 K_grow 可通过令 λ(k²) = 0 即 N(τ) = 0 解 k² 得到:
若 (5.28) 中常系数近似在小但有限的 τ 区间内合理有效,则在连续的 Δτ 区间内计算的解 F̃ 序列可产生一个解,先增大至最大值然后对所有之后 τ 衰减。增长阶段发生在 λ_+ 为正时。对方程作数值求解证实了这一预期行为;图 5.4 是参数值给定时的一次模拟。
F(τ) 最大值 F_max 的真正位置 τ = τ_crit 可能接近解析给出的 τ̃_crit:
Tyson 等(1999)将其与数值得到的 τ_crit 比较:两者很接近,且 k²(正比于 m²)和 α 越大吻合越好。τ̃_crit 与 τ_crit 的差异反映 τ_crit 处 d²F/dτ² 的大小。按定义 τ_crit 为 dF/dτ = 0 的时刻,因此 (5.28) 简化为
由于函数在极大处二阶导数为负,知 N(τ_crit) 为正。因此 τ 已经超过 N(τ) 改变符号的点,τ̃_crit 给出 τ_crit 的下界估计。τ̃_crit 与 τ_crit 相当接近,提示 N(τ_crit) 可能接近零,反过来表明 d²F/dτ² 在极大 F_max 处数值上较小。
由此鼓励省略二阶导数项后求解 (5.28)。经一些直接代数运算后,所得到的关于振幅 F(τ) 的解满足 F(τ_0) = 1 (= f(0)):
其一阶与二阶方程解 F(τ) 和 F₁(τ) 的曲线示于图 5.5 和图 5.6。
首先注意到两个函数在高度上差异显著。除此之外两者有许多相似之处:峰值大致出现在 τ 相同的位置,峰区间(定义为 F(τ) > F(τ_0) 的时间)也大致相同,尤其是较低频时;两曲线都向左偏斜。α 增大导致 F_max 和 F₁_max 都大幅增大。两解对所考察不同模式的行为也相似。k² 越大,τ_crit 越早到达,F 或 F₁ 大于初始值 F_0 的区间越短。
若将 F(τ) 和 F₁(τ) 数据归一化到区间 [0, 1],则两解几乎一一对应(图 5.6)。近似解与数值解之间的主要差别仅是缩放因子。该缩放因子很大,提示二阶导数项在极大邻域外并不小。
至此对 (5.28) 的解 F(τ) 的行为已获得直观理解。也获得了近似解析解 F₁(τ)((5.33)),可用来预测改变各种参数的效果。
5.4 解析结果的解读与数值解(Interpretation of the Analytical Results and Numerical Solutions)
特别关注模型对将形成聚集群数目以及它们可见多长时间(即 F(τ) 足够大)的预测。若非线性效应不太强,可预期聚集群数将由对应于各模式的解的联合效应决定。
图 5.7 给出一些数值结果。注意存在一个波数(模式)比其他都达到更高振幅,称该波数为 k_max;此处 k_max = 2.20。也注意到每个大于 k_max 的波数初始时比 k_max 有稍高增长率,但这些高频很快开始衰减,而 k_max 图样的振幅仍在快速增长。可以推测 k_max 模式可能是支配非线性偏微分方程完整系统的第一个模式。
一旦 k_max 解开始衰减,小波数的解按降序依次变得最大。每个 k < k_max 解的振幅一旦超过下一个较高模式即处于衰减状态。因此随时间 t 增加,所观察图样的波数应连续减小,伴随振幅减小。这与生物学上观察到的聚集群合并以及图样的最终消散一致。对这些图而言最大增长率的波数为 K_grow = 5.47,因此 k_max 比 K_grow 小 2 倍。在 Tyson(1996)和 Tyson 等(1999)研究中观察到的所有数值解里,k_max 都显著小于 K_grow。
一维数值模拟结果。比较线性理论的预测与偏微分方程的实际解行为。所有模拟均使用零通量边界条件。化学吸引剂浓度的初始条件在全域上为零,细胞密度初始为关于 u_0 = 1 的随机扰动。除其他数值检查外,所有解与守恒律 (5.35) 比较,因为液体模型中既无生长也无死亡。
图 5.8 显示 α = 80 时代表性的时间序列。序列在细胞密度峰数变化小、仅图样振幅减少的时间点截断。每幅图左列为各时刻 τ 的细胞密度剖面,右列为对应的功率谱密度。后者的密度轴被限制在 τ = τ_0 初始功率谱密度均值之上,从而突出增长的图样模式。
如 (5.32) 所预测,功率谱密度图显示高于 K = 10.9 的模式不增长。同时非零模式的分布随时间减少。在实际细胞密度分布中,最初观察到的图样有许多峰,峰数随时间减少。预测 k_max 是将支配解的空间图样模式,与图中所示差两倍:实际支配解的模式为 k_max ≈ 1.1,而预测值为 2.2。
二维数值模拟结果。二维情形下得到与一维情形相似的行为。同样以均匀分布附近的小随机扰动作为初始条件,生成的图样由随机排列的点组成(图 5.8)。图 5.9 的曲面图清楚显示聚集群与聚集群之间区域的密度差异。最初点很多,随时间邻近聚集群合并而减少,最终所有点都消失。
这与细菌实验中观察到的现象完全一致。最初在含均匀琥珀酸浓度的 petri 皿中加入细菌,将混合物充分搅拌后静置。此时皿内溶液状态由模型初始条件模拟,即均匀细胞与琥珀酸分布的小扰动。约 20 分钟后活细菌聚集成许多明显分离的小团——对应模型解中观察到的随机排列的点状聚集群,被近零细胞密度区域分隔开。实验中可见细菌聚集群彼此合并形成较少但更大的团。数学模型也是如此,将解作为帧间隔很小的电影显示时尤为清楚。模型与实验中点最终消失,且不能再诱发重新形成。这再次由数学模型解释为化学趋性响应饱和——化学吸引剂的产生持续增加,化学趋性不再起作用。
非均匀刺激物分布。迄今为止考虑的是刺激物 w 空间均匀分布下出现的细胞密度图样;上述分析在定性上捕捉了此情形下实验观察的行为。但实验中当刺激物以局部滴加方式加入培养基时还观察到其他图样。若模型本质正确,也应再现这些图样。本节作简要考察。
对非均匀刺激物分布,需在模型中加入刺激物的扩散。因此模型 (5.20)(5.21) 必须扩展为含刺激物方程,即 (5.14)–(5.16) 的无量纲形式。得到更一般的模型:
其中无量纲 d_s = D_s/D_c(参见 (5.19))。该形式要求 u(x, t) 和 w(x, t) 的均值在所有时刻都为 1。
图 5.10 显示该模型一次模拟的数值结果。如图 5.1 实验观察到的,刺激物加入处周围发育出高细胞密度的环。环内也形成一些聚集群(图 5.1 实验中也观察到),但远不如环那么密集。原因是环从其周界外招募细胞,因此环可用的细胞数远多于中心聚集群。
模拟显示随时间环的半径减小。最终半径小到环本质上只是一个点。实验上是否观察到该行为会很有意思。
在讨论更复杂的半固体实验图样建模前,简要回顾以上几节内容:相对简单但具启发性的分析解释了化学趋性模型中如何出现随机或环状排列点状的瞬态演化图样,匹配实验安排。核心想法是考察小时间区间内各模式的增长率,并外推到所有扰动频率的联合行为。低模式数扰动不稳定且振幅增大,但最终稳定且衰减,较高模式数先稳定。这不仅在定性上,而且在相当程度上定量地与实验与数值观察一致:聚集团形成、合并成更大聚集群、最终消失。
5.5 S. typhimurium 半固体相模型机制(Semi-Solid Phase Model Mechanism for S. typhimurium)
5.1 节讨论过半固体实验,Budsrene 和 Berg(1991)观察到两种截然不同的图样形成机制。S. typhimurium 情形下:薄薄的菌苔从接种点向外铺展,更浓的细菌环在菌苔边缘之后形成。每个环最终可能断裂为点,但通常在更大半径处形成更多环之后才发生(图 5.2)。第二种图样形成机制由 E. coli 展示:首先形成扩张的高密度细菌环(swarm ring);随 swarm ring 扩张,留下较小的细菌聚集群,构成图 5.3 所示的显著图样。
S. typhimurium 图样中菌苔在图样之前的事实提示:先建立了一个空间与时间均匀的稳态,然后图样在其上形成。如前所示,半固体实验模型没有非零稳态。但若营养消耗足够慢,可忽略之,模型简化为允许所需稳态的两方程。因此假设在第一种图样形成机制中营养消耗可忽略。该假设得到实验进一步支持:实验中食物供应量远超细胞饱和水平。第二种图样机制(E. coli)则涉及较强的营养消耗速率(Budsrene 和 Berg 1995);swarm ring 扩张后留下的区域内食物所剩无几,这种情形下半固体模型的所有三个方程都重要。
对 S. typhimurium 情形,半固体实验的模型系统 (5.11)–(5.13) 通过下述设置无量纲化:
其中 t^* 与空间标度都通过乘积 k_7 n_0 进入,扩散比 d_u、d_w 把细菌与营养的扩散系数相对于化学吸引剂标度化;化学趋性参数 α 与 5.3 节有相同形式,但增长率 ρ、承载能力 δ、产生系数 β、消耗系数 κ 都分别除以 k_7 或乘以 √k_9 等。得到无量纲模型(参数估计见表 5.3,为代数方便省去星号):
细胞方程中增殖项采用饱和形式 δ w²/(1 + w²) − u(营养供应充足时趋向 δ),化学吸引剂方程中产生项饱和于 u²/(μ + u²),消耗为线性项 u v;刺激物方程中消耗采用与增殖相同形式但带系数 κ。注意三个方程分别对应细菌、化学吸引剂、刺激物的守恒律,且化学吸引剂的扩散时间尺度与刺激物的不同。
回想 S. typhimurium 实验中图样形成过程分两步:第一步形成薄的圆盘状菌苔,提示暂时存在空间与时间均匀的稳态;第二步在菌苔前缘之后形成高密度环状聚集群,提示环图样形成于中间稳态之上。由于图样在菌苔建立之后很长时间仍在形成,菌苔内的营养消耗必须可忽略。因此可通过忽略食物消耗动力学来近似完整系统。同时由于营养浓度大但细胞浓度小,可研究简化的两方程模型:
此处 w 在两方程中作为参数进入(不再有时间演化),细胞增殖项中营养供应简化为 δ w/(1 + w)。方程 (5.41)(5.42) 的所有无量纲参数(表 5.3)均可由表 5.1 的维度值计算:d_u = 0.2–0.5(细菌相对化学吸引剂的扩散比)、d_w = 0.8–1.0(营养相对化学吸引剂的扩散比)、α = 87(化学趋性参数)、δ = 3.5(承载能力)。表中其余参数尚无实验估计,但因分析按无量纲形式进行,只需对参数的特定组合有估计即可。
对方程的分析比 5.3 节简单得多,因为存在 u 和 v 都非零的均匀稳态。可以用通常的线性分析确定该稳态是否不稳定以及空间图样是否可能形成。下文还讨论了完整的非线性分析部分结果,因其高度相关于具体形成的图样。
5.6 基本半固体模型的线性分析(Linear Analysis of the Basic Semi-Solid Model)
线性分析与第 2 章详细讨论过的方法相同,现在直接明了。在非零稳态 (u, v) 附近线性化 (5.41) 和 (5.42):
下述分析中为代数简单起见,对模型方程 (5.41)(5.42) 中的项使用一般形式:w 在此实际上是另一个参数。因此考虑
对比 (5.41)(5.42) 定义
用通常方法在稳态附近线性化,设
其中 0 < ε ≪ 1。代入 (5.44)(5.45) 得到线性化方程
由于 f = 0 且 g = 0(* 表示在稳态处求值)。将线性系统写成向量形式
其中矩阵 A 和 D 定义为
用通常方法寻找如下形式的解
其中 k 为波矢,c 为常数,色散关系 λ(k) 即增长率待确定。代入矩阵方程得
当且仅当系数矩阵行列式为零时有非平凡解:|λI + D|k|² − A| = 0(回忆第 2 章)。色散关系 λ(k²)(k² = |k|²)由特征方程给出
解记为 λ_+ 和 λ_−。
关注至少有一个正增长率的图样模式——即至少一个解有 Re λ > 0——因此聚焦于两个解中较大的 λ_+:
其中
回忆第 2 章的讨论:若 Re λ_+ 为正(负),则稳态附近的扰动将增长(衰减)。寻找对纯时间扰动(即 k² = 0)稳定(λ_+ < 0)但对至少一个空间模式不稳定(λ_+ > 0)的解,这就像纯扩散驱动失稳情形(针对某些非零 k₁² 和 k₂²)。数学上寻找满足以下条件的参数范围:
为满足 (5.55) 第一条必须有
注意这些条件意味着 λ_− 始终为负,因此只需聚焦 λ_+。作图而言,这些条件给出 λ(k²) 的反向抛物曲线,最大值在 k² = 0 的右侧;这就是扩散-化学趋性驱动失稳的最基本色散关系,第 2 章已详细讨论。图 5.11 给出典型色散关系,示例性地展示了它们随参数 μ 的变化:μ 足够小时不存在 k₁²、k₂² 范围。
关注曲线与 λ = 0 的交点,给出 k₁² 和 k₂²。这些对应参数空间的边界,由 (5.54) 通过下式两个 k² 解给出
一般而言,对每组参数值有两个、一个或零个 k² 满足 (5.57)。在分岔处(一个值),有
和
其中 k_c 是临界波数。对 α 解最后一个方程并代入 (5.58) 得临界值
给出模型系统中扩散-化学趋性驱动空间失稳的参数空间,通过函数 χ、f、g(由 (5.46) 定义)及其在稳态处的导数表示,依赖于其他模型参数。这里选取了二次方程的正根。(5.60) 定义参数空间中分隔正负 λ 区域的一个临界边界集。(5.59) 是 λ_+ = 0 的分岔曲线,在此波图样既不增长也不衰减。图 5.12 给出对应不同 μ 值的分岔曲线序列。每条曲线上方 λ_+ > 0,所有不稳定图样模式都增长,即对 (5.57) 所得 0 < k₁² < k² < k₂² 的所有 k² 成立。此刻尚不知可期什么图样,只知道空间图样可能。进一步确定哪些图样将出现是非线性问题。至今唯一解析确定这些的方法是所谓弱非线性分析——即在解从空间均匀分岔为空间异质附近的分析。第 2 章给出了部分反应扩散方程这样做的参考文献。该方法同样适用于这类方程,只是稍微复杂。Zhu 和 Murray(1995)作了非线性分析,并解析地评估了反应扩散和扩散-化学趋性系统的参数空间,比较它们生成空间图样的潜力。如前所述,Tyson(1996)用该方法分析本章讨论的具体模型方程。
这里不进行(相当复杂的)非线性分析(完整讨论参见 Zhu 和 Murray 1995 与 Tyson 1996)。仅勾勒程序并给出结果。首先由线性分析确定在哪些参数区域图样可能。非线性分析实际上给出将形成的图样类型。分析基于假设参数接近参数空间中的分岔曲线。以某个参数(如感兴趣的化学趋性参数 α)作关于其分岔值 α_c 接近、但其值使系统进入图样形成空间的渐近分析。从边界值问题的线性解开始,例如一维情形含 e^{λt} cos kx 的解,其中 k 在不稳定波数范围内。线性情形下该解以速率 e^{λ(k)t} 随时间指数增长。直观上考察系统方程 (5.38)–(5.40) 中未微分项可知,解不可能无限增长。在弱非线性分析中,先考虑线性解为 k² = k²_c 时的线性边界值问题解——即参数接近从无图样到图样的分岔曲线。然后用渐近扰动方法以解的振幅作为慢变时间函数的形式研究解。确定振幅在 t → ∞ 时有界的条件。该过程确定各种可能解中哪些将演化为稳定解。下文在讨论线性边界值问题的可能线性解形式时给出更多细节。
线性边界值问题。非线性分析的必要先决条件是求解相关的线性边界值问题。可能的图样类型依赖于边界条件允许的不同波矢 k 的数量。本章关注实验域由重复图样铺满。数值模拟中为简单选方形域,因此感兴趣的是方形或矩形瓦片产生的条纹与点。每个方形或矩形单元满足周期边界条件。一般而言,对规则铺砌,可以有方形、条纹、六边形等,正如第 2 章所讨论。
考虑矩形域 S,定义为 0 ≤ x ≤ l_x,0 ≤ y ≤ l_y,边记为 S₁: x = 0, 0 ≤ y ≤ l_y;S₂: x = l_x, 0 ≤ y ≤ l_y;S₃: y = 0, 0 ≤ x ≤ l_x;S₄: y = l_y, 0 ≤ x ≤ l_x。带周期边界条件的空间本征值问题为
偏微分方程 (5.61) 的可能本征解为
其中 k²_n 为允许的本征矢,下文讨论。将 (5.61) 的边界条件代入解 (5.62) 可得
其中 k^T_n = (k_x, k_y)。(5.63) 的一般解为
其中 m, n, p, q 均为整数。关注参数接近分岔曲线的情况,因此尤其感兴趣满足 |k_n|² = k²_c 的解 k_n。(5.63) 允许的这种解向量数依赖于 l_x 与 l_y 的关系。假设 k²_c = (2Mπ/l)² 且 S 是方形 l_x = l_y = l。若 M = 1,则有两种可能解向量:
若 M = 5,则有四种可能解向量:
这些解向量在非线性分析中很重要。
5.7 非线性分析纲要与结果(Brief Outline and Results of the Nonlinear Analysis)
通过线性分析只能确定当均匀稳态被驱动至空间异质扰动不稳定时 u 和 v 在该稳态附近的小振幅初始行为。这些空间不均匀解初期指数增长,显然对所有时间都无效。对本章这一类问题可以作非线性渐近分析,得到 O(ε)(原则上还可到更高阶,但代数过于繁重)的对所有时间有效的解。如前所述,过程细节由 Zhu 和 Murray(1995)针对包括扩散-化学趋性的多种图样形成机制给出。对更复杂的化学趋性机制,由 Tyson(1996)完成。本节勾勒解析过程——即确定系统方程 (5.41)(5.42)(一般形式为 (5.44)(5.45))对所有时间有效的小扰动解——的多尺度渐近分析。首先写出
其中 (u, v) 是空间均匀稳态,依赖于模型参数,由 5.6 节可知参数越过其分岔值时该稳态被驱动不稳定,导致空间不稳定解。通过下述方式重新标度时间
其中 ω_i (i = 1, 2, …) 需确定。
考虑方程 (5.44)。将展开式 (5.66)(5.67) 代入各项得
前两个表达式含 û 的线性项,后两个表达式含 û 和 v̂ 的线性、二次、三次及更高阶项;f* = 0 由稳态定义。所以 (5.44) 化为另一个含 û、v̂ 的线性、二次和更高阶项的方程。(5.45) 类似地变换。一般而言,(5.44)(5.45) 与 (5.67)(5.68) 取如下形式
其中 * 表示在稳态 (u, v) 处求值。Aû、Q(û)、C(û) 分别表示化学趋性函数和反应函数在稳态附近展开的线性、二次和三次项。矩阵 A 和 D 已由前述线性分析确定。关注模型参数取特定值使 λ = 0 的情形。由 (5.60) 定义该参数集,称为临界参数集。即参数集位于增长的空间异质解与空间均匀解之间的边界上。
现在扰动模型参数中的某一个 a(如 (5.44)(5.45) 中任意一个参数),关于其在临界集中的值,则时间增长的特征值为
因 λ(a_c) = 0 由 a_c 定义。扰动取为使
因此扰动对解 e^{λt + i k·x} 的影响限制为时间增长率 λ 的变化。若 Re(λ(a)) 因此变为正,按线性理论对应的 k²_c 图样模式预测会增长。依赖于参数,结果可以是稳定或不稳定的空间异质解。若增长足够慢,可以预测它是否会发展成时间稳定图样,并进一步预测图样特征(如点状或条纹)。先在临界集附近扰动稳态模型 (5.69)。为简化分析,仅扰动一个参数,并称其为 a。Tyson(1996)用模型方程的实际参数作了分析,下面给出她的结果。
考虑如下展开
代入 (5.69) 得系统
其中上标 c 表示在临界集处求值。临界参数 a 的变化只发生在 â 中,因此其效应可在分析中分离。代入所有小变量 ^() 的展开,收集并等同 ε 的同类项,对每个 ε 阶得方程组。为符号简洁,省略上标 c,分析的其余部分所有参数值均取自临界集。为展示这些方程的形式,给出 O(ε) 和 O(ε²) 系统;为完成非线性分析,还需要 O(ε³) 系统(代数极其复杂)。无需它们来勾勒过程。O(ε) 方程为
是线性的,并定义系数矩阵为线性算符 L。O(ε²) 阶方程为
需要 O(ε³) 方程的分析需要这些线性方程组的解。Tyson(1996)完成了解答(Zhu 和 Murray 1995 也作了等价分析)。代数庞大但必要,以得到对所有时间一致的 O(ε) 解并确定哪些具体图样将稳定。
为给出结果并解释非线性分析所得,需要 O(ε) 系统 (5.72) 的解。寻找如下形式的解
其中
是解的正弦部分。将该形式与 (5.62) 联系,ā_l(T) + a_l(T) ∝ B。代入 (5.74) 到 (5.72) 得 V_{1l} 的表达式直到任意常数倍。常选取使向量大小为 1,故
由于仅考虑临界参数集,故 |k_l|² = k²_c 对所有 l 成立,V_{1l} = V_{1} 对所有 l。此时并不知道 a_l(T) 和 ā_l(T),即 O(ε) 解的复振幅(是 (5.67) 定义的慢变时间 T 的函数);解的振幅为 |a_l(T)|。非线性分析的关键是确定振幅。所以需要解 O(ε²) 方程。由于这些线性方程(具有相同算符 L)右端含未微分项,可能有包含 x sin x 类的永真项(secular terms),对大 x 变得无界。
简单方程 u'' + u = −εu'(其中撇表示关于 x 的导数,0 < ε ≪ 1,且要求 u(0) = 1, u'(0) = 0)容易看出永真项如何出现。令 u = u_0 + ε u_1 + …,假设所有 u_i 都是 O(1)。代入并按 ε 收集同类项,得 u_0(x) = cos x,u_1 方程及边界条件为
解为 u_1(x) = (1/2)(sin x − x cos x)。所以 u_0 + ε u_1 + … 不是一致有效解,因为 u_1(x) 因 x cos x 项对所有 x 不是 O(1):这就是永真项。Murray(1984)关于渐近分析的教材对此类方程获得一致有效解的渐近过程作了详尽教学性描述。
回到 O(ε²) 方程讨论,结果表明这些方程不产生永真项,所以振幅函数 a_l(T) 和 ā_l(T) 在该阶仍未确定。但在 O(ε³) 处永真项出现。此时确定振幅方程:即使不实际求解该阶,也选择这些方程使永真项不出现在 O(ε³) 解中。得到这些方程(称为 Landau 方程)的代数极复杂。这些方程关键涉及 (5.74) 中的数 N——解中具有 |k_l|² = k²_c 的模式数。上节关于边界值问题的讨论中已看到解与本征矢如何随该数变化。
举例如,假设 N = 2;Tyson(1996)给出振幅(或 Landau)方程
其中 X_A、X_B 和 Y 是原系统 (5.44)(5.45)(即 (5.41)(5.42))参数的复杂函数。|a_1| 和 |a_2| 直接对应于 (5.62) 的 A 和 B,只是此处它们是时间的函数。稳定空间异质解是否存在取决于这些振幅方程在 t → ∞ 时的解。它们只是常系数常微分方程。它们有下述可能的稳态解,其存在性依赖于系数的符号:
(1) |a_1|² = 0, |a_2|² = 0 (2) |a_1|² = 0, |a_2|² = −Y/X_A (3) |a_1|² = −Y/X_A, |a_2|² = 0 (4) |a_1|² = −Y/(X_A + X_B), |a_2|² = −Y/(X_A + X_B) (5.79)
第一个稳态对应零振幅图样,即没有图样。第二、第三个对应一个方向零振幅、另一个方向非零振幅,产生条纹。第四个稳态在每个方向都有非零振幅,因此产生点状图样。若这些稳态都不稳定,分析不能确定形成的图样类型。这在参数图中被称为未确定区域,下文展示。这些稳态的稳定性条件汇总于表 5.4。通过在给定参数空间上计算 X_A、X_B 和 Y,可以用表 5.4 将该空间划分为发生点、条纹、无图样或未确定图样的区域。Tyson(1996)对本章研究的系统、Zhu 和 Murray(1995)对更简单的反应扩散-化学趋性系统计算了这些参数空间。
5.8 模拟结果、参数空间与基本图样(Simulation Results, Parameter Spaces and Basic Patterns)
本节给出两个模型——方程 (5.41)(5.42) 及其适度简化版本——的模拟。模拟所用初始条件严格适用于前述分析。即模拟从空间和时间均匀稳态解的小(ε 阶)扰动开始,并对其中一个参数作更小(ε² 阶)的扰动。回想边界条件为周期的。
简化模型。先考虑下述 (5.41)(5.42) 的简化模型:
该模型只有五个参数,因此参数空间比 (5.41)(5.42) 稍易探索。若固定 d_u 和 ρ,则可变化 δ 和 β,并由 (5.61) 第一式确定 α。δ、β 平面上的每一点将对应一个临界参数集,用于确定 X_A、X_B 和 Y。该方法产生图 5.13 所示图样区域。图中常 α 曲线为双曲线
对方程 (5.14)–(5.16) 取两套参数(一来自上条纹区域、一来自点区域)作数值解。模拟在方形域上执行,刚好能容纳一个完整图样周期。这要求选择 l_x = l_y = 2π / √(k²_c)。所以预期一套条纹参数产生一条条纹、一套点参数产生一个点。结果确实如此,图 5.14 显示细胞密度的典型结果。每次模拟有两幅图显示初始条件和稳态图样。
注意由于边界条件为周期,点状或条纹图样的最大值可出现在域的任何位置。另外对条纹图样,方形域上分析不决定条纹方向。所以同一组随机初始条件可能产生垂直条纹,另一组则产生水平条纹。具体的参数取值为:点参数 d_u = 0.25, α = 1.50, β = 0.1, δ = 15.0, ρ = 0.01, w = 1.0, u₀ = u = 15.0, v₀ = v = 1.50, k²_c = 3.0, l_x = l_y = l = 3.6276(能容纳一个 2π 振荡的域);条纹参数 d_u = 0.25, α = 2.25, β = 0.2, δ = 20.0, ρ = 0.01, w = 1.0, u₀ = u = 20.0, v₀ = v = 4.0, k²_c = 4.0, l_x = l_y = l = 3.1416。扰动大小为 ε = 0.1,时间为无量纲单位。白色表示高细胞密度,黑色反之。图 5.14(a)(b) 使用同一随机种子得到条纹(垂直方向);(c) 用不同种子得到水平条纹;振幅和大小相同。
完整模型。现在考虑更具生物学准确性的模型
含七个参数(含 w),对 α、d_u 和 δ 有实验估计。剩余四个参数之一可由分岔条件 (5.58) 确定。在简化模型中对 α 求临界值并固定其他参数;现在 α 已知,希望固定其值并解一个未知参数。最简单的解是 ρ:
剩余参数为 μ、β 和 w。需要探索 (β, μ)、(μ, w) 和 (β, w) 参数空间。系统为四维参数空间,但对给定 ρ 是三维。此时可通过设置
将模型重写为
这简化了代数,但仍可通过 (5.81) 用 w 增大时把 (β, δ) 点映射到 (β̄, δ̄) 的方式确定增减 w 的效应。该映射对 (β̄, δ̄) 平面上的每个 (β, δ) 对为 sigmoid 曲线。
用该形式未知参数为 β、δ 和 μ。若固定 μ 比 u² 大得多,则恢复刚讨论的简化模型 (5.80)。令人惊讶的是,β̄ 不必比 β 大很多就能恢复该图样域。对较小 μ,图样域不同。图 5.15 给出一个例子,参数 d_u = 0.25, α = 90.0, μ = 1000.0;随 β̄ 和 δ̄ 变化,由 (5.82) 确定对应的 ρ。扰动参数为 β。图样区域给出条纹(I)、点(II)和未确定图样(III)。该结果与简化模型相比,条纹-点-未确定区域的边界明显依赖于 μ 的大小。
Tyson(1996)对条纹、点与未确定区域的参数作了模拟,确认了解析预测——哪种图样类型将在小周期域上的均匀稳态的小随机扰动下出现。具体的点参数取 d_u = 0.25, α = 90.0, β̄ = 10.0, μ = 100.0, δ̄ = 5.1, ρ = 7.797, u₀ = u = 5.1, v₀ = v = 0.4747, l = 1.177;条纹参数 δ̄ = 4.6, ρ = 8.5133, u₀ = u = 4.6, v₀ = v = 0.3797, l = 1.213;未确定图样参数 δ̄ = 5.5, ρ = 7.139, u₀ = u = 5.5, v₀ = v = 0.4233, l = 1.159。扰动为 ε = 0.1。有趣的是她还发现在至少部分未确定区域存在稳态图样,而该处的分析目前尚不能预测将出现什么图样。这后者的解释是个有趣的解析问题。图 5.16 显示这些图样的例子。
5.9 实验初始条件下的数值结果(Numerical Results with Initial Conditions from the Experiments)
5.8 节讨论的半固体模型的非线性分析只适用于均匀正稳态的小随机扰动。但实验的初始条件(如 5.1 节所述)完全不同。最初没有化学吸引剂,只有小的局部接种——是均匀零稳态的较大扰动。在大域上也是零通量边界条件,而非小域上的周期边界条件。在图样之前的菌苔建立了与非线性行为相关的条件。因此用实验初始条件预期在图样域的条纹区域得到环、在点区域得到断裂(点状)环。Tyson(1996)在条纹与点区域各参数下作了一维和二维广泛模拟。这里给出她的一些结果及 Tyson 等(1999)给出的结果。
先用一维简化半固体模型 (5.80)。用与图 5.14 相同的参数值(曾产生条纹图样),再次得到一系列同心环。用曾产生点状图样的参数值,仅在初始扰动附近出现少量小脉冲,且缓慢衰减。似乎实验上得到的点状环图样要么来自非线性分析图样域的条纹区域,要么来自完全在非线性分析预测之外的某个图样区域。无论如何都不能从非线性分析得到完整的图样化场景。
即使是参数空间的粗略调查也显示不必把自己限制在非线性分析参数下寻找有趣图样及其随参数变化的规律。例如增大化学趋性系数 α,振幅增大且条纹波长也增大。稍微增大营养浓度 w,传播速度增大并导致波后脉冲衰减。稍微更大程度增大 w,脉冲消失,细胞密度快速逼近均匀正稳态。w 较低时仍形成脉冲,但图样传播较慢。增长率 ρ 对图样的影响与营养浓度完全相同。化学吸引剂产生 β 减小两倍使脉冲振幅和频率都增大。承载能力 δ 直接与图样振幅以及形成图样所需时间长度相关。排列组合无穷无尽。
二维情形下,用非线性分析给出条纹的参数作模拟,模拟给出同心环。径向图样的波长二维比一维小,振幅大致相同。环间距不随半径增大而变化(图 5.17)。
二维情形下把化学趋性参数从 α = 2.25 增大到 α = 5,得到由点组成的同心环序列。图 5.18 与 t = 70 时解的曲面图一起给出图样的时间序列。与一维情形相同,增大化学趋性系数使图样波长与振幅都增大。重要的是,营养消耗既不把点图样变为条纹图样,也不改变波长。
现在考虑完整模型系统 (5.84),它多引入一个参数 μ(衡量化学吸引剂产生的饱和水平)。简单模型中发现变化参数 α、β、ρ、δ 和 κ 可以改变同心环图样。基本上完整模型给出与简单模型相同的图样,但某些定量差异——如预期。平凡地,对足够大的化学吸引剂产生的饱和水平 μ 和可忽略的营养消耗,完整模型基本简化为简单模型并完全相同地表现。更重要的是,对允许天冬氨酸产生曲线中间及饱和部分参与的参数,也得到连续环和点状环(图 5.19(a)(b))。加入食物消耗不改变环图样性质,只是让它逐渐从中心向外消失。
模拟与实验的关系。模拟结果与 S. typhimurium 实验结果(5.1 节)吻合得极好。图样之前是低细胞密度的菌苔。每个环在前一个环之外某个离散径向距离形成,然后保持静止。点状环最初形成连续环、随后断裂成点。模型机制所得的这些特征都是 S. typhimurium 图样的典型特征。
若取图 5.19(b) 的参数值,得
在距离 x* = 10 上形成四个环,对应 x = 1.4 cm,接近实验观察值 x ≈ 1 cm(Woodward 等 1995)。
5.10 半固体相模型的 Swarm Ring 图样(Swarm Ring Patterns with the Semi-Solid Phase Model Mechanism)
Budsrene 和 Berg(1991)观察到的最显著的图样来自一束扩张的高密度细菌环——swarm ring。这些图样已在 5.1 节描述并图示。初始条件是 petri 皿上一处局部接种的细胞,皿中含有均匀分布的食物、无化学吸引剂,边界条件为零通量。
由于图样起始于连续的 swarm ring,仅在之后留下带角向变化的图样,自然先研究一维情形的行波或脉冲。初始构型为均匀食物浓度和一处致密的细菌接种。没有化学吸引剂。
最简单的 swarm ring 可在没有化学吸引剂产生的情况下形成。随时间增加,细菌消耗食物并向外扩散。留在中间的细菌变为不可动。在扩散细胞团的外缘,细胞密度低且食物浓度高,结果细胞增殖,提高局部细胞密度。同时在接种位置中心,食物已被消耗、浓度降到细胞死亡主导的程度。结果是接种位置细胞数减少、扩散前沿细胞数增加。这种情形演化产生行波脉冲。
加入化学吸引剂可使光滑 swarm ring 不稳定。可合理假设该不稳定能在二维(甚至更高维)中引发更复杂的几何,特别可产生 swarm ring 后实验观察到的点。
对该现象作了解析上的详细研究。从寻找简化模型的一维行波脉冲解开始,寻找以行波坐标 z = x − c t 表示的解,其中 c 为脉冲传播速度需确定。半固体模型方程 (5.38)–(5.40) 写作 u(x, t) = U(z), v(x, t) = V(z), w(x, t) = W(z),行波波形变为:
其中撇表示关于 z 的导数。(5.85) 可写成关于 U, U', V, V', W, W' 的 6 阶一阶系统,稳态为
和
第一个是脉冲后方存在的稳态,第二个是脉冲前方存在的稳态。实际解 U、V、W 必须非负有界,所以当接近两个稳态时特征值必须实,确保无振荡。在 (5.86) 和 (5.87) 附近线性化 (5.85) 并求解特征值,发现第一个稳态始终是焦点。要让第二个也是焦点,须有
对 Fisher-Kolmogorov 方程,卷 I 第 13 章详细研究过,在合适初始条件(即紧支撑)下稳定行波解以速度 c_min 演化。实验实际隐含这样的初始条件,所以假设 (5.85) 的行波脉冲解也将以 c_min 或接近 c_min 的速度传播。有趣的是,营养消耗速率 κ 和化学趋性系数 α 对 c_min 都无影响;按线性分析纯粹是扩散与增殖细菌的动力学决定脉冲速度。
Swarm Ring 的稳定性。一旦 swarm ring 形成并开始在 petri 皿上扩张,实验上观察到环周期性地断裂,留下点状图样。这提示数学上应寻找对径向方向局部稳定、对角向方向局部不稳定的 swarm ring 解。
以下概述如何可能获得 swarm ring 稳定性的某些信息。为使解析建议模仿实验安排,考虑矩形域(为解析方便),初始接种沿长度 l 的一条边放置。垂直于该边,预期解为行波脉冲。将二维模型在 (x, y, t) 坐标下用一个行波坐标 z、平行于波前的坐标 y 重写,方程 (5.38)–(5.40) 变为:
现在假设 z 方向有行波解 U(z, y) = U(z), V(z, y) = V(z), W(z, y) = W(z) 对所有 y ∈ [0, l]。则按定义
若该解对 z 方向扰动稳定(如对某些参数范围很可能是这样),则可通过只考虑涉及空间坐标 y 的时空扰动得到垂直于波传播方向扰动效应的某些概念。作为这种稳定性分析的首次尝试,可考虑如下形式的解 u、v、w:
其中 ū、v̄、w̄ 为小常数。代入 (5.89) 并收集 O(ū)、O(v̄)、O(w̄) 项得 O(1) 线性方程组:
其中矩阵 A 为
其中
按通常方式,扰动 ū、v̄、w̄ 的非零解当且仅当 |A| = 0 成立,由此得到色散关系 λ 的特征方程。令 |A| = 0 得到关于 λ 的三次方程:
其中 A、B、C 是参数、波数 k 以及通过行波解 U、V、W 及其导数作为 z 的函数的函数。若暂时假设 A、B、C 为常数,保证 Re(λ) < 0 的 Routh-Hurwitz 条件(参见卷 I 附录 A)为:
要确保不稳定(即 Re(λ) > 0 对某个 k² ≠ 0),至少一个条件须被违反。将 (5.93) 的系数矩阵 A 重写为
其中始终为负的项前置负号,明显为正的项前置正号。该记法下 (5.95) 有:
考虑稳定性的三个必要条件。对 A,唯一可能为负的项是 c_1,且仅当 (5.94) 中 H_1 足够大时才发生。但这依赖于 z,因此依赖于在行波前缘上的位置。对其他稳定性条件类似论证。显然行波前缘(swarm ring)是否对横向扰动不稳定依赖于行波变量 z。对至少某些 z 若不能违反条件 (5.96),从而对非零波数有 Re(λ) > 0 并引起横向空间不稳定(即 swarm ring 断裂为点),将会令人惊讶。由于化学趋性参数 α 出现在 H_1 和 H_2 中,化学趋性再次扮演关键角色。这些 swarm ring 稳定性的完整分析是一个有挑战性的未解决问题。
Tyson(1996)求解了一维方程,对各种参数值从完整模型获得单一的行波脉冲或两到四个行波脉冲组成的脉冲串,且营养消耗不可忽略。模型参数 α、β、μ、δ、ρ 和 w 对行波图样的影响类似于它们对前述静止图样的影响。脉冲串中脉冲的波长主要受化学趋性系数 α 影响。
预测与数值计算的波速吻合得很好。计算波速始终比预测值 c_min 大(约 5–10%),如预期那样。由于有限域上必须施加边界条件,所以模拟中的波速不会如 c_min 那样小。
二维 Swarm Ring 的数值结果。二维情形下发现行波脉冲串变成一个 swarm ring 后跟着一两圈点。这些点来自内圈——内圈发展出角向不稳定性并随后断裂为点。发展良好的 swarm ring 图样的图像图示于图 5.20,展示了其发展中的多个时刻。
该模型发展和分析是第一个针对 E. coli 和 S. typhimurium 的数学模型,无须假设任何外部生物学活动即可生成 swarm ring 产生点。Woodward 等(1995)通过假设存在化学吸引营养物获得了有趣图样。S. typhimurium 的实验确实在这种情形下进行过,但 E. coli 实验中不存在这样的营养物。本模型显示朝向细胞产生的天冬氨酸的化学趋性加上食物消耗足以产生实验观察到的行为。这与实验者的直觉一致。
5.11 Bacillus subtilis 的分支图样(Branching Patterns in Bacillus subtilis)
前述图样虽然复杂,但属于点和环形成的相当规则的图样。当 Bacillus subtilis 接种到营养很少的琼脂培养基上时,可展现出相当不同的分形状图样,与扩散限制聚集(diffusion-limited aggregation)中看到的类似(参见 Matsuyama 和 Matsushita 1993)。但琼脂半固体时,B. subtilis 形成的菌落是密集分支图样,外包光滑包络。培养基的硬度影响形成的图样。Shigesada 及其同事(Kawasaki 等 1997)研究了这种特定细菌并构造了相对简单的反应扩散模型,捕捉了实验发现的许多图样特征:他们将结果与实验比较。本节简要描述他们的模型并展示部分结果。虽然他们的模型是反应扩散型的,但它原始且与本章此前研究的反应扩散系统有本质不同。它再次凸显了这类相对简单系统图样形成的丰富性。
他们提出的模型由针对细菌细胞和营养的守恒方程组成:
其中 n 和 b 分别为营养浓度和细菌细胞密度。函数 k n b / (1 + γ n)(k、γ 为常数)是细菌对营养的消耗率;θ (k n b / (1 + γ n)) 是细胞的增长率,θ 为转化率因子。D_n 和 D_b 分别为营养与细胞的扩散系数。现在论证 D_b 形式背后的原因。D_b = σ n b 的形式基于 Ohgiwara 等(1992)的观察,他们详细观测了细菌运动,发现细胞在菌落内部(营养水平低)不怎么运动,但在菌落外围(营养水平高)运动剧烈。他们还指出在最外层前沿(细胞密度相当低)细胞再次相对不活跃。Kawasaki 等(1997)由此论证:细菌在营养 n 或细菌密度 b 小时不活跃。他们通过将细菌扩散取为正比于 n b 来建模这些效应,比例因子为 σ。他们还观察到虽然每个细胞以典型随机方式运动,但其中一些表现出随机涨落。他们通过设置 σ = 1 + ξ 量化了这一点,其中参数 ξ 衡量偏离通常随机扩散的随机涨落。
Kawasaki 等(1997)研究了这些模型方程在二维情形下受初始条件
的图样形成潜力,其中 n_0 为初始均匀营养浓度,b_0(x) 为初始细菌接种。由于实验中营养浓度相对较低,γ n 项体现的饱和效应可忽略,所以营养消耗可近似取为 k n b,下文即用此函数形式。
通过下述设置对模型作无量纲化:
应用上述 γ = 0 的消耗近似并省略星号作代数简化,模型机制变为:
含唯一参数 σ,初始条件为
Kawasaki 等(1997)在多种不同情形下求解该系统,发现解展现出惊人范围的复杂图样。图 5.21 给出一些例子。
Kawasaki 等(1997)还研究了当随机参数 ξ = 0 时形成的图样。图样仍形成并给出分支状图样,但由于没有各向异性,图样规则且对称得多。当然实际中营养与细菌密度中的小随机扰动不会产生如此规则的图样,正如 Kawasaki 等(1997)所指出的。
虽然图样在二维中演化,每个尖端本质上在一维中生长——除非发生分支。这使得可以用模型方程的一维版本获得一些关于尖端生长的近似解析结果,因此也能获得菌落生长的结果:
这些一维方程的数值模拟给出与二维方程所得吻合良好的增长率。他们通过将 σ₀ ν₀ b 和 ν₀ b (1 − b/K) 替换 (5.108) 中的 σ n b 和 n b 进一步近似模型,得到关于 b 的标量方程:
此处细菌增长受营养限制,按典型 logistic 增长,K 为细菌的饱和水平。若在无扩散情形下考虑 (5.107)(5.108),可对方程相加并积分,且 n + b = ν_0,得细菌的 logistic 形式 db/dt = ν₀ b (1 − b/ν₀),由此将 K 与 ν₀ 关联。(5.109) 的形式与卷 I 第 13 章 4 节详细讨论的方程相同,该方程有精确行波解,波速(菌落生长速度)为:
该速度是从 (5.107)(5.108) 获得的菌落生长速度的很好近似(稍高估)。
营养水平不低时须使用完整模型 (5.99)(5.100) 并令 γ ≠ 0。Kawasaki 等(1997)考虑了这种情形,发现 γ 增大时分支宽度增大、复杂度降低。他们包含随机性的方法有趣且重要,因为它允许一些随机性而不带有此类研究中通常涉及的复杂性。这一概念显然有更广泛的应用,例如本书中研究的许多模型。
本章个人批注
本章给我最深的印象是作者把一个具体的实验系统(Berg 实验室的细菌图样)当作"机制验证场",并以此示范如何把抽象的反应扩散-化学趋性模型一步步落地到真实生物学。本章与第 4 章(生长域)形成强烈对比——第 4 章是从生物学出发展开建模,本章则是先有数学模型再去验证生物学。两章风格互补,Murray 在第 5 章的开篇就坦白说:"Howard Berg initially got in touch with me",整个建模工作的起点不是数学上的优雅,而是实验现象得不到直觉解释。这是本章叙事的一条隐线,我作为读者更喜欢这种带故事感的写法。
5.2 节中关于扩散系数的讨论令人想起第 1 卷第 13 章 Fisher-Kolmogorov 等的处理方式——把"扩散系数"作为可测量的运动学量而非物理常数。Phillips 等(1994)汇编的 D_n 在 10⁻⁴ 到 10⁻⁷ cm²/s 之间相差三个数量级,作者把这种不确定性归因于测量维度差异(capillary 阵列的尺寸)。这是一个诚实的、引人思考的细节:建模时若直接采用某个引用值而不问"该值是在什么几何条件下测的",可能引入系统偏差。
5.3 节是我个人最喜欢的部分。Tyson 等(1999)的"启发性解析方法"——把二阶常微分方程的系数视为慢变,然后在小 τ 区间内逐段当常数解、再拼接——这种思路在数学物理中其实是古老技巧(WKB 类的多尺度分析),但作者把它用得极为自然,并把解析结果与数值 F(τ)、F₁(τ) 的对比坦诚地摆在图 5.5–5.6 里。结果显示二者形状几乎一致但缩放因子很大,作者把这归因于"二阶导数项在极大邻域外并不小"。这种坦白很可贵。
但 5.4 节给出的预测与数值之间存在一个明显的数量级偏差:预测 k_max = 2.2,实际支配解的是 k_max ≈ 1.1,差了 2 倍。Murray 用 Kgrow = 5.47 解释了这种差异——线性增长率最大的波数 Kgrow 与最终支配的非线性模式 k_max 之间没有简单对应。我作为读者会对这种"线性预测只决定初始增长率、非线性筛选决定最终模式"感到信服,但同时也怀疑:在更广的参数扫描中,这种 2 倍偏差是否稳定?还是说仅在某些参数下出现?书中没给系统性的误差分析。
5.6 节的扩散-化学趋性色散关系 (5.53)–(5.54) 与第 2 章 2.4 节的纯反应扩散色散关系形式上极为相似——都是关于 λ 的二次方程,b(k²)、c(k²) 的二次与四次结构相同。差异只在 D 矩阵的 −α u χ 项,它让"扩散"部分从标量变为 2×2 矩阵。这种结构相似性意味着第 2 章发展的所有 Turing 失稳的判据(零扰动稳定、有限波数带不稳定)可几乎逐字搬到本章。但作者没有显式强调这一点——也许是因为矩阵 D 的非对角形式需要额外小心。
5.7 节的非线性分析是我读得最吃力的部分。Landau 方程 (5.77)(5.78) 的系数 X_A、X_B、Y "是原系统参数的复杂函数"——作者没有给出具体形式,仅说"代数极其复杂"。这有点令人失望:作为读者我会希望看到这些系数的最终表达式(哪怕是 Tyson 1996 中的),以便能亲手验证点图样 vs 条纹图样的分界。但 Murray 显然选择了"勾勒方法 + 引参考文献"的教学策略,这对教学有效但对希望复现的读者不友好。
5.10 节的 swarm ring 稳定性分析是我在数学上感觉最不确定的部分。Routh-Hurwitz 条件 (5.96) 仅在 A、B、C 为常数时严格成立——但作者把它们写成 z 的函数(通过 U、V、W 及其导数),然后说"对于至少某些 z 若不能违反条件即可"——这种"局部不稳定性证明全局不稳定"的论据我不太放心。Murray 自己也说"完整分析是一个有挑战性的未解决问题"。这是科学写作中少见的诚实。
关于 Bacillus subtilis(第 5.11 节):D_b = σ n b 的形式很巧妙——把"细菌只在营养充足且密度适中时才活跃运动"的观察转化为一个简单乘积。但 σ = 1 + ξ 引入随机涨落的做法让我觉得有点临时——为什么不用 Itô 或 Stratonovich 的随机微分方程版本?作者给出"避免了通常此类研究的复杂性"的理由,但我怀疑这种简化在某些参数下可能改变分形维数。
我两个待解答的疑问:(a) 5.10 节基于 Routh-Hurwitz 的 swarm ring 稳定性分析——是否已有完整数值结果证明对某些参数条件 (5.96) 至少一条被违反?Murray 提到"令人惊讶"如果不能违反,但没有给出存在性证明。(b) 5.3 节中 Tyson 等(1999)的多尺度分析与 (5.32) 给出的 τ̃_crit 解析表达式——这个公式在多大参数范围内与数值 τ_crit 吻合?书里只说"很接近",没给误差图。
与上下章的衔接(一段话)
第 4 章(生长域上图样形成)处理的是发育生物学中"形态发生素场在生长胚胎上如何铺设图样"这一核心问题——所有讨论都围绕胚胎、生长曲线、组织层次。本章则把视角转向种群动力学与化学通信:Berg 实验室的细菌图样是发育生物学之外、但仍涉及"局部相互作用产生全局图样"这一普适框架的具体实例。Murray 选择细菌系统作为第 5 章的核心不是偶然:细菌实验数据完备(半固体 petri 皿、接种位置、刺激物浓度全部可控)、时间尺度可测量、且图样极其规则(同心环、点状环、螺旋),使得数学模型能从参数估计直接走到数值验证。Murray 在本章的具体方法论选择——先用 Keller-Segel(1970)这种经典化学趋性框架写出三项守恒方程,再以文献值定参数,最后用色散关系 + Landau 方程作弱非线性分析——正是第 2 章和第 3 章发展的反应扩散分析技术在"扩散张量含非对角项"(即含化学趋性)情形下的直接推广。第 6 章(机械化学图样形成)和第 7 章(演化与形态发生约束)会继续延伸这一思想:第 6 章把机械力纳入守恒方程(细胞不再仅扩散和化学趋化,还受机械应力),第 7 章则把"为什么是这样的图样而不是其他"的问题放到演化选择的框架下。本章因此在抽象数学工具(第 2、3 章)和后续更复杂的生物学应用(第 6、7 章及伤口愈合、血管生成)之间占据承上启下的位置——既展示了反应-扩散-化学趋性数学工具在最干净的实验系统(细菌图样)中的应用,又为后续章节处理更具生物学复杂度的系统提供了基础参考案例。