http:∥xuebao.gxnu.edu.cn

数是使用 Tanimoto 系数及其化学结构的乘积[45]计算得到的);

⑥ 蛋白质序列相似性网络(基于基因组序列使用成对的史密斯-沃特曼得分[46]获得蛋白质相似性网

络,用于训练卷积网络以更新药物靶标特征)。

2.2 模型性能

现有的药物靶标预测主要将药物靶标的已知相互作用作为阳性实例,未知相互作用作为阴性实例。

实验中采用 10 倍交叉验证,并随机选择 10%的数据作为测试集,其余 90%的数据作为训练集。 实验将

NGDTI 与 NeoDTI

[31]

、DTINet

[22]

、BLMNI

[12]

、NetLapRLS

[47]和 HNM

[48] 等 5 种方法进行比较。 在实验中使

用 AUPR(精确召回曲线下的面积)来衡量 NGDTI 的预测效果。 从图 3( a)可知,NGDTI 比其他方法效果

更好,其中 AUPR 比最佳方法高 0.01。 与 DTINet 相比,尽管使用的都是降维之后的网络扩散特征,但

NGDTI 方法在此基础上进一步使用图卷积模型优化特征,从而获得更好的结果。 在正负样本比例设置为

1 ∶ 1 的情况下,实验结果显示 NeoDTI 比 DTINet 更糟糕。 为了验证 NGDTI 在稀疏阳性样本下的性能,实

验对样本数进行修改,并为阳性和阴性样本指定了 1 ∶ 10 的比例。 从图 3(b)可知,每种算法的性能都下

降了。 相比之下,NGDTI 仍然取得了最佳的预测性能,且比对比方法中最好的 AUPR 值多 0.015。 这表

明,即使在标记稀疏的情况下,对比方法的预测性能仍然不如 NGDTI 方法。

由于数据可能是冗余的,例如,数据集中存在一种蛋白质的多种同源蛋白质或一种药物的多种高度相

似的药物,这可能会对预测性能产生负面影响。 因此,本文采用与 Luo 等[22]相同的策略,通过在药物靶标

矩阵中删除那些具有相似药物或靶标的药物-靶标关联来减少数据冗余的影响。 实验还在消除药物-靶标

关联情况下进行测试,其中关联网络中的 Jaccard 相似度大于 0.6,药物化学相似性网络中的结构相似性得

分超过 0.6,蛋白质-蛋白质序列相似性网络中的同一性得分超过 0.4。 在这些实验中,阴性和阳性样本的

比例保持为 1 ∶ 1。 从图 3(c)(d)(e)(f)的实验结果来看,虽然冗余的药物靶标关联删除后,NGDTI 性能

下降,但其仍优于其他预测方法。

图 3 NGDTI 方法与其他方法的比较

Fig. 3 Comparison of NGDTI with other baseline methods

97

广西师范大学学报(自然科学版),2022,40(2)

2.3 模型参数的影响

本文提出的多网络集成的药物靶标预测算法 NGDTI 的核心是使用变分图自编码器(VGAE)来更新

药物和靶点特征。 与文献[49]中使用图卷积神经网络(GCN)的目的不同,NGDTI 是通过 GCN 学习药物

和蛋白质的网络特征,使得药物和蛋白质的特征编码含更丰富的生物信息,而不是通过 GCN 来学习每种

药物的分子结构。 VGAE 能够聚合邻域特征来进一步提高特征的可用性,它使用基于谱的图卷积网络

(GCN)方法从图信号处理的角度引入滤波器来定义图卷积,其中的图卷积操作可以被认为是从图信号中

去除噪声。 为了验证 VGAE 部件的有效性,实验中实现了一个不包含 VGAE 部件的多网络集成框架来预

测药物-靶标相互作用。 实验首先对 NGDTI 中是否有 VGAE 部件、药物特征维数、蛋白质特征维数进行不

同设置,并比较在不同条件下的预测效果,实验结果如表 1 所示。 从表 1 可以看出,在药物和蛋白质特征

维数相同的情况下,有 VGAE 部件的 NGDTI 方法预测效果更好,而且在药物特征维数为 100、蛋白质特征

维数为 400 时预测效果最好。

表 1 NGDTI 在不同设置下的预测性能(正负样本比例为 1 ∶ 1)

Tab. 1 Prediction performance of NGDTI under different settings (ratio of positive and negative samples is 1 ∶ 1)

是否有 VGAE 部件 药物特征维数 蛋白质特征维数 AUPR AUROC

否 100 200 0.889 0.863

是 100 200 0.901 0.880

否 200 200 0.894 0.875

是 200 200 0.914 0.895

否 100 400 0.924 0.904

是 100 400 0.943 0.901

否 200 400 0.921 0.900

是 200 400 0.928 0.910

在之后的实验中,主要评估参数的影响和 NGDTI 的鲁棒性。 这些实验通过更改与药物或靶标相关的

网络数量以及 NGDTI 的超参数来测试 NGDTI 的鲁棒性。 所有实验结果均通过多次实验取平均获得。

首先,通过实验验证聚合多个异构网络对预测结果的影响。 在实验中只使用部分网络的情况下进行

性能评估,同时将预测结果与使用所有网络的预测结果进行对比,结果如图 4(a)所示。 通过结果可以观

察到随着多数据源数据的整合,预测性能显著提高。 在添加了疾病和副作用相关信息的网络之后,模型的

预测效果也得到改进,这也表明整合多数据源数据的有效性(NGDTI 可以整合多种蛋白质或药物相关数

据来改善预测性能)。

图 4 整合更多与药物或靶标相关的信息的效果和重启随机游走概率 p 的影响

Fig. 4 Effect of integrating more information related to the drug or target and the effect of restarting the random walk probability

98

http:∥xuebao.gxnu.edu.cn

此外,本文还探索模型超参数对实验性能的影响。 在这里,主要研究随机游走重启概率 p 对实验结果

的影响。 如图 4(b)所示,在测试中改变重新启动概率值为 0.4 到 0.7,以观察不同概率下的性能稳定性。

在图 4(b)中可以看出,当重启概率为 0.4 至 0.7 时,NGDTI 实现了稳定的性能。 从以上实验可以得出,模

型的参数对实验性能的影响较小。

2.4 NGDTI 预测的药物靶标相互作用

最终的预测结果选取可信度排名前 10 位的药物-靶标相互作用,其中有 4 个药物-靶标相互作用有相

关文献研究的支持。 例如,nifedipine 是一种被批准用于辅助治疗高原肺水肿的药物,而 NR3C1 的多态性

与高原肺水肿有着重要的关联。 这一预测可以被先前的一项研究支持,该研究表明 NR3C1 多态性与高原

肺水肿的易感性有关[50]

。 此前有研究表明硝苯地平是一种可以抑制自发性心律失常的药物,而 SCN5A

在心律失常中起重要作用[51]

,这一关系也被 NGDTI 所预测。 此外,sorafenib 被批准用于治疗晚期肾细胞

癌,NGDTI 预测 sorafenib 与集落刺激因子受体(CSF1R)存在相互作用,已有研究也证实了 CSF1R 确实在

乳腺的发展和乳腺癌变中起着重要作用[52-53]

。 最后,NGDTI 预测 rivastigmine 和 CES1 的相互作用也在文

献[54]中得到支持。 简而言之,NGDTI 预测的药物-靶标相互关系有一部分是有文献研究支持的,这进一

步表明了 NGDTI 优秀的预测能力。

3 结语

本文提出一种名为 NGDTI 的模型,用于集成来自不同异构网络的信息来预测新的药物-靶标相互作

用。 NGDTI 可以通过网络扩散过程从异构网络中提取低维隐藏特征信息,为了对节点特征进行平滑和降

噪,该模型添加了图卷积编码来获得更加有效的节点特征。 从实验结果看,NGDTI 获得了比其他基准方

法更好的预测性能,而且 NGDTI 具有很强的鲁棒性。 此外,NGDTI 是可扩展的框架,其他有关药物和靶标

的更多信息也可以轻松地纳入到当前框架中。 因此,NGDTI 可以为加强药物开发和药物靶标预测提供有

用的工具。 后续将进一步优化 NGDTI 模型,整合更多异构信息,并改善模型的预测结果。 在本研究中,

NGDTI 模型仅用于预测未知的药物-靶标相互作用,但 NGDTI 模型也可以扩展应用于其他研究领域。

参 考 文 献

[1] 徐国保, 陈媛晓, 王骥. 基于图卷积网络的药物靶标关联预测算法[J]. 计算机应用, 2021, 41(5): 1522-1526. DOI:

10.11772 / j.issn.1001-9081.2020081186.

[2] 余冬华, 郭茂祖, 刘晓燕, 等. 药物靶标作用关系预测结果评价及查询验证[J]. 计算机研究与发展, 2019, 56(9):

1881-1888. DOI: 10.7544 / issn1000-1239.2019.20180830.

[3] 宫莉, 康经武. 基于化学蛋白质组学的药物靶标鉴定[J]. 色谱, 2020, 38(8): 877-879.

[4] 李梢. 网络靶标: 中药方剂网络药理学研究的一个切入点[J]. 中国中药杂志, 2011, 36(15): 2017-2020.

[5] PAUL S M, MYTELKA D S, DUNWIDDIE C T, et al. How to improve R&D productivity: the pharmaceutical industry’ s

grand challenge[J]. Nature Reviews Drug Discovery, 2010, 9(3): 203-214. DOI: 10.1038 / nrd3078.

[6] 彭利红, 刘海燕, 任日丽, 等. 基于多标记学习预测药物-靶标相互作用[ J]. 计算机工程与应用, 2017, 53(15):

260-265.

[7] MORRIS G M, HUEY R, LINDSTROM W, et al. AutoDock4 and AutoDockTools4: automated docking with selective

receptor flexibility[J]. Journal of Computational Chemistry, 2009, 30(16): 2785-2791. DOI: 10.1002 / jcc.21256.

[8] KEISER M J, ROTH B L, ARMBRUSTER B N, et al. Relating protein pharmacology by ligand chemistry[ J]. Nature

Biotechnology, 2007, 25(2): 197-206. DOI: 10.1038 / nbt1284.

99

广西师范大学学报(自然科学版),2022,40(2)

[9] KLIPP E, WADE R C, KUMMER U. Biochemical network-based drug-target prediction[J]. Current opinion in

biotechnology, 2010, 21(4): 511-516. DOI: 10.1016 / j.copbio.2010.05.004.

[10] YAMANISHI Y, ARAKI M, GUTTERIDGE A, et al. Prediction of drug-target interaction networks from the integration of

chemical and genomic spaces[J]. Bioinformatics, 2008, 24(13): i232-i240. DOI: 10.1093 / bioinformatics/ btn162.

[11] BLEAKLEY K, YAMANISHI Y. Supervised prediction of drug-target interactions using bipartite local models [ J ].

Bioinformatics, 2009, 25(18): 2397-2403. DOI: 10.1093 / bioinformatics/ btp433.

[12] MEI J P, KWOH C K, YANG P, et al. Drug-target interaction prediction by learning from local information and neighbors

[J]. Bioinformatics, 2013, 29(2): 238-245. DOI: 10.1093 / bioinformatics/ bts670.

[13] VAN LAARHOVEN T, NABUURS S B, MARCHIORI E. Gaussian interaction profile kernels for predicting drug-target

interaction[J]. Bioinformatics, 2011, 27(21): 3036-3043. DOI: 10.1093 / bioinformatics/ btr500.

[14] PAHIKKALA T, AIROLA A, PIETILÄ S, et al. Toward more realistic drug-target interaction predictions[ J]. Briefings in

Bioinformatics, 2015, 16(2): 325-337. DOI: 10.1093 / bib / bbu010.

[15] HAO M, WANG Y L, BRYANT S H. Improved prediction of drug-target interactions using regularized least squares

integrating with kernel fusion technique[J]. Analytica Chimica Acta, 2016, 909: 41-50. DOI: 10.1016 / j.aca.2016.01.014.

[16] NASCIMENTO A C A, PRUDÊNCIO R B C, COSTA I G. A multiple kernel learning algorithm for drug-target interaction

prediction[J]. BMC Bioinformatics, 2016, 17: 46. DOI: 10.1186 / s12859-016-0890-3.

[17] LIU Y, WU M, MIAO C Y, et al. Neighborhood regularized logistic matrix factorization for drug-target interaction prediction

[J]. PLOS Computational Biology, 2016, 12(2): e1004760. DOI: 10.1371 / journal.pcbi.1004760.

[18] GÖNEN M. Predicting drug-target interactions from chemical and genomic kernels using Bayesian matrix factorization[ J].

Bioinformatics, 2012, 28(18): 2304-2310. DOI: 10.1093 / bioinformatics/ bts360.

[19] ZHENG X D, DING H, MAMITSUKA H, et al. Collaborative matrix factorization with multiple similarities for predicting

drug-target interactions[C]∥ Proceedings of the 19th ACM SIGKDD International Conference on Knowledge Discovery and

Data Mining. New York, NY: ACM Press, 2013: 1025-1033. DOI: 10.1145 / 2487575.2487670.

[20] HAO M, BRYANT S H, WANG Y L. Predicting drug-target interactions by dual-network integrated logistic matrix

factorization[J]. Scientific Reports, 2017, 7: 40376. DOI: 10.1038 / srep40376.

[21] MIZUTANI S, PAUWELS E, STOVEN V, et al. Relating drug-protein interaction network with drug side effects [ J].

Bioinformatics, 2012, 28(18): i522-i528. DOI: 10.1093 / bioinformatics/ bts383.

[22] LUO Y N, ZHAO X B, ZHOU J T, et al. A network integration approach for drug-target interaction prediction and

computational drug repositioning from heterogeneous information [ J]. Nature Communications, 2017, 8: 537. DOI: 10.

1038 / s41467-017-00680-8.

[23] NATARAJAN N, DHILLON I S. Inductive matrix completion for predicting gene-disease associations[ J]. Bioinformatics,

2014, 30(12): i60-i68. DOI: 10.1093 / bioinformatics/ btu269.

[24] CHEN Q F, LAI D H, LAN W, et al. ILDMSF: inferring associations between long non-coding RNA and disease based on

multi-similarity fusion[J]. IEEE/ ACM Transactions on Computational Biology and Bioinformatics, 2021, 18( 3): 1106-

1112. DOI: 10.1109 / TCBB.2019.2936476.

[25] LAN W, LAI D H, CHEN Q F, et al. LDICDL: LncRNA-disease association identification based on collaborative deep

learning[J/ OL]. IEEE/ ACM Transactions on Computational Biology and Bioinformatics, 2020 [ 2021-07-23]. https:∥

ieeexplore.ieee.org / document / 9246263. DOI: 10.1109 / TCBB.2020.3034910.

[26] 李琳, 梁永全, 刘广明. 基于重启随机游走的图自编码器[ J]. 计算机应用研究, 2021, 38(10): 3009-3013. DOI:

10.19734 / j.issn.1001-3695.2021.03.0083.

[27] 翟正利, 梁振明, 周炜, 等. 变分自编码器模型综述[ J]. 计算机工程与应用, 2019, 55(3): 1-9. DOI: 10.3778 / j.

issn.1002-8331.1810-0284.

[28] 郭景峰, 董慧, 张庭玮, 等. 主题关注网络的表示学习[J]. 计算机应用, 2020, 40(2): 441-447. DOI: 10.11772 / j.

issn.1001-9081.2019081529.

[29] 王佩琪, 高原, 刘振宇, 等. 深度卷积神经网络的数据表示方法分析与实践[J]. 计算机研究与发展, 2017, 54(6):

100

http:∥xuebao.gxnu.edu.cn

1348-1356. DOI: 10.7544 / issn1000-1239.2017.20170098.

[30] 刘思, 刘海, 陈启买, 等. 基于网络表示学习与随机游走的链路预测算法[ J]. 计算机应用, 2017, 37(8): 2234-

2239. DOI: 10.11772 / j.issn.1001-9081.2017.08.2234.

[31] WAN F P, HONG L X, XIAO A, et al. NeoDTI: neural integration of neighbor information from a heterogeneous network for

discovering new drug-target interactions[J]. Bioinformatics, 2019, 35(1): 104-111. DOI: 10.1093/ bioinformatics/ bty543.

[32] GILMER J, SCHOENHOLZ S S, RILEY P F, et al. Neural message passing for quantum chemistry[C]∥ Proceedings of

the 34th International Conference on Machine Learning: Volume 70. Sydney: PMLR, 2017: 1263-1272.

[33] ROSSI R A, ZHOU R, AHMED N K. Deep inductive graph representation learning[ J]. IEEE Transactions on Knowledge

and Data Engineering, 2020, 32(3): 438-452. DOI: 10.1109 / TKDE.2018.2878247.

[34] WU Z H, PAN S R, CHEN F W, et al. A comprehensive survey on graph neural networks[ J]. IEEE Transactions on

Neural Networks and Learning Systems, 2021, 32(1): 4-24. DOI: 10.1109 / TNNLS.2020.2978386.

[35] 焦李成, 杨淑媛, 刘芳, 等. 神经网络七十年: 回顾与展望[J]. 计算机学报, 2016, 39(8): 1697-1716.

[36] SUN C, CAO Y K, WEI J M, et al. Autoencoder-based drug-target interaction prediction by preserving the consistency of

chemical properties and functions of drugs[J]. Bioinformatics, 2021, 37(20): 3618-3625. DOI: 10.1093 / bioinformatics/

btab384.

[37] XUAN P, ZHANG Y, CUI H, et al. Integrating multi-scale neighbouring topologies and cross-modal similarities for drugprotein interaction prediction[J]. Briefings in Bioinformatics, 2021, 22(5): bbab119. DOI: 10.1093 / bib / bbab119.

[38] GAO K Y, FOKOUE A, LUO H, et al. Interpretable drug target prediction using deep neural representation [ C]∥

Proceedings of the Twenty-Seventh International Joint Conference on Artificial Intelligence. Marina del Rey, CA: IJCAI,

2018: 3371-3377. DOI: 10.24963 / ijcai.2018 / 468.

[39] WANG S, CHO H H, ZHAI C X, et al. Exploiting ontology graph for predicting sparsely annotated gene function[ J].

Bioinformatics, 2015, 31(12): i357-i364. DOI: 10.1093 / bioinformatics/ btv260.

[40] KIPF T N, WELLING M. Variational graph auto-encoders[EB/ OL]. (2016-11-21)[2021-07-23]. https:∥arxiv.org / abs/

1611.07308.

[41] KNOX C, LAW V, JEWISON T, et al. DrugBank 3.0: a comprehensive resource for ‘Omics’ research on drugs[ J].

Nucleic Acids Research, 2011, 39(suppl_1): D1035-D1041. DOI: 10.1093 / nar/ gkq1126.

[42] PRASAD T S K, GOEL R, KANDASAMY K, et al. Human protein reference database: 2009 update[ J]. Nucleic Acids

Research, 2009, 37(suppl_1): D767-D772. DOI: 10.1093 / nar/ gkn892.

[43] DAVISA P, MURPHY C G, JOHNSON R, et al. The comparative toxicogenomics database: update 2013[ J]. Nucleic

Acids Research, 2013, 41(D1): D1104-D1114. DOI: 10.1093 / nar/ gks994.

[44] KUHN M, CAMPILLOS M, LETUNIC I, et al. A side effect resource to capture phenotypic effects of drugs[J]. Molecular

Systems Biology, 2010, 6(1): 343. DOI: 10.1038 / msb.2009.98.

[45] ROGERS D, HAHN M. Extended-connectivity fingerprints[ J]. Journal ofChemical Information and Modeling, 2010, 50

(5): 742-754. DOI: 10.1021 / ci100050t.

[46] SMITH T F, WATERMAN M S. Identification of common molecular subsequences[J]. Journal of Molecular Biology, 1981,

147(1): 195-197. DOI: 10.1016 / 0022-2836(81)90087-5.

[47] XIA Z, WU L Y, ZHOU X B, et al. Semi-supervised drug-protein interaction prediction from heterogeneous biological

spaces[J]. BMC Systems Biology, 2010, 4(Suppl 2): S6. DOI: 10.1186 / 1752-0509-4-S2-S6.

[48] WANG W H, YANG S, ZHANG X, et al. Drug repositioning by integrating target information through a heterogeneous

network model[J]. Bioinformatics, 2014, 30(20): 2923-2930. DOI: 10.1093 / bioinformatics/ btu403.

[49] DUVENAUD D, MACLAURIN D, AGUILERA-IPARRAGUIRRE J, et al. Convolutional networks on graphs for learning

molecular fingerprints[EB/ OL]. (2015-11-03)[2021-07-23]. https:∥arxiv.org / abs/ 1509.09292v2.

[50] YANG Y Z, DU H, LI Y H, et al. NR3C1 gene polymorphisms are associated with high-altitude pulmonary edema in Han

Chinese[J]. Journal of Physiological Anthropology, 2019, 38: 4. DOI: 10.1186 / s40101-019-0194-1.

[51] LIU G X, REMME C A, BOUKENS B J, et al. Overexpression of SCN5A in mouse heart mimics human syndrome of

101

广西师范大学学报(自然科学版),2022,40(2)

enhanced atrioventricular nodal conduction[ J]. Heart Rhythm, 2015, 12( 5): 1036-1045. DOI: 10.1016 / j. hrthm.2015.

01.029.

[52] THOMAS G, GURUNG I S, KILLEEN M J, et al. Effects of L-type Ca

2+

channel antagonism on ventricular

arrhythmogenesis in murine hearts containing a modification in the Scn5a gene modelling human long QT syndrome 3[ J].

The Journal of Physiology, 2007, 578(1): 85-97. DOI: 10.1113 / jphysiol.2006.121921.

[53] ULLRICH K, WURSTER K D, LAMPRECHT B, et al. BAY 43-9006 / Sorafenib blocks CSF1R activity and induces

