跳转至

第 11 章:反应扩散、趋化性与非局部机制(Reaction Diffusion, Chemotaxis, and Nonlocal Mechanisms)

11.1 简单随机游走与扩散方程的推导(Simple Random Walk and Derivation of the Diffusion Equation)

本节从最简单的一维随机游走出发推导出宏观扩散方程。作者开篇指出:在一群粒子中(细胞、细菌、化学物、动物等),每个粒子通常做随机运动,这种个体不规则运动让粒子群扩散开来。当这种微观不规则运动在群体上引起某种宏观或整体有规律的运动时,我们可以把它视为扩散过程。当然,粒子之间可能有相互作用,或环境可能给出某种偏向(bias),此时整体运动不再是简单扩散。从个体微观行为得到宏观行为"太难",所以作者们用连续介质模型方程,以粒子密度或浓度为变量描述整体行为。先从随机过程出发用概率方法初等处理,再推导确定性模型是"有教益的"。为简单起见先考虑一维运动与最简随机游走;高维推广从一维方程可"直观地"得到。

考虑粒子沿直线以固定步长 \(\Delta x\)、固定时间步 \(\Delta t\) 做无偏随机行走(每步向右或向左等概率)。经过 \(N\) 步后粒子可在 \(-N\Delta x\)\(N\Delta x\) 之间任一点——若以粒子起点为原点。若在 \(x=0\) 附近释放一群粒子,则空间分布明显不会是均匀的:经过 \(N\) 步后到达 \(x=N\Delta x\) 的概率远小于到达 \(x\) 接近 \(0\) 之处的概率。作者先给出二项式概率公式 \(p(m,n)=(1/2^n)(n!/(a!(n-a)!)))\)\(a=(n+m)/2\)(要求 \(n+m\) 为偶数),表示 \(n\) 步后到达第 \(m\) 个空间位置的概率;这就是 (11.1) 给出的二项分布。具体推导:设粒子向右走 \(a\) 步、向左走 \(b\) 步到达 \(m\Delta x\),则 \(m=a-b\)\(a+b=n\),故 \(a=(n+m)/2\)\(b=n-a\);到达该点 \(x=m\Delta x\) 的可能路径数为 \(n!/(a!b!)=n!/(a!(n-a)!)=C_a^n\)(二项式系数,由 \((x+y)^n=\sum_{a=0}^n C_a^n x^{n-a}y^a\) 给出)。\(n\) 步路径总数为 \(2^n\),故概率 \(p(m,n)=(\text{有利情形})/(\text{总情形})=1/2^n \cdot n!/(a!(n-a)!)\)\(a=(n+m)/2\)\(n+m\) 为偶数。

可用 \(\sum_{m=-n}^n p(m,n)=\sum_{a=0}^n C_a^n (1/2)^{n-a}(1/2)^a=(1/2+1/2)^n=1\) 验证——总概率必须为 1,作者明文写出这一恒等式作为校验。用 Stirling 公式 \(n!\sim(2\pi n)^{1/2}n^n e^{-n}\)\(n,m\gg 1\) 下做渐近展开(Stirling 公式可由 \(n!=\int_0^\infty e^{-t}t^n\,dt\)(gamma 函数定义)配合 Laplace 方法对大 \(n\) 渐近逼近得到;见 Murray 1984 渐近分析教材),得到正态(高斯)概率分布 \(p(m,n)\sim(2/(\pi n))^{1/2}\exp[-m^2/(2n)]\)。作者提醒这一近似在 \(n\) 并不很大时已很准确——例如 \(n=8,m=6\) 时偏差 < 5%,\(n=10,m=4\) 时偏差 < 1%;\(n>6\) 时即可放心使用。作者评论:"渐近近似常在比想象更广的范围上准确得令人惊讶。"

接着令 \(m\Delta x=x, n\Delta t=t\),但作者强调不能直接把 \(p\) 当成最终变量(因为 \(\Delta x\to 0\) 时概率本身也趋于零),更合适的依赖变量是 \(u=p/(2\Delta x)\),即 \(2u\Delta x\) 是粒子落在区间 \((x, x+\Delta x)\) 内的概率。代入并取 \(\lim_{\Delta x\to 0,\Delta t\to 0}(\Delta x)^2/(2\Delta t)=D\neq 0\) 的极限,得 \(u(x,t)=(1/(4\pi Dt))^{1/2}\exp[-x^2/(4Dt)]\)\(D\) 即扩散系数(或扩散率),量纲为 (长度)\(^2\)/(时间);它度量粒子从高密度向低密度扩散的效率。作者给出两个具体数值例子:血液中血红蛋白分子 \(D\sim 10^{-7}\,\mathrm{cm^2\,sec^{-1}}\),而血液中氧 \(D\sim 10^{-5}\,\mathrm{cm^2\,sec^{-1}}\)

下一段把这一结果与 Fick 扩散定律联系起来。Fick 定律说物质通量 \(J\)(可以是细胞、化学物、动物等)与浓度梯度成正比:一维情况下 \(J\propto -\partial c/\partial x\),即 \(J=-D\,\partial c/\partial x\)\(c(x,t)\) 是浓度,\(D\) 是扩散率),其中负号表示扩散将物质从高浓度向低浓度输运。对区间 \([x_0, x_1]\) 写出守恒方程——区域内物质变化率 = 边界流入率 + 区域内生成率(无生成时仅边界流)——得到 \(\partial/\partial t \int_{x_0}^{x_1} c\,dx = J(x_0,t)-J(x_1,t)\)。令 \(x_1-x_0\to 0\) 并用 Fick 定律,给出经典扩散方程 \(\partial c/\partial t=-\partial J/\partial x=\partial(D\,\partial c/\partial x)/\partial x\);若 \(D\) 为常数则简化为 \(\partial c/\partial t=D\,\partial^2 c/\partial x^2\)。若在 \(t=0\)\(x=0\) 释放每单位面积 \(Q\) 个粒子(初值 \(c(x,0)=Q\delta(x)\)\(\delta(x)\) 是 Dirac 三角函数),则 (11.8) 的解为 \(c(x,t)=Q/(2(\pi Dt)^{1/2})\exp[-x^2/(4Dt)]\)\(t>0\)(见 Crank 1975);这与从随机游走方法所得的 (11.4) 在 \(Q=1\) 时一致(图 11.1 定性展示)。

