具有Michaelis-Menten饱和函数的趋化模型的Turing分岔分析

Turing Bifurcation Analysis of Chemotaxis Model with Michaelis-Menten Saturation Function

摘要: 本文研究了一类具有Michaelis-Menten饱和函数的趋化模型，以趋化性系数为参数，分析了系统在齐次Neumann边界条件下的Turing分岔行为。首先通过对系统在正平衡点处的特征方程进行讨论，得到了正平衡点的稳定性和Turing分岔的存在性。其次利用中心流形和正则形式理论，得到了Turing分岔的稳定性和分岔方向。最后，通过数值模拟验证了理论分析结果。

Abstract: In this paper, we investigate a chemotaxis models with Michaelis-Menten saturation functions subject to the homogeneous Neumann boundary condition. And discussed the Turing bifurcation by choosing the chemotaxis coefficient as the bifurcation parameter. The stability of the positive equilibrium and the existence of Turing bifurcation are obtained by the analysis of the corresponding characteristic equation. Moreover, we derive the stability and direction of the Turing bifurcation by using center manifold and normal form theory. Some numerical simulations are also carried out to illustrate the theoretical results.

1. 引言

趋化性指的是生物体对化学物质的空间梯度的响应。具体来说，生物体或细胞会向着浓度高的区域移动(正趋化性)，或远离浓度高的区域(负趋化性)。为了研究这种趋化性运动，Keller和Segel [1]-[3]通过建立如下抛物型的偏微分系统，描述了菌群密度与菌群生存环境中化学物质溶液浓度的动态关系，