apoptosis in various classical Hodgkin lymphoma cell lines[ J]. British Journal of Haematology, 2011, 155(3): 398-402.

DOI: 10.1111 / j.1365-2141.2011.08685.x.

[54] TSURKAN L G, HATFIELD M J, EDWARDS C C, et al. Inhibition of human carboxylesterases hCE1 and hiCE by

cholinesterase inhibitors[J]. Chemico-biological interactions, 2013, 203(1): 226-230. DOI: 10.1016 / j.cbi.2012.10.018.

Study on Multi-information Integration for Drug Target Prediction

TAN Kai

1

, LI Yongjie

1

, PAN Haiming

1

, HUANG Kexin

2

, QIU Jie

2

, CHEN Qingfeng

1∗

(1. School of Computer, Electronics and Information, Guangxi University, Nanning Guangxi 530004, China;

2. Guangxi Medical University, Nanning Guangxi 530021, China;

3. School of Computer Science and Engineering, Yulin Normal University, Yulin Guangxi 537000, China)

Abstract: Accurate determination of drug-target interactions is crucial in drug discovery process and

repositioning. Traditional methods for DTI prediction are either time-consuming ( simulation-based methods) or

heavily dependent on domain expertise ( similarity-based and feature-based methods). Existing computationbased methods using single data information or sparse data, always suffer from high false positive rates. Although

integrating multiple heterogeneous networks has been prevalent for drug target prediction, how to retain as much

structural information as possible is still a big challenge. This paper proposes a novel framework NGDTI, which

extracts relevant biological properties and association information from the network while maintaining the topology

information. Further, the graph neural network is applied to update the extracted feature information. The learned

topology-preserving representations of drugs and targets promote DTI prediction. Compared with the state-of-theart methods, NGDTI increases the AUPR value by nearly 0.01. The results demonstrate that NGDTI is promising

for drug development and repositioning.

Keywords: drug target association prediction; network embedding; network integration; matrix decomposition;

graph neural network

(责任编辑 黄 勇)

102

第 40 卷 第 2 期

2022 年 3 月

广西师范大学学报(自然科学版)

Journal of Guangxi Normal University (Natural Science Edition)

Vol. 40 No. 2

Mar. 2022

DOI: 10.16088 / j.issn.1001-6600.2021052801 http: xuebao.gxnu.edu.cn

徐王军, 曹进德, 伍代勇, 等.一类具有迁移和 Allee 效应的食饵-捕食者系统稳定性[J]. 广西师范大学学报(自然科学版), 2022, 40(2):

103-115. XU W J, CAO J D, WU D Y, et al. Stability of a prey-predator model with migration and Allee effects [ J]. Journal of Guangxi Normal

University(Natural Science Edition), 2022, 40(2): 103-115.

一类具有迁移和 Allee 效应的食饵-捕食者

系统稳定性

徐王军1

, 曹进德2

, 伍代勇1

, 申传胜1∗

(1. 安庆师范大学 数理学院, 安徽 安庆 246133; 2. 东南大学 数学学院, 江苏 南京 211189)

摘 要: 研究一类食饵具有 Allee 效应且捕食者具有人工控制迁移的食饵-捕食者系统, 该系统具有平方根项的功能性反

应函数。 首先通过定性分析, 证明解的有界性, 分析平衡点的存在性, 得到系统平衡点的局部稳定性的充分条件。 接着

讨论平衡点的 Hopf 分岔存在性, 并通过计算第一李雅普诺夫系数, 研究平衡点 Hopf 分岔的稳定性和方向。 最后通过数

值模拟验证所得结论的正确性, 结果表明 Allee 效应和人工控制迁移率对食饵种群和捕食者种群的生存与灭绝具有重要

意义。

关键词: 食饵-捕食者模型; Allee 效应; 迁移率; Hopf 分岔; 稳定性

中图分类号: Q141; O175 文献标志码: A 文章编号: 1001-6600(2022)02-0103-13

在生态系统中,任何种群都不是独立存在的,而是与生物群落中其他种群相互联系、互相共存。 同时,

种群之间一般存在着种内竞争、种间竞争、捕食与被捕食及互惠等关系。 其中,关于捕食与被捕食系统,人

们研究了多种功能性反应函数,如 Holling 型功能性反应函数[1-6]

、Beddington-DeAngelis 型功能性反应函

数[7]

,以及带有平方根项的功能性反应函数[8]

。 这里平方根项是基于食饵集体防御行为,即食饵种群常

常抱团形成一个圆形防御区,处于圆形区域最外围的个体数 R(x)将成为被捕食对象。 注意到 R(x)与食

饵种群圆形防御区域的周长 l 成比例,即 R(x)≈l。 假设该圆形防御区域的面积为 S, 则有 l = 2 πS 。一

般地,食饵种群密度 x, 即单位表面积内食饵的个体数与其所占据的区域面积 S 成正比, 即 x∝S。 因此,

有 R(x)∝ x 。 然而目前关于平方根项功能性反应函数的研究鲜有报道。 英国人口学家 Malthus(1766—

1834)根据人口统计资料,于 1798 年提出人口指数增长模型[9]

,该模型假设单位时间内人口的增长量与

当年的人口数量成正比。 后来生态学家们又提出了密度制约的 Logistic 增长模型[10]

。 实际上,当种群数

量低于某一个阈值时,由于环境的作用,种群密度往往呈现负增长,例如自然界中一些密度很小的珍稀动

物(大熊猫、金丝猴、华南虎),如果人类不予保护,顺其自然,即使不加捕杀也会灭绝,这种现象被称为

Allee 效应。 近年来,Allee 效应[11-16]受到广泛关注,例如 Petrovskii 等[17] 研究了具有 Allee 效应的生物入

侵机制的食饵-捕食者系统,Rao 等[18]讨论了具有 Allee 效应的食饵-捕食者扩散模型。

近年来,农作物病虫害绿色防控方法与措施受到高度关注,尤其是生物防治技术受到生态学家们的普

遍欢迎,如利用释放禽类等害虫天敌来防控病虫害。 Barclay

[19] 提出利用捕食者释放、栖息地管理和农药

释放组合防治害虫的模型。 Tang 等[20]研究了虫害综合治理的最佳时机:天敌释放和杀虫剂施用率的模

拟。 成定平[21]建立鼠类与天敌系统渐近稳定性的数学分析。 Chen 等[22]讨论了具有捕食者迁移的时滞捕

食者-食饵模型,该模型中捕食者的迁移率依赖于有效食饵,当捕食者捕获率较低时,不能有效控制食饵的

数量,从而导致食饵急剧泛滥,故通过人为控制迁移捕食者达到控制捕食者的种群密度。 最近,文献[23]

收稿日期: 2021-05-28 修回日期: 2021-06-28

基金项目: 国家自然科学基金(11975025); 国家自然科学基金委员会与英国皇家学会合作交流项目(12011530158)

通信作者: 申传胜(1975—), 男, 安徽六安人, 安庆师范大学教授, 博士。 E-mail: csshen@mail.ustc.edu.cn

广西师范大学学报(自然科学版),2022,40(2)

在三种群的食饵-捕食者模型中引入食饵的时变迁移速率,发现迁移诱导出现 2 个 Hopf 分叉和 2 个极限

环,有趣的是,第 2 个 Hopf 分叉后的灭绝点出现类似于不连续的一级相变。 关于其他密度依赖的迁移对

食饵-捕食者模型稳定性的影响研究参见文献[24-27]。

文献[22]建立的模型没有考虑 Allee 效应以及群防御行为。 实际上,当食饵的数量较大时,由于被捕

食者群防御行为,捕食者捕获到的食饵往往不能呈线性增长。 其次,捕获率还会受到其他因素影响,如时

间、环境因素等。 为了更贴近实际生态意义,本文基于文献[22]建立一类同时具有 Allee 效应和人工控制

迁移率的 2 种群食饵-捕食者模型,分析了模型解的有界性、平衡点的存在性、局部稳定性以及 Hopf 分岔。

此外,通过数值模拟验证了以上动力学性质的正确性。

1 模型的建立

经典两物种食饵-捕食者模型为

dx

dt

= h0(x)-R(x)y,

dy

dt

= -uy+eR(x)y。

ì

î

í

ï

ï

ï

ï

(1)

式中:x 和 y 分别表示食饵种群的密度和捕食者种群的密度;h0(x)是遵循 Logistic 增长的函数,即 h0(x)=

rx(1-x / K),这里 r 和 K 分别表示食饵的内禀自然增长率和环境容纳量;u 和 e 分别表示捕食者的死亡率

和捕食者捕获食饵后的转化率;R(x)= α x ,α 表示捕食者的捕获率,R(x)描述了在食饵种群中观察到的

群体防御行为[8]

。

对模型(1)中第 1 个方程引入 Allee 效应项, 得到

h(x)= rx 1-

x

K

( ) 1-

A+c

x+c

( ) ,

式中 A 和 c 分别表示 Allee 效应阈值和辅助参数,这里 0<A<K。 显然,当食饵种群密度小于 Allee 阈值时,

食饵呈现负增长,可见 Allee 效应对食饵种群的增长具有很大的影响。

对模型(1)中第 2 个方程引入人工控制迁移函数

P(x,y)= ε(x-ρy),

式中:ε 表示捕食者的迁移率,ρ 表示单位时间内捕食者捕获食饵的消耗率。 因此,得到如下模型:

dx

dt

= h(x)-α x y,

dy

dt

= -uy+eα x y+P(x,y)。

ì

î

í

ï

ï

ï

ï

(2)

为了简便起见,对模型(2)引入无量纲变换:

?t = rt,x?=

x

K

,y?=

y

K

,α?=

α K

r

,c?=

c

K

,?e = e,u?=

u

r

,ε?=

ε

r

,ρ?= ρ,θ =

A

K

。

式中 0<θ<1。 于是,去掉“短横”,得到系统

dx

dt

=

x

x+c

(1-x)(x-θ)-α x y,

dy

dt

= -uy+eα x y+ε(x-ρy)。

ì

î

í

ï

ï

ï

ï

(3)

考虑到生物学的实际意义,系统(3)的可行域为 E ={(x,y) x≥0,y≥0} 。 系统中的所有参数都是正常数。

104

http:∥xuebao.gxnu.edu.cn

2 系统解的有界性

定理 1 系统(3)初值满足(x0 ,y0 )∈E 的解非负且有界。

证明 对于任何初始值(x0 ,y0 )∈E,并且对任意 t>0,由系统(3)可知:

x(t) = x(0)exp ∫

t

0

(1 - x(s))(x(s) - θ)

x(s) + c

-

αy(s)

x(s)

é

ë

ê

ê

ù

û

ú

ú

ds ( ) ≥ 0,

y(t) = y(0)exp ∫

t

0

- u + eα x(s) +

εx(s)

y(s)

- ερ

é

ë

ê

ê

ù

û

ú

ú ( ds) ≥ 0。

ì

î

í

ï

ï

ï

ï

因此,对任意 t>0,都有 x(t)≥0,y(t)≥0。 系统(3)的解非负。

当 x∈(0,θ)∪(1,+∞ )时,恒有

dx

dt

<0,即有lim

t→∞

x(t)= 0,而当 x∈(θ,1)时,有

dx

dt

≤

x(1-x)(x-θ)

x+c

≤(1-x)(x-θ)。 (4)

考虑方程

dz

dt

= (1-z)(z-θ),

由方程解得lim sup

t→∞

z(t)≤1。 通过比较定理,得lim sup

t→∞

x(t)≤1。

定义一个关于解的和函数 L(t)= x(t)+

y(t)

e

,其中 e>0,那么

dL

dt

=

x

x+c

(1-x)(x-θ)-

u

e

y+

ε

e

(x-ρy)≤(1-x)(x-θ)-

u+ερ

e

( ) y+

ε

e

x =

-x

2+(θ+1)x-θ+

ε

e

x-(u+ερ)(L-x)= -x

2+(θ+

ε

e

+u+ερ+1)x-θ-(u+ερ)L。 (5)

记

g(x)= -x

2+(θ+

ε

e

+u+ερ+1)x-θ,

则有

dg

dx

= -2x+θ+

ε

e

+u+ερ+1。 (6)

令

dg

dx

= 0,解得

xa

=

(θ+u+ρε+1)e+ε

2e

,

且

d

2

g

dx

2

= -2<0。

(ⅰ)当 xa<1 时,g(x)在 x = xa 取得最大值。

g(xa )=

(θe+ue+ερe+e+ε)

2

4e

2

- θ,

即有

dL

dt

≤g(xa )-(u+ερ)L,

从而

105

广西师范大学学报(自然科学版),2022,40(2)

lim sup

t→∞

L(t)≤

g(xa )

u+ερ

。 (7)

(ⅱ)当 xa≥1 时,g(x)在 x = 1 取得最大值

g(1)=

ε+ue+ρεe

e

,

即有

dL

dt

≤g(1)-(u+ερ)L,

从而

lim sup

t→∞

L(t)≤

g(1)

u+ερ

。 (8)

综合式(7)和式(8),可得

lim sup

t→∞

x(t)+

y(t)

e

( ) ≤

g(xa )

u+ερ

, xa<1,

g(1)

u+ερ

, xa≥1。

ì

î

í

ï

ï

ï

ï

(9)

于是结论成立。 证毕。

3 平衡点分析

3.1 平衡点的存在性

确定系统(3)的非负平衡点,令(x

∗

,y

∗

)为其平衡点,则:

x

∗

(1-x

∗

)(x

∗ -θ)-α(x

∗ +c) x

∗

y

∗ = 0,

-uy

∗ +eα x

∗

y

∗ +ε(x

∗ -ρy

∗

)= 0。 { (10)

1)系统(3)恒存在灭绝平衡点 E0(0,0)。

2)当 ε = 0,系统(3)有 2 个轴向平衡点 E11(1,0)和 E12( θ,0)以及正平衡点 E1( x

∗

1 ,y

∗

1 ),记 x

∗

1

= η,

y

∗

1

= ξ。 当 η∈(θ,1),由方程组(10)可得

η =

u

eα

( )

2

,ξ =

u(1-η)(η-θ)

eα

2

(η+c)

。 (11)

3)当 ε≠0,系统(3)有正平衡点 E2(x

∗

,y

∗

),即:

y

∗ =

x

∗

(1-x

∗

)(x

∗ -θ)

α(x

∗ +c)

,

(eα x

∗ -u-ρε)y

∗ +εx

∗ = 0。

ì

î

í

ï

ï

ï

ï

(12)

从而

x

∗

(1-x

∗

)(x

∗ -θ)

α(x

∗ +c)

=

-εx

∗

eα x

∗ -u-ρε

。 (13)

记 δ = x

∗

,由式(13)得

eαδ

5-(u+ρε)δ

4-α(e(θ+1)+ε)δ

3+(u+ρε)(θ+1)δ

2-α(εc-eθ)δ-(u+ρε)θ = 0。 (14)

再令

M(δ)= a0

δ

5+a1

δ

4+a2

δ

3+a3

δ

2+a4

δ+a5 , (15)

式中:

a0

= eα>0,a1

= -(u+ρε)<0,a2

= -α(e(θ+1)+ε)<0,a3

= (u+ρε)(θ+1)>0,

106

http:∥xuebao.gxnu.edu.cn

a4

= -α(εc-eθ),a5

= -(u+ρε)θ<0。

若 M(δ)= 0 的 5 个根 δ 是实数,由笛卡尔符号法则可知,则 M(δ)= 0 正根的个数等于它的系数序列的变

号数,即 M(δ)= 0 有 5 个正实数根,但是它的系数序列的变号数为 3,这与结论相矛盾。 因此,M(δ)= 0 至

多存在 3 个正实数根满足 δ = x

∗

>0 且 x

∗ ∈( θ,1),即系统(3)至多存在 3 个正平衡点 E2( x

∗

,y

∗

),同

时,因为 a0>0 且 a5<0,所以 M(δ)的系数序列的变号数至少存在一个变号数,即系统(3)至少存在一个正

平衡点 E2(x

∗

,y

∗

)。

3.2 稳定性分析

对系统(3)的所有非负平衡点进行局部稳定性分析,还要注意,对任意 x∈(0,θ],恒有 dx / dt<0,随着

时间的不断演化,食饵最终会趋于灭绝(x = 0),由于食饵灭绝,捕食者的增长率呈现负增长,即有dy / dt<0,

捕食者也趋于灭绝(y = 0),因此,灭绝平衡点 E0 是稳定的。 至于系统(3)的其他平衡点的局部稳定性,可

以通过计算其 Jacobian 矩阵的特征值,并运用 Lyapunov 第一种方法和 Routh-Hurwitz 准则确定各平衡点局

部渐近稳定的充分条件。

令

dx

dt

=

x

x+c

(1-x)(x-θ)-α x y:=Q(x,y),

dy

dt

= -uy+eα x y+ε(x-ρy):=W(x,y)。

ì

î

í

ï

ï

ï

ï

(16)

于是,有

∂Q

∂x

=

(2θ-3c+1)x

2-2x

3+2c(θ+1)x-cθ

(x+c)

2

-

αy

2 x

,

∂Q

∂y

= -α x ,

∂W

∂x

=

eαy

2 x

+ε,

∂W

∂y

= -u+eα x -ερ。

ì

î

í

ï

ïï

ï

ï

(17)

(ⅰ)系统(3)在轴向平衡点 E11处的 Jacobian 矩阵为

J11

=

θ-1

c+1

-α

0 eα-u

( ) 。 (18)

因此,轴向平衡点 E11对应的特征方程为

λθ-1

c+1

( ) (λ-eα+u)= 0, (19)

它的特征值为 λ1

=

θ-1

c+1

,λ2

= eα-u。 显然,λ1<0。 若 u>eα,即 λ2<0,则 E11局部渐近稳定,此时 E11为稳定结

点;若 u<eα,即 λ2>0,则 E11为鞍点。

(ⅱ)系统(3)在轴向平衡点 E12处的 Jacobian 矩阵为

J12

=

θ(1-θ)

θ+c

-α θ

0 eα θ -u

( )

。 (20)

因此,轴向平衡点 E12对应的特征方程为

λθ(1-θ)

θ+c

( ) (λ-eα θ +u)= 0, (21)

它的特征值为 λ1

=

θ(1-θ)

c+1

,λ2

= eα θ -u。 显然,λ1>0。 若 θ<

u

eα

( )

2

,即 λ2<0,则 E12是鞍点;若 θ>

u

eα

( )

2

,即

λ2>0,则 E12是不稳定结点。

(ⅲ)系统(3)在正平衡点 E1 处的 Jacobian 矩阵为

107

广西师范大学学报(自然科学版),2022,40(2)

J1

=

a11 a12

a21 0 ( ) 。 (22)

式中:

a11

=

(θ-3c+1)η

2-2η

3+2c(θ+1)η-cθ

(η+c)

2

-

αξ

2 η

,a12

= -

u

e

,a21

=

e

2

α

2

ξ

2u

。 (23)

因此,E1 对应的特征方程为

λ

2-a11λ-a12 a21

= 0, (24)

设其特征值为 λi(i = 1,2),若

θ<

3η

3+(5c-1)η

2-3cη

η

2+(3c+1)η-c

, (25)

则有

λ1

+λ2

= a11<0,

λ1λ2

