理论生物学最基本的问题之一即是揭示现实世界模式的生成机制.通过在未搅拌的凝胶反应器中进行绿泥石-碘-丙二酸 (CIMA) 反应试验, Castets和Kepper等[1-2]第一次发现了Turing模式的实验证据. Turing[3]发现:耦合的反应扩散方程组可以用来描述生态系统的模式生成, 即在某些条件下, 反应扩散系统对应的常微分系统的稳定平衡态会随着扩散的产生而失稳, 并在空间形成一些有规律的斑图结构.这种不稳定性, 称之为Turing不稳定性 (Turing instability) 或扩散导致的不稳定性.
随后, Turing的发现引起了越来越多研究者的关注.为揭示现实世界多样的模式生成机制, 各种生化模型被相继提出, 比如: CIMA反应模型[4-5]、活化-抑制 (Activator-Inhibitor) 模型[6-9]、Schnakenberg模型[12-13]等.关于模式生成相关研究工作的进展, 详见文献[14].
通过方程逼近, Gierer和Meinhardt[8]推导出一个标准模型, 说明均匀分布的形态因子可导致模式生成.他们通过选择合适的参数区域, 构建了一些 (空间一维情形) 分子模型.文献[8]中的活化-抑制模型 (15a)-(15b) 可重新改写为如下模型
其中$u(x, t)$、$v(x, t)$分别表示活化剂和抑制剂的浓度. $k_1, \cdots, k_5, D_{u}, D_{v}$都是正常数. $k_1$表示活化剂的常数输入率, $k_{3}$ ($k_4$) 表示活化剂 (抑制剂) 的分布.活化剂和抑制剂的移除函数可用$k_2u, k_5v$表示, 或被酶降解, 或泄漏, 或由源重新吸收.抑制剂比活化剂扩散的快得多[8] (活化剂和抑制剂的扩散系数分别为$D_{u}, D_{v}$, $D_{u}<D_{v}$).活化剂可以激活活化剂源和抑制剂源, 但抑制剂只可抑制活化剂源.
令$L=1$.对模型 (1.1) 中参数进行无量纲化, Page等[9]提出了下述模型
文献[9]在参数空间不均匀的情形下, 讨论了空间模式的生成问题 (先前研究的Turing反应扩散模型, 大多在参数空间均匀情形下讨论), 并通过数值模拟与分析结果进行了比较.
对系统 (1.1) 中参数进行特殊的无量纲化.令$k_3=ak_1, k_2=bk_1, k_5=ck_4, k_4=rk_1, \overline{t}=k_1t, D_u=k_1D_1, D_v=k_1D_2$, 并略去$t$的上划线.于是, 系统 (1.1) 等价于
对应初边值条件
本文结构安排如下:第2节证明模型 (1.3) 对应ODE系统内部平衡态的线性稳定性, Hopf分歧 (极限环) 的存在性、稳定性及走向; 第3节说明模型 (1.3) 内部常数平衡态的Turing不稳定性, 及其附近Turing分歧、Hopf分歧的存在性; 第4节采用Matlab软件中"ode45、pdepe、bvp4c"函数对模型 (1.3) 进行数值模拟, 验证理论分析结果.
本节主要讨论ODE系统
内部平衡态的线性稳定性, 及其附近Hopf分歧的存在性、稳定性.
容易验证, 系统 (2.1) 没有边界平衡态, 且近存在惟一内部平衡态$U_{*}=(u_{*}, v_{*})$, 其中$u_{*}=\frac{1+ac}{b}, v_{*}=\frac{(1+ac)^{2}}{b^{2}c}$.简单计算知, 系统 (2.1) 在$U_{*}=(u_{*}, v_{*})$处的Jacobian矩阵为
方便起见, 令$h_{0}=\frac{2au_{*}}{v_{*}}-b$, $h=cr$, $h$作为主要参数.实际上, $r $是正真的参数, 因为$h$是$r$的线性函数.由$(u_{*}, v_{*})$定义知
因此, 系统 (2.1) 在$U_{*}$处的线性化特征方程为
显然, 方程 (2.3)有两根, 它们是
已知$\det J>0$.因此, 当${\rm trace}J<0$ (即$h>h_0$) 时, 方程 (2.3) 的两根的实部均小于零, 从而$U_{*}=(u_{*}, v_{*})$线性稳定.反之, 若$h<h_0$, 则$U_{*}$不稳定.
现讨论参数$h$跨过临界值$h_0$时, 系统 (2.1) 在$U_{*}=(u_{*}, v_{*})$处是否产生Hopf分歧.令$\beta=\sqrt{\det J}$.简单计算知, 当$h=h_0$时, 方程 (2.3) 有一对纯虚根$\pm {\rm i} \beta$.由Hopf分歧定理知, 要想系统 (2.1) 从$U_{*}=(u_{*}, v_{*})$分歧出小的周期解, 只需要横截条件 (transversality condition) 成立.将$h$作为分歧参数, 验证Hopf分歧定理中的横截条件.当$|h-h_0|$充分小时, 设$\lambda=p+{\rm i}q (p, q\in\mathbb{R})$是方程 (2.3) 的一个根.易知, 当$h=h_{0}$时, $\lambda={\rm i}\beta$.将$\lambda=p+{\rm i}q$代入方程 (2.3), 然后分离实部和虚部得
易知, 当$h=h_0$时, $p=0$.方程 (2.4) 两边对变量$h$求导得, $p'_{h}(h)|_{h=h_{0}}=-\frac{1}{2}<0$.根据Andronov-Hopf分歧定理[10-11], 我们得到下面的结论:
定理2.1 固定参数$a, b, c$, 设$ac>1$.
1) 若$h>h_{0}$, 则系统 (2.1) 内部平衡态$U_{*}=(u_{*}, v_{*})$渐近稳定; 若$h<h_{0}$, 则不稳定.
2) 当$h=h_{0}$时, 系统 (2.1) 在$U_{*}$处产生Hopf分歧.
注2.1 若$ac>1$, 则$f_{u}(u_{*}, v_{*})=\frac{2au_{*}}{v_{*}}-b\equiv h_{0}>0$.实际上, $\frac{2au_{*}}{v_{*}}-b=\frac{b(ac-1)}{1+ac}$.
现在分析定理2.1 ($ac>1$) 中Hopf分歧产生的周期解稳定性及走向.令$\overline{u}=u-u_{*}, \overline{v}=v-v_{*}$, 将内部平衡态$U_{*}$平移到原点.方便起见, 略去$u$和$v$的上划线, 则系统 (2.1) 转化为
将系统 (2.5) 在原点 (0, 0) 处进行三次泰勒展开, 系统 (2.5) 转化为
其中
且
已知$h=cr$.系统 (2.6) 在$U_{*}=(0, 0)$处线性化Jacobian矩阵为
由前文$(u_{*}, v_{*})$及$h_{0}$定义知线性化特征值为$\lambda_{1, 2}=\frac{{\rm trace}L\pm {\rm i}\sqrt{4\det L-({\rm trace}L)^{2}}}{2}\equiv p(h)\pm {\rm i}q(h)$, 其中
由注2.1知, $p(h_{0})=0$, $q(h_{0})=\sqrt{bh_{0}}>0$.令
其中$M=\frac{p(h)-a_1+b}{a_2}, N=-\frac{q(h)}{a_2}$.当$|h-h_{0}|$充分小时, $M, N>0$.令$(u, v)^{{\rm T}}=K(x, y)^{{\rm T}}$, 则系统 (2.6) 转化为
将系统 (2.7) 改写成以下极坐标形式
将系统 (2.8) 在$h=h_{0}$处泰勒展开
为确定前述周期解的稳定性, 需要计算系数$a(h_0)$的符号.
由$p(h)$表达式知, $\mu_{2}\equiv-\frac{a(h_{0})}{p'(h_{0})}=2a(h_{0})$.根据Andronov-Hopf分歧定理 (其中$a(h_{0})$决定分歧周期解的稳定性, $\mu_{2}$决定Hopf分歧的方向), 我们得到以下结论.
定理2.2 设$ac>1$.当$a(h_{0})<0$时, 定理2.1中的Hopf分歧是上临界的、渐近稳定的; 当$a(h_{0})>0$时, Hopf分歧是下临界的、不稳定的.
由上一节知, 当$h>h_{0}$时, 系统 (2.1) 的内部平衡态$U_{*}$渐近稳定, 而且当$h$跨越临界值$h_{0}$时, 内部平衡态$U_{*}$附近产生了小振幅的周期解.本节在$h>h_{0}$条件下, 讨论"扩散"效应对系统 (1.3)-(1.4) 内部平衡态稳定性的影响.
将空间区域限制在$(0, l\pi)$, 系统 (1.3)-(1.4) 等价于
令$(u, v)^{{\rm T}}=(\zeta, \tau)^{{\rm T}}{\rm e}^{\lambda t}\cos(kx)$.定义$J_{k}\equiv J-k^{2}{\rm diag}(D_{1}, D_{2})$, 其中$J$在第二节给出.系统 (3.1) 在内部常数平衡态$U_{*}$处的线性化特征方程为
当$h>h_{0}$时, ${\rm trace}J<0$.因此, ${\rm trace}J_{k}<0$, $k>0$. $U_{*}$不稳定, 只需找到某个$k>0$, 使得方程 (3.2) 至少存在一个正实根, 或只需$D_{1}h-D_{2}h_{0}<0$.假设$D_{1}h-D_{2}h_{0}<0$成立.方便起见, 记$H(k^{2})\equiv\det J_{k}$, 则$H(k^2)$在$k^2=k^2_{\min}=-\frac{D_{1}h-D_{2}h_{0}}{2D_1d_2}$处取得最小值.事实上, 存在$k>0$使得$H(k^2)<0$等价于$H(k^2_{\min})<0$.简单计算知
不等式$(3.4)\Rightarrow\frac{D_1}{D_2}>\frac{h_0}{h}$, 这与假设条件$D_{1}h-D_{2}h_{0}<0$矛盾.于是, 我们得到以下结论:
定理3.1 固定参数$a, b, c$, 设$h>h_{0}$ (此时, (2.1) 内部平衡态$U_{*}$渐近稳定), $D_{1}h-D_{2}h_{0}<0$.若
则系统 (3.1) 内部平衡态$U_{*}$不稳定, 即Turing不稳定发生.
注3.1 由定理3.1前分析过程知, 条件$D_{1}h-D_{2}h_{0}<0$, $h>h_{0}$是系统 (3.1) 内部平衡态$U_{*}$ Turing不稳定的必要条件.该条件隐含着$D_1<D_2$, 即抑制剂$v$的扩散速度较活化剂$u$快.
类似文献[15-16], 将$D_{2}$作为分歧参数, 讨论系统 (3.1) 内部平衡态$U_{*}$附近的Turing和Hopf分歧.已知${\rm trace}J=h_{0}-h, \det J=bh$, $k\equiv k(n)=n/l, n=0, 1, 2, \cdots$.方程 (3.2) 可改写为
或者
参数$h_{0}, h$, $k$已经在前文中给出, 并在后文中简记$k(n)=n/l$为$k$.
定理3.2 当$ac>1$时, 下面结论成立:
1) 设$h_{0}<\min\{h, \frac{D_{1}}{l^{2}}\}$.对任意的$D_{2}>0$, 系统 (3.1) 的内部平衡态$U_{*}$局部渐近稳定.
2) 设$\min\{h, \frac{D_{1}}{l^{2}}\}<h_{0}<\max\{h, \frac{D_{1}}{l^{2}}\}$.
a) 当$h<\frac{D_{1}}{l^{2}}$时, 对任意的$D_{2}>0$, 系统 (3.1) 的内部平衡态$U_{*}$不稳定.
b) (H1):对任意的$1\leq n\leq N_{1}$, 当$m\neq n, 1\leq m, n\leq N_{1}$时, $D_{2, m}^{S}\neq D_{2, n}^{S}$, 其中$D_{2, n}^{S}=\frac{h(D_{1}k^{2}+b)}{k^{2}(h_{0}-D_{1}k^{2})}$, $N_{1}$为满足$h_{0}-\frac{D_{1}n^{2}}{l^{2}}>0$的最大正整数.若$h>\frac{D_{1}}{l^{2}}$且条件 (H1) 成立, 则对任意的$1\leq n\leq N_{1}$, 系统 (3.1) 在$(D_{2, n}^{S}, U_{*})$处产生平衡态分歧.若进一步假设$D_{2}<\overline{D}_{2}\equiv\min_{1\leq n\leq N_{1}}D_{2, n}^{S}$, 则系统 (3.1) 的内部平衡态$U_{*}$局部渐近稳定.
3) 设$\max\{h, \frac{D_{1}}{l^{2}}\}<h_{0}<h+\frac{D_{1}}{l^{2}}$.若2) b) 中条件 (H1) 成立, 则对任意的$1\leq n\leq N_{1}$, 系统 (3.1) 在$(D_{2, n}^{S}, U_{*})$处产生平衡态分歧.若进一步假设$D_{2}<\overline{D}_{2}\equiv\min_{1\leq n\leq N_{1}}D_{2, n}^{S}$, 则系统 (3.1) 的内部平衡态$U_{*}$不稳定.
4) 设$h_{0}>h+\frac{D_{1}}{l^{2}}$.令$N_{2}$为满足$h_{0}-h-\frac{D_{1}n^{2}}{l^{2}}>0$的最大正整数.若$l^{2}<\frac{D_{1}(h_{0}+h)}{h_{0}(h_{0}-h)}$, 则对任意的$1\leq n\leq N_{2}$, 系统 (3.1) 在$(D_{2, n}^{H}, U_{*})$处产生Hopf分歧, 其中$D_{2, n}^{H}=\frac{h_{0}-h-D_{1}k^{2}}{k^{2}}$.
证 由注2.1及条件$ac>1$知, $h_{0}>0$.
1) 设$h_{0}<\min\{h, \frac{D_{1}}{l^{2}}\}$, 则$h_{0}<h$, $h_{0}<\frac{D_{1}}{l^{2}}$.对任意的$D_{2}>0$, 由 (3.7)-(3.8) 式知: ${\rm trace}J_{0}=h_{0}-h<0, \det J_{0}=bh>0$, 且${\rm trace}J_{n}<0, \det J_{n}>0$, $n\geq 1$.因此, $J_{k}$或$J_{n}$的所有特征值均具有负实部.即对任意的$D_{2}>0$, 系统 (3.1) 的内部平衡态$U_{*}$局部渐近稳定.
2) 假设$\min\{h, \frac{D_{1}}{l^{2}}\}<h_{0}<\max\{h, \frac{D_{1}}{l^{2}}\}$.
a) 若$h<\frac{D_{1}}{l^{2}}$, 则$h<h_{0}<\frac{D_{1}}{l^{2}}$.类似1), 对任意的$D_{2}>0$, ${\rm trace}J_{n}<0$, $n\geq 1$, 且$\det J_{n}<0$, $n\geq 0$,但${\rm trace}J_{0}=h_{0}-h>0$.因此, $J_{k}$或$J_{n}$存在一对具有正实部的特征值, 其余特征值具有负实部.即对任意的$D_{2}>0$, 系统 (3.1) 的内部平衡态$U_{*}$不稳定.
b) 若$h>\frac{D_{1}}{l^{2}}$, 则$\frac{D_{1}}{l^{2}}<h_{0}<h$.类似1), 对任意的$D_{2}>0$, ${\rm trace}J_{n}<0$, $n\geq 0$, $\det J_{0}>0$.由$N_{1}$和$D_{2, n}^{S}$定义知, $\det J_{n}(D_{2.n}^{S})=0$, $1\leq n\leq N_{1}$.假设$1\leq n\leq N_{1}$时, 由系统 (3.7) 知:当$1\leq m\leq N_{1}, m\neq n$时, $D_{2, m}^{S}\neq D_{2, n}^{S}$, 进而$\det J_{m}(D_{2.n}^{S})\neq0$; 当$m>N_{1}$时, $D_{1}k^{2}-h_{0}>0$, 进而$\det J_{m}(D_{2.n}^{S})>0$.已知$\frac{\rm d}{{\rm d}D_{2}}\det J_n(D_{2, n}^{S})=\left(\frac{n^{2}}{l^{2}}\right)\left[D_{1}\left(\frac{n^{2}}{l^{2}}\right)-h_{0}\right]<0$, $1\leq n\leq N_{1}$.于是, 当$1\leq n\leq N_{1}$时, 系统 (3.1) 在$(D_{2, n}^{S}, U_{*})$产生平衡态分歧.同上, 对任意$n\geq 1$及$D_{2}<\overline{D}_{2}$, $\det J_{n}>0$.因此, 当$D_{2}<\overline{D}_{2}$时, $J_{k}$或$J_{n}$的所有特征值均具有负实部.即对任意的$D_{2}>0$, 系统 (3.1) 的内部平衡态$U_{*}$局部渐近稳定.
3) 已知$\max\{h, \frac{D_{1}}{l^{2}}\}<h_{0}<h+\frac{D_{1}}{l^{2}}$, ${\rm trace}J_{n}=-D_{2}k^{2}+h_{0}-h-D_{1}k^{2}$.类似2) b), 系统 (3.1) 在$(D_{2, n}^{S}, U_{*})$处产生平衡态分歧.这里, 我们略去分歧存在性的证明过程, 只说明其不稳定性.类似2) b), $\det J_{n}>0$, $n\geq 0$, 且${\rm trace}J_{0}>0$, 但${\rm trace}J_{n}<0$, $n\geq 1$.于是, $J_{k}$或$J_{n}$存在一对具有正实部的特征值, 其余特征值具有负实部.即对任意的$D_{2}>0$, 系统 (3.1) 的内部平衡态$U_{*}$不稳定.
4) 设$h_{0}>h+\frac{D_{1}}{l^{2}}$成立.同上, 对任意的$D_{2}>0$, ${\rm trace}J_{0}>0, \det J_{0}=bh>0$.由$N_{2}$及$D_{2, n}^{H}$定义知, $D_{2, m}^{H}\neq D_{2, n}^{H}$, $1\leq m, n \leq N_{2}$, $m\neq n$.假设$1\leq m, n \leq N_{2}$.若$m=n$, 则${\rm trace}J_{m}(D_{2, n}^{H})=0$; 若$m\neq n$, 则${\rm trace}J_{m}(D_{2, n}^{H})\neq 0$.假设$1\leq n \leq N_{2}$.对任意$D_{2}>0$及$m>N_{2}$, ${\rm trace}J_{m}(D_{2, n}^{H})\neq 0$.现说明:对任意$1\leq n\leq N_{2}$及$m\geq 1$, $\det J_{m}(D_{2, n}^{H})>0$.事实上, 由$\frac{D_{1}}{h_{0}-h}<l^{2}<\frac{D_{1}(h_{0}+h)}{h_{0}(h_{0}-h)}$知, $D_{2, 1}^{H}h_{0}<D_{1}h$.于是, $D_{2, n}^{H}h_{0}\leq D_{2, 1}h_{0}<D_{1}h$, $1\leq n\leq N_{2}$.因此, 对任意的$m\geq 1$及$1\leq n\leq N_{2}$,
因此, 系统 (3.1) 在$(D_{2.n}^{H}, U_{*})$ ($1\leq n\leq N_{2}$) 处产生Hopf分歧, 只需要横截条件成立.实际上, $\frac{\rm d}{{\rm d}D_{2}}{\rm trace}J_n(D_{2, n}^{H})=-\frac{n^{2}}{2l^{2}}<0$, $1\leq n\leq N_{2}$, 横截条件成立.定理证毕.
本节借助Matlab软件和"ode45、pdepe、bvp4c"函数进行数值模拟, 分析和验证前几节的理论分析结果.首先考虑如下ODE系统
其中$r=0.1973$或$r=0.1132$, 初值$(u(0), v(0))=(1.4, 3.5)$.简单计算知, 系统 (4.1) 具有内部平衡态$U_{*}=(u_{*}, v_{*})=(1.9024, 3.9750)$, $h_{0}=0.1037$, $ac=1.1973>1$.
由引理2.1知, 若$h>h_{0}$ ($h<h_{0}$), 则$U_{*}$渐近稳定 (不稳定).另外, 当$h=h_{0}$时, 系统 (2.1) 在$U_{*}$处产生极限环.取$r=0.1973$, 图 1 (左边) 表明, 当$0.1796=h>h_{0}=0.1037$时, 内部平衡态$U_{*}$渐近稳定.
当$r=0.1132$时, 简单计算知, 系统 (2.10) 中$a(h_{0})=-2.3040<0$, $0.1031=h<h_{0}=0.1037$. 图 1 (右边) 表明, 在不稳定内部平衡态$U_{*}$附近, 系统 (4.1) 产生Hopf分歧 (极限环存在), 并且该Hopf分歧是上临界的、渐近稳定的.
ODE系统 (4.1) 对应的PDE系统为
初值选取$u(x, 0)=1.9024+0.02\sin(x), v(x, 0)=3.9750+0.02\cos(x)$, 其他参数选择如下
令$l\pi=50$.在条件 (4.3) 和 (4.4) 下, 不等式 (3.5) 成立.此时, $0.1796=h>h_{0}=0.1037$, $\frac{h_{0}}{h}=0.5774$.简单计算知$\frac{hh_0+2\det J-\sqrt{(hh_0+2\det J)^2-h^2h_0^2}}{h^2}=0.1024$.由定理3.1知:当$h>h_{0}$, $\frac{D_1}{D_2}<0.1024$时, $U_{*}$不稳定, 即Turing不稳定产生.
由条件 (4.5)-(4.6) 知, 作为系统 (2.1) 内部平衡态, $U_{*}$不稳定.从而, 作为系统 (3.1) 的内部平衡态, $U_{*}$亦不稳定.取$D_{2, 1}^{H}=0.0665$, 则定理3.2中条件4) 成立.简单计算知$N_{2}=1$.根据定理3.2, 系统 (3.1) 在$D_{2, 1}^{H}=0.0665$处产生Hopf分歧.由模拟图 2-3可知, 当$D_{2}>D_{2, 1}^{H}$或$D_{2}<D_{2, 1}^{H}$时, 系统 (3.1) 存在稳定的Hopf分歧正周期解.
由条件 (4.7), 取$D_{2, 1}^{S}=4.1855$, 则定理3.2中条件3) 成立.简单计算知, $N_{1}=1$.根据定理3.2, 系统 (3.1) 在$D_{2, 1}^{H}=375.4966$处产生Turing分歧.即当$D_{2}<D_{2, 1}^{S}$时, 系统 (3.1) 存在稳定的非常数平衡态正解.