$\{\begin{array}{l}{u}_t=\nabla \cdot \left(d_1\nabla u-\chi u\nabla v\right)+f\left(u,v\right),\\ v_t=d_2\Delta v+g\left(u,v\right).\end{array}$ (1)

其中，u是细胞的浓度，v是化学物质浓度，$d_1$ 和$d_2$ 分别表示细胞和化学物质的随机扩散率，$\chi$ 是趋化系数，它表示化学物质浓度的空间变化对细菌定向运动产生的刺激。$f\left(u,v\right)$ 和$g\left(u,v\right)$ 分别表示细菌繁殖与死亡，化学物质产生与消耗。

对于系统(1)，Chen [6]通过调整趋化参数的取值范围，得到了Turing不稳定性的充分条件，并以趋化系数为分岔参数，给出了Turing-Hopf分岔的发生条件。Kong [4]分析了Keller-Segel趋化模型的局部一维尖峰模式的存在性、线性稳定性。Wang [5]研究了具有齐次Neumman边界条件的一维密度抑制运动模型，通过将化学扩散速率视为分岔参数来表示解的单调性，从而得到了全局分岔图。Mi [9]建立了具有均匀时间界的经典解的存在性，分析了空间齐次共存稳态在一定参数条件下的局部稳定性和全局稳定性。

饱和效应是指当某种物质的浓度达到一定的阈值后，进一步增加该物质的浓度将不会显著提高反应速率。目前郭飞燕[7]研究了具有饱和效应的自催化反应扩散模型Hopf分支的存在性及稳定性。刘晓惠[8]研究了具有饱和效应的生化反应模型Hopf分支的存在性和稳定性并建立了稳态分支的存在性。Wang [10]研究了具有饱和效应的扩散Sel'kov模型的空间齐次Hopf分岔周期解的Turing不稳定性。以上研究结果表明，趋化机制和饱和效应都能够影响系统的动力学行为。但是当物质浓度持续增长，直至达到饱和阈值的过程中，趋化机制对系统动力学行为的影响是否会有所不同呢？本文在系统(1)的基础之上，建立了一个具有Michaelis-Menten饱和函数的趋化模型：

$\{\begin{array}{l}{u}_t=D\Delta u-\chi \nabla \cdot \left(u\nabla v\right)+\mu u\left(1-u\right),x\in \left(0,l\pi \right),t>0,\\ v_t=\Delta v+\frac{u}{u+\gamma }-v,x\in \left(0,l\pi \right),t>0,\\ u_x\left(0,t\right)=v_x\left(0,t\right)=u_x\left(l\pi ,t\right)=v_x\left(l\pi ,t\right)=0,t>0,\\ u\left(x,0\right)=u_0\left(x\right)\ge 0,v\left(x,0\right)=v_0\left(x\right)\ge 0,x\in \left(0,l\pi \right).\end{array}$ (2)

其中，$D,\chi >0$ 分别表示扩散速率与趋化速率，$\mu >0$ 是Logistic增长率，$\frac{u}{u+\gamma }$ 是饱和函数，$\gamma >0$ 是控

制化学引诱剂产生和降解速率的常数。本文将通过选择趋化性系数作为分岔参数，对Turing分岔的特征方程进行分析，来讨论趋化与饱和效应对系统Turing分岔存在性、稳定性等动力学行为的影响。

2. Turing分岔的存在性

不难得到模型(2)有两个常数解，$\left(0,0\right)$ 和$\left(1,\frac{1}{1+\gamma }\right)$ ，很容易可以判定$\left(0,0\right)$ 总是不稳定的。系统(2)在$\left(u^{*},v^{*}\right):=\left(1,\frac{1}{1+\gamma }\right)$ 处的线性化系统为

$\{\begin{array}{l}{u}_t=D\Delta u-\chi \Delta v-\mu u,\text{}\\ v_t=\Delta v+\frac{\gamma u}{{\left(1+\gamma \right)}^2}-v.\end{array}$ (3)

记算子$-\Delta$ 在齐次Neumann边界条件下的特征值为 ${\lambda }_k=\frac{k^2}{l^2}\left(k=0,1,2,\cdots \right)$ ，其对应的特征函数为

${\phi }_k\left(x\right)$ 。则系统(3)所对应的线性算子L为

$L=\left(\begin{array}{cc}-D\frac{k^2}{l^2}-\mu & \chi \frac{k^2}{l^2}\\ \frac{\gamma }{{\left(1+\gamma \right)}^2}& -\frac{k^2}{l^2}-1\end{array}\right).$

很容易求出其特征方程为${\lambda }^2\text{+}{T}_k\lambda +Q_k=0$ ，其中

$\begin{array}{l}{T}_k\text{=}\frac{k^2}{l^2}\left(D+1\right)+1+\mu >0,\\ Q_k=D{\left(\frac{k^2}{l^2}\right)}^2+\frac{k^2}{l^2}\left[D+\mu -\frac{\gamma \chi }{{\left(1+\gamma \right)}^2}\right]+\mu ,\end{array}$ (4)

不难得到，对于对应的常微分系统，即无扩散及趋化时，$\left(u^{*},v^{*}\right)$ 总是稳定的。然而在引入扩散项和趋化项后，$\left(u^{*},v^{*}\right)$ 变得不稳定了。这说明了趋化在Turing不稳定的出现上，起到了关键作用。将$\chi$ 作为分岔参数，不难得到发生Turing不稳定的临界值为

${\chi }_T\left(k\right)=\frac{{\left(1+\gamma \right)}^2\left[D+\mu +D\frac{k^2}{l^2}+\frac{l^2\mu }{k^2}\right]}{\gamma },$ (5)

记(5)为Turing分岔曲线$L_k$ 。根据特征方程可以计算出

${\frac{\text{d}\mathrm{Re}\lambda \left(\chi \right)}{\text{d}\chi }|}_{\chi ={\chi }_T\left(k\right)}=\frac{\frac{\gamma k^2}{l^2{\left(1+\gamma \right)}^2}}{\frac{k^2}{l^2}\left(D+1\right)+1+\mu }>0.$ (6)

因此，我们有如下结果：

定理2.1：存在$k_0$ ，使得${\chi }_{T\left(k_0\right)}=\underset{k\in Z^{+}}{\min }{\chi }_{T\left(k\right)}:={\chi }^{*}$ ，

1) 当$\chi <{\chi }^{*}$ 时，$\left(u^{*},v^{*}\right):=\left(1,\frac{1}{1+\gamma }\right)$ 总是稳定的。

2) 当$\chi ={\chi }^{*}$ 时，系统在$\left(u^{*},v^{*}\right)$ 处产生Turing分岔。

3. Turing不稳定的规范型

由(5)和(6)可知，当$\chi$ 横截穿过Turing分岔曲线$L_k,k\in Z^{+}$ 时，系统(2)在正平衡点$\left(u^{*},v^{*}\right)$ 处会发生Turing分岔。接下来，我们将通过中心流形理论和正则形式理论来探讨系统(2)的Turing分岔范式的新表达式，该表达式可以确定Turing分岔的性质。

选择$\chi$ 作为分岔参数，并定义以下Sobolev空间

${\rm X}=\left\{{\left(\phi \left(x\right),\psi \left(x\right)\right)}^{\text{T}}\in W^{2,2}\left(0,l\pi \right)|{\frac{\partial \phi }{\partial x}|}_{x=0,l\pi }={\frac{\partial \psi }{\partial x}|}_{x=0,l\pi }=0\right\}$ .

和内积

$\left[U_1,U_2\right]={\displaystyle {\int }_0^{l\pi }\left({\phi }_1{\phi }_2+{\psi }_1{\psi }_2\right)\text{d}x},$

其中$U_1={\left({\phi }_1,{\psi }_1\right)}^{\text{T}},U_2={\left({\phi }_2,{\psi }_2\right)}^{\text{T}}\in {\rm X}$ 。

令$\widehat{u}\left(\cdot ,t\right)=u\left(\cdot ,t\right)-u^{*},\widehat{v}\left(\cdot ,t\right)=v\left(\cdot ,t\right)-v^{*},\widehat{U}\left(t\right)=\left(\widehat{u}\left(\cdot ,t\right),\widehat{v}\left(\cdot ,t\right)\right),g\left(u\right)=\frac{u}{u+\gamma }$，然后去掉^，我们可以将系统重写为以下形式：

$\frac{\text{d}U\left(t\right)}{\text{d}t}=\left(\begin{array}{cc}D& -\chi \left(u+u^{*}\right)\\ 0& 1\end{array}\right)\Delta U+\left(\begin{array}{c}-\chi u_x{v}_x+\mu \left(u+u^{*}\right)\left(1-\left(u+u^{*}\right)\right)\\ g\left(u+u^{*}\right)-\left(v+v^{*}\right)\end{array}\right)$ . (7)

进一步，令$\epsilon =\chi -{\chi }^{\ast }$ ，这意味着$\epsilon =0$ 是分岔值，此时系统(2)等价于

$\frac{\text{d}U\left(t\right)}{\text{d}t}=d\Delta U+L_0\left(U\right)+f\left(U,\epsilon \right)$ (8)

其中

$\begin{array}{l}d=\left(\begin{array}{cc}D& -\chi u^{*}\\ 0& 1\end{array}\right),L_0\left(U\right)=\left(\begin{array}{c}\mu \left(1-2u^{*}\right)u\\ u{g}^{\prime }\left(u^{*}\right)-v\end{array}\right),\\ f\left(U,\epsilon \right)={\displaystyle \sum _{i+j\ge 2}\frac{1}{i!j!}{f}_{ij}{u}^i{v}^j+\frac{1}{2!}{f}_2^{\chi }}+{\displaystyle \sum _{i\ge 3}\frac{1}{i!}{f}_i^{\chi }},\end{array}$ (9)

且$f_{ij}=\left(\begin{array}{c}{f}_{ij}^{\left(1\right)}\\ f_{ij}^{\left(2\right)}\end{array}\right)$ ，$f_{ij}^{\left(s\right)}=\frac{{\partial }^{i+j}\tilde{f}\left(s\right)\left(0,0\right)}{\partial u^i\partial v^j},s=1,2$

$\begin{array}{l}\left(\begin{array}{c}\tilde{f}\left(1\right)\left(u,v\right)\\ \tilde{f}\left(2\right)\left(u,v\right)\end{array}\right)=\left(\begin{array}{c}0\\ g\left(u+u^{*}\right)-v^{*}-u{g}^{\prime }\left(u^{*}\right)\end{array}\right),\\ f_2^{\chi }=\left(\begin{array}{c}-\left(2{\chi }^{*}u{v}_{xx}+2u^{*}\epsilon v_{xx}+2{\chi }^{*}{u}_x{v}_x\right)\\ 0\end{array}\right),\\ f_2={\displaystyle \sum _{i+j=2}\frac{1}{i!j!}{f}_{ij}{u}^i{v}^j,}\\ f_i^{\chi }=\left(\begin{array}{c}A\\ 0\end{array}\right),i\ge 3,\end{array}$

其中

$\begin{array}{l}A=-[i{\chi }^{*}{\left(u^{*}\right)}^{i-1}{u}^{i-1}{v}_{xx}+i\left(i-1\right){\left(u^{*}\right)}^{i-2}\epsilon u^{i-2}{v}_{xx}\\ +i\left(i-1\right)\left(i-2\right){\left(u^{*}\right)}^{i-2}\epsilon u^{i-3}{u}_x{v}_x+i\left(i-1\right){\chi }^{*}{\left(u^{*}\right)}^{i-1}{u}^{i-2}{u}_x{v}_x].\end{array}$

令

$M_k=\left(\begin{array}{cc}-D\frac{k^2}{l^2}+\mu \left(1-2u^{*}\right)& {\chi }^{*}{u}^{*}\frac{k^2}{l^2}\\ g^{\prime }\left(u^{*}\right)& -\frac{k^2}{l^2}-1\end{array}\right)$ . (10)

对于奇点，我们有${\Lambda }_k=\left\{0\right\},B_k=0$ 。其中${\Lambda }_k$ 表示系统的特征值，$B_k$ 表示涉及到反应和扩散的线性算子。通过令$M_k{p}_k={\left(0,0\right)}^{\text{T}},M_k^{\text{T}}{q}_k={\left(0,0\right)}^{\text{T}}$ ，我们可以选择

$p_k=\left(\begin{array}{c}1\\ \frac{D\frac{k^2}{l^2}-\mu \left(1-2u^{*}\right)}{{\chi }^{*}{u}^{*}\frac{k^2}{l^2}}\end{array}\right),q_k=\left(\begin{array}{c}\frac{\frac{k^2}{l^2}+1}{T_k}\\ \frac{{\chi }^{*}{u}^{*}\frac{k^2}{l^2}}{T_k}\end{array}\right)$

使得$\langle q_k^{\text{T}},p_k\rangle =1$ 。

由

$\left(\begin{array}{c}u\\ v\end{array}\right)=p_{k}z{\phi }_k\left(x\right)+\omega ,z\in R,\omega ={\left({\omega }_1,{\omega }_2\right)}^{\text{T}}$ (11)

结合(9)和(11)，可以得到

$\begin{array}{c}\frac{1}{2}{f}_2^1\left(z,0,\epsilon \right)=\frac{1}{2}{f}_2^1\left(p_{k}z{\phi }_k\left(x\right),\epsilon \right)\\ =\frac{1}{2}{q}_k\left(\begin{array}{c}0\\ g^{″}\left(u^{*}\right)p_{k1}^2{z}^2{\phi }_k^2\left(x\right)\end{array}\right)\\ +\frac{1}{2}{q}_k\left(\begin{array}{c}2{\chi }^{*}\frac{k^2}{l^2}{p}_{k1}{p}_{k2}{\phi }_k^2\left(x\right)z^2+2\frac{k^2}{l^2}{p}_{k2}{u}^{*}{\phi }_k\left(z\right)\epsilon z-2{\chi }^{*}\frac{k^2}{l^2}{p}_{k1}{p}_{k2}{\xi }_k^2\left(x\right)z^2\\ 0\end{array}\right)\\ =\frac{1}{2}{q}_k\left(\begin{array}{c}2{\chi }^{*}\frac{k^2}{l^2}{p}_{k1}{p}_{k2}{\phi }_k^2\left(x\right)z^2+2\frac{k^2}{l^2}{p}_{k2}{u}^{*}{\phi }_k\left(z\right)\epsilon z-2{\chi }^{*}\frac{k^2}{l^2}{p}_{k1}{p}_{k2}{\xi }_k^2\left(x\right)z^2\\ g^{″}\left(u^{*}\right)p_{k1}^2{z}^2{\phi }_k^2\left(x\right)\end{array}\right),\end{array}$ (12)

其中${\xi }_k^2\left(x\right)=\frac{\sin \frac{k}{l}x}{{\Vert \cos \frac{k}{l}x\Vert }_{2,2}}$ 。

由于$B_k=0$ ，易得

$I_m{\left(M_2^1\right)}^c=span\left\{z^2,z\epsilon ,{\epsilon }^2\right\},I_m{\left(M_3^1\right)}^c=span\left\{z^3,z^2\epsilon ,z{\epsilon }^2,{\epsilon }^3\right\}$ (13)

因此，根据文献[11]，我们有

$\frac{1}{2}{g}_2^1\left(z,0,\epsilon \right)=\frac{1}{2}Pro{j}_{I_m{\left(M_2^1\right)}^c}{f}_2^1\left(z,0,\epsilon \right)=S_{k11}z\epsilon +S_{k21}{z}^2$ . (14)

由于${\displaystyle {\int }_0^{l\pi }{\phi }_k^2\left(x\right)\text{d}x}={\displaystyle {\int }_0^{l\pi }{\xi }_k^2\left(x\right){\phi }_k\left(x\right)\text{d}x}=0$ ，我们可以得到

$S_{k11}=u^{*}\frac{k^2}{l^2}{p}_{k2}{q}_{k1}{\displaystyle {\int }_0^{l\pi }{\phi }_k^2\left(x\right)\text{d}x}=u^{*}\frac{\frac{k^2}{l^2}\left(\frac{k^2}{l^2}+1\right)\left(D\frac{k^2}{l^2}-\mu \left(1-2u^{*}\right)\right)}{{\chi }^{*}{u}^{*}\frac{k^2}{l^2}{T}_k},$

$\begin{array}{l}{S}_{k21}=\left(u^{*}\frac{k^2}{l^2}{p}_{k1}{p}_{k2}{q}_{k1}+\frac{1}{2}{g}^{″}\left(u^{*}\right)p_{k1}^2{q}_{k2}\right){\displaystyle {\int }_0^{l\pi }{\phi }_k^3\left(x\right)\text{d}x}\\ -{\chi }^{*}{p}_{k1}\frac{k^2}{l^2}{p}_{k2}{q}_{k1}{\displaystyle {\int }_0^{l\pi }{\xi }_k^2\left(x\right){\phi }_k\left(x\right)\text{d}x}=0.\end{array}$ (15)

接下来，我们继续计算第三项：$\frac{1}{3!}{g}_3^1\left(z,0,\epsilon \right)$ 。因为 $B_k=0$ ，很容易得到${\left(Ker\left({\left(M_3^1\right)}^c\right)\right)}^c=0$ ，这意味

着$U_2^1=0$ 。由(14)式，我们很容易推导出$g_2^1\left(z,0,0\right)=2S_{k21}{z}^2=0$ 。因此由文献[11]能得到

$\frac{1}{3!}{g}_3^1\left(z,0,\epsilon \right)=\frac{1}{3!}Pro{j}_S{\tilde{f}}_3^1\left(z,0,0\right)+o\left(\epsilon {\left|\left(z,\epsilon \right)\right|}^2\right),$

其中$S=span\left\{z^3\right\}$ ，且

$\begin{array}{l}{\tilde{f}}_3^1\left(z,0,0\right)=f_3^1\left(z,0,0\right)+\frac{3}{2}[\left(D_z{f}_2^1\right)\left(z,0,0\right)U_2^1\left(z,0\right)+\left(D_{\omega }{f}_2\right)\left(z,0,0\right)U_2^2\left(z,0\right)\\ +\left(D_{\omega ,{\omega }_x,{\omega }_{xx}}{f}_2^{\chi }\right)\left(z,0,0\right){\left(U_2^2\left(z,0\right),U_{2x}^2\left(z,0\right),U_{2xx}^2\left(z,0\right)\right)}^{\text{T}}].\end{array}$ (16)

由(9)和(11)，我们可以得到

$\frac{1}{3!}{f}_3^1\left(z,0,0\right)=\frac{1}{3!}{f}_3^1\left(p_{k}z{\phi }_k\left(x\right),0\right)=\frac{1}{3!}{q}_k\left(\begin{array}{c}0\\ g^{‴}\left(u^{*}\right)p_{k1}^3{\phi }_k^3\left(x\right)z^3\end{array}\right),$ (17)

结合

${\displaystyle {\int }_0^{l\pi }{\phi }_k^4\left(x\right)\text{d}x}=\frac{3}{2l\pi }$ ，${\displaystyle {\int }_0^{l\pi }{\phi }_k^4\left(x\right)\text{d}x}=\frac{3}{2l\pi }$

可得

$\frac{1}{3!}Pro{j}_S{\tilde{f}}_3^1\left(z,0,0\right)=\frac{1}{4l\pi }{g}^{‴}\left(u^{*}\right)p_{k1}^3{q}_{k2}{z}^3$ . (18)

令

$U_2^2\left(z,0\right)\triangleq h\left(z\right)={\displaystyle \sum _{j\ge 0}{h}_{kj}^{\text{T}}}{z}^2\left(\begin{array}{c}{\beta }_j^1\\ {\beta }_j^2\end{array}\right)\in X^S$

和

$h_{kj}=\left(\begin{array}{c}{h}_{kj}^{\left(1\right)}\\ h_{kj}^{\left(2\right)}\end{array}\right),h_{kj}^{\left(1\right)},h_{kj}^{\left(2\right)}\in R$ ，

通过(9)，我们得到

$\begin{array}{l}{q}_k^{\text{T}}\left(\begin{array}{c}\left[D_k{f}_2\left(p_{k}z{\phi }_k\left(x\right),0,0\right)\left(h_{kj}^{\text{T}}{z}^2\left(\begin{array}{c}{\beta }_j^1\\ {\beta }_j^2\end{array}\right)\right),{\beta }_k^1\right]\\ \left[D_k{f}_2\left(p_{k}z{\phi }_k\left(x\right),0,0\right)\left(h_{kj}^{\text{T}}{z}^2\left(\begin{array}{c}{\beta }_j^1\\ {\beta }_j^2\end{array}\right)\right),{\beta }_k^2\right]\end{array}\right)\\ =2g^{″}\left(u^{*}\right)p_{k1}{q}_{k2}{h}_{kj}^{\left(1\right)}{z}^3{\displaystyle {\int }_0^{l\pi }{\phi }_k^2\left(x\right){\phi }_j\left(x\right)\text{d}x},\end{array}$ (19)

且

$\begin{array}{l}{q}_k^{\text{T}}\left(\begin{array}{c}\left[D_{\omega ,{\omega }_x,{\omega }_{xx}}{f}_2^{\chi }\left(p_{k}z{\phi }_k\left(x\right),0,0\right)\left(\begin{array}{c}\left(h_{kj}^{\text{T}}{z}^2\left(\begin{array}{c}{\beta }_j^1\\ {\beta }_j^2\end{array}\right)\right)\\ {\left(h_{kj}^{\text{T}}{z}^2\left(\begin{array}{c}{\beta }_j^1\\ {\beta }_j^2\end{array}\right)\right)}_x\\ {\left(h_{kj}^{\text{T}}{z}^2\left(\begin{array}{c}{\beta }_j^1\\ {\beta }_j^2\end{array}\right)\right)}_{xx}\end{array}\right),{\beta }_k^1\right]\\ \left[D_{\omega ,{\omega }_x,{\omega }_{xx}}{f}_2^{\chi }\left(p_{k}z{\phi }_k\left(x\right),0,0\right)\left(\begin{array}{c}\left(h_{kj}^{\text{T}}{z}^2\left(\begin{array}{c}{\beta }_j^1\\ {\beta }_j^2\end{array}\right)\right)\\ {\left(h_{kj}^{\text{T}}{z}^2\left(\begin{array}{c}{\beta }_j^1\\ {\beta }_j^2\end{array}\right)\right)}_x\\ {\left(h_{kj}^{\text{T}}{z}^2\left(\begin{array}{c}{\beta }_j^1\\ {\beta }_j^2\end{array}\right)\right)}_{xx}\end{array}\right),{\beta }_k^2\right]\end{array}\right)\\ =2q_{k1}{\chi }^{*}[\left(\frac{k^2}{l^2}{\displaystyle {\int }_0^{l\pi }{\phi }_k^2{\phi }_j\left(x\right)\text{d}x}-\frac{kj}{l^2}{\displaystyle {\int }_0^{l\pi }{\xi }_k\left(x\right){\zeta }_j\left(x\right){\phi }_k\left(x\right)\text{d}x}\right)p_{k2}{h}_{kj}^{\left(1\right)}\\ +\left(\frac{j^2}{l^2}{\displaystyle {\int }_0^{l\pi }{\phi }_k^2{\phi }_j\left(x\right)\text{d}x}-\frac{kj}{l^2}{\displaystyle {\int }_0^{l\pi }{\xi }_k\left(x\right){\zeta }_j\left(x\right){\phi }_k\left(x\right)\text{d}x}\right)p_{k1}{h}_{kj}^{\left(2\right)}]z^3.\end{array}$ (20)

由于

$c_{kj}={\displaystyle {\int }_0^{l\pi }{\phi }_k^2{\phi }_j\left(x\right)\text{d}x}=\{\begin{array}{l}\frac{1}{\sqrt{l\pi }}\text{, }j=0,k\ne 0,\\ \frac{1}{\sqrt{2l\pi }},j=2k\ne 0,\\ 0,其他,\end{array}$

和

${\tilde{c}}_{kj}={\displaystyle {\int }_0^{l\pi }{\xi }_k^2{\xi }_j\left(x\right){\phi }_k\left(x\right)\text{d}x}=\{\begin{array}{l}\frac{1}{\sqrt{2l\pi }},j=2k\ne 0,\\ 0,j\ne 2k.\end{array}$

联合(19)和(20)，我们可以得到

$\begin{array}{l}\frac{1}{3!}Pro{j}_S\left(D_{\omega }{f}_2\right)\left(z,0,0\right)\left(h\right)=\frac{1}{3}{g}^{″}\left(u^{*}\right)p_{k1}{q}_{k2}\left[\frac{1}{\sqrt{l\pi }}{h}_{k0}^{\left(1\right)}+\frac{1}{\sqrt{2l\pi }}{h}_{k2k}^{\left(1\right)}\right]z^3,\\ \frac{1}{3!}Pro{j}_S\left(D_{\omega ,{\omega }_x,{\omega }_{xx}}{f}_2^{\chi }\right)\left(z,0,0\right)\left(h\right)=\frac{1}{3!}{q}_{k1}{\chi }^{*}\left\{\frac{1}{\sqrt{l\pi }}\frac{k^2}{l^2}{p}_{k2}{h}_{k0}^{\left(1\right)}+\frac{1}{\sqrt{2l\pi }}\left[-\frac{k^2}{l^2}{p}_{k2}{h}_{k2k}^{\left(1\right)}+\frac{2k^2}{l^2}{p}_{k1}{h}_{k2k}^{\left(2\right)}\right]z^3\right\},\end{array}$

接下来，我们继续计算$h_{kj}^{\left(1\right)},h_{kj}^{\left(2\right)}$ 。结合文献[11]，我们有

$\left(\begin{array}{c}\left[f_2^2,{\beta }_j^1\right]\\ \left[f_2^2,{\beta }_j^2\right]\end{array}\right)=\left(\begin{array}{c}\left[M_2^2\left(h_{kj}^{\text{T}}{z}^2\left(\begin{array}{c}{\beta }_j^1\\ {\beta }_j^2\end{array}\right)\right),{\beta }_j^1\right]\\ \left[M_2^2\left(h_{kj}^{\text{T}}{z}^2\left(\begin{array}{c}{\beta }_j^1\\ {\beta }_j^2\end{array}\right)\right),{\beta }_j^2\right]\end{array}\right)=\frac{j^2}{l^2}\left(\begin{array}{cc}D& -{\chi }^{*}{u}^{*}\\ 0& 1\end{array}\right)h_{kj}{z}^2-L_0{h}_{kj}{z}^2$ (21)

和

$\left(\begin{array}{c}\left[f_2^2,{\beta }_j^1\right]\\ \left[f_2^2,{\beta }_j^2\right]\end{array}\right)=c_{kj}\left(\begin{array}{c}4{\chi }^{*}\frac{k^2}{l^2}{p}_{k1}{p}_{k2}\\ g^{″}\left(u^{*}\right)p_{k1}^2\end{array}\right)z^2.$ (22)

从(21)和(22)，我们可以推导出

$\begin{array}{l}{h}_{kj}=-c_{kj}{M}_j^{-1}\left(\begin{array}{c}0\\ g^{″}\left(u^{*}\right)p_{k1}^2\end{array}\right),j=0,\\ h_{kj}=-c_{kj}{M}_j^{-1}\left(\begin{array}{c}4{\chi }^{*}\frac{k^2}{l^2}{p}_{k1}{p}_{k2}\\ g^{″}\left(u^{*}\right)p_{k1}^2\end{array}\right),j=2k.\end{array}$ (23)

令

$\begin{array}{l}{S}_k=g^{‴}\left(u^{*}\right)p_{{}_{k1}}^3{q}_{k2},\\ S_{k,0}=\left(g^{″}\left(u^{*}\right)p_{k1}{q}_{k2}+q_{k1}{\chi }^{*}\frac{k^2}{l^2}{p}_{k2}\right)h_{k0}^{\left(1\right)},\\ S_{k,2k}=\left[g^{″}\left(u^{*}\right)p_{k1}{q}_{k2}-q_{k1}{\chi }^{*}\frac{k^2}{l^2}{p}_{k2}\right]h_{k2k}^{\left(1\right)}+q_{k1}{\chi }^{*}\frac{2k^2}{l^2}{p}_{k1}{h}_{k2k}^{\left(2\right)}.\end{array}$ (24)

结合上式，可以导出截断为三项的范式如下：

$\stackrel{\cdot }{z}=S_{k11}z\epsilon +S_{k30}{z}^3$ (25)

其中

$S_{k30}=\frac{1}{4l\pi }{S}_k+\frac{1}{\sqrt{2l\pi }}{S}_{k,0}+\frac{1}{2\sqrt{2l\pi }}{S}_{k,2k}$ . (26)

结合文献[12]可知，$S_{k30}$ 的符号决定了Turing分岔稳定性和分岔方向，即：

定理3.1：如果参数满足定理2.1中的条件(2)，则系统在$\chi ={\chi }^{\ast }$ 处产生Turing分岔且

1) 当$S_{k30}<0$ 时，在正平衡点$\left(u^{*},v^{*}\right)$ 附近发生的Turing分岔是超临界的，且分岔出的空间非齐次稳态解是稳定的。

2) 当$S_{k30}>0$ 时，在正平衡点$\left(u^{*},v^{*}\right)$ 附近发生的Turing分岔是次临界的，且分岔出的空间非齐次稳态解是不稳定的。

4. 数值模拟

在本节中，我们将通过数值模拟验证前几节中的理论结果。首先，我们选择一组参数值为

$\gamma =0.5,\mu =0.6,l=3,g^{\prime }\left(u^{*}\right)=\frac{\gamma }{{\left(u^{*}+\gamma \right)}^2},g^{″}\left(u^{*}\right)=-\frac{2\gamma }{{\left(u^{*}+\gamma \right)}^2},g^{‴}\left(u^{*}\right)=\frac{6\gamma }{{\left(u^{*}+\gamma \right)}^4}$ 。通过直接计算得到

$\left(u^{*},v^{*}\right)=\left(1,0.6667\right)$ 。此外，我们可以绘制出模型(2)的Turing分岔图，并找出$D-\chi$ 平面上的稳定区域，见图1。从图1中，我们可以看到，当$\chi$ 穿过Turing分岔曲线$L_k$ 时，系统(2)在正平衡点$\left(u^{*},v^{*}\right)$ 附近发生Turing分岔。

进一步，当$k=3$ 时，令$D=0.55$ 可以得到${\chi }_T\left(3\right)=10.3500$ 。结合(25)和(26)，可以得到$S_{311}=2.3000>0$ ，$S_{330}=0.1656>0$ 。根据定理3.1的理论分析结果可知，当$\chi ={\chi }_T\left(3\right)$ 时，系统在正平衡点$\left(u^{*},v^{*}\right)=\left(1,0.6667\right)$ 处发生次临界Turing分岔，见图2。

Figure 1. The Turing bifurcation diagram

图1. 图灵分岔图

Figure 2. Spatially inhomogeneous steady state solutions

图2. 空间非齐次稳态解

接下来，我们通过选择参数$\gamma =0.5,\mu =0.6,D=1.5,l=1$ ，来研究特征值$\lambda$ 与趋化系数$\chi$ 的关系。如图3所示，随着$\chi$ 的逐渐增大，特征值$\lambda$ 相应由负向正增大。这表明系统从稳定状态进入了不稳定状态。这种变化是Turing斑图形成的关键机制，系统将在不均匀性条件下发展出各种空间图案。因此我们可以通过调整趋化系数，控制图案的形成和性质，从而为相关应用研究提供理论依据。

Figure 3. Relationships between eigenvalues $\lambda$ and chemotaxis coefficients $\chi$

图3. 特征值$\lambda$ 的关系和趋化性系数$\chi$ 的关系图

5. 结论

本文研究了在齐次Neumann边界条件下具有Michaelis- Menten饱和函数的趋化模型。首先将趋化系数作为分岔参数，分析了模型(2)的平衡点的稳定性和Turing分岔的存在性。然后根据Turing分岔范式的新表达式，确定了Turing分岔的稳定性和方向。最后通过数值模拟发现系统在一定的参数范围内出现了不稳定的次临界解。这意味着系统的动力学行为在临界点附近发生了变化，但这种变化不是完全对称的，尽管系统产生了一个略微偏离平衡点的分岔，但并不会发生显著的改变或崩溃。随着参数的继续变化，系统可能会产生稳定的斑图。对Turing分岔和次临界解的深入研究，可以帮助我们更好地对不同生长条件下菌落模型的形态演变进行预测，从而为污染治理、生态恢复等问题提供理论依据。

NOTES

^*通讯作者。