之后作者从随机游走出发再做一次"准严格"推导:设在时间 \(t\) 粒子到达 \(x\) 的概率为 \(p(x,t)\)——从 (11.4) 出发,\(p(x,t)\) 是粒子在 \(t=0\)\(x=0\) 出发、于时间 \(t\) 到达 \(x\) 的概率。设 \(\alpha,\beta\) 表示左右移动概率(\(\alpha+\beta=1\),无偏各向同性时 \(\alpha=\beta=1/2\)),则 \(p(x,t)=\alpha p(x-\Delta x, t-\Delta t)+\beta p(x+\Delta x, t-\Delta t)\)。对右端做 Taylor 展开得 \(\partial p/\partial t=[(\Delta x)^2/(2\Delta t)]\,\partial^2 p/\partial x^2+(\Delta t/2)\,\partial^2 p/\partial t^2+\cdots\);取 \(\lim (\Delta x)^2/(2\Delta t)=D\)(与 (11.3) 取极限的方式相同)后得 \(\partial p/\partial t=D\,\partial^2 p/\partial x^2\)。若释放总数 \(Q\) 个粒子,则浓度 \(c(x,t)=Qp(x,t)\),于是得到 (11.8)。作者承认这一推导仍依赖于步长与时间步以"相当特定的方式"趋于零(使得 \(D\) 存在)——更严谨的方式是借助 Fokker–Planck 方程结合 Markov 过程:所谓 Markov 过程即时刻 \(t\) 的状态只依赖时刻 \(t-\Delta t\) 的状态,"一世代时间依赖"过程;作者指出 Skellam(1973)与 Okubo(1980, 1986)给出了更精细的推导,并讨论了多种扩散方程的建立——Okubo 1980 专著为上述极限过程给出"相当多的证明",Okubo 1986 综述也讨论了多种扩散方程的推导。

11.2 反应扩散方程(Reaction Diffusion Equations)

本节将上一节的守恒-扩散框架推广到三维并加入反应项。设 \(S\) 为任意包围体积 \(V\) 的曲面,守恒方程说 \(V\) 内物质量变化率 = 通过 \(S\) 流入 \(V\) 的通量 + \(V\) 内生成量。写为 \(\partial/\partial t\int_V c\,dv=-\int_S J\cdot ds+\int_V f\,dv\),其中 \(J\) 是通量,\(f\) 是源项(可为 \(c,x,t\) 的函数)。用散度定理化简曲面积分,假设 \(c\) 连续,得到积分形式的 \(\int_V[\partial c/\partial t+\nabla\cdot J-f]\,dv=0\);由于 \(V\) 任意,被积函数必须为零,于是有 \(\partial c/\partial t+\nabla\cdot J=f\)。该方程对任意通量机制(扩散或其他)都成立。

对经典 Fick 扩散,\(J=-D\nabla c\),方程变为 \(\partial c/\partial t=f+\nabla\cdot(D\nabla c)\);这里 \(D\) 可为 \(x\)\(c\) 的函数,\(f\) 可为 \(c,x,t\) 的函数。作者指出 \(D\) 依赖于空间坐标的情形在生物医学建模中日益重要——典型例子是转基因生物在异质环境中的扩散,以及大脑白质/灰质对脑肿瘤生长与扩散的影响。

从生态学角度看,\(f\) 可表示出生-死亡过程,\(c\) 是种群密度 \(n\)。对 logistic 增长 \(f=rn(1-n/K)\)\(r\) 为线性增长率,\(K\) 为环境承载量)配合常数 \(D\) 给出 \(\partial n/\partial t=rn(1-n/K)+D\nabla^2 n\)——这是著名的 Fisher–Kolmogorov 方程。Fisher(1937)首次以一维形式把它作为有利基因在种群中扩散的模型提出,Kolmogoroff 等(1937)深入研究并给出主要解析结果;第 13 章会详细讨论这一方程。

若把 (11.16) 进一步推广到多个相互作用的物种或化学物质,得到向量 \(u(x,t)\)\(i=1,\ldots,m\),各分量以自身扩散系数 \(D_i\) 扩散、向量源项 \(f\) 相互作用,则有 \(\partial u/\partial t=f+\nabla\cdot(D\nabla u)\)\(D\) 是一般扩散系数矩阵(若物种间无交叉扩散,\(D\) 即对角阵);\(\nabla u\) 是张量,故 \(\nabla\cdot D\nabla u\) 是向量。方程 (11.18) 称为反应扩散系统。这一机制是 Turing(1952)作为形态发生的化学基础提出的,被作者誉为"20 世纪理论生物学最重要的论文之一";自 1970 年前后起被广泛研究。本卷主要讨论 \(D\) 对角且常值、\(f\) 仅依赖 \(u\) 的情形;进一步推广可在 \(f\) 中加入反映种群历史的积分项。在某些癌症模型中(如脑胶质母细胞瘤及其他涉及突变癌细胞的肿瘤)会出现交叉扩散项与不相等的对角扩散项——这些情形的数学推广似乎是"无穷无尽"的。作者提醒读者:对真实世界的大多数实际模型而言,在简化版本已与实验/观察数据相比较之前,就花大量时间于复杂推广是"至少说过分乐观"——脚注引用 de Tocqueville 语:"推广无意义,因为上帝知道所有特例"(As de Tocqueville remarked, there is no point in generalising since God knows all the special cases.)。