= -a12 a21>0。 { (26)

从而它的特征值 λi<0(i = 1,2),于是可得 E1 是局部渐近稳定的。

(ⅳ)系统(3)在正平衡点 E2 的 Jacobian 矩阵为

J2

=

b11

b12

b21

b22

( ) , (27)

式中:

b11

=

(2θ-3c+1)x

∗ 2-2x

∗ 3+2c(θ+1)x

∗ -cθ

(x

∗ +c)

2

-

αy

∗

2 x

∗

,

b12

= -α x

∗

,

b21

=

e1αy

∗

2 x

∗

+ε,

b22

= eα x

∗ -u-ερ。

ì

î

í

ï

ï

ï

ï

ï

ï

ï

ï

ïï

(28)

于是,它所对应的特征方程为

λ

2+ϑ1λ+ϑ2

= 0, (29)

式中:

ϑ1

= -(b11

+b22 ),

ϑ2

= b11

b22

-b12

b21 。 { (30)

根据 Routh-Hurwitz 判据准则,有 Δ1

=ϑ1>0,Δ2

=

ϑ1 1

0 ϑ2

=ϑ1ϑ2>0。 如果 ϑ1>0 和 ϑ2>0 成立,那么 E2 是

局部渐近稳定的。

3.3 Hopf 分岔

以 Allee 效应阈值 θ 为分岔参数,讨论在正平衡点 E1 处发生 Hopf 分岔的条件。 设 λ( θ) = λR(θ) +

iλI(θ)为特征方程(24)的特征值,代入特征方程(24),然后将其实部和虚部分离,得到:

λ

2

R

-λ

2

I

-a11λR

-a12 a21

= 0,

2λRλI

-a11λI

= 0。 { (31)

当特征值出现实部为零的复共轭时,正平衡点 E1 发生 Hopf 分岔失去其稳定性。 此时,Hopf 分岔点定义为

θH ,则有 λR(θH )= 0,将其代入式(31),得到:

-λ

2

I

-a12 a21

= 0,

a11λI

= 0(λI≠0)。 { (32)

108

http:∥xuebao.gxnu.edu.cn

从方程(32)解得 a11(θH )= 0,而 λI(θH )= -a12(θH )a21(θH ) >0,即有

det(J1 ) θ=θH

= -a12(θH )a21(θH )>0。

联合式(23),得到

θH

=

αξ η (η+c)

2+4η

4-4cη

2+2(3c-1)η

3

2η(η

2+2cη-c)

。 (33)

接着对方程(31)两边关于 θ 求导,注意 λR(θ)= 0,于是有

-a11

dλR

dθ

-2λI

dλI

dθ

=

d(a12 a21 )

dθ

,

2λI

dλR

dθ

-a11

dλI

dθ

=λI

d(a11 )

dθ

。

ì

î

í

ï

ïï

ï

ï

(34)

解得

dλR

dθ θ=θH

=

2λ

2

I da11

dθ

-

a11 d(a12 a21 )

dθ

a

2

11

+4λ

2

I θ=θH

≠0 ⇔

2λ

2

I da11

dθ θ=θH

≠

a11 d(a12 a21 )

dθ θ=θH

。 (35)

定理 2 如果det(J1 ) θ=θH

= -a12(θH )a21(θH ) >0 且

2λ

2

I da11

dθ θ=θH

≠

a11 d(a12 a21 )

dθ θ=θH

成立,那么当 θ<θH

时,系统(3)的正平衡点 E1 局部渐近稳定;并且 E1 在 θ = θH 处发生 Hopf 分岔。

定理 3 定义 l

1 为

l

1 := 3sxxx

s

2

y

-(sxy

+2sxxy

syη)⌊sxx

sy

+3sxxx

s

2

yη- -syhx

sy

sxxy

ξ」, (36)

式中:

sy

= a12

= -α η ,sxx

=

-2η

3-6cη

2-6c

2

η+2c

2

θ+2c

2+2cθ

(η+c)

3

+

αξ

4η η

,sxy

=

-α

2 η

,

sxxx

=

6c

2

η-6c

2

θ-6c

2-6cθ

(η+c)

4

-

3αξ

8η

2

η

,sxxy

=

α

4η η

,hx

= a21

=

e

2

α

2

ξ

2u

。

ì

î

í

ï

ïï

ï

ïï

(37)

若 l

1>0,则正平衡点 E1 的 Hopf 分岔是亚临界的;若 l

1<0,则正平衡点 E1 的 Hopf 分岔是超临界的。

证明 为了研究 Hopf 分岔的稳定性和方向,下面计算第一李雅普诺夫系数,令 u = x-η,v = y-ξ,那么

系统(3)转变为

du

dt

=

u+η

u+η+c

(1-u-η)(u+η-θ)-α u+η (ξ+v):= s(u,v),

dv

dt

= -u(ξ+v)+eα u+η (ξ+v)+ε(u+η-ρ(ξ+v)):= h(u,v)。

ì

î

í

ï

ï

ï

ï

(38)

现在对方程(38)在(u,v)= (0,0)考虑其三阶泰勒展开式,得到:

du

dt

= a11 u+a12

v+s1(u,v),

dv

dt

= a21 u+a22

v+h1(u,v)。

ì

î

í

ï

ï

ï

ï

(39)

这里 s1(u,v)和 h1(u,v)是关于 u 和 v 的高阶项,则有

s1(u,v)= suu u

2+suvuv+svv

v

2+suuu u

3+suuvu

2

v+suvvuv

2+svvv

v

3

,

h1(u,v)= huu u

2+huvuv+hvv

v

2+huuu u

3+huuvu

2

v+huvvuv

2+hvvv

v

3

。 { (40)

式中:

109

广西师范大学学报(自然科学版),2022,40(2)

su

= a11

=

(θ-3c+1)η

2-2η

3+2c(θ+1)η-cθ

(η+c)

2

-

αξ

2 η

,sv

= a12

= -α η ,

suu

=

-2η

3-6cη

2-6c

2

η+2c

2

θ+2c

2+2cθ

(η+c)

3

+

αξ

4η η

,suv

=

-α

2 η

,

svv

= 0,suuu

=

6c

2

η-6c

2

θ-6c

2-6cθ

(η+c)

4

-

3αξ

8η

2

η

,suuv

=

α

4η η

,suvv

= 0,svvv

= 0,

hu

= a21

=

e

2

α

2

ξ

2u

,hv

= a22

= 0,huu

= -

eαξ

4η η

,huv

=

eα

2 η

,

hvv

= 0,huuu

=

3eαξ

8η

2

η

,huuv

=

-eα

4η η

,huvv

= 0,hvvv

= 0。

ì

î

í

ï

ï

ï

ï

ï

ï

ïï

ï

ï

ï

ï

ï

ï

ï

(41)

方程(40)中所有偏导数都是在分岔点计算的,即(u,v)= (0,0)。 因此系统(39)可以写成

Ṗ = JE1

P+D(P), (42)

式中:

P = (u,v)

T

,D = (s1(u,v),h1(u,v))

T

,s1(u,v)= suu u

2+suvuv+suuu u

3+suuvu

2

v,

h1(u,v)= huu u

2+huvuv+huuu u

3+huuvu

2

v。

当 su

= a11

= 0 时,正平衡点 E1 在 Hopf 分岔点的特征值是纯虚数,即 λ = i -a12 a21

= i -svhu 。 再令该

特征值对应的特征向量为 u?= ( u1 ,u2 )

T

。 记正平衡点 E1 对应的 Jacobian 矩阵为 A =

0 sv

hu 0 ( ) ,则有 Au?=

λu?,解得 u?= (sv,i -svhu )

T

。

定义

F= (Re(u?)-Im(u?))=

sv 0

0 - -svhu

é

ë

ê

ê

ù

û

ú

ú

,

令 P =FZ,其中 Z = (z1 ,z2 )

T

,得

u

v

( ) =

sv 0

0 - -svhu

é

ë

ê

ê

ù

û

ú

ú

z1

z2

é

ë

ê

ê

ù

û

ú

ú

,

从而

u = sv

z1 ,

v = - -svhu

z2 。 {

因此,通过变换得到如下系统

Ż = (F

-1

J1F)Z+F

-1D(FZ),

即

ż

1

ż

2

é

ë

ê

ê

ù

û

ú

ú

=

0 - -svhu

-svhu 0

é

ë

ê

ê

ê

ù

û

ú

ú

ú

z1

z2

é

ë

ê

ê

ù

û

ú

ú

+

D1(z1 ,z2 )

D2(z1 ,z2 )

é

ë

ê

ê

ù

û

ú

ú

, (43)

式中:

D1(z1 ,z2 )=

1

sv

[suu (sv

z1 )

2- -svhu

suv

sv

z1

z2

+suuu(sv

z1 )

3-suuv(sv

z1 )

2 -svhu

z2 ],

D2(z1 ,z2 )= -

1

-svhu

huu(sv

z1 )

2- -svhu huv

sv

z1

z2

+huuu(sv

z1 )

3- -svhu huuv(sv

z1 )

2

z2

[ ] 。

ì

î

í

ï

ïï

ï

ï

(44)

Hopf 分岔的方向由第一李雅普诺夫系数的符号决定,于是有

l

1 :=

1

16

∂

3D1

∂z

3

1

+

∂

3D2

∂z

2

1 ∂z2

é

ë

ê

ê

ù

û

ú

ú

+

1

16 -svhu

∂

2D1

∂z1 ∂z2

∂

2D1

∂z

2

1

-

∂

2D2

∂z1 ∂z2

∂

2D2

∂z

2

1

-

∂

2D1

∂z

2

1

∂

2D2

∂z

2

1

é

ë

ê

ê

ù

û

ú

ú

, (45)

110

http:∥xuebao.gxnu.edu.cn

式中:

∂

3D1

∂z

3

1

= 6suuu

s

2

v ,

∂

3D1

∂z1 ∂z

2

2

= 0,

∂

3D2

∂z

2

1 ∂z2

= 2huuv

s

2

v ,

∂

3D2

∂z

3

2

= 0,

∂

2D1

∂z1 ∂z2

= - -svhu

suv

-2 -svhu

suuvu,

∂

2D2

∂z1 ∂z2

= svhuv

+2huuv

svu,

∂

2D1

∂z

2

1

= 2suu

sv

+6suuu

svu+2vsv

suuv,

∂

2D1

∂z

2

2

= 0,

∂

2D2

∂z

2

1

= -

2huu

s

2

v

-svhu

-

6huuu

s

2

v u

-svhu

-

2huuv

s

2

v

v

-svhu

,

∂

2D2

∂z

2

2

= 0。

ì

î

í

ï

ï

ï

ï

ï

ï

ï

ï

ï

ï

ï

ï

(46)

通过简化得到 l

1 的表达式为

l

1 := 3suuu

s

2

v

-(suv

+2suuv

svη)[suu

sv

+3suuu

s

2

v η- -svhu

sv

suuv

ξ]。 (47)

若 l

1>0,则正平衡点 E1 的 Hopf 分岔是亚临界的;若 l

1<0,则正平衡点 E1 的 Hopf 分岔是超临界的。

4 数值模拟

下面,通过选择不同的 Allee 效应阈值 θ 和人工控制迁移参数 ε,采用数值模拟来验证各个平衡点存

在条件及其稳定性。

在图 1 中,固定参数 r = 0.9,c = 1.5。 红色曲线和蓝色曲线指食饵种群的环境容纳量为 15 和 25。 当食

饵种群不受 Allee 效应的影响(θ = 0)且环境容纳量较大时,h0(x)的最大值也较大(见图 1(a));当食饵种

群受到 Allee 效应的影响时,食饵种群在 0<x<5 呈现负增长,在 5<x<15 呈现正增长。 与没有 Allee 效应相

比,当食饵种群的环境容纳量相同时,食饵种群受到 Allee 效应影响后的增长率会减小(见图 1( b))。 因

此,Allee 效应对种群的存活和灭绝具有非常重要的作用。

(a)没有 Allee 效应时 (b)具有 Allee 效应时

图 1 参数 r = 0.9,c = 1.5 时食饵种群的增长率 h0(x)和 h(x)

Fig. 1 Growth rate h0(x) and h(x) of the prey population when r = 0.9,c = 1.5

在图 2 中,固定参数 α= 0.24,u = 0.05,e = 0.32,c = 0.01,根据方程(11)可得 η = 0.423 85。 当 θ = 0.04 满

足条件(25)时,食饵种群和捕食者种群的相图曲线趋于一个稳定点 E1(0.423 85,1.382 7),此时,E1 是局

部渐近稳定的(如图 2(a)); 当 θ = 0.075不满足条件(25)时,食饵种群和捕食者种群的相图曲线趋于一个

稳定的极限环(如图 2(b))。 因此,从图 2(a)和图 2(b)可以看出,当 Allee 阈值 θ 较小时,则系统趋于正

平衡点的可能性较大。

111

广西师范大学学报(自然科学版),2022,40(2)

图 2 系统(3)在 ε= 0 时的相图

Fig. 2 Phase diagram of system (3) at ε= 0

在图 3 中,假设系统(3)没有人为控制迁移。 固定参数 α = 0.24,u = 0.05,e = 0.32,c = 0.01, 根据方程

(11)解得 η = 0.423 85,再代入方程(33),得到 θH≈0.06,由方程(36)解得 l

1

= -0.031 8<0。 图 3(a)绘制了

参数为 θ 的食饵种群分岔图,其中插图是临界值(θc

= 0.075 6)附近的放大。 由图 3( a)可知,随着 θ 的增

大,系统经历从稳定状态到 Hopf 分岔(HB),以及极限环的不稳定状态和急速下降,直到食饵种群的密度

低于 θc 时,食饵进入灭绝状态。 可见,模拟得到的分岔点与理论计算阈值(θH≈0.06)一致。 若 θ>θc,则灭

绝平衡点 E0 是局部渐近稳定的。 系统的正平衡点 E1 在 θ = θH 发生 Hopf 分岔,当 0<θ<θH 时,正平衡点 E1

是局部渐近稳定的。 图 3(b)绘制了参数为 θ 的捕食者种群分岔图。 由图 3(b)可知,捕食者种群经历的

状态和食饵种群经历的状态相同,但是捕食者种群在稳定状态(θ∈(0,0.06))的密度是单调递减的,而食

饵种群的密度一直保持稳定不变,即 x = 0.423 85。 因此,食饵种群受到 Allee 效应的影响越小,对保护一

些密度较小的珍稀动物的种群多样性是有益的。

图 3 系统(3)中参数为 θ 且 ε= 0 的食饵和捕食者分岔图

Fig. 3 Bifurcation diagram of prey population and predator population in system (3) with parameter θ and ε= 0

在图 4 中,假设系统(3)存在人为控制迁移。 固定参数 α= 0.24,u = 0.05,e = 0.32,c = 0.01。 图 4(a)绘

制了参数为 θ 的食饵种群分岔图,其中图 4(a)插图分别是分岔点(eθh

= 0.125)附近的放大和临界值(θ1

=

0.135)附近的放大。 由图 4(a)可知,随着 θ 的增大,系统经历从稳定状态到 Hopf 分岔(HB),以及极限环

的不稳定状态和急速下降,直到食饵种群的密度低于 θ1 时,食饵进入灭绝状态。 若 θ>θ1 ,则灭绝平衡点

E0 是局部渐近稳定的。 系统的正平衡点 E2 在 θ = θh 发生 Hopf 分岔,当 0<θ<θh 时,随着 θ 的增大,食饵种

群的密度在不断递减,此时 E2 在稳定状态是局部渐近稳定的。 图 4(b)绘制了参数为 θ 的捕食者种群分

岔图, 其中图 4(b)插图是临界值(θ2

= 0.156)附近的放大。 由图 4(b)可知,系统经历稳定状态和急剧下

降,直到捕食者种群灭绝,并且捕食者种群在稳定状态的密度也在逐渐减小。 因此, 生物控制迁移对维护

生态系统的动态平衡有着非常重要的意义。

112

http:∥xuebao.gxnu.edu.cn

图 4 系统(3)中参数为 θ 且 ε= 0.2 的食饵和捕食者分岔

Fig. 4 Bifurcation diagram of prey population and predator population in system (3) with parameter θ and ε= 0.2

在图 5 中,固定参数 α= 0.3,u = 0.2,e = 1,c = 0.001,ρ = 0.5,并解出 θ = 0.103 不满足条件(25)。 当系统

(3)没有人为控制迁移时,食饵种群和捕食者种群的相图曲线趋于一个稳定的极限环(如图 5( a));当系

统(3)存在人为控制迁移时,食饵种群和捕食者种群的相图曲线趋于一个稳定点 E2

= (0.46,0.95)(如图 5

(b))。 因此,人为控制迁移有利于生态系统的稳定性。

图 5 系统(3)在 θ= 0.103 时的相图

Fig.5 Phase diagram of system (3) at θ= 0.103

在图 6 中,固定参数 α= 0.24,u = 0.05,c = 0.01,ρ = 0.5,θ = 0.1。 绘制参数为 e 的捕食者种群分岔图,其

中图 6(a)插图分别是第 1 个分岔点(eH1

= 0.16)附近的放大、第 2 个分岔点(eH2

= 0.306)附近的放大和临

界值(ec

= 0.312)附近的放大,图 6(b)插图是临界值(e0

= 0.443)附近的放大。 当系统(3)没有人为控制迁

移时,由图 6(a)可知,e 从 0 增加到 0.4,系统(3)经历从灭绝状态到 Hopf 分岔(HB),以及极限环的不稳定

状态和急剧上升,再到经历微小的 Hopf 分岔(HB)和急剧下降,直到捕食者种群趋于灭绝,即 y = 0。 若捕

食者转化率高于临界值 ec,则灭绝平衡点 E0 是局部渐近稳定的;若捕食者转化率低于分岔点 eH1

,则轴向

平衡点 E11是局部渐近稳定的;若捕食者转化率 e∈(0.24,0.306),则正平衡点 E1 是局部渐近稳定的,如图

6(a)所示。 其次,如果系统(3)受到人为控制迁移时,系统经历稳定状态和急剧下降,直到捕食者种群灭

绝。 若捕食者转化率 e∈(0,0.443),则捕食者种群的密度随着 e 的增大而逐渐增加,并且正平衡点 E2 是

局部渐近稳定的,如图 6(b)所示。 因此,适当的人为控制迁移不仅提高了捕食者的转化率,而且有利于持

久维护物种的多样性。

113

广西师范大学学报(自然科学版),2022,40(2)

图 6 系统(3)中参数为 e 的捕食者种群分岔

Fig. 6 Bifurcation diagram of predator population with parameterein system (3)

5 结论

数值模拟结果和分析表明,食饵-捕食者系统在 Allee 效应和人工控制迁移作用下,出现了稳定的极限

环、稳定点以及 Hopf 分岔等复杂的动力学行为。 食饵种群受到 Allee 效应的影响越小越有利于维持其物

种多样性,从而对保护一些密度较小的珍稀动物种群是有益的,例如北极熊、金丝猴等;生态系统中有时需

要适当的人为控制迁移来维护系统的稳定性,如自然界中一些鸟类在捕获害虫时, 因为鸟类自身捕获率

较低,不能抑制害虫数量的增长,所以需要人为控制迁移一些鸟类种群来捕获害虫,从而使得鸟类和害虫

构成的生态系统保持动态平衡。

参 考 文 献

[1] MA Z H, WANG S F, WANG T T, et al. Stability analysis of prey-predator system with Holling type functional response

and prey refuge[J]. Advances in Difference Equations, 2017, 2017: 243.

[2] LIU B, ZHANG Y J, CHEN L S. Dynamic complexities of a Holling I predator-prey model concerning periodic biological

and chemical control[J]. Chaos, Solitons & Fractals, 2004, 22(1): 123-134.

[3] ZHANG Y J, XU Z L, LIU B, et al. Dynamic analysis of a Holling I predator-prey system with mutual interference

concerning pest control[J]. Journal of Biological Systems, 2005, 13(1): 45-58.

[4] ZU J, MIMURA M. The impact of Allee effect on a predator-prey system with Holling type II functional response [ J].

Applied Mathematics and Computation, 2010, 217(7): 3542-3556.

[5] SUN G, SARWARDI S, PAL P J, et al. The spatial patterns through diffusion driven instability in modified Leslie-Gower

and Holling-type II predator-prey model[J]. Journal of Biological Systems, 2010, 18(3): 593-603.

[6] YI F Q, WEI J J, SHI J P. Bifurcation and spatiotemporal patterns in a homogeneous diffusive predator-prey system[ J].

Journal of Differential Equations, 2009, 246(5): 1944-1977.

[7] FAN M, KUANG Y. Dynamics of a nonautonomous predator-prey system with the Beddington-DeAngelis functional response

[J]. Journal of Mathematical Analysis and Applications, 2004, 295(1): 15-39.

[8] AJRALDI V, PITTAVINO M, VENTURINO E. Modeling herd behavior in population systems[J]. Nonlinear Analysis: Real

World Applications, 2011, 12(4): 2319-2338.

[9] KOSOBUD R F, O'NEILL W D. On the dependence of population growth on income: New results in a ricardian-Malthus

model[J]. De Economist, 1981, 129(2): 206-223.

[10] TSOULARIS A, WALLACE J. Analysis of logistic growth models[J]. Mathematical Biosciences, 2002, 179(1): 21-55.

[11] BOUKAL D S, SABELIS M W, BEREC L. How predator functional responses and Allee effects in prey affect the paradox of

enrichment and population collapses[J]. Theoretical Population Biology, 2007, 72(1): 136-147.

114

http:∥xuebao.gxnu.edu.cn

[12] KRAMER A M, DENNIS B, LIEBHOLD A M, et al. The evidence for Allee effects[ J]. Population Ecology, 2009, 51

(3): 341-354.

[13] BEREC L, ANGULO E, COURCHAMP F. Multiple Allee effects and population management[ J]. Trends in Ecology &

Evolution, 2007, 22(4): 185-191.

[14] BISWAS S, SAIFUDDIN M, SASMAL S K, et al. A delayed prey-predator system with prey subject to the strong Allee

effect and disease[J]. Nonlinear Dynamics, 2016, 84(3): 1569-1594.

[15] COURCHAMP F, BEREC L, GASCOIGNE J. Allee effects in ecology and conservation[ J]. Environmental Conservation,

2008, 36(1): 80-85.

[16] WU D Y, ZHAO H Y. Complex dynamics of a discrete predator-prey model with the prey subject to the Allee effect[ J].

Journal of Difference Equations Applications, 2017, 23(11): 1765-1806.

[17] PETROVSKII S, MOROZOV A, LI B L. Regimes of biological invasion in a predator-prey system with the Allee effect[J].

Bulletin of Mathematical Biology, 2005, 67(3): 637-661.

[18] RAO F, KANG Y. The complex dynamics of a diffusive prey-predator model with an Allee effect in prey[ J]. Ecological

Complexity, 2016, 28: 123-144.

[19] BARCLAY H J. Models for pest control using predator release, habitat management and pesticide release in combination

[J]. Journal of Applied Ecology, 1982, 19(2): 337-348.

[20] TANG S Y, TANG G Y, CHEKE R A. Optimum timing for integrated pest management: Modelling rates of pesticide

application and natural enemy releases[J]. Journal of Theoretical Biology, 2010, 264(2): 623-638.

[21] 成定平. 鼠类-天敌系统渐近稳定性的数学分析[J]. 生物数学学报, 2003, 18(3): 283-286.

[22] CHEN Y M, ZHANG F Q. Dynamics of a delayed predator-prey model with predator migration[ J]. Applied Mathematical

Modelling, 2013, 37(3): 1400-1412.

[23] JIN M, XU F, SHEN C S, et al. Nontrivial effect of time varying migration on the three species prey-predator system[ J].

Communications in Theoretical Physics, 2019, 71(1): 127-131.

[24] ABDLLAOUI A E, AUGER P, KOOI B W, et al. Effects of density-dependent migrations on stability of a two-patch

predator-prey model[J]. Mathematical Biosciences, 2007, 210(1): 335-354.

[25] HUANG Y, DIEKMANN O. Predator migration in response to prey density: What are the consequences? [ J]. Journal of

Mathematical Biology, 2001, 43(6): 561-581.

[26] SUN G Q, JIN Z, LIU Q X, et al. Dynamical complexity of a spatial predator-prey model with migration[ J]. Ecological

Modelling, 2008, 219(1 / 2): 248-255.

[27] PETROVSKII S, LI B L. An exactly solvable model of population dynamics with density-dependent migrations and the Allee

effect[J]. Mathematical Biosciences, 2003, 186(1): 79-91.

Stability of a Prey-predator Model with Migration and Allee Effects

XU Wangjun

1

, CAO Jinde

2

, WU Daiyong

1

, SHEN Chuansheng

1∗

(1. School of Mathematics and Physics, Anqing Normal University, Anqing Anhui 246133, China;

2. School of Mathematics, Southeast University, Nanjing Jiangsu 211189, China)

Abstract: A kind of prey-predator system with Allee effect and artificially controlled migration of predators is

studied. The system has a square root functional response function. Firstly, by qualitative analysis of the model,

the boundedness of the solution is proved, and the existence of the equilibrium point is analyzed. Sufficient

conditions for the local stability of the equilibrium point of the system are obtained. Then, the existence of the

Hopf-bifurcation of the equilibrium point is discussed, and the stability and direction of the equilibrium Hopfbifurcation are studied by calculating the first Lyapunov coefficient. Finally, the correctness of the conclusion is

verified by numerical simulation. The results indicate that the Allee effect and artificially controlled migration rate

are important for the survival and extinction of prey and predator populations.

Keywords: prey-predator model; Allee effect; migration rate; Hopf bifurcation; stability

(责任编辑 苏凯敏)

115

第 40 卷 第 2 期

2022 年 3 月

广西师范大学学报(自然科学版)

Journal of Guangxi Normal University (Natural Science Edition)

Vol. 40 No. 2

Mar. 2022

DOI: 10.16088 / j.issn.1001-6600.2021052401 http: xuebao.gxnu.edu.cn

蒋群群, 王林峰.一类非线性 p-Laplace 方程的 Liouville 定理[ J]. 广西师范大学学报(自然科学版), 2022, 40(2): 116-124. JIANG Q Q,

WANG L F. Liouville theorems for a nonlinear p-Laplace equation[J]. Journal of Guangxi Normal University (Natural Science Edition), 2022, 40

(2): 116-124.

一类非线性 p-Laplace 方程的 Liouville 定理

蒋群群, 王林峰∗

(南通大学 理学院, 江苏 南通 226019)

摘 要: 在带有适当曲率条件的完备流形上研究非线性 p-Laplace 方程 Δpu+au

p-1

ln u+λu

p-1 = 0, 式中 a、 λ 和 p>1 为给定

常数。 通过考虑几何量沿 p-Laplace 方程的演化, 在 Ricci 曲率有下界的紧致流形上建立上述方程的微分不等式。 借助截

断函数及 Hessian 比较定理, 在截面曲率有下界的非紧流形上也建立类似不等式。 作为应用得到了 Liouville 定理。

关键词: 微分不等式; 非线性; p-Laplace 方程; Liouville 定理; Ricci 孤立子

中图分类号: O186.1 文献标志码: A 文章编号: 1001-6600(2022)02-0116-09

假设 M 是一个带度量 g 的 n 维完备流形。 本文研究如下非线性 p-Laplace 方程

Δpu+au

p-1

ln u+λu

p-1 = 0, (1)

式中 p>1,且 Δpu = div( | u |

p-2

u)。

如果 p = 2, 式(1)为

Δu+auln u+λu = 0。 (2)

文献[1]得到 Ricci 曲率为负下界的完备流形上正调和函数的一个最优微分不等式。 λ 为正时方程 Δu+

λu = 0 的微分不等式的证明方法类似于文献[1],也可参看文献[2]。 文献[1]不等式中的等号能够取到,

因此这个微分不等式是最优的,这也意味着著名的 Yau 的 Liouville 定理[3]

, 即在具有非负 Ricci 曲率的完

备流形上不存在非常数正的调和函数。 光滑度量测量空间上类似的最优微分不等式在文献[4] 中被

建立。

当 a<0 时,Ma

[5]得到式(2)的正解的一个微分不等式,还通过观察式(2)与膨胀梯度 Ricci 孤立子之

间的关系,说明该微分不等式最优。 膨胀梯度 Ricci 孤立子由下式定义[6]

Ric+Hess f =

a

2

g, (3)

式中: f 为势函数;a<0。 Ma 观察到如果 f 是式(3)的势函数,则对某个常数 λ,u = e

-f满足式 (2)。

注意式(2)也与 Perelman 的 W 函数有关。 定义

λ = inf W(u) u > 0,u ∈ C

∞

(M),∫

M

u

2 { dx = 1} ,

式中

W(u) = ∫

M

(4 u

2 - u

2

ln u

2

)dx。

文献[7]证明在紧流形上 λ 可以由满足

4Δu+2uln u+λu = 0

的正光滑函数 u 取到。 文献[6]利用这个最小化结果来排除紧致流形上非平凡收缩的梯度 Ricci 孤立子

的存在性。 在其他情况下 W 泛函的研究可以参看文献[8-9]。 由文献[7]可知,式(2)也与著名的 Gross

收稿日期: 2021-05-24 修回日期: 2021-07-08

基金项目: 国家自然科学基金(11771223)

通信作者: 王林峰(1973—), 男, 江苏南通人, 南通大学教授, 博士。 E-mail: wlf711178@ntu.edu.cn

http:∥xuebao.gxnu.edu.cn

对数 Sobolev 不等式密切相关。

对 p>1,Kotschwar 等[10]在截面曲率有下界的假设下建立 p-调和函数的局部微分不等式。 由这个微分

不等式能推出 Liouville 定理,该定理表明在截面曲率非负的完备流形上不存在非常数正的 p-调和函数。

Wang 等[11]在截面曲率有下界的非紧流形上证明方程 Δpu+λu

p-1 = 0 正解的一个最优微分不等式。 与 pLaplace 相关的 Liouville 定理的研究还可以参看文献[12-13]及其中的参考文献。

对于 p>1 和某个给定常数 a,定义

λa,p

= inf Wa,p(u) u > 0,u ∈ C

∞

(M),∫

M

u

p { dx = 1} ,

式中

Wa,p(u) = ∫

M

u

p -

a

p

u

p

lnu

p

( ) dx。

在紧致流形上能使 λa,p取到的函数 u 满足 Euler-Lagrange 方程(1),式中 λ =λa,p。 关于 p-Laplace 方程的微

分不等式已有不少研究[14-15]

。 本文将建立非线性 p-Laplace 方程(1)正解的微分不等式。

1 本文主要结果及预备知识

为方便起见,设

? =

2

n+2

, a>0,

6(p-1)

6(p-1)+(4(p-1)+1)n

, a<0,

ì

î

í

ï

ï

ï

ï

(4)

以及

A =

(p-1)

1

4

(p-1)+

1

2

( ) n,a>0,

4(p-1)

2+8(p-1)+1

12

n, a<0。

ì

î

í

ï

ï

ï

ï

(5)

下面结果是紧流形上的微分不等式。

定理 1 设 M 为 Ricci 曲率以-(n-1)K 为下界的 n 维紧致流形,式中 K≥0。 假设 u 是式(1)的正解,

令 u

? = v,并且 M= sup

x∈M

{u(x)} ,则

① 对 a>0,如果 p≥2,则

v ≤ p

1

p a

1

p (2A)

-

1

p ?

p-2

p +A

-

1

2 ((n-1)K)

1

[ 2 ] M

?

, (6)

如果 1<p≤2,又有

v ≤ p

1

2 (p-1)

-

1

2 a

1

p (2A)

-

1

p ?

p-2

p +(p-1)

-

1

2 A

-

1

2 ((n-1)K)

1

[ 2 ] M

?

。 (7)

② 对 a<0,如果 p≥2,则

v ≤max pa(2A)

-1

?

p-2+A

-

p

2 ((n-1)K)

p

{ 2 ,0}

1

p M

?

, (8)

如果 1<p≤2,又有

v ≤A

-

1

2 ((n-1)K)

1

2 M

?

。 (9)

定理 1 能推出如下 Liouville 定理。

定理 2 设 M 为 Ricci 曲率以-(n-1)K 为下界的 n 维紧致流形,式中 K≥0。 假设 u 是 a<0 时式(1)

的正解。 对 p≥2,如果

pa(2A)

-1

?

p-2+A

-

p

2 ((n-1)K)

p

2 ≤0, (10)

或者当 1<p<2 时 K = 0,则 u≡e

-

λ

a 是常数。

117

广西师范大学学报(自然科学版),2022,40(2)

流形非紧时也可以建立式(1)的微分不等式。 引入截断函数时计算中包含距离函数的 Hessian,所以

假设截面曲率有下界[10]

。

定理 3 设 M 为 n 维完备非紧流形,截断曲率以-K 为下界,式中 K≥0。 假设 1<p≤2,u 是式(1)的正

解,设 M= sup

x∈M

{u(x)} ,如果 a>0,那么 v = u

? 满足

v ≤ p

1

2 (p-1)

-

1

2 a

1

p (2A)

-

1

p ?

p-2

p +(p-1)

-

1

2 A

-

1

2 ((n-1)K)

1

[ 2 ] M

?

, (11)

如果 a<0,则

v ≤A

-

1

2 ((n-1)K)

1

2 M

?

。 (12)

如第 3 章推论 2 所述,由定理 3 能得到关于截面曲率为负的非紧流形的 Liouville 定理。

假设 u>0 是式(1)的连续弱解[10,15]

(以下简称解),并且 M= sup

x∈M

{u(x)} 。 对任意 δ>0,存在只依赖于 p

和 δ 的常数 C(p,δ),使得

u

p-1

ln u ≤C(p,δ)max 1,u

p-1+δ

{ }≤C(p,δ)max 1,M

p-1+δ

{ } ,

则文献[16]中的梯度估计意味对某些 α>0, u

2∈C

α

,且对某些 β>1 有 u

2∈W

1,β

loc 。 实际上,u 在

{ u≠0}处光滑。 由 v = u

?

,得 v = ?u

?-1

u,ln v = ?ln u,且式(1)可改写为

Δp

v =

(p-1)(?-1)

?

v

p

v

-a?

p-2

v

p-1

ln v-λ?

p-1

v

p-1

。 (13)

对给定的满足式(13)的函数 v,定义算子 L 为

L(w)= div( v

p-2 w)+(p-2)div( v

p-4

( v· w) v)。

下面引理可以看作是 Bochner 公式的推广(见文献[17]式(2.5))。

引理 1 在任意 v≠0 点处,有

L( v

p

)= p v

2p-4 ‖Hess v‖2

g?

( +Ric( v, v) ) +p v

p-2 Δp

v· v, (14)

式中

g?= g+(p-2)

dv?dv

v

2

,

‖Hess v‖2

?g

= Hess v

2

g

+

p-2

2

v

2 2

v

2

+

(p-2)

2

4

v· v

2 2

v

4

。

基于引理 1 可以推导出以下估计。

引理 2 设 M 为 Ricci 曲率以-( n-1)K 为下界的 n 维完备非紧流形,式中 K≥0。 假设 v = u

? 是式

(13)的解,则

L( v

p

)+(n-1)Kp v

2p-2≥

pA

v

2p

v

2

-ap?

p-2

v

p-2

v

p+

p(p-1)(?-1)

?

υ

p-2

v

p· v

v

(15)

在所有 v≠0 的点处成立,这里 ? 和 A 分别见式(4)、(5)。

证明 由式(13)得

Δp

v· v =

(p-1)(?-1)

?

v

p· v

v

-

v

p-2

v

( 2 ) -

(p-1)a?

p-2

ln v+a?

p-2+(p-1)λ?

p-1

[ ] v

p-2

v

2

。 (16)

易知,

‖Hess v‖2

g?≥

1

n

(Trg?(Hess v))

2 =

1

n

Δv+

p-2

2

v· v

2

v

( 2 )

2

=

1

n

( v

2-pΔp

v)

2 =

1

n

v

4-2p

(Δp

v)

2 =

1

n

v

4-2p (p-1)(?-1)

?

v

p

v

- a?

p-2

v

p-1

ln v-λ?

p-1

v

é p-1

ë

ê

ê

ù

û

ú

ú

2

。 (17)

把式(16)、(17)代入式(14),得

L( v

p

)≥

p(p-1)

2

(?-1)

2

n?

2

-

p(p-1)(?-1)

? ( )

v

2p

v

2

-

118

http:∥xuebao.gxnu.edu.cn

p(p-1) 1+

2(?-1)

n?

( )

v

p

v

(a?

p-2

v

p-1

ln v+λ?

p-1

v

p-1

)+

p

n

(a?

p-2

v

p-1

ln v+λ?

p-1

v

p-1

)

2-(n-1)Kp v

2p-2-

ap?

p-2

v

p-2

v

p+

p(p-1)(?-1)

?

v

p-2

v

p· v

v

。

? 的定义意味着

a?

p-2

1+

2(?-1)

n?

( ) ≥0。 (18)

假设 δ 满足

a

n

?

p-2-

p-1

δ

1+

2(?-1)

n?

( ) = 0, (19)

在点 x∈M,如果 a?

p-2

v

p-1

ln v+λ?

p-1

v

p-1≤

a?

p-2

δ v

p

v

,则

L( v

p

)+(n-1)Kp v

2p-2≥

p(p-1)

(p-1)(?-1)

2

n?

2

-

?-1

?

- a?

p-2

δ 1+

2(?-1)

n?

( )

é

ë

ê

ê

ù

û

ú

ú

v

2p

v

2

-

ap?

p-2

v

p-2

v

p+

p(p-1)(?-1)

?

v

p-2

v

p· v

v

。 (20)

在点 x 处,如果 a?

p-2

v

p-1

ln v+λ?

p-1

v

p-1≥

a?

p-2

δ v

p

v

,则

L( v

p

)+(n-1)Kp v

2p-2≥

p(p-1)

2

(?-1)

2

n?

2

-

p(p-1)(?-1)

? ( )

v

2p

v

2

+

p(p-1)(?-1)

?

v

p-2

v

p· v

v

+

p

n

-

p(p-1)

a?

p-2

δ

1+

2(?-1)

n?

( )

é

ë

ê

ê

ù

û

ú

ú

(a?

p-2

v

p-1

ln v+λ?

p-1

ν

p-1

)

2-ap?

p-2

v

p-2

v

p =

p(p-1)

2

(?-1)

2

n?

2

-

p(p-1)(?-1)

? ( )

v

2p

v

2

- ap?

p-2

v

p-2

v

p+

p(p-1)(?-1)

?

v

p-2

v

p· v

v

,

这里最后一个等式来自式(19)。 因此式(20)在所有点都成立。 由式(19),

L( v

p

)+(n-1)Kp v

2p-2≥

p(p-1) -

3(p-1)

n

?-1

?

( )

2

-(4(p-1)+1)

?-1

?

- n(p-1)

é

ë

ê

ê

ù

û

ú

ú

v

2p

v

2

-

ap?

p-2

v

p-2

v

p+

p(p-1)(?-1)

?

v

p-2

v

p· v

v

=

pA

v

2p

v

2

- ap?

p-2

v

p-2

v

p+

p(p-1)(?-1)

?

v

p-2

v

p· v

v

,

从而得到式(15)。 证毕。

注 1 式(4)中选择的 ? 是最优的,因为

-

3(p-1)

n

?-1

?

( )

2

-(4(p-1)+1)

?-1

?

- n(p-1)

在式(18)的条件下达到最大值。

119

广西师范大学学报(自然科学版),2022,40(2)

2 紧致流形

本章证明定理 1 和定理 2。

定理 1 的证明 假设 v

p 在 x0∈M 取到最大值,并且 v

p

(x0 )>0。 因此,在 x0 处,有 v

p = 0,

L( v

p

)≤0。 由式(15)得

pA

v

2p

v

2 ≤(n-1)Kp v

2p-2+ap?

p-2

v

p-2

v

p

。 (21)

首先考虑 a>0 的情况。 如果 p≥2,容易看出

(n-1)K v

p-2≤

(p-2)A

pv

2

v

p+2((n-1)K)

p

2 p

-

p

2 (p-2)

p-2

2

(p-2)A

pv

( 2 )

-

p

2

+1

。 (22)

把式(22)代入式(21),得

v

p

(x0 )≤ pa(2A)

-1

?

p-2+A

-

p

2 ((n-1)K)

p

[ 2 ] v

p

(x0 )≤

pa(2A)

-1

?

p-2+A

-

p

2 ((n-1)K)

p

[ 2 ] M

?p

。

因此,对任意 x∈M,

v

p

(x)≤ pa(2A)

-1

?

p-2+A

-

p

2 ((n-1)K)

p

[ 2 ] M

?p

。

从而得式(6)。

如果 1<p<2,容易看出

a?

p-2

v

p-2

v

2-p≤(2-p)Av

-2

v

2+p2

-

2

p (a?

p-2

)

2

p A

1-

2

p 。 (23)

将式(23)代入式(21),并与 p≥2 的情况类似进行讨论,得到式(7)。

对 a<0 的情况,当 p≥2 时,将式(22)代入式(21),得

v

p

(x0 )≤ pa(2A)

-1

?

p-2+A

-

p

2 ((n-1)K)

p

[ 2 ] v

p

。

当 1<p<2 时,由式(21) 直接得 pA

v

2p

v

2 ≤( n-1) Kp v

2p-2在 x0 处成立,于是可以得到式(8)、(9)。

证毕。

现在证明定理 2。

定理 2 的证明 当 a<0 时,若对 p≥2,K 满足式(10),或者对 1<p<2,K = 0,从定理 1 可得 v≡0,因

此,u 是常数。 由式(1)知 u≡e

-

λ

a 。 证毕。

在定理 2 中令 p = 2 得到对式(2)的 Liouville 定理。

推论 1 设 M 为 Ricci 曲率以-(n-1)K 为下界的 n 维紧致流形,式中 K≥0。 假设 a<0 时 u 是式(2)

的正解。 如果 a+(n-1)K≤0,那么 u≡e

-

λ

a 是常数。

3 非紧流形

本章将利用截断函数给出非紧致流形上式(1) 正解的微分不等式。 由于计算涉及到距离函数的

Hessian,需要假设截断曲率有下界。 首先证明定理 3。

定理 3 的证明 设 O∈M 为定点,r( x) = dist(O,x) 为由 O 决定的距离函数。 对足够大的 R>0,用

B(O,R)表示中心为 O 半径为 R 的测地球。

120

http:∥xuebao.gxnu.edu.cn

考虑如下函数[15]

θ(t):[0,+∞ )→[0,1],θ(t)=

1,0≤t≤1,

0,t≥2, {

使得

-10θ

1

2 ≤θ′≤0,θ″>-10θ。 (24)

定义截断函数 φ:M→R 为 φ(x,t)= θ

r(x)

R

( ) ,计算表明:

L(φ v

p

)= div( v

p-2

(φ v

p

))+(p-2)div( v

p-4

( v· (φ v

p

)) v)=

φL( v

p

)+2(p-1) v

2p-4

v

2· φ+(p-2) v

p-2

v· φΔp

v +

v

2p-2 Δφ+

(p-2)Hess φ( v, v)

v

( 2 ) +

(p-1)(p-2) v

2p-6

( v· v

2

)( v· φ)。 (25)

设 φ v

p =G,假设 G 在 x0∈M 时取到最大值,不妨假设 G(x0 )>0,那么在 x0 点,φ>0, v≠0, G= 0,即

v

2 = -

2

pφ

v

2

φ。 (26)

在 x0 点,还有 L(G)≤0。 对任意 σ∈(0,1),有

‖Hess v‖2

g?≥

1

n

( v

2-pΔp

v)

2 =

σ

n

v

4-2p

(Δp

v)

2+

1-σ

n

v

4-2p

(Δp

v)

2 =

σ

n

v

4-2p (p-1)(?-1)

?

v

p

v

- a?

p-2

v

p-1

ln v-λ?

p-1

v

é p-1

ë

ê

ê

ù

û

ú

ú

2

+

1-σ

n

v

4-2p

(Δp

v)

2

。 (27)

令

?σ

=

2σ

n+2σ

, a≥0,

6(p-1)σ

6(p-1)σ+(4(p-1)+1)n

, a<0,

ì

î

í

ï

ï

ï

ï

以及

Aσ

=

(p-1)

1

4

(p-1)+

1

2

( )

n

σ

, a≥0,

4(p-1)

2+8(p-1)+1

12

n

σ

, a<0。

ì

î

í

ï

ï

ï

ï

假设 δσ 满足

aσ

n

?

p-2

σ

-

p-1

δσ

1+

2σ(?σ

-1)

n?σ

( ) = 0。

与式(15)类似,有

L( v

p

)+(n-1)Kp v

2p-2-

p(1-σ)

n

(Δp

v)

2≥

pAσ

v

2p

v

2

- ap?

p-2

σ

v

p-2

v

p+

p(p-1)(?σ

-1)

?σ

v

p-2

v

p· v

v

。

由式(25)、(26)、 (27),以及在 x0 处 L(G)≤0,得在 x0 处有

0≥pAσ φ

v

2p

v

2

- apφ?

p-2

σ

v

p-2

v

p-(n-1)Kpφ v

2p-2+

121

广西师范大学学报(自然科学版),2022,40(2)

p(1-σ)φ

n

(Δp

v)

2-

p(p-1)(?σ

-1)

?σ

v

2p-2

v· φ

v

-

4(p-1)

p

v

2p-2 φ

2

φ

+(p-2) v

p-2

v· φΔp

v+

v

2p-2 Δφ+

(p-2)Hessφ( v, v)

v

( 2 ) -

2(p-1)(p-2)

p

v

2p-4

v· φ

2

φ

。

易见

p(1-σ)φ

n

(Δp

v)

2+(p-2) v

p-2

v· φΔp

v≥-

(p-2)

2

n

4p(1-σ)

v

2p-4

v· φ

2

φ

,

以及 v· φ ≤ v φ ,则在 x0 处有

0≥pAσ φ

v

2p

v

2

- apφ?

p-2

σ

v

p-2

v

p- (n-1)Kpφ v

2p-2+

p(1-σ)φ

n

(Δp

v)

2-

p(p-1)(?σ

-1)

?σ

v

2p-2

v· φ

v

-

4(p-1)

p

v

2p-2 φ

2

φ

+(p-2) v

p-2

v· φΔp

v+

v

2p-2 Δφ+

(p-2)Hess φ( v, v)

v

( 2 ) -

2(p-1)(p-2)

p

v

2p-4

v· φ

2

φ

,

及

Δφ+(p-2) v

-2Hess φ( v, v)=

θ″

R

2

g+

p-2

v

2

dv?dv ( ) ( r, r)+

θ′

R

<g+

p-2

v

2

dv?dv,Hess r>。

由式(24)、 Hessian 比较定理[18]

,Hess r≤

1+Kr

r

g,以及 g?= g+

p-2

v

2

dv?dv 正定,得

Δφ+

(p-2)Hess φ( v, v)

v

2 ≥-

10φ

R

2

-

10(n+p-2) φ

R

1+KR

R

。

易见 φ =

θ′

R

≤

10 φ

R

。 那么在 x0 处有

0≥pAσ

G

2

v

2

φ

- ap?

p-2

σ

v

p-2G-(n-1)Kpφ

2

p

- 1G

2 -

2

p +

10p(p-1)(?σ

-1)

?σ

φ

1

p

-

3

2

Rv

G

2 -

1

p -

400(p-1)

pR

2

φ

2

p

-2G

2 -

2

p -

10φ

R

2

+

10(n+p-2) φ

R

1+ K R

R

é

ë

ê

ê

ù

û

ú

ú

φ

2

p

- 2G

2 -

2

p -

2(p-1) p-2

p

+

(p-2)

2

n

4p(1-σ)

é

ë

ê

ê

ù

û

ú

ú

100

R

2

φ

2

p

-2G

2 -

2

p 。

由 0≤φ≤1,得

pAσ

G

2

v

2 ≤ap?

p-2

σ

v

p-2G+(n-1)KpG

2 -

2

p +

C1φ

1

p

-

1

2

Rv

G

2 -

1

p +

C2

+C3R

R

2

φ

2

p

-1G

2 -

2

p , (28)

式中 Ci,i = 1,2,3,…是依赖于 n、p、σ、K、?σ 的常数。 由 1<p<2,0≤φ≤1,得

pAσ

G

2

v

2 ≤ap?

p-2

σ

v

p-2G+(n-1)KpG

2-

2

p +

C1

Rv

G

2-

1

p +

C2

+C3R

R

2

G

2-

2

p 。 (29)

首先考虑 a>0 的情形。 易知

a?

p-2

σ

v

p-2G≤(2-p)Aσ

v

-2G

2+2

-

2

p p(a?

p-2

σ )

2

p A

1-

2

p

σ G

2-

2

p 。 (30)

122

http:∥xuebao.gxnu.edu.cn

容易看出,对充分小的 δ?>0,

C1

Rv

G

2-

1

p ≤δ?Aσ

G

2

v

2

+

C

2

1

4δ?Aσ R

2

G

2-

2

p 。 (31)

将式(30)、(31)代入式(29),得

(p(p-1)-δ?)Aσ G

2

p (x0 )≤ p

2

2

-

2

p (a?

p-2

σ )

2

p A

1-

2

p

σ

+(n-1)Kp+

C

2

1

4δ?Aσ R

2

+

C2

+C3R

R

2

é

ë

ê

ê

ù

û

ú

ú

M

2

。

因此,对所有 x∈B(O,R),

(p(p-1)-δ?)Aσ

v

2

(x)= (p(p-1)-δ?)Aσ G

2

p (x)≤(p(p-1)-δ?)Aσ G

2

p (x0 )≤

p

2

2

-

2

p (a?

p-2

σ )

2

p A

1-

2

p

σ

+(n-1)Kp+

C

2

1

4δ?Aσ R

2

+

C2

+C3R

R

2

é

ë

ê

ê

ù

û

ú

ú

M

2

。

令 R→∞ ,对所有 x∈M,

(p(p-1)-δ?)Aσ

v

2

(x)≤ p

2

2

-

2

p (a?

p-2

σ )

2

p A

1-

2

p

σ

[ +(n-1)Kp] M

2

。

令 δ?↘0,得

(p-1)Aσ

v

2

(x)≤ p2

-

2

p (a?

p-2

σ )

2

p A

1-

2

p

σ

[ +(n-1)K] M

2

。 (32)

注意当 σ↗1 时,?σ→?,Aσ→A,在式(32)中令 σ↗1,得

(p-1)A v

2

(x)≤ p2

-

2

p (a?

p-2

)

2

p A

1-

2

[ p +(n-1)K] M

2

。

这个不等式意味着式(11),通过类似讨论可以得到式(12)。 证毕。

推论 2 在具有非负截面曲率的 n 维非紧流形上,当 a<0,且 1<p≤2 时,式(1)无非平凡的正有界解。

参 考 文 献

[1] LI P, WANG J P. Complete manifolds with positive spectrum, II[J]. Journal of Differential Geometry, 2002, 62: 143-162.

[2] LI P. Geometric analysis[M]. Cambridge: Cambridge University Press, 2012.

[3] YAU S T. Harmonic functions on complete Riemannian manifolds[J]. Communications on Pure and Applied Mathematics,

1975, 28(2): 201-228.

[4] WANG L F. The upper bound of the L

2

μ

spectrum[J]. Annals of Global Analysis and Geometry, 2010, 37(4) : 393-402.

[5] MA L. Gradient estimates for a simple elliptic equation on complete non-compact Riemannian manifolds[ J]. Journal of

Functional Analysis, 2006, 241(1): 374-382.

[6] PERELMAN G. The entropy formula for the Ricci flow and its geometric applications[EB/ OL]. (2002-11-11) [2021-05-

24]. https:∥arxiv.org / abs/ math / 0211159.

[7] ROTHAUS S. Logarithmic Sobolev inequalities and the spectrum of Schrödinger operators [ J ]. Journal of Functional

Analysis, 1981, 42(1): 110-120.

[8] NI L. The entropy formula for linear heat equation[J]. Journal of Geometric Analysis, 2004, 14(1): 87-100.

[9] WANG L F. A new entropy formula for the linear heat equation[J]. Archiv der Mathematik, 2011, 96(5): 473-481.

[10] KOTSCHWAR B, NI L. Local gradient estimates of p-harmonic functions, 1 / H-flow, and an entropy formula[ J]. Annales

Scientifiques de l′École Normale Supérieure, 2009, 42(1): 1-36.

[11] WANG L F, ZHU Y P. A sharp gradient estimate for the weighted p-Laplacian[ J]. Applied Mathematics-A Journal of

Chinese Universities, 2012, 27(4): 462-474.

[12] 周俞洁, 张泽宇, 王林峰. 黎曼流形上 p-Laplace 算子的 Liouville 定理[ J]. 西南大学学报(自然科学版), 2017,

123

广西师范大学学报(自然科学版),2022,40(2)

39(10): 62-68.

[13] 魏公明. 具奇系数发展型 p-Laplace 不等方程整体解的不存在性[J]. 数学年刊 A 辑, 2007, 28(3): 387-394.

[14] 汪悦. Riemann 流形上改进的 p-Laplace 方程的梯度估计[J]. 中国科学(数学), 2014, 44(3): 288-294.

[15] 王林峰. p-Laplace Schrödinger 热方程的椭圆型梯度估计[J]. 数学学报, 2010, 53(4): 643-654.

[16] TOLKSDORF P. Regularity for a more general class of quasilinear elliptic equations[ J]. Journal of Differential Equations,

1984, 51(1): 126-150.

[17] WANG L F. Gradient estimates on the weighted p-Laplace heat equation [ J]. Journal of Differential Equations, 2018,

264(1): 506-524.

[18] SCHEON R, YAU S T. Lectures on differential geometry[M]. Cambridge, MA: International Press, 1994.

Liouville Theorems for a Nonlinear p-Laplace Equation

JIANG Qunqun, WANG Linfeng

∗

(School of Sciences, Nantong University, Nantong Jiangsu 226019, China)

Abstract: In this paper the nonlinear p-Laplace equation Δpu + au

p-1

lnu + λu

p-1 = 0 is studied on complete

manifolds with some suitable curvature condition, where a, λ and p > 1 are some given constants. Differential

inequalities for the p-Laplace equation on compact manifolds with Ricci curvature bounded from below are

established, based on the evolution of the geometric quantity along the p-Laplace equation. Similar inequalities

can also be established on a noncompact manifold whose sectional curvature is bounded from below, based on the

skills of cut off function and the Hessian comparison theorem. As an application, Liouville theorems are obtained.

Keywords: differential inequality; nonlinear; p-Laplace equation; Liouville theorem; Ricci soliton

(责任编辑 吴佃华)

124

第 40 卷 第 2 期

2022 年 3 月

广西师范大学学报(自然科学版)

Journal of Guangxi Normal University (Natural Science Edition)

Vol. 40 No. 2

Mar. 2022

DOI: 10.16088 / j.issn.1001-6600.2021061201 http: xuebao.gxnu.edu.cn

王涵, 张映辉. 模拟趋化现象的三维双曲-抛物系统的最优衰减率[ J]. 广西师范大学学报(自然科学版), 2022, 40( 2): 125-131. WANG

H, ZHANG Y H. Optimal time-decay rates of the hyperbolic-parabolic system modeling chemotaxis in R

3

[ J]. Journal of Guangxi Normal University

(Natural Science Edition), 2022, 40(2): 125-131.

模拟趋化现象的三维双曲-抛物系统的最优衰减率

王 涵, 张映辉∗

(广西师范大学 数学与统计学院, 广西 桂林 541006)

摘 要: 本文研究一个模拟趋化现象的三维双曲-抛物系统的 Cauchy 问题解的大时间行为, 得到其解及其各阶空间导数

的最优时间衰减率。 跟已有结果相比, 本文主要创新在于给出解的最高阶空间导数的最优衰减率, 且该衰减率与热方程

的衰减率一样。 研究方法主要基于高频-低频分解和精细的能量估计。

关键词: 双曲-抛物系统; 最优衰减速率; 高频-低频分解; 大时间行为; 趋化现象

中图分类号: O29 文献标志码: A 文章编号: 1001-6600(2022)02-0125-07

1 本文主要成果

本文研究一个模拟趋化现象的三维双曲-抛物系统解的时间衰减率,该系统可表示为如下形式:

vt

- v = 0,x∈R

3

, t>0,

ut

- ·(uv)= DΔu,x∈R

3

, t>0。 { (1)

初值为

(v, u)(x,0)= (v0 ,u0 )(x)→(0,u?), | x |→∞ , (2)

式中常数 u?>0。 式(1) 与以下系统密切相关,

∂p

∂t

=D · p ln

p

Φ(w)

( ( ) ) ,

∂w

∂t

= βpw。

ì

î

í

ï

ï

ï

ï

(3)

该模型由 Othmer 等[1]以及 Levine 等[2]通过大量生物学估计和数值计算得到。 式中:p( x,t)表示粒子密

度;w(x,t)表示化学物质浓度;D 表示粒子扩散速率且 D>0;Φ 表示化学势能。

类似文献[3-9],假设 Φ(w)= w

-α

,α 是大于零的常数,系统(3)可以改写为

pt

=DΔp+Dα · p

w

w

( ) ,

wt

= βpw,

ì

î

í

ïï

ï

(4)

进一步,假设

q = (ln w)=

w

w

,

系统(4)可改写为

pt

=DΔp+Dα ·(pq),

qt

= β p。 { (5)

收稿日期: 2021-06-12 修回日期: 2021-07-02

基金项目: 国家 自 然 科 学 基 金 ( 11771150, 11571280 ); 广 西 自 然 科 学 基 金 ( 2019JJG110003, 2019AC20214,

2019JJA110071)

通信作者: 张映辉(1981—), 男, 湖南祁阳人, 广西师范大学教授, 博士。 E-mail: yinghuizhang@gxnu.edu.cn

广西师范大学学报(自然科学版),2022,40(2)

最后,正常数 A、B 和 c1 由以下式子给出。 令 τ = At,ξ =Bx,u = p,v = c1

q, 则将系统(5)改写为

vτ

=

βBc1

A

ξu,

uτ

=

DB

2

A

Δξu+

DαB

Ac1

ξ·(uv)。

ì

î

í

ï

ï

ï

ï

(6)

若在式(6)中取

βBc1

A

= 1,

B

2

A

= 1,

DαB

Ac1

= 1,

ì

î

í

ï

ï

ï

ï

ï

ï

ï

ï

即,

A =Dαβ>0,B = Dαβ >0,c1

=

Dα

β

>0,

那么易知 u 和 v 满足

vτ

- ξu = 0,

uτ

- ξ·(uv)= DΔξu。 { (7)

用变量(x,t)替换(τ,ξ),则式(7)恰好是式(1)。

对系统(1)的一维情形, 文献[5,10]和[3]分别考虑了初边值问题和 Cauchy 问题光滑解的存在性和

渐近性;文献[5]考虑了系统(1)的初边值问题, 当 u0

-1

2

H2+ v0

2

H2足够小时,证明系统(1)光滑解的全局

存在性;文献[3]证明在初始值充分大的条件下,系统(1)Cauchy 问题光滑解的全局存在性;文献[11]考

虑带流量限制的趋化模型的初边值问题, 证明经典解整体存在且关于时间一致有界;文献[12]研究非线

性扩散-趋化聚集模型,证明生物趋化模型弱解的全局有界性及唯一性;文献[13] 得到具有不同分数阶扩

散趋化模型的衰减估计;文献[14] 研究一个耦合双曲-抛物系统的全局光滑解;对于高维情形,文献

[9,15-17]分别研究系统(1)初边值问题和 Cauchy 问题的光滑小解全局存在性。 但是, 对高维趋化模型

全局解的大时间行为方面的研究结果很少。 文献[4]首次详细证明在定常状态下, 初始扰动的 H

2 范数足

够小,并且 L

1 范数有界时,三维问题强解的全局存在性和最优时间衰减速率

(v,u)(t) L

2≤C0(1+t)

-

3

4 和 (v,u)(t) H1≤C0(1+t)

-

5

4 。 (8)

最近, 文献[14]得到当(v0 ,u0 )∈H

N∩L

p

(p∈[1,2]), 且 (v0 ,u0 ) H2≤δ0 时,L

p -L

2 型时间衰减速率,其

主要结果为

ℓ(v, u)(t) HN-ℓ≤C0(1+t)

-

3

2

1

p

-

1

2 ( ) +

ℓ

2 ,ℓ = 0,…,N-1。 (9)

一方面,容易发现式(8)中解(v,u)的 2 阶(即最高阶)空间导数的 L

2 衰减速率是(1+t)

-

5

4 ,与它的 1

阶空间导数具有相同的衰减速率,并且比热方程的 L

2 衰减速率(1+t)

-

7

4 慢。 另一方面,如果在式(9)中取

p = 1,ℓ =N-1,那么(v,u)的 N 阶(即最高阶)空间导数的 L

2 衰减速率为(1+ t)

-

3+2(N-1)

4 ,与它的 N-1 阶空间

导数的衰减速率相同,并且比热方程的 L

2 衰减速率(1+t)

-

3+2N

4 慢。 因此,在这种意义下,式(8)、(9)中解的

最高阶空间导数的衰减率都不是最优。

本文旨在对上述问题给出一个明确答案,更精确地说,本文得到 Cauchy 问题(1)、(2)的解及其各阶

空间导数(1 阶到最高阶 N 阶)的最优时间衰减速率。 值得一提的是,本文得到解的任一阶空间导数的衰

减速率都与热方程的衰减速率一样,证明方法是基于高频-低频分解和精细的能量估计。

类似文献[4],重新改写 Cauchy 问题(1)、(2)。 假设

126

http:∥xuebao.gxnu.edu.cn

λ =

1

u?

,λ1

= u?。

作变量替换(v,u)→(λv,u+u?),在(0,u?)附近将系统(1)线性化,Cauchy 问题(1)、(2)改写为

vt

-λ1 u = 0,

ut

-λ1 ·v-DΔu =λ ·(uv),

(v,u)(x,0)= (v0 ,u0 )(x)→(0,0), | x |→∞ 。

ì

î

í

ï

ï

ïï

(10)

在本处以及后续工作中,为了简单起见,仍用(v,u)表示重新构造的变量。 为方便起见,本文采用以下记

号:H

s

(R

3

)(s∈R)表示通常 Sobolev 空间,具有范数 · Hs;L

p

(R

3

) (1≤p≤∞ )表示通常的 L

p 空间,且范

数为 · L

p;

ℓ 表示任意 ℓ 阶空间导数,式中 ℓ≥0;当 ℓ<0 或 ℓ 不是正整数时,

ℓ 表示 Λ

ℓ 并且 Λ

ℓ

f: =

F

-1

( | ξ |

iF f),式中 F 表示通常的 Fourier 变换算子,F

-1表示其逆算子。

本文假设 C 为大于零的常数,符号 a≲b 表示 a≤Cb,式中 C 只依赖于问题参数,且 C>0。 设存在径向

函数 φ∈C

∞

0 (R

3

ξ),当| ξ |≤1 时满足 φ(ξ)= 1,当 | ξ | ≥2 时满足 φ(ξ)= 0。 定义 f 的低频部分和高频部分

分别为

f

l =F

-1

[φ(ξ)

^

f], f

h=F

-1

[(1-φ(ξ))

^

f]。

为了简单起见, 记 (A,B) X := A X

+ B X 。

本文主要结果如下。

定理 1 假设 ×v0

= 0,(v0 ,u0 )∈H

N∩ L

1

,整数 N≥2,v0 、u0∈L

1

,存在常数 δ0 ,使得

(v0 ,u0 ) H2≤δ0 , (11)

则 Cauchy 问题(10)存在一个唯一全局解(v,u),使得对所有 t≥0,以下衰减估计成立:

k(v,u)(t) HN-k≲(1+t)

-

3+2k

4 ,k = 0,…,N。 (12)

注 1 与文献[18]的结果(式(9))相比, 本文主要结果(式(11))中关于解(v,u)的最高阶(即 N)空

间导数的 L

2 衰减速率的结果是全新的。 特别地,本文结果蕴含解(v,u)的最高阶(N)空间导数的 L

2 衰减

速率为(1+t)

-

3+2N

4 ,比文献[18]的结果(式(9))中的(1+t)

-

3+2(N-1)

4 要快, 并且该衰减速率与热方程的衰减速

率相同。 因此,本文的衰减速率是最优的。

现在简单给出本文分析要点。 基于文献[18]结果,只需要证明式(12)中关于解( v,u)的最高阶(N)

空间导数的 L

2 衰减速率。 证明主要包括如下 3 个步骤。

首先得到解的能量估计:

d

dt

N(v,u)

2

L

2+C

N+1 u

2

L

2≲δ(

N

v

2

L

2+ N+1 u

2

L

2 )。 (13)

其次,注意到上述能量估计中不包含 v 的耗散。 因此,为了揭示 v 的耗散,可以利用文献[14]方法,

即通过构造 v 和 u 之间的相互作用能量泛函,得

-

d

dt

∫

R3

N-1 u

N

vdx + C

N

v

2

L

2 ≲ Nu

2

H1 。 (14)

注意到式(14)中时间导数的积分项包含 N-1 u,所以,通过式(13)、(14)似乎不可能封闭解的 N 阶能量估

计。 为此将充分利用高频-低频分解好的性质来克服这一困难。 具体地说,代替构造 v 和 u 之间的相互作

用能量泛函,本文构造 v

h 和 u

h 之间新的相互作用能量泛函,得

-

d

dt

∫

R3

N-1 u

h N

v

h

dx + C

N

v

h 2

L

2 ≲ Nu

h 2

L

2 + N+1 u

h 2

L

2 + (1 + t)

-

3+2N

2 。 (15)

于是,取足够大的 T0 和 D0 ,定义能量泛函

E(t) = D0

N(v,u)

2

L

2 - ∫

R3

N-1 u

h N

v

h

dx,

当 t≥T0 时,由于 D0 充分大,则 E(t)等价于 N(v,u)

2

L

2 。 因此,结合式(13)、(15),可得

d

dt

E(t)+CE(t)≲(1+t)

-

3+2N

2 + N(v

l

,u

l

)

2

L

2 。 (16)

127

广西师范大学学报(自然科学版),2022,40(2)

最后, 利用 Plancherel 定理、Hausdorff-Young 不等式和 Duhamel 原理可得到解的最高阶空间导数的低

频衰减速率。

N(v

l

,u

l

)

2

L

2≲(1+t)

-

3+2N

2 。 (17)

结合式(16)、(17),并利用 Gronwall 不等式可得到解的最高阶空间导数的最优衰减速率,从而完成定理 1

的证明。

2 非线性能量估计

本节将推导式(10)的先验非线性能量估计。 假设对足够小的 δ>0,有先验估计

(v,u)(t) H2≤δ。

本节将广泛使用 Gagliardo-Nirenberg Sobolev 插值不等式。

以下引理是文献[19]第 125 页定理的特例。

引理 1 假设 0≤i,j≤k,则

i

f

L

p≲ Δ

j

f

1-θ

L

q

k

f

θ

L

r,

式中 θ 满足

i

3

-

1

p

=

j

3

-

1

q

( ) (1-θ)+

k

3

-

1

r

( ) θ。

引理 2 假设整数 k≥1,则

k(fg) L

p≲ f

L

p

1

kg L

p

2

+ k

f

L

p

3

g L

p

4

,

式中 p、p2 、p3∈[1,+∞ ],且

1

p

=

1

p1

+

1

p2

=

1

p3

+

1

p4

。

证明 对 p = p2

= p3

= 2,可以利用引理 1 证明,详情可参看文献[20]第 98 页。 对于一般情形,参看文

献[21]中引理 3.1。

引理 3 若对任意 2≤p≤∞ , f∈L

p

(R

3

),则

f

l

L

p+ f

h

L

p≲ f

L

p。

证明 对 2≤p≤∞ , 根据卷积的 Young 不等式,对于低频可得

f

l

L

p≲ F

-1

φ L

1 f

L

p≲ f

L

p,

因此,

f

h

L

p≲ f

L

p+ f

l

L

p≲ f

L

p。

证毕。

接下来给出一类包含 u 耗散估计的能量估计。

引理 4 假设 0≤k≤N-1, 则

d

dt

k+1(v,u)

2

L

2+2D k+2 u

2

L

2≲δ(

k+1

v

2

L

2+ k+2 u

2

L

2 )。 (18)

证明 对 0≤k≤N-1,将 k+1作用于式(10)的第 1、2 式子,且将所得等式分别乘以 k+1

v 和 k+1 u,再将

其结果求和,然后在 R

3 上进行分部积分,可得

1

2

d

dt

∫

R3

|

k+1(v,u) |

2

dx + D k+2 u

2

L

2 = ∫

R3

k+1(λ ·(uv))

k+1 udx。 (19)

下面估计等式(19)的右边项,由分部积分、Hölder 不等式、Cauchy-Schwarz 不等式和引理 1、2,可得

∫

R3

k+1(λ ·(uv))

k+1 udx = - λ∫

R3

k+1(uv)

k+2 udx ≲ k+1(vu) L

2

k+2 u L

2 ≲ (

k+1 u L

6 v L

3 +

128

http:∥xuebao.gxnu.edu.cn

k+1

v L

2 u L∞ )

k+2 u L

2 ≲ ( v H1

k+2 u L

2 + u H2

k+1

v L

2 )

k+2 u L

2 ≲ δ(

k+1

v

2

L

2 + k+2 u

2

L

2 )。

(20)

结合式(19)、(20)即可证得式(18)成立。 证毕。

最后,考虑以下线性系统:

vt

-λ1 u = 0,

ut

-λ1 ·v-DΔu = 0,

(v,u)(x,0)= (v0 ,u0 )(x)→(0,0), | x |→∞ 。

ì

î

í

ï

ï

ïï

(21)

给出线性系统(21)的 L

2 时间衰减速率。

以下引理 5 见文献[4]。

引理 5 假设(v?,u?)是系统(21)的解,初值(v?0 ,u?0 )∈H

N∩L

1

,则对 0≤k≤N,有

k(v?

l

,u?

l

)(t) L

2≲(1+t)

-

3

4

-

k

2 ( (v?0 ,u?0 ) L

1+ k(v?0 ,u?0 ) L

2 )。

引理 6 在定理 1 的假设条件下,非线性系统(10)的解(v,u)满足以下衰减估计

k(v

l

,u

l

)(t) L

2≲(1+t)

-

3

4

-

k

2 ,k = 0,…,N。 (22)

证明 定义 S = ( S

1

,S

2

)

t = (0,λ ·( uv))

t

,根据引理 5、引理 2、引理 3、Plancherel 定理、HausdorffYoung 不等式、Duhamel 原理和 Hölder 不等式,可得

k(v

l

,u

l

)(t) L

2 ≲ (1 + t)

-

3

4

-

k

2 (v,u)(0) L

1 + ∫

1

2

0

(1 + t - τ)

-

3

4

-

k

2 S(τ) L

1 dτ +

∫

t

t

2

(1 + t - τ)

-

5

4

| ξ |

k-1

S^

l

(τ) L∞ dτ ≲ (1 + t)

-

3

4

-

k

2 + ∫

t

t

2

(1 + t - τ)

-

5

4

| ξ |

k-1

S^

l

(τ) L∞ dτ。 (23)

另一方面,估计 | ξ |

k-1

S^

l

(τ) L∞ 。

| ξ |

k-1

S^

l

(τ) L∞ ≲ k-2(λ ·(uv)(τ) L

1 =λ

k-1(uv)(τ) L

1≲λ( u(τ) L

2

k-1

v(τ) L

2+

v(τ) L

2

k-1 u(τ) L

2 )≲λ( (v,u)(τ) L

2

k-1(v,u)(τ) L

2 )≲(1+τ)

-

3

4 (1+τ)

-

3

4

-

k-1

2 ≲(1+τ)

-1-

k

2 。 (24)

将式(24)代入式(23)即可证式(22)。 证毕。

3 定理 1 的证明

本章将利用低频和高频分解完成定理 1 的证明。 定理 1 中 l≤N-1 的情形在文献[14]已证明,只需

证明当 l = N 时

N(v,u)(t) L

2≲(1+t)

-

3+2N

4 (25)

成立,即可完成定理 1 的证明。

步骤 1

N

v

h 的耗散。 将 N-1F

-1

(1-φ(ξ))作用到式(10)的第 2 个式子,所得结果与 N

v

h 相乘,并

在 R

3 上积分,得

∫

R3

N-1 u

h

t

N

v

h

dx - λ1 ∫

R3

N

v

h N

v

h

dx - D∫

R3

N+1 u

h N

v

h

dx = λ∫

R3

N(uv)

h N

v

h

dx。

又因为

∫

R3

N-1 u

h

t

N

v

h

dx =

d

dt

∫

R3

N-1 u

h N

v

h

dx - ∫

R3

N-1 u

h N

v

h

t dx =

d

dt

∫

R3

N-1 u

h N

v

h

dx - λ1 ∫

R3

N-1 u

h N+1 u

h

dx。

于是,

-

d

dt

∫

R3

N-1 u

h N

v

h

dx + λ1

N

v

h 2

L

2 = - λ1 ∫

R3

N-1 u

h N+1 u

h

dx - D∫

R3

N+1 u

h N

v

h

dx -

λ∫

R3

N(uv)

h N

v

h

dx: = I1

+ I2

+ I3 。 (26)

下面对等式(25) 右边进行估计。 首先估计 I1 。 由分部积分得

129

广西师范大学学报(自然科学版),2022,40(2)

| I1

| = λ1

| ∫

R3

N-1 u

h N+1 u

h

dx | ≲ Nu

h 2

L

2 。

其次估计 I2 。 根据 Hölder 不等式和 Young 不等式,得

| I2

| = D | ∫

R3

N+1 u

h N

v

h

dx | ≲ N+1 u

h

L

2

N

v

h

L

2 ≲ N+1 u

h 2

L

2 +

λ1

4

N

v

h 2

L

2 。

最后估计 I3 。 根据 Hölder 不等式、引理 3、引理 2、引理 1、式(8)、(9)和 Young 不等式,得

|I3

| =λ | 〈

N(uv)

h

,

N

v

h

〉 |≲ N(uv)

h

L

2

N

v

h

L

2≲ (v,u) L∞

N(v,u) L

2

N

v

h

L

2≲

(v,u)

1

2

L

2

2(v,u)

1

2

L

2

N(v,u) L

2

N

v

h

L

2≲(1+t)

-

5

4

×

1

2 (1+t)

-

7

4

×

1

2 (1+t)

-

3+2(N-1)

4 N

v

h

L

2≲

(1+t)

-

3+2N

4 (1+t)

-1 N

v

h

L

2≲(1+t)

-

3+2N

2 +(1+t)

-2 N

v

h 2

L

2 。

因此,

-

d

dt

∫

R3

N-1u

h N

v

h

dx + λ1

N

v

h 2

L

2 ≲

λ1

4

+ (1 + t)

-2

( )

N

v

h 2

L

2 + Nu

h 2

L

2 + N+1u

h 2

L

2 + (1 + t)

-

3+2N

2 。 (27)

步骤 2 式(12)的证明。 在式(18)中取 k =N-1,得

d

dt

N(v,u)

2

L2

+2D N+1 u

2

L

2≲δ(

N

v

2

L

2+ N+1 u

2

L

2 。 (28)

取足够大的 T1 和 D1 ,定义能量泛函

E(t) = D1

N(v,u)

2

L

2 - ∫

R3

N-1 u

h N

v

h

dx。

当 t≥T1 时,由于 D1 充分大,则 E(t)等价于 N(v,u)

2

L

2 。 将式(28)乘以 D1 ,然后将所得结果与式(27)相

加,注意到当 t≥T1 时,由于 δ 充分小,再结合引理 3 的高频性质和引理 6,整理可得

d

dt

E(t)+C1E(t)≲(1+t)

-

3+2N

2 + N(v

l

,u

l

)

2

L

2 。

结合式(22)与 Gronwall 不等式,得

E(t)≲(1+t)

-

3+2N

2 。

故式(25)成立。

综上所述,定理 1 得证。

参 考 文 献

[1] OTHMER H G, STEVENS A. Aggregation, blowup, and collapse: the ABCs of taxis in reinforced random walks[J]. SIAM

Journal on Applied Mathematics, 1997, 57(4): 1044-1081.

[2] LEVINE H A, SLEEMAN B D. A system of reaction diffusion equations arising in the theory of reinforced random walks[J].

SIAM Journal on Applied Mathematics, 1997, 57(3): 683-730.

[3] GUO J, XIAO J X, ZHAO H J, et al. Global solutions to a hyperbolic-parabolic coupled system with large initial data[J].

Acta Mathematica Scientia, 2009, 29(3): 629-641.

[4] XIE W J, ZHANG Y H, XIAO Y D, et al. Global existence and convergence rates for thestrong solutions in H

2

to the 3D

chemotaxis model[J]. Journal of Applied Mathematics, 2013, 2013(1): 1-9.

[5] ZHANG M, ZHU C J. Global existence of solutions to a hyperbolicparabolic system[ J]. Proceedings of the American

Mathematical Society, 2007, 135(4): 1017-1027.

[6] 张映辉, 谭忠, 赖柏顺.一个模拟趋化现象的广义双曲-抛物系统的光滑解的全局分析[ J]. 数学年刊 A 辑, 2012,

33(1): 27-38.

[7] ZHANG Y H, TAN Z, SUN M B. Global existence and asymptotic behavior of smooth solutions to a coupled hyperbolicparabolic system[J]. Nonlinear Analysis: Real World Applications, 2013, 14 (1): 465-482.

[8] ZHANG Y H, DENG H Y, SUN M B. Global analysis of smooth solutions to a hyperbolic-parabolic coupled system[ J].

130

http:∥xuebao.gxnu.edu.cn

Frontiers of Mathematics in China, 2013, 8(6): 1437-1460.

[9] ZHANG Y H, XIE W J. Global existence and exponential stability for the strong solutions in H

2

to the 3-D chemotaxis model

[J]. Boundary Value Problems, 2015, 116: 1-13.

[10] LI D, PAN R H, ZHAO K. Quantitative decay of a one-dimensional hybrid chemotaxis model with large data [ J].

Nonlinearity, 2015, 28(7): 2181-2210.

[11] 涂馨予, 穆春来, 郑攀. 带有流晕限制的拟线性趋化模型解的整体有界性[ J]. 中国科学: 数学, 2021, 51(6):

1003-1012.

[12] 张婕燕, 辛巧, 穆春来. 具有非线性扩散的趋化模型弱解的有界性[ J]. 西南大学学报(自然科学版), 2021,

43(5): 88-95.

[13] 王曦, 刘祖汉, 周玲. 具不同分数阶扩散趋化模型的衰减估计[J]. 数学年刊 A 辑, 2020, 41(2): 175-200.

[14] 张映辉, 谭忠, 孙明保.一个耦合双曲-抛物系统的全局光滑解[J]. 数学年刊 A 辑, 2013, 34(1): 29-46.

[15] LI D, PAN R H, ZHAO K. On a hybrid type chemotaxis model on bounded domains with large data[J]. SIAM Journal on

Applied Mathematics, 2012, 72(7): 417-443.

[16] HAO C C. Global well-posedness for a multidimensional chemotaxis model in critical Besov spaces [ J]. Zeitschrift für

Angewandte Mathematik und Physik, 2012, 63: 825-834.

[17] LI D, LI T, ZHAO K. On a hyperbolic-parabolic system modeling chemotaxis[ J]. Mathematical Models and Methods in

Applied Sciences, 2011, 21(8): 1631-1650.

[18] ZHANG Y H, LI C, XIE W J. Decay of a 3-D hyperbolic parabolic system modeling chemotaxis[J]. Journal of Information

and Optimization Sciences, 2018, 39(7): 1505-1525.

[19] NIRENBERG L. On elliptic partial differential equations [ J]. Annali della Scuola Normale Superiore di Pisa-Classe di

Scienze, 1959, 13(2): 115-162.

[20] MAJDA A J, BERTOZZI A L. Vorticity and incompressible flow[M]. Cambridge: Cambridge University Press, 2002.

[21] JU N. Existence and uniqueness of the solution to the dissipative 2D quasi-geostrophic equations in the Sobolev space[ J].

Communications in Mathematical Physics, 2004, 251: 365-376.

Optimal Time-decay Rates of the Hyperbolic-parabolic

System Modeling Chemotaxis in R

3

WANG Han, ZHANG Yinghui

∗

(School of Mathematics and Statistics, Guangxi Normal University, Guilin Guangxi 541006, China)

Abstract: The large-time behavior of solutions to the Cauchy problem of a 3D hyperbolic-parabolic system

modeling chemotaxis is investigated. The optimal time decay rates of the higher-order spatial derivatives of the

solutions are obtained. Compared with previous results, the main innovation of this paper is to give the highest

order spatial derivative of the solutions which is the same as that of the heat equation. The proof is mainly based

on low-frequency and high-frequency decomposition and delicate energy estimates.

Keywords: hyperbolic-parabolic system; optimal decay rates; low-frequency and high-frequency decomposition;

large-time behavior; chemotaxis

(责任编辑 吴佃华)

131

第 40 卷 第 2 期

2022 年 3 月

广西师范大学学报(自然科学版)

Journal of Guangxi Normal University (Natural Science Edition)

Vol. 40 No. 2

Mar. 2022

DOI: 10.16088 / j.issn.1001-6600. 2021070902 http: xuebao.gxnu.edu.cn

张琬婧, 林支桂. 增长区域上一类寄生虫-宿主模型的 Turing 不稳定[J]. 广西师范大学学报(自然科学版), 2022, 40(2): 132-139. ZHANG

W J, LIN Z G. Turing instability of a parasite-host model on growing domains[J]. Journal of Guangxi Normal University (Natural Science Edition),

2022, 40(2): 132-139.

增长区域上一类寄生虫-宿主模型的 Turing 不稳定

张琬婧, 林支桂∗

(扬州大学 数学科学学院, 江苏 扬州 225002)

摘 要: 为理解区域演化对寄生虫-宿主 Turing 不稳定的影响, 本文以寄生虫-宿主传染病模型为主体, 研究增长区域上

的反应扩散问题, 通过线性化和谱分析给出模型产生 Turing 不稳定的条件, 再利用数值模拟验证理论结果。 结果表明扩

散系数的增加有利于 Turing 斑图的形成, 而区域增长对 Turing 斑图形成起破坏作用。

关键词: 寄生虫-宿主模型; 增长区域; Turing 不稳定; 数值模拟; 传染病

中图分类号: R181; O175 文献标志码: A 文章编号: 1001-6600(2022)02-0132-08

近期 H1N1

[1]

、H7N9

[2]

、登革热[3]等传染病反复发作,新型冠状病毒肺炎[4] 不断流行,传染病的传播

引起流行病学专家以及社会广泛关注。 传染病受季节、温度等因素影响,其传播范围随时间变化,意味着

受传染病影响的区域与时间有关。 考虑到种群在空间上的迁徙,在模型中研究种群密度的时空演化能更

好地描述传染病的变化规律,故而在常微分方程的基础上引入扩散作用,从而得到反应扩散方程。 此类方

程描绘了生物种群随时空变化而产生的一系列规律[5]

。 但是引入扩散作用后,原本模型的一系列稳定性

质就会被打乱,其中 Turing 不稳定就是一个例子:由均一的平衡状态变为不均一状态,这一工作得到了人

们的极大关注[6-9]

。 由此,为理解在增长的栖息地区域中,扩散对物种生存的影响,本文在增长区域下讨

论一类寄生虫-宿主传染病模型的 Turing 不稳定。

1 扩散的寄生虫-宿主传染病模型

2004 年,在 Kermack 等[10]建立的经典仓室模型基础上,Berezovsky 等[11]提出 1 个包含可变人口、接触

传播、因病死亡等因素的寄生虫-宿主模型。

dS

dt

= rN 1-

N

K

( ) -β

SI

N

-(μ+m)S,

dI

dt

= β

SI

N

-(μ+d)I,

ì

î

í

ï

ï

ï

ï

(1)

式中:N 表示物种总数,它包含易感人群 S 和染病人群 I,且满足 N= S+I;r 表示易感人群 S 的内禀增长率;

K 表示易感人群 S 的承载能力;β 表示染病人群 I 的接触传播率;μ 表示自然死亡率;d 表示病死率;m 表示

非染病者的人均移民率。

Wang 等[12]将模型(1)进行无量纲化,并引入扩散

收稿日期: 2021-07-09 修回日期: 2021-08-16

基金项目: 国家自然科学基金(11771381, 11911540464)

通信作者: 林支桂(1965—), 男, 江苏兴化人, 扬州大学教授, 博士。 E-mail: zglin@yzu.edu.cn

http:∥xuebao.gxnu.edu.cn

?S

?t

-d1ΔS = νRd(S+I)(1-S-I)-R0

SI

S+I

-νS, x∈Ω,t>0,

?I

?t

-d2ΔI =R0

SI

S+I

-I, x∈Ω,t>0,

?S

?υ

= 0,

?I

?υ

= 0, x∈?Ω,t>0,

S(x,0)= S0(x)≥0,I(x,0)= I0(x)≥0, x∈Ω,

ì

î

í

ï

ï

ï

ïï

ï

ï

ï

ï

(2)

式中:R0

=

β

μ+d

为基本再生数;Rd

=

r

μ+m

为人口再生数;ν =

μ+m

μ+d

表示易感染者与感染者平均寿命的比值;Ω

是具有光滑边界的有界开区域;υ 表示边界上的单位外法向量;正常数 d1 、d2 为扩散系数;S0(x)、I0(x)为

非负光滑且不恒为零的函数。

由偏微分方程理论,易知系统(2)的解存在唯一,且当 Rd >1 时,系统(2)存在一个无病平衡点 E0

=

(S

Δ

,I

Δ

)= 1-

1

Rd

( ,0) ;当 Rd >R

∗

d : =

R0

+v-1

vR0

,且 R0 >1 时,系统(2) 存在一个染病平衡点 E

∗ = ( S

∗

,I

∗

),

式中:

S

∗ =

vR0Rd

-R0

+1-v

vR

2

0Rd

; I

∗ = (R0

-1)

vR0Rd

-R0

+1-v

vR

2

0Rd

= (R0

-1)S

∗

。

其对应的常微分方程的正平衡解所对应的 Jacobi 矩阵为

J =

-

R

2

0

+νR0Rd

+νR0

-4R0

-2ν+3

R0

2R0

+2ν-νR0Rd

-3

R0

(R0

-1)

2

R0

1

R0

-1

æ

è

ç

ç

ç

ç

ç

ö

ø

÷

÷

÷

÷

÷

△=

a11 a12

a21 a22

( ) 。

下面讨论系统(2)解的 Turing 失稳条件。

定理 1 当 1<R0<2 时,如果

Rd>

1

R0(2-R0 )

,

R0

-1

R0Rd

-1

<ν<

(R0

-1)(3-R0 )

R0

-2+R0Rd

,

(d1 a22

+d2 a11 )

2-4d1 d2(a11 a22

-a12 a21 )>0,

且存在某个 δi>0,使得不等式 0<k1<δi<k2 成立,那么系统(2)在 E

∗ 处 Turing 不稳定,式中:δi 是 Ω 中带有

齐次 Neumann 边界条件的-Δ 算子的特征值;

k1

=

d1 a22

+d2 a11

- (d1 a22

+d2 a11 )

2-4d1 d2 det(J)

2d1 d2

;

k2

=

d1 a22

+d2 a11

+ (d1 a22

+d2 a11 )

2-4d1 d2 det(J)

2d1 d2

。

证明 易知,当 R0>1,Rd>

R0

-1+ν

vR0

时,E

∗对系统(2)所对应的常微分系统稳定。

下面设 0 = δ0<δ1<δ2<…是 Ω 中带有齐次 Neumann 边界条件的-Δ 算子的特征值,在 C

1

(?Ω)中 δi 对应

的特征空间为 E(δi),设

X = {U= (φ,ψ)∈C

1

(?Ω)×C

1

(?Ω) | ∂ηφ= ∂ηψ= 0,x∈∂Ω},

{φij,j = 1,2,…,dim E(δi)}是 E(δi)的一组标准正交基,

Xij

= c·φij

| c∈R

2

{ } ,Xi

=?dimE(δI

)

j = 1 Xij,X=?+∞

i = 1Xi。

μ 是 Xi 上的特征值当且仅当 μ 是矩阵

Ai

=

-d1

δi

+a11 a12

a21

-d2

δi

+a22

( )

133

广西师范大学学报(自然科学版),2022,40(2)

的特征值。 此时特征方程为

μ

2+[(d1

+d2 )δi

-a11

-a22 ]μ+[d1 d2

δ

2

i

-(d1 a22

+d2 a11 )δi

+a11 a22

-a12 a21 ] = 0。

显然,要证明 E

∗对于系统(2)不稳定,只需要 det(Ai)<0,也就是特征方程至少有 1 个正解。 因为 δi >0,所

以不稳定第一个必要条件就是 d1 a22

+d2 a11>0,这等价于 a11>-d1 a22

/ d2>0。

要使 a11

= -

R

2

0

+vR0Rd

+vR0

-4R0

-2v+3

R0

>0,即要 ν<

(R0

-1)(3-R0 )

R0

-2+R0Rd

,R0 <2,Rd >

1

R0(2-R0 )

。 下面证明要

使特征方程存在 1 个正解,则必须满足

(d1 a22

+d2 a11 )

2-4d1 d2(a11 a22

-a12 a21 )>0。

综上所述,要使系统(2)在平衡解 E

∗处 Turing 不稳定,那么 Rd 、R0 、ν 必须满足

1<R0<2,Rd>

1

R0(2-R0 )

,

R0

-1

R0Rd

-1

<ν<

(R0

-1)(3-R0 )

R0

-2+R0Rd

,

(d1 a22

+d2 a11 )

2-4d1 d2(a11 a22

-a12 a21 )>0。

事实上,此情形还不能保证一定存在 δi>0,使得

h(δi)△=d1 d2

δ

2

i

-(d1 a22

+d2 a11 )δi

+a11 a22

-a12 a21<0。 (3)

所以直接计算

h(δi)= δi(d1

δi

-a11 ) d2

-

M-d1 a22

δi

δi(a11

-d1

δi)

( ) ,

式中 M= det J = a11 a22

-a12 a21 。

进一步,假设 a11>0,且 d2 >Di: =

M-d1 a22

δi

δi(a11

-d1

δi)

成立,h(δi) <0。 于是特征方程至少有 1 个正根,故平衡

解 E

∗对于系统(2)不稳定。 如果对所有 δi>0 都有 di<Di,那么平衡解 E

∗对于系统(2)稳定。

在满足上述条件的情况下,方程 h(δi)= 0 存在 2 个正解:

k1

=

d1 a22

+d2 a11

- (d1 a22

+d2 a11 )

2-4d1 d2(a11 a22

-a12 a21 )

2d1 d2

,

k2

=

d1 a22

+d2 a11

+ (d1 a22

+d2 a11 )

2-4d1 d2(a11 a22

-a12 a21 )

2d1 d2

。

从而有结论:如果存在 δi>0,使不等式 0<k1<δi <k2 成立,那么系统(2)在染病平衡点 E

∗ 处 Turing 不稳定。

证毕。

注 1 假设 Ω = (0,lπ),容易计算出( 0,lπ) 上具有齐次 Neumann 边界条件的-Δ 算子特征值是

i-1

l

( )

2

,特征函数是 cos

i-1

l

x,i = 1,2,3,…。 由此得 k1<

i-1

l

( )

2

<k2 ,从而,

l k1 <i-1<l k2 ,

最后只需要满足(k1

+k2 )-2 k1

k2 >1 / l

2

。

定理 2 假设 E

∗是系统(2)的染病平衡点,如果 Ω= (0,lπ),满足

Rd>

1

R0(2-R0 )

,

R0

-1

R0Rd

-1

<ν<

(R0

-1)(3-R0 )

R0

-2+R0Rd

,

d1(1-R0 )-d2(R

2

0

+vR0Rd

+vR0

-4R0

-2v+3)-2 d1 d2(vR0Rd

-R0

+1-v)(R0

-1) >

R0 d1 d2

l

2

,

那么系统(2)在平衡解 E

∗处 Turing 不稳定。

故而当 i→∞ 时,δi→∞ ,并且特征值 δi

{ }

+∞

i = 1会随区域 Ω 的大小变化而相应连续发生变化。 当区域 Ω

134

http:∥xuebao.gxnu.edu.cn

足够大时,不等式 0<k1<δi<k2 自然满足,从而定理 1 自然成立;当区域 Ω 充分小时,不等式 0<k1 <δi <k2 不

能满足,即定理 1 不成立。 因此,空间大小会影响平衡解的稳定性。

2 增长区域上的寄生虫-宿主传染病模型

第 1 章通过线性化方法给出了固定区域下寄生虫-宿主传染病模型的 Turing 不稳定的条件,结果表

明,扩散会导致此模型的 Turing 不稳定。 本章将考虑增长区域上具 Neumann 边界条件的寄生虫-宿主模

型。 记

f(S,I)= νRd(S+I)(1-S-I)-R0

SI

S+I

-νS,

g(S,I)= R0

SI

S+I

-I。

ì

î

í

ï

ï

ï

ï

运用 Plaza 等[13]建立的框架,将区域的增长率和曲率融入传染病动力学模型,即设二维平面 X 嵌入三维空

间 R

3 中,表示为

X(ξ,η,t)= (x(ξ,η,t),y(ξ,η,t),z(ξ,η,t)),

并定义

h1

= Xξ , h2

= Xη , (4)

式中 · 是 Euclidean 范数,下标表示偏导数。 由此,系统(2)转化为

?S

?t

=

d1

h1 h2

h2

h1

( Sξ )

ξ

+

h1

h2

( Sη )

η

é

ë

ê

ê

ù

û

ú

ú

-?t(ln(h1 h2 ))S+f(S,I), X(t)∈Ω(t), t>0,

?I

?t

=

d2

h1 h2

h2

h1

I ( ξ )

ξ

+

h1

h2

I ( η )

η

é

ë

ê

ê

ù

û

ú

ú

-?t(ln(h1 h2 ))I+g(S,I), X(t)∈Ω(t), t>0,

?S

?υ

= 0,

?I

?υ

= 0, X(t)∈?Ω(t), t>0,

S(X(0),t)= S0(X(0))≥0,I(X(0),t)= I0(X(0))≥0, X(t)∈Ω(0)。

ì

î

í

ï

ï

ï

ï

ï

ï

ï

ï

ï

ï

(5)

式中:S0(X(0))和 I0(X(0))是正有界函数;Ω(0)为初始区域。 由此得到增长区域下具 Neumann 边界条

件下的寄生虫-宿主传染病模型。

为刻画区域演化率的影响,考虑栖息地的演化是各向同性的。 首先探索限制在平面区域上增长时的

情形,数学上表达式为

X(ξ,η,t)= ρ(t)Y(ξ,η)= ρ(t)[ξ,η,0]

T

,

式中 ρ(t)是增长函数,满足 ρ(t)在[0,+∞ )上连续可微,

ρ(0)= 1,

̇ρ(t)≥0,lim

t→∞

ρ(t)= ρ∞ >1,lim

t→∞

̇ρ(t)= 0。

经计算,得

h

2

1

= Xξ

2 = ρ

2

(t),h

2

2

= Xη

2 = ρ

2

(t),h1

/ h2

= 1。

由此,式(5)转化为平面增长区域下的寄生虫-宿主传染病扩散模型

?S

?t

=

d1

ρ

2

(t)

(Sξξ

+Sηη )+f(S,I)-

2

̇ρ(t)

ρ(t)

S, (ξ,η)∈Ω(0),t>0,

?I

?t

=

d2

ρ

2

(t)

(Iξξ

+Iηη )+g(S,I)-

2

̇ρ(t)

ρ(t)

I, (ξ,η)∈Ω(0),t>0,

?S

?υ

= 0,

?I

?υ

= 0, (ξ,η)∈?Ω(0),t>0,

S(ξ,η,0)= S0(ξ,η)≥0,I(ξ,η,0)= I0(ξ,η)≥0, (ξ,η)∈Ω(0)。

ì

î

í

ï

ï

ï

ï

ï

ï

ï

ï

ïï

(6)

135

广西师范大学学报(自然科学版),2022,40(2)

3 增长区域上的 Turing 不稳定

本章将通过线性化和谱分析得出正常数平衡点的稳定性,从而给出此情形下所需要的 Turing 不稳定

条件。 为具体展现区域演化对 Turing 不稳定的影响,取增长函数 ρ( t) = e

kt

,则式(6) 中的稀释项即为

̇ρ(t) / ρ(t)= k。 为简化起见,记

?f(S,I)= f(S,I)-2kS, g?(S,I)= g(S,I)-2kI, (7)

从而式(6)变为

?S

?t

=

d1

ρ

2

(t)

(Sξξ

+Sηη )+?f(S,I), (ξ,η)∈Ω(0),t>0,

?I

?t

=

d2

ρ

2

(t)

(Iξξ

+Iηη )+g?(S,I), (ξ,η)∈Ω(0),t>0,

?S

?υ

= 0,

?I

?υ

= 0, (ξ,η)∈?Ω(0),t>0,

S(ξ,η,0)= S0(ξ,η)≥0,I(ξ,η,0)= I0(ξ,η)≥0, (ξ,η)∈Ω(0)。

ì

î

í

ï

ï

ï

ï

ï

ï

ï

ï

ïï

(8)

式(8)对应的动力系统为

dS

dt

=?f(S,I),

dI

dt

= g?(S,I),

S(0)= S0≥0, I(0)= I0≥0。

ì

î

í

ïï

ï

(9)

定理 3 当 R0>1+2k,Rd >

(R0

+ν-1)(1+2k)

νR0

时,式(9) 对应的正平衡点为 E

Δ = ( S

Δ

,I

Δ

),式中:S

Δ =

(1+2k) vR0Rd

-(R0

[ -1+v)(1+2k) ]

vR

2

0Rd

;I

Δ =

(R0

-1-2k)S

Δ

1+2k

,并且正平衡点 E

Δ 局部渐近稳定。

证明 为了得到具有每种分支类型(Hopf、Turing 和 Turing-Hopf)下的参数值的条件,根据文献[14],

先做下列计算。 式(7)对应的 ODE 系统的零斜率线为

I1(S)= -

νRdR

2

0

(1+2k)

3

S

2+

νRdR0

(1+2k)

2

-

v+2k

1+2k ( ) S 和 I2(S)=

R0 S

1+2k

-S。 (10)

I1 和 I2 相交所满足的等式为

νRdR

2

0

(1+2k)

3

S

2-

νRdR0

(1+2k)

2

-

ν+2k+R0

1+2k

+1 ( ) S = 0。 (11)

因此,给出条件

νRdR0

1+2k

-v-R0

+1>0 和 R0>1+2k, (12)

从而可以保证式(11)中 S

Δ

>0,I

Δ

>0。 故 E

Δ = (S

Δ

,I

Δ

)是正平衡解。

其次,记(S

Δ

,I

Δ

)处的 Jacobi 矩阵为 J(E

Δ

),这里

J(E

Δ

)=

?f

S?f

I

g?S g?I

( ) , (13)

式中:

?f

S

= -

R

2

0

+νR0Rd

+νR0

-(4+6k)R0

-2(1+2k)v+(1+2k)(2k+3)

R0

;

?f

I

=

-νR0Rd

+(2R0

+2ν-3-2k)(1+2k)

R0

;

g?S

=

(R0

-1-2k)

2

R0

, g?I

=

(1+2k)

R0

2

-1-2k<0。

136

http:∥xuebao.gxnu.edu.cn

考虑其特征方程

μE-J(E

Δ

) = μ

2-(?f

S

+g?I)μ+?f

S g?I

-?f

Ig?S

= 0, (14)

经计算,

?f

S

+g?I

= -

(R0

-1-2k)(R0

-1+ν)+νR0Rd

-(R0

-1+ν)(1+2k)

R0

,

?f

S g?I

-?f

Ig?S

= detJ(E

Δ

)=

(R0

-1-2k) (1+2k)(-R0

-ν+1)+vR0Rd

[ ]

R0

。

由条件 R0>1+2k,Rd>

(R0

+v-1)(1+2k)

vR0

,得 ?f

S

+g?I<0 且?f

S g?I

-?f

Ig?S >0,故特征值 μ 具有负实部。 根据常微分

方程稳定性理论,平衡解(S

Δ

,I

Δ

)对式(9)局部渐近稳定。 证毕。

类似定理 1 的证明,通过线性化得特征矩阵,讨论特征方程,给出其至少存在 1 个正根的必要条件,即

得系统产生 Turing 不稳定的必要条件。

定理 4 式(8)产生 Turing 不稳定的必要条件为

f

S

+gI

-4k<0,

f

S gI

-f

IgS

-2k(f

S

+gI)+4k

2

>0,

(d2

f

S

+d1 gI)-2k(d1

+d2 )>0,

- (d2

f

S

+d1 gI)-2k(d1

[ +d2 ) ]

2+4d1 d2

f

S gI

-f

IgS

-2k(f

S

+gI)+4k

2

[ ] <0。

把 f、g 的表达式代入计算可得

1+2k<R0<2+4k, (15)

Rd>

(1+2k)

2

R0(2-R0

+2k)

, (16)

(R0

-1)(1+2k)

R0Rd

-1-2k

<v<

(R0

-1-2k)(1+2k)+(-R0

+4+4k)(R0

-1)

R0

-2-4k+R0Rd

, (17)

- -

d2 [R

2

0

+νR0Rd

+νR0

-(4+4k)R0

-2(1+2k)ν+(1+2k)(2k+3)]

R0

+d1

(1+2k)

2

R0

( -1-4k)

é

ë

ê

ê

ù

û

ú

ú

2

+

4d1 d2

(R0

-1-2k) [R0((1+2k)(4k+3-ν-3R0 )+νR0Rd )+(1+2k)8kν]

R

2

0

é

ë

ê

ê

ù

û

ú

ú

<0。 (18)

由此,通过分析发生 Turing 不稳定的必要条件,可以发现:扩散系数 d2 越大,模型更容易出现斑图模

式;而区域增长快,即 k 大时,上述必要条件不容易满足,所以对 Turing 斑图起破坏作用。

4 数值模拟及解释

本章将通过数值模拟来研究固定区域上 Turing 不稳定在此类寄生虫-宿主模型中的动力学行为。 在

实际流行病学中,空间上的扩散影响物种初始活动空间范围[15]

。 为研究模型非平凡解的动态分布,选取

参数 R0

= 1.3,Rd

= 1.7,ν = 0.3,d1

= 0.1。

对时间和空间都采用有限差分方法,且取时间步长 Δt = 0.01,空间步长 Δx = 1.5,最终时间 t = 1 000,初

始条件是种群初值的随机分布。

由图 1 可以得出,当扩散系数 d2 越小,即(d1 a22

+d2 a11 )

2-4d1 d2(a11 a22

-a12 a21 )<0 时,种群最后会收敛

到一个平衡解,不发生 Turing 不稳定;当扩散系数 d2 较大时,Turing 不稳定的充分条件满足,因此种群在

不同位置的值不同,于是常数解不稳定。 图 2 则是产生 Turing 不稳定时种群 S 和 I 所产生的斑图。

137

广西师范大学学报(自然科学版),2022,40(2)

图 1 不同扩散系数下种群 S 在位置(25,25)、(50,50)、(75,75)时的密度

Fig. 1 Densities of population S at locations (25, 25), (50, 50) and (75, 75)

with different diffusion coefficients

图 2 d2

= 5 时种群 S 和 I 所形成的斑图

Fig. 2 Spatial patterns of population S and I when d2

= 5

参 考 文 献

[1] NISHIURA H, WILSON N, BAKER M. Estimating the reproduction number of the novel influenza A virus(H1N1) in a

Southern Hemisphere setting: preliminary estimate in New Zealand [ J ]. The New Zealand Medical Journal, 2009,

122(1299): 73-77.

[2] ZHOU J F, WANG D Y, GAO R B, et al. Biological features of novel avian influenza A (H7N9) virus[J]. Nature, 2013,

499(7459): 500-503.

[3] GUZMAN M G, HALSTEAD S B, ARTSOB H, et al. Dengue: a continuing global threat[J]. Nature Reviews Microbiology,

2010, 8(12 Suppl): S7-S16.

[4] 周涛, 刘权辉, 杨紫陌, 等. 新型冠状病毒肺炎基本再生数的初步预测[ J]. 中国循证医学杂志, 2020, 20( 3):

359-364.

[5] KANG H, RUAN S G. Nonlinear age-structured population models with nonlocal diffusion and nonlocal boundary conditions

[J]. Journal of Differential Equations, 2021, 278: 430-462.

[6] 丁亚君, 张存华. Lengyel-Epstein 反应扩散系统的稳定性和 Turing 不稳定性[ J]. 高校应用数学学报 A 辑, 2018,

33(3): 272-280.

[7] 林支桂. 数学生态学导引[M]. 北京: 科学出版社, 2013: 94-107.

[8] ZUÑIGA-GALINDO W A. Reaction-diffusion equations on complex networks and Turing patterns, via p-adic analysis[ J].

138

http:∥xuebao.gxnu.edu.cn

Journal of Mathematical Analysis and Applications, 2020, 491(1): 124239.

[9] 张道祥, 赵李鲜, 胡伟.一类三种群食物链模型中交错扩散引起的 Turing 不稳定[ J]. 山东大学学报(理学版),

2017, 52(1): 88-97.

[10] KERMACK W O, McKendrick A G. A contribution to the mathematical theory of epidemics[J]. Proceedings of the Royal

Society A: Mathematical, Physical and Engineering Sciences, 1927, 115(772): 700-721.

[11] BEREZOVSKY F, KAREV G, SONG B, et al. A simple epidemic model with surprising dynamics [ J]. Mathematical

Biosciences and Engineering, 2005, 2(1): 133-152.

[12] WANG W M, CAI Y L, WU M J, et al. Complex dynamics of a reaction-diffusion epidemic model[J]. Nonlinear Analysis:

Real World Applications, 2012, 13(5): 2240-2258.

[13] PLAZA R G, SÁNCHEZ-GARDUÑO F, PADILLA P, et al. The effect of growth and curvature on pattern formation [ J].

Journal of Dynamics and Differential Equations, 2004, 16(4): 1093-1121.

[14] SANCHEZ-GARDUNO F, KRAUSE A L, CASTILLO J A, et al. Turing-Hopf patterns on growing domains: The torus and

the sphere[J]. Journal of Theoretical Biology, 2019, 481: 136-150.

[15] 王玮明, 蔡永丽. 生物数学模型斑图动力学[M]. 北京: 科学出版社, 2020: 8-11.

Turing Instability of a Parasite-host Model on Growing Domains

ZHANG Wanjing, LIN Zhigui

∗

(School of Mathematical Science, Yangzhou University, Yangzhou Jiangsu 225002, China)

Abstract: In order to study the influence of growing domain on Turing instability of the parasite-host model, this

paper takes the parasite-host infectious disease model as the main body and analyzes its Turing instability

conditions through linearization and spectral analysis. Numerical simulations are used to verify the theoretical

results. The results show that the increase of the diffusion coefficient is benefitial for formation of Turing pattern,

but regional growth has a destructive effect on Turing pattern formation.

Keywords: parasite-host model; growing domain; Turing instability; numerical simulation; infectious diseases

(责任编辑 吴佃华)

139

第 40 卷 第 2 期

2022 年 3 月

广西师范大学学报(自然科学版)

Journal of Guangxi Normal University (Natural Science Edition)

Vol. 40 No. 2

Mar. 2022

DOI: 10.16088 / j.issn.1001-6600.2021041901 http: xuebao.gxnu.edu.cn

刘奇文,李丹,黄小芳,等. 纳米金催化甲酸还原磷钼酸耦合电置换反应-共振瑞利散射测定痕量汞[ J]. 广西师范大学学报(自然科学版),

2022, 40(2): 140-148. LIU Q W, LI D, HUANG X F, et al. A new strategy for the determination of trace mercury by resonance Rayleigh scattering

method based on nano-gold catalytic amplification and galvanic replacement reaction-phosphomolybdic acid[J]. Journal of Guangxi Normal University

(Natural Science Edition), 2022, 40(2): 140-148.

纳米金催化甲酸还原磷钼酸耦合电置换反应共振瑞利散射测定痕量汞

刘奇文1,2

, 李 丹1,2

, 黄小芳1,2

, 梁爱惠1,2∗

, 蒋治良1,2∗

(1. 珍稀濒危动植物生态与环境保护教育部重点实验室(广西师范大学), 广西 桂林 541006;

2. 广西环境污染控制理论与技术重点实验室(广西师范大学), 广西 桂林 541006)

摘 要: 在 pH= 3.1 的 HCOOH-HCOONa 缓冲液中, 磷钼酸粒子在 450 nm 处产生一个共振瑞利散射(RRS)峰。 金纳米粒

子(AuNPs)可催化磷钼酸-甲酸反应生成磷钼蓝, 使得 450 nm 处磷钼酸的 RRS 强度线性降低。 Hg

2+可与 AuNPs 发生电置

换反应, 从而抑制 AuNPs 的催化作用, RRS 峰增强。 在 2.5×10

-4

~ 3.5 μmol / L, 随着 Hg

2+浓度的增加, AuNPs 的催化作

用逐渐减弱, 反应液的颜色逐渐从蓝色变为无色, 体系在 450 nm 处的 RRS 峰值(ΔI) 线性增高, 其线性方程为 ΔI =

0.32C+46.1, 检出限为 0.18 nmol / L。 该法用于废水中 Hg

2+的检测, 结果令人满意。

关键词: 汞; 纳米金催化; 静电取代; 磷钼蓝; 共振瑞利散射

中图分类号: O657.3 文献标志码: A 文章编号: 1001-6600(2022)02-0140-09

电置换反应(GR)指一种较活泼金属被一种较不活泼的金属侵蚀,反应较快,常在几分钟内发生[1]

。

由于纳米表面的高催化活性,纳米粒子表面的 GR 反应更快。 GR 反应由于具有高度的可调性而引起人们

的极大兴趣,可用于研究金属纳米结构表面错综复杂的合金化和去合金化反应[2]

;GR 反应还提供一种非

常简易、多功能路径制备可控空心和多孔壁纳米结构材料方法[3]

。 Sun 等[4] 采用 Ag 作为模板金属,和

Au

3+发生氧化还原反应,制备了 Au 中空结构,并对其反应机理进行了解释。 Liu 等[5] 利用超小 ssDNA-模

板化银纳米簇与 Cu(II)发生拮抗电置换反应制备了 Ag-Cu 合金,并用光散射技术检测反应的发生。 Bi

等[6]通过 Ag / AgCl 核壳纳米线和 H2PtCl

6发生电置换反应得到较高电催化活性的纳米管。 近年来,GR 反

应在纳米分析检测技术中的应用也较活跃[7-10]

。 Netzer 等[7] 将 GR 反应与银纳米线 Langmuir-Blodgett 膜

技术有机结合成功地制备了高活性 SERS 基底,基于 1 394 cm

-1处的灵敏拉曼峰可用于检测 8 nmol / L 4-氨

基硫酚。 共振瑞利散射(RRS)不仅是一种简便灵敏的分子光谱分析技术,也是一种研究纳米微粒反应的

灵敏光谱技术[11-13]

,纳米催化放大是提高分析方法灵敏度的重要途径之一[14]

。 RRS 与纳米催化技术结合

构建了重金属环境污染物 RRS 分析新方法。 Wang 等[15]以掺钯共价有机骨架(TpPaPd)为催化剂,催化次

磷酸钠还原 Ni(II)形成 Ni-P 合金,该产物具有 RRS 信号,将纳米催化反应与基于 Pb

2+的 DNA 酶反应相

结合,建立了检测 0.001~0.100 nmol / L Pb

2+的 RRS 分析方法,其检出限为 0.4 pmol / L。 Zhang 等[16]利用掺

金碳点(CDAu )对 AgNO3

-葡萄糖反应具有强催化,其产物银纳米粒子(AgNPs)具有 SPR 效应。 结合 Apt

特异性反应与纳米催化放大信号,建立了 RRS 法检测水中痕量 As

3+

,线性范围为 0.025~0.750 μg / L ,检出

限为 0.01 μg / L。 但尚未见基于纳米金催化甲酸还原磷钼酸耦合金纳米粒子(AuNPs)-汞离子电置换反

应,RRS 测定痕量汞的研究报道。

收稿日期: 2021-04-19 修回日期: 2021-04-29

基金项目: 国家自然科学基金(21767004)

通信作者: 梁爱惠(1965—), 女, 广西岑溪人, 广西师范大学研究员。 E-mail: ahliang2008@163.com

蒋治良(1965—), 男, 广西全州人, 广西师范大学教授, 博导。 E-mail: zljiang@gxnu.edu.cn

http:∥xuebao.gxnu.edu.cn

汞离子是重金属中最重要的污染物之一,其不可生物降解,容易通过食物链在鱼类和贝类中积累[17]

。

此外,汞离子在自然环境中可以转化为高毒性甲基汞,比 Hg

2+更容易穿过细胞膜,对脑组织、肾脏、神经系

统、泌尿系统和内分泌系统造成严重损害[18]

。 鉴于 Hg

2+在环境和生物系统中的毒性,因此,开发一种快

速、高效、灵敏检测痕量 Hg

2+的方法至关重要。 最近,已有学者给出了一些用于检测 Hg

2+的传感平台,包

括比色法[19]

、电化学法[20]

、荧光法[21]和 RRS 法[22]

。 在这些方法中,RRS 法因其高灵敏度、高选择性、操

作简单和响应速度快等独特优势而备受关注。 Zhang 等[23] 合成了二维纳米复合物 rGO/ PEI/ Pd,借助其

过氧化物酶活性,开发了一种快速、高度选择性和超痕量肉眼比色法检测水溶液中 Hg

2+的通用策略。 主

要基于在汞离子存在下,rGO/ PEI/ Pd 可以促进 3,3’,5,5’-四甲基联苯胺(TMB)的有效氧化,使其颜色变

为肉眼和吸收光谱法可检测到的深蓝色,最低检测浓度为 0.39 nmol / L Hg

2+

。 Tan 等[24]制备并表征了电化

学衍生的还原型氧化石墨烯化学阻抗传感器,该传感器对 Hg

2+表现出选择性响应,将其用于水样中 Hg

2+

的检测,最低检测浓度为 0.5 nmol / L。 Yu 等[25]以金华佛手柑为碳源,水热法制备了具有较高光致发光性

水溶性荧光碳点,基于 Hg

2+会导致碳点发生荧光猝灭,建立一个荧光检测 0.01~100 μmol / L Hg

2+的分析方

法,其检出限为 5.5 nmol / L。 Ngernpimai 等[26]利用 Hg

2+可诱导单链 DNA 的构象变化,进一步导致金纳米

棒(GNR)聚集,引起 RRS 强度的增强,建立共振瑞利散射法检测 Hg

2+的策略,其检出限为 0.23 nmol / L。

但上述方法有的灵敏度欠佳,有的分析过程复杂。

本文利用 AuNPs 催化磷钼酸-甲酸反应,结合 AuNPs 与 Hg

2+的 GR 反应,以磷钼酸颗粒为 RRS 指示

剂,建立一种简便、灵敏的检测 Hg

2+的 RRS 分析新方法,检出限为 0.18 nmol / L。

1 实验部分

1.1 主要仪器和试剂

日立 F-7000 荧光分光光度计(日立高新技术公司);TU-1901 型双光束紫外可见分光光度计(北京普

析通用仪器有限责任公司);SYZ-550 型石英亚沸蒸馏水器(江苏晶玻仪器厂);KQ3200DB 型数控超声波

清洗器(昆山市超声仪器有限公司,功率 150 W,工作频率 40 kHz);pH 计(梅特勒-托利多仪器上海有限公

司);Nano-23 型纳米粒度与电位分析仪(英国 Malvern 公司);FEI Quanta 200-场发射扫描电子显微镜(FEI

公司)。

5.0 mmol / L 磷钼酸(中国菱湖化工试剂厂,分析纯);0.1 mol / L 甲酸、质量分数 0.05% NaBH4(汕头市

西陇化工厂);0.1 mol / L NaOH(天津市百世化工有限公司);质量分数 1%柠檬酸三钠(汕头市西陇化工

厂);0.01 mol / L HAuCl

4(国药集团化学试剂公司);0.01 mol / L HgSO4(汕头市西陇化工厂)储备液;pH3.1

HCOOH-HCOONa 缓冲液:10 mL 离心管中加入 2.0 mL 0.1 mol / L NaOH 溶液和 8.0 mL 0.1 mol / L HCOOH

溶液,充分混合,得浓度为 5.0 mmol / L 的 HCOOH-HCOONa 缓冲溶液(浓度以醋酸根计)。

金纳米粒子(AuNPs):常温下,量取 40.0 mL 二次水于锥形瓶,在搅拌下加入 0.5 mL 质量分数 1%

HAuCl

4 溶液和 3.5 mL 质量分数 1%柠檬酸钠水溶液,充分混合后,在搅拌下缓慢滴加 4.0 mL 质量分数

0.05%的 NaBH4 溶液,加完后继续搅拌 10 min,最后定容至 50.0 mL,其浓度为 58 mg / L AuNPs。

1.2 实验方法

在 5.0 mL 的刻度试管中,依次移取 170 μL 5 mmol / L 磷钼酸、500 μL pH = 3.1 的 HCOOH-HCOONa 缓

冲液、50 μL 58 mg / L 的 AuNPs 及一定浓度的 Hg

2+

,定容至 2.0 mL,混匀。 80 ℃ 水浴 35 min 后,冰水终止

反应。 取溶液于石英皿内,在 400 V,激发和发射狭缝 10 nm,用荧光分光光度计扫描,获得共振瑞利散射

光谱;不加 Hg

2+做空白,测定溶液 450 nm 处的 RRS 峰值,计算 ΔI = I450 nm

-(I450 nm )0 。

2 结果与讨论

2.1 分析原理

在实验条件下,磷钼酸具有共振瑞利散射信号,可作为该分析反应的指示组分。 AuNPs 可以催化磷

钼酸与甲酸反应,生成磷钼蓝而吸收瑞利散射光,使得体系的 RRS 信号降低。 Hg

2+与 AuNP 发生电置换反

141

广西师范大学学报(自然科学版),2022,40(2)

应生成稳定的金汞齐(AuHg)而覆盖在其表面,抑制了 AuNPs 的催化活性,使得生成的磷钼蓝降低,导致

体系的 RRS 信号增强,据此可以建立测定 Hg

2+的共振瑞利散射新方法(图 1)。

图 1 纳米金催化甲酸还原磷钼酸耦合电置换反应-RRS 测定汞原理

Fig. 1 Analytical principle of the RRS measurement of Hg

2+

by coupling AuNP-phosphomolybdic

acid-formic acid catalytic amplification with GR reaction of AuNP-Hg(Ⅱ)

2.2 RRS 光谱

在 pH = 3.1 的 HCOOH-HCOONa 缓冲溶液及常温下,磷钼酸与甲酸不反应,即使在沸水加热条件下反

应也极难进行。 AuNPs 可催化磷钼酸-甲酸反应生成磷钼蓝,随着 AuNPs 浓度的增加,反应液的颜色逐渐

从无色变为蓝色,体系在 450 nm 处的共振瑞利散射峰逐渐降低(图 2A),系生成蓝色的磷钼蓝对 RRS 光

吸收所致。 Hg

2+可与 Au 发生电置换反应,生成的金汞齐覆盖在 AuNPs 表面抑制 AuNPs 的催化作用,在

2.5×10

-4

~3.5 μmol / L 随着 Hg

2+浓度的增加,AuNPs 的催化作用逐渐减弱,反应液的颜色逐渐从蓝色变为

无色,体系在 450 nm 处的共振瑞利散射峰线性升高(图 2B)。

a: 0.425 mmol / L 磷钼酸+5.0 mmol / L HCOOHHCOONa; b: a+0.29 mg / L AuNPs; c: a+0.58

mg / L AuNPs; d: a+ 1.16 mg / L AuNPs; e: a+

1.45 mg / L AuNPs; f: a+2.03 mg / L AuNPs; g:

a+2.61 mg / L AuNPs

a: 0. 425 mmol / L 磷钼酸 + 5. 0 mmol / L HCOOHHCOONa+ 1. 45 mg / L AuNPs; b: a + 0. 25 nmol / L

Hg

2+

; c: a + 5. 0 nmol / L Hg

2+

; d: a + 1. 0 μmol / L

Hg

2+

; e: a + 1. 5 μmol / L Hg

2+

; f: a + 2. 5 μmol / L

Hg

2+

; g: a+3.5 μmol / L Hg

2+

图 2 磷钼酸-HCOOH-HCOONa-AuNPs-Hg

2+体系 RRS 光谱

Fig. 2 RRS spectra of phosphomolybdic acid-HCOOH-HCOONa-AuNPs-Hg

2+

system

142

http:∥xuebao.gxnu.edu.cn

2.3 紫外吸收光谱

由图 3A 可知,AuNPs 可催化磷钼酸-甲酸反应,其吸光度随着 AuNPs 浓度的增加线性增强。 随着

Hg

2+浓度的增加,AuNPs 对磷钼酸-甲酸反应的催化作用逐渐减弱,生成的磷钼蓝逐渐减少,在 1.0×10

-2

~

3.5 μmol / L 范围内,随着 Hg

2+浓度的增大,体系在 700 nm 处的 Abs 信号线性降低(图 3B),此结论与 RRS

光谱结果一致。 虽然紫外吸收光谱也可用于汞离子测定,但不及 RRS 灵敏。

2.4 扫描电镜(SEM)

图 4A 是磷钼酸的扫描电镜图,由图可见其颗粒均匀分散,这是由于磷钼酸分子具有较强的疏水性并

聚集形成粒径约为 8 μm 的颗粒。 随着 Hg

2+的加入,AuNP-磷钼酸-甲酸钠催化反应逐渐被抑制,但磷钼酸

颗粒的形状变化不大(图 4B)。

a: 0. 425 mmol / L 磷钼酸+ 5. 0 mmol / L HCOOHHCOONa; b: a+0.29 mg / L AuNPs; c: a+0.58 mg /

L AuNPs; d: a+0.87 mg / L AuNPs; e: a+1.16 mg /

L AuNPs; f: a+1.45 mg / L AuNPs; g: a+2.03 mg /

L AuNPs; h: a+2.32 mg / L AuNPs; i: a+2.61 mg /

L AuNPs

a: 0.425 mmol / L 磷钼酸+ 5. 0 mmol / L HCOOHHCOONa+ 1. 45 mg / L AuNPs; b: a + 10 nmol / L

Hg

2+

; c: a+ 0.2 μmol / L Hg

2+

; d: a + 1. 0 mmol / L

Hg

2+

; e: a + 1. 5 mmol / L Hg

2+

; f: a + 2. 5 mmol / L

Hg

2+

; g: a+3.5 mmol / L Hg

2+

图 3 磷钼酸-AuNPs-HCOOH-HCOONa-Hg

2+体系吸收光谱

Fig. 3 Absorption spectra of phosphomolybdic acid-AuNPs-HCOOH-HCOONa-Hg

2+

system

2.5 激光散射

按实验方法制得反应液,用纳米粒度与 Zeta 电位分析仪记录其粒度。 从实验结果可以看出,无 Hg

2+

存在时,磷钼酸-甲酸钠反应在 AuNPs 的催化作用下反应充分,体系中的磷钼酸少量聚集,颗粒粒度较小

(图 5a),其粒径分布在 100~500 nm;随着 Hg

2+的加入,磷钼酸-甲酸钠反应逐渐被抑制,体系的颗粒分布

在 140~400 nm(图 5b)。 此激光散射粒径与 SEM 粒径不一致的主要原因是电镜制备样品需要脱水,此过

程会导致颗粒聚集。

2.6 分析条件优化

考察磷钼酸浓度对体系 RRS 信号的影响,加入 5.0 mmol / L HCOOH-HCOONa 缓冲溶液,1.45 mg / L

AuNPs,1.0 μmol / L HgSO4 和不同浓度的磷钼酸溶液,其余按照 1. 2 节的实验步骤操作,当加入 0. 425

143

广西师范大学学报(自然科学版),2022,40(2)

图 4 体系扫描电镜

Fig. 4 Scanning electron micrograph of the system

a: 0.425 mmol / L 磷钼酸+ 5.0 mmol / L HCOOH-HCOONa+1.45 mg / L AuNPs; b: a+0.35 μmol / L Hg

2+

图 5 粒度分析

Fig. 5 Particle size analysis

mmol / L 磷钼酸时 ΔI 达到最大,故选用 0.425 mmol / L 磷钼酸(图 6A)。 随着 AuNPs 浓度的增加,其催化作

用逐渐增强,当加入 1. 45 mg / L AuNPs 时 ΔI 达到最大,故选用 1. 45 mg / L AuNPs ( 图 6B)。 考察了

HCOOH-HCOONa 缓冲溶液 pH 对体系 RRS 信号的影响,当加入 pH = 3.1 的 HCOOH-HCOONa 缓冲溶液时

ΔI 达到最大,故选用 pH = 3.1 的 HCOOH-HCOONa 缓冲溶液(图 6C)。 HCOOH-HCOONa 缓冲溶液的浓度

同样会对体系 RRS 信号产生影响,当加入 5.0 mmol / L HCOOH-HCOONa 缓冲溶液时 ΔI 达到最大,故选用

5.0 mmol / L 的 HCOOH-HCOONa 缓冲溶液(图 6D)。 考察反应温度对体系 RRS 信号的影响,当反应温度

为 80 ℃时 ΔI 达到最大,故选用 80 ℃作为反应温度(图 6E)。 考察反应时间对体系 RRS 信号的影响,当

反应时间为 35 min 时 ΔI 达到最大,故选用 35 min 作为反应时间(图 6 F)。

2.7 工作曲线

在最佳实验条件下,AuNPs 表现出对于磷钼酸-甲酸钠 RRS 体系的催化作用,在 0.29 ~ 2. 61 mg / L

AuNPs,450 nm 处的 RRS 强度变化 ΔI 与 AuNPs 浓度呈线性关系, AuNPs 的催化线性方程为 y =

288.39CAuNPs

+26.9,线性相关系数 R

2 为 0.994(图 2A)。 基于磷钼酸-AuNPs-甲酸钠体系表现出对 Hg

2+浓

度的响应,建立一个测定 Hg

2+的共振瑞利散射方法,在 2.5×10

-4

~3.5 μmol / L Hg

2+

,450 nm 处的 RRS 强度

变化 ΔI 与 Hg

2+浓度呈线性关系,线性方程为 y = 0.32CHg

2+ +46.1,线性相关系数 R

2 为 0.994 6(图 2B)。 对

144

http:∥xuebao.gxnu.edu.cn

图 6 磷钼酸-AuNPs-HCOOH-HCOONa-Hg

2+体系反应条件优化

Fig. 6 Optimization of the conditions of phosphomolybdic acid-AuNPs-HCOOH-HCOONa-Hg

2+

system

于磷钼酸-甲酸钠 UV 体系,在 0.29~2.61 mg / L AuNP,700 nm 处的 UV 强度变化 ΔA 与 AuNPs 浓度呈线性

关系,线性方程为 y = 0.241CAuNPs

+0.003,线性相关系数 R

2 为 0.994 8(图 3A)。 已知磷钼酸-AuNPs-甲酸钠

体系可用于检测 Hg

2+

,在 1.0×10

-2

~3.5 μmol / L Hg

2+

,700 nm 处的 UV 强度变化 ΔA 与 Hg

2+浓度呈线性关

系,线性方程为 y = 0.000 1CHg

2+ +0.028,线性相关系数 R

2 为 0.986 9(图 3B)。 表 1 是本方法与部分文献报

道方法的对比。

表 1 本法与已报道方法分析特性比较

Tab. 1 Comparison of analysis characteristics between this method and the reported method

方法 线性范围/ (nmol·L

-1

) 检出限/ (nmol·L

-1

) 注解 参考文献

荧光法 50~ 20 000 13.2 高灵敏但底物合成繁琐 [27]

伏安法 2.4~ 220 0.8 选择性好但电极修饰繁杂 [28]

比色法 10

4

~ 10

5

3.6 稳定性好但检测范围不宽 [29]

比色法 5 000~ 35 000 2 370 快速方便但检测灵敏度低 [30]

本法 0.25~ 3 500 0.18 灵敏、简便、快速 本文

2.8 干扰试验

按实验方法考察了共存离子对磷钼酸-AuNPs-HCOOH-HCOONa-Hg

2+ RRS 测定 0.5 μmol / L Hg

2+的干

扰情况。 由表 2 表明,5 μmol / L Mg

2+

、SO

2-

4 、Fe

3+

、Pb

2+

;25 μmol / L Ca

2+

、Br

-

、Cr

3+

、Al

3+

、HSA、血红蛋白、血

浆蛋白;50 μmol / L IO

-

3 、Na

+

、Zn

2+

、Ba

2+

、NH

+

4 、NO

-

3 、Co

2+

、Cu

2+

、Mn

2+

、BSA、鲱鱼精 DNA、牛白蛋白、核糖核酸

酶不干扰测定,表明本法具有较好的选择性。

2.9 样品分析

从污水厂取得 3 份工业废水,过滤处理后分别取样品 200 μL,按 1.2 节试验方法进行测定,检测结果

见表 3。 水样中 Hg

2+含量为 1.9 ~ 5.1 μmol / L,相对标准偏差(RSD) 为 2.6% ~ 7.0%,回收率为 93.7% ~

98.4%,具有较好的回收率和重现性。

145

广西师范大学学报(自然科学版),2022,40(2)

表 2 干扰离子对体系的影响

Tab. 2 Influence of interfering ions on the system

干扰离子 相对倍数 相对误差/ % 干扰离子 相对倍数 相对误差/ %

Mg

2+

10 5.4 SO

2-

4 10 -7.0

Fe

3+

10 1.6 Pb

2+

10 3.5

Ca

2+

50 9.5 Br

-

50 4.5

Cr

3+

50 -8.8 Al

3+

50 2.8

HSA 50 9.4 血红蛋白 50 -5.2

血浆蛋白 50 -3.4 IO

-

3 100 8.8

Na

+

100 5.6 Zn

2+

100 4.6

Ba

2+

100 1.9 NH

+

4 100 9.1

NO

-

3 100 2.9 Co

2+

100 4.4

Cu

2+

100 -8.8 Mn

2+

100 -8.7

BSA 100 -3.1 鲱鱼精 DNA 100 3.3

牛白蛋白 100 -0.2 核糖核酸酶 100 3.1

表 3 样品的 RRS 测定结果

Tab. 3 RRS measurement results of the samples

样品

平均值

/ (nmol·L

-1

)(n = 5)

加入 Hg

2+

/ (nmol·L

-1

)

测得值/

(nmol·L

-1

)

回收率/ % RSD/ %

Hg

2+含量/

(μmol·L

-1

)

水样 1 192.6 200 386.2 96.8 7.0 1.9

水样 2 261.7 200 449.1 93.7 4.9 2.6

水样 3 506.9 200 703.7 98.4 2.6 5.1

3 结论

本文以磷钼酸颗粒具有 RRS 效应作为分析反应的指示剂,借助 AuNPs 的催化放大作用以及 AuNP 与

Hg

2+的电置换反应,建立了一个简单、灵敏检测痕量 Hg

2+ 的共振瑞利散射分析方法,线性范围为 0.25 ~

3 500 nmol / L,检测限为 0.18 nmol / L。

参 考 文 献

[1] ZHANG C L, LUO L, LUO J, et al. A process-analysis microsystem based on density gradient centrifugation and its application in

the study of the galvanic replacement mechanism of Ag nanoplates with HAuCl

4[ J]. Chemical Communications, 2012, 48

(58): 7241-7243.

[2] JAIN P K, HUANG X H, EL-SAYED I H, et al. Noble metals on the nanoscale: optical and photothermal properties and

some applications in imaging, sensing, biology, and medicine[J]. Accounts of Chemical Research, 2008, 41: 1578-1586.

[3] LU X M, CHEN J Y, SKRABALAK S E, et al. Galvanic replacement reaction: a simple and powerful route to hollow and

porous metal nanostructures[J]. Proceedings of the Institution of Mechanical Engineers, Part N: Journal of Nanomaterials,

Nanoengineering and Nanosystems, 2007, 221(1): 1-16.

146