本节再附注一点:交叉扩散在真正实用的模型中并不常见——一个例子是第 II 卷第 1.2 节将要讨论的。交叉扩散系统在良态性(well-posedness)等问题上可引出有趣的数学问题。

本节余下部分介绍与反应扩散方程相关但本章不深入的两个重要方向。第一个是促进扩散(facilitated diffusion),它与肌肉中氧的运输密切相关。O₂ 与血红蛋白、肌红蛋白可逆结合的生化动力学使 O₂ 在血液-肌肉界面间传递——例如氧可与血红蛋白可逆结合、与肌红蛋白可逆结合;肌红蛋白在肌肉中至关重要,但作为促进剂效率不如血红蛋白。Wyman(1966)启发式论证:若无促进扩散,肌组织无法存活。该领域在实验上由 Wittenberg(1970 综述;Wittenberg et al. 1975)及其同事做了深入工作;蛋白促进质子方面由 Gros 等(1976, 1984)做了实验,揭示出"质子由血红蛋白与其他蛋白通过旋转扩散方式促进跨组织扩散"的现象。Wyman(1966)、Murray(1971, 1974)、Murray 与 Wyman(1971)建立了氧与一氧化碳的促进扩散数学理论。促进有效的条件是组织内存在反应平衡区——这意味着近表面必存在非平衡边界层(Murray 1971, Mitchell 与 Murray 1973, Rubinow 与 Dembo 1977)。这一非平衡边界层的存在条件解释了为何血红蛋白比肌红蛋白是更好的 O₂ 促进剂、为何一氧化碳不被肌红蛋白促进。促进扩散在 CO 中毒与解毒的困难中也起关键作用(Britton 与 Murray 1977)。第二个方向是质子促进扩散,它由 Gros 等(1976, 1984,及其早期引用文献)实验揭示:质子使血红蛋白分子旋转,从而提高跨组织扩散。Murray 与 Smith(1986)给出了旋转扩散的数学理论,复杂得多。

11.3 动物扩散模型(Models for Animal Dispersal)

本节把扩散模型应用于昆虫与动物的扩散、入侵。Okubo(1980, 1986)、Shigesada(1980)、Lewis(1997)详细讨论了动物种群模型的各方面。相互作用物种的扩散由 Shigesada 等(1979)讨论,竞争物种的扩散由 Shigesada 与 Roughgarden(1982)讨论。Kareiva(1983)证明许多物种的扩散可用常扩散系数的反应扩散模型描述,并给出了他本人从多种昆虫实验测得的扩散系数实测值;Kot 等(1996)综合了真实数据研究生物扩散(含 Kot 2001)。昆虫种群的一个常见特征是离散时间的种群增长——按直觉这会对空间扩散有重要影响。模型方程涉及离散时间与连续空间的耦合,Kot(1992)与 Neubert 等(1995)研究过这一主题。Tilman 与 Kareiva(1998)编辑的论文集是该方向空间作用的实用资料来源,文章涉及濒危物种持续性、生物多样性、疾病动力学、多物种竞争等议题。Renshaw(1991)与 Williamson(1996)专著是研究物种入侵的优秀教材,含大量实例。Shigesada 与 Kawasaki(1997)专著更具数学与建模导向,讨论了哺乳动物、鸟类、昆虫、植物的各种生物入侵形式,扩散只是其中一种机制——对认真建模这些现象的人而言,这些书是必读。

昆虫扩散有多种形式,其中一个对昆虫尤其相关的是扩散随种群压力增大——即 \(D\) 是种群密度 \(n\) 的增函数:\(J=-D(n)\nabla n\)\(dD/dn>0\)。典型形式 \(D(n)=D_0(n/n_0)^m\)\(m>0\)\(D_0,n_0>0\))。无增长项的扩散方程为 \(\partial n/\partial t=D_0\nabla\cdot[(n/n_0)^m\nabla n]\),一维形式 \(\partial n/\partial t=D_0\,\partial/\partial x[(n/n_0)^m \partial n/\partial x]\),它有解析解

\[ n(x,t)=\frac{n_0}{\lambda(t)}\left[1-\left(\frac{|x|}{r_0\lambda(t)}\right)^{2}\right]^{1/m}, \quad |x|\le r_0\lambda(t); \qquad 0, \quad |x|>r_0\lambda(t) \]

其中 \(\lambda(t)=(t/t_0)^{1/(2+m)}\)\(r_0=[Q\Gamma(1/m+3/2)/\{\pi^{1/2}n_0\Gamma(1/m+1)\}]\)\(t_0=r_0^2 m/(2D_0(m+2))\)\(\Gamma\) 是 gamma 函数,\(Q\) 是释放于原点的初始昆虫数。容易验证 (11.21) 对所有 \(r_0\) 都是 (11.20) 的解(\(r_0\) 的值由要求 \(\int n\,dx=Q\) 给出)。在另一语境下 (11.20) 称为多孔介质方程。该解在结构上与经典扩散解 (11.10) 截然不同——因为 \(D(0)=0\)。解代表一种波,在 \(x=x_f=r_0\lambda(t)\) 处有前沿;此处 \(n\) 的导数不连续。前沿——即 \(n=0\) 之处——传播速度 \(dx_f/dt=r_0 d\lambda/dt\),由 (11.22) 知对所有 \(m\) 该速度随时间减小(图 11.2 示意)。蝗虫的扩散模式表现出与该模型相似的行为(Aikman 与 Hewitt 1972)。无源项时,\(t\to\infty\)\(n\to 0\)。Shigesada(1980)提出 \(D(n)\propto n\) 的线性密度依赖模型;见 Shigesada 与 Kawasaki(1997)。

径向对称的平面扩散(\(Q\) 昆虫释于 \(r=0\))满足 \(\partial n/\partial t=(D_0/r)\partial/\partial r[r(n/n_0)^m \partial n/\partial r]\),解为

\[ n(r,t)=\frac{n_0}{\lambda^2(t)}\left[1-\left(\frac{r}{r_0\lambda(t)}\right)^2\right]^{1/m}, \quad r\le r_0\lambda(t); \quad 0, \quad r>r_0\lambda(t) \]

其中 \(\lambda(t)=(t/t_0)^{1/(2(m+1))}\)\(t_0=r_0^2 m/(4D_0(m+1))\)\(r_0^2=(Q/(\pi n_0))(1+1/m)\)\(m\to 0\) 时这些解趋于常扩散解(用 \(\exp[s]=\lim_{m\to 0}(1+ms)^{1/m}\) 证明 (11.21)、(11.24) 趋于 (11.10) 等常扩散解;这一极限需要"一些代数与指数定义")。

昆虫在低密度下倾向于聚集。一种一维模型取通量 \(J=Un-D(n)\partial n/\partial x\),其中 \(U\) 是输运速度;若吸引中心是原点、吸引速度为常数,Shigesada 等(1979)取 \(U=-U_0\,\mathrm{sgn}(x)\),扩散方程变为

\[ \frac{\partial n}{\partial t}=U_0\frac{\partial}{\partial x}[n\,\mathrm{sgn}(x)]+D_0\frac{\partial}{\partial x}\left[\left(\frac{n}{n_0}\right)^m\frac{\partial n}{\partial x}\right] \]

这一方程"并不容易求解"。但可分区间近似:\(t\approx 0\)\(x=0\) 附近 \(n\) 梯度大,扩散项主导,对流项小,解近似为 (11.21);\(t\to\infty\) 时种群达到稳态空间非均匀分布,对流与扩散效应平衡。把 (11.25) 中 \(\partial n/\partial t=0\) 积分两次(用 \(n\to 0,\partial n/\partial x\to 0\)\(|x|\to\infty\))得稳态空间分布

\[ \lim_{t\to\infty}n(x,t)\to n(x)=n_0\left(1-\frac{mU_0|x|}{D_0}\right)^{1/m}, \quad |x|\le D_0/(mU_0); \quad 0, \quad |x|>D_0/(mU_0) \]

这一推导留作练习。解 (11.26) 表明扩散在 \(x\) 上是有限的;\(m=1/2\) 时的形式与 Okubo 与 Chiang(1974)观察到的特定蚊群分布相似(见 Okubo 1980 图 9.6;图 11.3 示意)。

作者补充:昆虫扩散是一个非常重要但仍未被完全理解的课题;上述模型虽简单,但仍能为可能的昆虫扩散行为给出一些指示。若在 (11.25) 右端加入种群增长/死亡项(练习 3 与 4)则昆虫种群如预期那样灭绝(无出生只有死亡),但有趣的是昆虫仅从原点移动有限距离。扩散模型在动物与昆虫扩散研究中正被越来越广泛地应用(Okubo 1980 与 Shigesada & Kawasaki 1997 给出大量例子)。Okubo(1986)综述讨论了多种模型并专门处理动物群组、昆虫群与聚集;Mogilner 与 Edelstein-Keshet(1999)讨论基于非局部相互作用的群组模型——含长程吸引与排斥,证明当排斥项的密度依赖阶数高于吸引项时,群有恒定内部密度与锐利边缘(即"看起来像一群")。第 II 卷第 14 章专门讨论狼的扩散与领地形成,以及狼-鹿(狼的主食)的生存;其中将考虑一种更现实的吸引中心项形式(对狼来说吸引中心是夏季窝),不产生图 11.3 那样的梯度间断。

11.4 趋化性(Chemotaxis)

本节研究由化学梯度引导的定向运动。许多昆虫与动物(包括人类)依赖敏锐嗅觉在同种个体间传递信息;涉及的化学物质称为信息素(pheromone)。例如雌性蚕蛾 Bombyx mori 释放称为 bombykol 的性信息素,雄蚕蛾有极其灵敏的触角滤波器测量其浓度,并沿浓度上升方向移动;建模涉及流体力学与过滤理论跨尺度耦合(Murray 1977)。许多深海鱼类的敏锐嗅觉对交流与捕食至关重要。本节对定向运动(即趋化性)建立模型——它与扩散不同,是沿浓度梯度向上的运动。

趋化性不仅在动物与昆虫生态中重要,也在生物学过程里起关键作用——例子之一是细菌感染被白细胞向感染源靠趋化性移动所抵抗;有说服力的证据表明白细胞向细菌感染区域靠趋化性移动(由感染造成的化学梯度驱动)以对抗感染(Lauffenburger 与 Keller 1979, Tranquillo 与 Lauffenburger 1986, 1988, Alt 与 Lauffenburger 1987)。第 II 卷第 5 章将详细讨论细菌趋化性。一个被广泛研究的趋化性现象是粘菌 Dictyostelium discoideum——单细胞变形虫向其自身产生的化学物质 cAMP(环腺苷酸)相对高浓度区域移动;实验观察到有趣的波状运动与空间斑图(第 II 卷第 1 章)。Segel(1984)专著给出了早期数学模型;随着生物系统认识加深,模型必然改变,Martiel 与 Goldbeter(1987)、Goldbeter(1996)、Monk 与 Othmer(1989)提出了更复杂、更生物上现实的模型,都表现振荡行为。Spiro 等(1997)给出了细菌趋化性中激发与适应的模型,纳入详细生化数据并用于阐释实验过程。

粘菌的空间斑图模型大多基于化学引诱物与细胞的连续介质模型。Dallon 与 Othmer(1997)提出了一个新模型,细胞视为离散实体、化学引诱物浓度视为连续;结果与许多已有实验结果吻合,可用于研究不同细胞运动规则对聚集模式与波动的影响(包括普遍的螺旋波起源)。第 13 章与第 II 卷第 1 章会详细讨论波动传播(含螺旋波)。本节将构建一组基础反应-扩散-趋化方程。

设引诱物 \(a(x,t)\) 存在梯度会驱动细胞(如)向梯度方向运动;细胞通量随细胞数 \(n(x,t)\) 增加,故可合理取趋化通量 \(J=n\chi(a)\nabla a\),其中 \(\chi(a)\) 是引诱物浓度的函数。在 \(n\) 的一般守恒方程 \(\partial n/\partial t+\nabla\cdot J=f(n)\) 中,\(J=J_{\text{diffusion}}+J_{\text{chemotaxis}}\)——扩散项用 (11.15),趋化项用 (11.27)——得

\[ \frac{\partial n}{\partial t}=f(n)-\nabla\cdot(n\chi(a)\nabla a)+\nabla\cdot(D\nabla n) \]

其中 \(D\) 是细胞扩散系数。引诱物 \(a\) 也是化学物质,会扩散并被变形虫产生,故需附加方程 \(\partial a/\partial t=g(a,n)+\nabla\cdot(D_a\nabla a)\)\(D_a\)\(a\) 的扩散系数,\(g\) 是动力学/源项(可依赖 \(n\)\(a\)),一般 \(D_a>D\)。若多个物种或细胞类型都对引诱物响应,可把 (11.28) 推广为向量形式,每物种有自身的 \(\chi(a)\);第 II 卷第 5 章会展示 \(f(n),g(a,n),\chi(a)\) 的具体形式对生成模式的决定性影响。

Keller 与 Segel(1971)的开创性粘菌模型取 \(g(a,n)=hn-ka\)\(h,k>0\)),\(hn\) 表示引诱物的自发产生正比于变形虫数 \(-ka\) 表示引诱物活性的指数衰减。简化版本可取 \(f(n)=0\)(变形虫产生率可忽略,对应模式形成阶段),\(\chi(a)=\chi_0\)(正常数);与常扩散系数结合,\(g\) 用上述线性形式,得到一维非线性系统

\[ \frac{\partial n}{\partial t}=D\frac{\partial^2 n}{\partial x^2}-\chi_0\frac{\partial}{\partial x}\left(n\frac{\partial a}{\partial x}\right),\quad \frac{\partial a}{\partial t}=hn-ka+D_a\frac{\partial^2 a}{\partial x^2} \]

第 II 卷第 1 章会研究此系统,其中 \(n\) 视为细菌种群,\(a\) 为其消耗的食物。其他趋化因子 \(\chi(a)\) 的形式包括对数律 \(\chi(a)=\chi_0/a\) 与受体律 \(\chi(a)=\chi_0 K/(K+a)^2\)\(\chi_0,K>0\))——\(a\) 减小时趋化效应增大;第 II 卷第 5 章讨论 Salmonella 与 E. coli,并给出实验推导的 \(f(n),g(a,n),\chi(a)\) 形式。

实际可测的趋化性指数 \(I\) 有多种定义方式。以细胞向位于 \(x_s\) 的化学引诱物源做平面运动为例:细胞起始于 \(x_A\),距源 \(D_1\);无趋化时细胞运动纯随机,在给定时间 \(T\) 内朝 \(x_s\) 方向移动的平均距离 \(D_2\) 为零;有趋化时随机运动被修正,细胞一般有向源运动的趋势,相同时间 \(T\)\(D_2>0\)。定义 \(I=D_2/D_1\)\(I\) 越大表示趋化越强。Tranquillo 与 Lauffenburger(1988)分析了白细胞详细的化学感受运动以确定其趋化参数;Woodward 等(1995)、Tyson(1996)、Murray 等(1998)、Tyson 等(1998, 1999)给出了 Salmonella 与 E. coli 的实验趋化参数。

有些细胞会受外加电场影响,沿场方向运动,称为电趋性(galvanotaxis);其强度可类似定义。若 \(V\) 是电势,电趋性通量可取 \(J\propto nG(V)\nabla V\)\(G\) 可为电压 \(V\) 的函数。作者最后提醒:(11.28) 与 (11.30) 中扩散项与趋化项的符号不同:扩散是稳定化力(一般),趋化是去稳定化力(一般),类似"负扩散"。可合理推测 (11.30) 模型中稳定与去稳定力的平衡可产生 \(n\)\(a\) 的稳态空间斑图或非稳态波状非均匀结构(即细胞密度的非均匀空间模式;见第 II 卷第 1 章与第 5 章)。若趋化效应足够强,可能发生解爆破(solution blow-up);Jäger 与 Luckhaus(1992)等给出了相关参考文献。

11.5 非局部效应与长程扩散(Nonlocal Effects and Long Range Diffusion)

本节说明经典扩散框架的限制并引入长程扩散。经典 Fick 扩散(\(J=-D\nabla c\)\(J=-D(n)\nabla c\))只严格适用于稀薄系统;但实际应用更广。然而这些形式本质上是"局部"或"短程"扩散——可从 Laplacian 算符 \(\nabla^2 n\) 看出:\(\nabla^2 n\propto (\langle n(x,t)\rangle-n(x,t))/R^2\)\(R\to 0\),其中 \(\langle n\rangle\) 是以 \(x\) 为球心、半径 \(R\) 的球内平均密度 \(n_{av}=\langle n(x,t)\rangle\equiv (3/(4\pi R^3))\int_V n(x+r,t)\,dr\)\(V\) 是半径 \(R\) 的球);这一 Laplacian 的解释最早由 Maxwell 在 1871 年提出(见 Maxwell 1952,论文集)。因 \(R\to 0\) 可将 \(n(x+r,t)\)\(x\) 处 Taylor 展开为 \(n(x,t)+(r\cdot\nabla)n+1/2(r\cdot\nabla)^2 n+\cdots\) 并代入球平均积分;由对称性一阶项积分为零;若忽略 \(O(r^3)\) 及更高项,积分得 \(n_{av}=n(x,t)+(3/10)R^2\nabla^2 n(x,t)\)(比例常数 10/3 即 (11.32) 的精确形式)。

但许多生物领域(如胚胎发育)涉及细胞密度较大、单纯梯度比例的局部扩散通量不足。当第 II 卷第 6 章讨论生物模式形成的力学理论时,会说明在某些情况下包含长程效应是直观合理甚至必要的。设通量 \(J=G_{r\in N(x)}[\nabla n(x+r,t)]\)\(N(x)\)\(x\) 的某邻域(在该邻域内 \(x\) 处"感受到"效应),\(G\) 是关于梯度的泛函;由对称性与各向同性假设可得 \(J\) 对简单线性 \(\nabla n\) 的首阶修正为 \(\nabla(\nabla^2 n)\) 项,即 \(J=-D_1\nabla n+\nabla(D_2\nabla^2 n)\)\(D_1>0\)\(D_2\) 常数);\(D_2\) 度量长程效应,其量级一般小于 \(D_1\)。这一方法由 Othmer(1969)给出,详尽讨论了一般泛函 \(G\) 的形式、推导与构造;本节及 11.6 节给出长程 \(D_2\) 项的不同动机。将 (11.36) 代入守恒方程 (11.14)(\(f\equiv 0\))得 \(\partial n/\partial t=\nabla\cdot D_1\nabla n-\nabla\cdot\nabla(D_2\nabla^2 n)\)——第一项代表最近邻平均,第二项(双调和项)代表最近邻平均的"平均"。

双调和项的稳定/去稳定性质:\(D_2>0\) 时稳定,\(D_2<0\) 时去稳定。设 (11.37) 的试探解 \(n(x,t)\propto\exp[\sigma t+ik\cdot x]\)\(k=|k|\),代入得色散关系 \(\sigma=-D_1 k^2-D_2 k^4\)。由 (11.40) 可见:\(k^2>D_1/|D_2|\) 时,若 \(D_2>0\)\(n\to 0\),若 \(D_2<0\)\(n\to\infty\)\(t\to\infty\))。经典 Fick 扩散 \(D_2\equiv 0\) 时所有 \(k\) 都有 \(n\to 0\)\(D_2>0\) 双调和贡献是稳定的,\(D_2<0\) 是去稳定的。色散关系在很多上下文中重要(第 II 卷第 2 章详细讨论)。

另一种重要的长程效应建模方法是积分方程形式。\(x\)\(n\) 的变化率依赖所有其他位置 \(x'\)\(n\) 的影响。一维例子为

\[ \frac{\partial n}{\partial t}=f(n)+\int_{-\infty}^{\infty} w(x-x')n(x',t)\,dx' \]

其中 \(w(x-x')\) 是核函数,量化邻近 \(n(x',t)\)\(n(x,t)\) 的影响;该形式假设影响只依赖 \(x\)\(x'\) 的距离。\(f(n)\) 是通常的源/动力学项(与应用到神经元细胞时称为发放率);假设 \(w\to 0\)\(|x-x'|\to\infty\),且 \(w\) 空间对称 \(w(x-x')=w(x'-x)\)。该模型直接包含长程效应:\(w\) 衰减快(如 \(\exp[-(x-x')^2/s]\)\(0<s\ll 1\))则长程效应弱,\(s\gg 1\) 则长程效应强。

要确定 (11.41) 解的时空性质,必须指定核 \(w\)。这涉及具体生物现象的建模。对神经元细胞(可自发发放的细胞,\(n\) 表示发放率):\(f(n)\) 是自主空间无关的发放率,无邻近细胞影响时发放率演化到 \(f(n)\) 的零决定的稳定稳态(与第 1 章相同)。例如 \(f(n)\) 如图 11.4(a) 时发放率演化到单一稳态 \(n_1\);图 11.4(b) 时存在临界发放率 \(n_c\)\(n<n_c\) 时发放率趋于零(灭绝),\(n>n_c\) 时趋于 \(n_1\)

若考虑空间效应,必须包含邻近细胞影响,即指定核 \(w\)。假设细胞同时接受兴奋与抑制输入,最强兴奋信号来自其自身(即高发放态倾向于增加,类似自催化——"如果一个细胞处于高发放态 \(n\),它倾向于增加")——图 11.5 是这种核的典型形式(具有兴奋-抑制核的典型形状)。

可将积分方程方法与长程扩散近似联系起来。令 \(y=x-x'\),则 \(\int_{-\infty}^\infty w(x-x')n(x',t)\,dx'=\int_{-\infty}^\infty w(y)n(x-y,t)\,dy\)。将 \(n(x-y)\)\(x\) 处 Taylor 展开(与对 (11.33) 的展开相同的方式),得 (11.43) 所示的逐项展开。由核 \(w(y)\) 的对称性,\(\int_{-\infty}^\infty y^{2m+1}w(y)\,dy=0\)\(m=0,1,2,\ldots\)),所有奇次项消失。定义核矩 \(w_{2m}=1/(2m)!\int_{-\infty}^\infty y^{2m}w(y)\,dy\)\(m=0,1,2,\ldots\),则 (11.41) 化为

\[ \frac{\partial n}{\partial t}=f(n)+w_0 n+w_2\frac{\partial^2 n}{\partial x^2}+w_4\frac{\partial^4 n}{\partial x^4}+\cdots \]

典型核的高阶矩迅速减小(从 (11.45) 的定义可直观看出)。若在四阶矩截断,得到与 (11.37) 形式相近的双调和方程。解的行为关键依赖于核矩的符号与核的具体形式(如图 11.5):\(w_2<0\) 时"短程扩散"项去稳定(对比 (11.39) 中 \(D_1<0\) 的情形),\(w_4<0\) 时"长程扩散"项稳定(与 \(D_2>0\) 类似)。这一积分方程方法更具描述性,更贴合生物实际:它以更"描述性"的方式反映生物上正在发生什么。第 II 卷第 12 章讨论产生空间斑图的神经模型时会深入研究这种模型。

11.6 细胞势能与扩散/长程效应的能量方法(Cell Potential and Energy Approach to Diffusion and Long Range Effects)

本节给出长程扩散项的另一推导方法,处理对象是细胞种群密度;先回顾经典 Fick 扩散推导再推广。处理方式遵循 Cohen 与 Murray(1981)。如果存在势 \(\mu\) 的梯度,它可驱动通量 \(J\) 经典地正比于 \(\nabla\mu\)。把势理解为改变状态所做功(即能量的变分导数):设 \(n(x,t)\) 为细胞密度,对空间分布 \(n\) 关联一能量密度 \(e(n)\),使体积 \(V\) 内总能量 \(E[n]=\int_V e(n)\,dx\)。改变量 \(\delta E\)(即改变状态 \(\delta n\) 做的功)是变分导数 \(\delta E/\delta n\),它定义一势 \(\mu(n)=e'(n)\)。势梯度产生通量 \(J=-D\nabla\mu(n)\)\(D\) 是比例参数(可依赖 \(x,t,n\))。连续性方程变为 \(\partial n/\partial t=-\nabla\cdot J=\nabla\cdot[D\nabla\mu(n)]=\nabla\cdot[D e''(n)\nabla n]\),即 \(\partial n/\partial t=\nabla\cdot[D^*(n)\nabla n]\)\(D^*(n)=De''(n)\)。经典常扩散情形 \(e(n)=n^2/2\)\(\mu(n)=n\),给出 \(\partial n/\partial t=D\nabla^2 n\)。若 \(D\) 依赖 \(x,t,n\),守恒方程为 \(\partial n/\partial t=\nabla\cdot[D^*(x,t,n)\nabla n]\)\(n\) 可为细胞物种向量。

此推导假设 \(e(n)\) 仅依赖 \(n\)。若细胞对环境的响应超出紧邻范围,则合理的假设是"维持空间异质性所需的能量"依赖于邻近密度的梯度。取一更现实的能量泛函(要求在反射 \(x_i\to -x_i\) 与旋转下不变):\(E[n]=\int_V[e(n)+k_1\nabla^2 n+k_2(\nabla n)^2+\cdots]\,dx\)\(k_i\) 可依赖 \(n\)。用 Green 定理(边界 \(\partial n/\partial N=0\) 零通量)化简 \(\int k_1\nabla^2 n\,dx=-\int k_1'(n)(\nabla n)^2\,dx\)(边界项忽略),得到空间异质情形下的能量泛函 \(E[n]=\int_V[e(n)+k_2'(\nabla n)^2+\cdots]\,dx\)\(k_2'=-k_1'(n)+k_2\)。其中 \(e(n)\) 是空间均匀下的能量密度,其他项代表依赖邻近空间密度变化的"梯度"能量密度。

再按 (11.48)→(11.53) 的同样步骤求势 \(\mu\)\(\mu=\mu(n,\nabla n)=\delta E[n]/\delta n=-k\nabla^2 n+e'(n)\)(用变分法,取 \(k\) 为常数)。通量 \(J=-D^*\nabla\mu(n,\nabla n)\),扩散方程为

\[ \frac{\partial n}{\partial t}=D^*\nabla^2[-\nabla^2 n+e'(n)]=-kD^*\nabla^4 n+D^*\nabla\cdot[e''(n)\nabla n] \]

\(D^*\)\(k\) 都取常数。对 \(e(n)\) 的基本假设:因能量密度不能依赖 \(n\) 的符号,故 \(e(n)\) 只含 \(n\) 的偶次幂。Landau–Ginzburg 自由能形式 \(e(n)=an^2/2+bn^4/4\) 代入 (11.58) 得 \(\partial n/\partial t=-D^*k\nabla^4 n+D^* a\nabla^2 n+D^* b\nabla^2 n^3\)。设 \(D_1=D^*a\)\(D_2=D^*k\)\(D_3=D^*b\),得广义扩散方程

\[ \frac{\partial n}{\partial t}=D_1\nabla^2 n-D_2\nabla^4 n+D_3\nabla^2 n^3 \]

注意额外的非线性项 \(D_3\) 的出现。若 \(e(n)\) 仅含通常的 \(n^2\)\(b=0\)),则 (11.59) 化为 (11.37)。把反应/动力学项 \(f(n)\) 加入 (11.59) 即得到等价于 (11.14) 的广义反应扩散方程。Cohen 与 Murray(1981)证明 (11.59) 一维标量形式配 logistic 增长 \(f(n)\) 后,可呈现稳态空间非均匀解。Lara-Ochoa(1984)分析了二维情形,证明该模型能反映多细胞系统形态发生的若干方面。

本章个人批注

本章是 Murray 第一卷中数学内容最丰富、跨领域最广的一章——从最基础的随机游走与扩散方程推导出发,一路推进到趋化性、非局部效应、长程扩散与能量方法,最后用 Landau–Ginzburg 自由能推导出带双调和项的广义反应扩散方程 (11.59)。它在前 10 章建立的 ODE 框架(种群动力学、流行病、振子)与第 II 卷将要深入讨论的 PDE 时空斑图之间扮演了承上启下的过渡角色。

让我印象最深的是 11.5 节"积分方程→核矩展开→双调和项"这一连串推导。把 (11.41) 中的卷积 \(w*n\) 用 Taylor 展开逐阶分离出 \(\int y^{2m}w(y)\,dy\) 这种核矩,然后把 \(w_2\) 与经典扩散联系、\(w_4\) 与长程扩散联系——这给 PDE 框架一个严格的"微观→宏观"基础。Othmer 1969 年的工作至今仍是这一思路的源头,但在生物 PDE 教材中很少被系统地展开。读到这里让我想到 Ogawa(1980s)关于非局部 PDE 的工作与 Murray 这里的处理在精神上是相通的——只是 Murray 选择了更"工程化"的色散关系语言,而 Ogawa 更偏泛函分析。

11.4 节的趋化性模型也有值得反复品味的地方。Keller-Segel 模型 (11.30) 看似简单,物理内容却极其丰富:扩散项(稳定)+ 趋化项(去稳定)+ 衰减项的动力学。Murray 在末尾点出 (11.28) 中扩散与趋化项的符号差异是理解整个系统行为的钥匙——这种"用符号判稳定性"的分析哲学与第 7 章中生物振子极限环稳定性分析(附录 B 的 Routh-Hurwitz 判据)形成对照。Jäger-Luckhaus(1992)的解爆破(blow-up)现象则是该方向几十年来最显著的发现之一,提示 Keller-Segel 模型在高趋化强度下的强非线性。

11.6 节"能量方法"则让我看到 1990 年代应用数学界向"梯度流"思想靠拢的趋势。Landau-Ginzburg 自由能 + 变分推导给出带 \(\nabla^4 n\)\(\nabla^2 n^3\) 项的 (11.59),这种推导与统计物理中 Allen-Cahn 方程、Cahn-Hilliard 方程的推导完全同构(作者也明确引用了 Cahn 1958, 1959, Cahn-Hilliard 1958, 1959 与 Huberman 1976)。这一节对做形态发生学(morphogenesis)研究的人来说是一座"语法金矿"——它告诉你 PDE 每一项都来自哪一项能量假设。但反过来,这种推导也带来自洽性问题:若 \(e(n)\) 形式改变,所有非线性项都会跟着改;模型对 \(e(n)\) 的具体形式相当敏感,这是能量方法的局限。

最后想提一下作者在 11.2 节末段关于"避免过早推广"的那段提醒(脚注 1 引 de Tocqueville:"推广无意义,因为上帝知道所有特例")。这与全书第 1 章反复强调的"模型服务于实验/观察"是一以贯之的——Murray 不止一次在章末提醒读者:把 \(\chi(a)\) 换成 log law 或 receptor law 之前,先用最简的 \(\chi(a)=\chi_0\) 形式与实验/数值对比;引入交叉扩散矩阵 \(D\) 之前,先确认对角常 \(D\) 的版本不充分。这种"先简单后复杂"的方法论纪律值得我记下来。

与上下章的衔接(一段话)

本章是第 I 卷"从 ODE 到 PDE 时空建模"的关键桥梁。其前第 10 章是 ODE 流行病学的集大成——SIR、HIV 病毒动力学、寄生虫免疫、牛结核元胞自动机等都仍属于集中参数或时间连续的常微分方程框架;本章则把舞台从时间扩展到空间:从最基础的随机游走 (11.1) 推导出扩散方程 (11.8),再把守恒-扩散框架与源项 \(f\) 整合得到反应扩散系统 (11.18),然后进一步引入趋化性 (11.28)、非局部效应 (11.41)、核矩展开与色散关系 (11.39)、以及能量方法导出的双调和广义扩散方程 (11.59)——这些是第 13 章(生物波:Fisher-KPP 方程、belousov-zhabotinsky 反应、螺旋波)的数学基础,也是第 II 卷第 1 章(连续介质形态发生、Turing 不稳定性)、第 5 章(粘菌与细菌趋化性的深入模型)、第 6 章(基于力学的形态发生)、第 12 章(神经场方程)的入门。Turing(1952)形态发生化学基础这一概念在本章 (11.18) 之后将延续到第 II 卷第 1 章深入讨论;扩散方程与 Fick 定律的等价性 (11.5)↔(11.7) 是后续所有空间模型的起点;从本章 (11.59) 中蕴含的双调和与 \(\nabla^2 n^3\) 非线性项出发,读者可以追溯到 Cahn-Hilliard 方程、相场模型(phase-field)等 1990s–2000s 凝聚态物理的活跃领域。下一章(第 12 章)将进入"振子产生的波现象"——把第 7 章的孤立振子放到二维空间,通过 FitzHugh-Nagumo 等可激发介质模型展示行波、靶波、螺旋波的几何与数学结构,与本章"扩散"主题在数学上共用 PDE 框架但物理上更强调"非平衡可激发性"这一新维度。