Stability and Bifurcation of a Pathogen-Immune Model with Delay and Diffusion Effects

Wang Jingnan,, Yang Dezhong

 基金资助: 国家自然科学基金.  11801122黑龙江省自然科学基金.  A2018008

 Fund supported: the NSFC.  11801122the NSF of Heilongjiang Province.  A2018008

Abstract

In order to understand the effects of diffusion and time-delay factors on the dynamics between pathogens and immune cells, a delayed pathogen-immune reaction diffusion model with homogeneous Neumann boundary condition is established. By using the diffusion ratio of pathogen-immune cells and immune delay as two parameters, the characteristic root distribution of the linearized system at the positive steady state is analyzed and the necessary and sufficient conditions for the positive steady state to undergo Turing instability and Hopf bifurcation are obtained by using the bifurcation theory of functional differential equations. In addition, the dynamic behavior close to the critical value of Turing instability and Hopf bifurcation is intuitively shown by Matlab numerical simulation. The biological and medicinal significance of corresponding dynamic behaviors are discussed. Furthermore, the obtained results provide certain theoretical support for controlling the growth of pathogen.

Keywords： Pathogen immunity ; Reaction diffusion ; Delay ; Turing bifurcation ; Hopf bifurcation

Wang Jingnan, Yang Dezhong. Stability and Bifurcation of a Pathogen-Immune Model with Delay and Diffusion Effects. Acta Mathematica Scientia[J], 2021, 41(4): 1204-1217 doi:

1 引言

$$$\left\{\begin{array}{l} {\frac{{\rm d}P}{{\rm d}t} = rP-\frac{c_{s}PB}{1+a_{s}P}-\frac{c_{u}P M}{1+a_{u}P}}, \\ {\frac{{\rm d}B}{{\rm d}t} = \frac{kPB}{1+k_{m}P}-\delta B+h}, \end{array}\right.$$$

$$$\left\{\begin{array}{l} {u^{\prime} = \alpha u-\frac{u v}{1+\beta_{s} u}-\frac{m u}{1+\beta_{u} u}}, \\ {v^{\prime} = \frac{u v}{1+\gamma u}-v+\eta}, \end{array}\right.$$$

(H1) 当病原体密度特别高时, 特异性免疫细胞的复制率要高于其衰减率, 即$k/k_m>\delta$.

(H2) 在高密度病原体的免疫反应过程中, 特异性免疫强于非特异性免疫, 即$a _{ u }> a _ { s }\geq0$.

$$$\left\{\begin{array}{l} {\frac{\partial u(x, t)}{\partial t} = d_{1} \Delta u(x, t)+\alpha u(x, t)-\frac{u(x, t) v(x, t)}{1+\beta_{s} u(x, t)}-\frac{m u(x, t)}{1+\beta_{u} u(x, t)}, \quad x \in (0, l \pi), t>0}, \\ {\frac{\partial v(x, t)}{\partial t} = d_{2} \Delta v(x, t)+\frac{u(x, t-\tau) v(x, t-\tau)}{1+\gamma u(x, t-\tau)}-v(x, t)+\eta, \quad x \in (0, l \pi), t>0}, \\ {u_{x}(0, t) = v_{x}(0, t) = 0, u_{x}(l \pi, t) = v_{x}(l \pi, t) = 0, \quad t>0}, \\ {u(x, t) = u_{0}(x, t) \geq 0, v(x, t) = v_{0}(x, t), \quad x \in [0, l \pi], t \in[-\tau, 0]}.\end{array}\right.$$$

2 扩散因素与时滞因素对稳态解稳定性的影响

$$$\rho_{B}\left(d_{2}\right) = \left\{\begin{array}{ll} {\rho_{1}\doteq \frac{(\sqrt{B C}-\sqrt{B C-A(1-D)})^{2}}{(1-D)^{2}}}, \ & {0< d_2 \leq \tilde{d}_{2}} , \\ {\rho_{2}(d_2)\doteq \frac{A}{d_{2}+(1-D)}}, & {d_2\geq \tilde{d}_{2}}, \end{array}\right.$$$

本部分证明参考了文献[24]的证明思路. 当$n\in N$时, (2.8)式中的$D(n, \rho) $$n = n_{\min} 处达到极小值 \min_{n\in R^+} D(n, \rho) , 其中 $$n_{\min} = l\sqrt{\frac{-(\rho(1-D)-A)}{2 \rho d_{2}}},$$ $$\min\limits_{n\in R^+} D(n, \rho) = B C+A(D-1)-\frac{\left(\rho (1-D)- A\right)^{2}}{4 \rho}.$$ 经计算可得, 如果 \rho>\rho_{1} , 则 \min_{n\in R^+} D(n, \rho)>0 . 再由 A+D<1 , 可知 T(n, \rho)>0 , 即对于任意的 n\in R^+ , 特征方程(2.6)不可能有零根. 因此, 若特征方程(2.6)存在零根, 则需 \rho<\rho_{1} . 另外, 当 \rho<\rho_{2}(d_2) 时, 可推出 n_{\min}>1/\sqrt{2} , 即存在大于1的正数 n , 使得当 \rho<\rho_{1} 时, D(n, \rho) = 0 , 即特征方程(2.6)有零根. 再由 \rho_2(d_2)$$ d_2$的严格递减函数, 且当$d_2 = 0$时, $\rho_2(0)>\rho_1$, 因此一定存在唯一的$\tilde{d}_{2} = \frac{A(1-D)^{2}}{(\sqrt{B C}-\sqrt{B C-A(1-D)})^{2}}+D-1$, 满足$\rho_2(\tilde{d}_2) = \rho_1$, 且当$d_2\in(0, \tilde{d}_2)$时, $\rho_2({d}_2)>\rho_1$; 当$d_2\in(\tilde{d}_2, +\infty)$时, $\rho_2({d}_2)<\rho_1$. 综上, 若方程(1.4)在$(\rho, d_2)$平面上存在Turing分支点, 那么$(\rho, d_2)\in S$. 证毕.

$$$\rho\left(d_{2}, n\right) = \frac{A(d_{2} \frac{n^{2}}{l^{2}}-Z_{1})}{d_{2} \frac{n^{2}}{l^{2}}\left(d_{2} \frac{n^{2}}{l^{2}}+Z_{2}\right)},$$$

$$$D\left(n, \rho\right) = \rho d_{2}^{2}\left(\frac{n^{2}}{l^{2}}\right)^{2}+\left((1-D) \rho- A\right)d_{2} \frac{n^{2}}{l^{2}}+B C+A(D-1) = 0 .$$$

$\rm(2)$对于$d_{2}>0$, 方程$\rho\left(d_{2}, n\right) = \rho\left(d_{2}, n+1\right)$有且仅有一个正根$d_{n, n+1} \in(d_{M}(n+1), $$d_{M}(n)) , 其表达式为 $$d_{n, n+1} = \frac{Z_{1}}{2}\left[\frac{l^{2}}{n^{2}}+\frac{l^{2}}{(n+1)^{2}}+\sqrt{\left(\frac{l^{2}}{n^{2}}+\frac{l^{2}}{(n+1)^{2}}\right)^{2}+\frac{4 Z_{2} l^{4}}{Z_{1} n^{2}(n+1)^{2}}}\ \right].$$ 而且当 d_{2}>d_{n, n+1} 时, 则 \rho\left(d_{2}, n\right)>\rho\left(d_{2}, n+1\right)>\rho\left(d_{2}, n+2\right)>\cdots 成立. \rm(3) 如果定义 $$\rho_{*}\left(d_{2}\right) = \rho\left(d_{2}, n\right), d_{2} \in\left[d_{n, n+1}, d_{n-1, n}\right), \quad d_{0, 1} = +\infty,$$ 那么, 对 0<d_{2}<\infty , 不等式 \rho_{*}\left(d_{2}\right) \leq \rho_{B}\left(d_{2}\right) 成立, 当且仅当 d_{2} = d_{M}(n), n \in N 等号成立, 此时 \rho_{*}\left(d_{2}\right) = \rho_{B}\left(d_{2}\right) , 其中 \rho_{B}\left(d_{2}\right) 如(2.9)式所定义. 由引理2.2和引理2.3可得如下引理. 引理2.4 对于给定的正整数 n_{T} \in N , 如果 \rho = \rho_{*}\left(d_{2}\right), d_{2} \in\left[d_{n_{T}, n_{T}+1}, d_{n_{T}-1, n_{T}}\right) , 那么, 当 \tau = 0 时, 特征方程(2.4)除了一个单零根 \lambda = 0 以外, 其他根都具有严格负实部. 令 \lambda = \lambda(n, \tau, \rho) 是特征方程(2.4) 的根, 且在 n = n_{T} 时, 满足 \lambda = \lambda(n_{T}, 0, \rho_*(d_2)) = 0 , 且当 (\tau, \rho)\in [0, +\infty)\times (\rho_{*}\left(d_{2}\right)-\epsilon, \rho_{*}\left(d_{2}\right)+\epsilon) 时, 其中 \epsilon 为充分小的正数, 则有 \left.\frac{\mathrm{d} \lambda\left(n_{T}, \tau, \rho\right)}{\mathrm{d} \rho}\right|_{\rho=0 \atop \rho_{*}\left(d_{2}\right)} <0. 当 n = n_{T}, \tau = 0 时, 由(2.12), (2.13)与(2.15)式可知, D\left(n_{T}, \rho_{*}\left(d_{2}\right)\right) = 0 , 即此时, 特征方程(2.4) 存在一个零根 \lambda = 0 . 为了进一步判断 \lambda = 0 是否为特征方程(2.4)的单根. 对特征方程(2.4)中的 \Delta_{n}(\lambda, \tau) 关于 \lambda 求导, 并将 C_{n_{T}} = -B_{n_{T}} 代入, 可得 其中, A_{n_{T}}, B_{n_{T}}$$ C_{n_{T}}$分别为(2.5)式中$A_{n}, B_{n} $$C_{n}$$ n = n_{T}$时的值. 可见此时, 当$\tau = 0$时, 如果$n = n_{T}$, 特征方程(2.4) 存在一个单零根$\lambda(n_T, 0, \rho) = 0$.$A+D-1<0, \quad B C+A(D-1)>0$, 那么当$n \neq n_{T}$时, $T(n, \rho)>0$, $D(n, \rho)>0$, 特征方程(2.4) 其他所有根都具有严格负实部.

$\rm(2) $$\rho = \rho_{*}\left(d_{2}\right), d_{2} \in\left[d_{n_{T}, n_{T}+1}, d_{n_{T}-1, n_{T}}\right) 时, 模型(1.4)在正稳态解 E_{+} 处经历波数为 n_{T} 的Turing分支, 其中 \rho_{*}\left(d_{2}\right) 如(2.15)式所示. 上面得到了扩散因素对正稳态解稳定性的影响, 下面讨论免疫时滞对正平衡解稳定性的影响. 在 A+D-1<0 , A(D-1)+B C>0$$ \rho>\rho_{*}\left(d_{2}\right), d_{2}>0$同时成立的条件下, 由前面的推导可知, 若$\tau = 0$时, 模型(1.4)在正稳态解$E_{+}$是渐近稳定的. 现研究随着免疫时滞$\tau$的增大(即$\tau>0$)模型(1.4)正稳态解$E_{+}$稳定性的变化情况, 即判断模型(1.4)在正稳态解$E_{+}$附近历Hopf分支的条件. 假设特征方程(2.4)存在纯虚根(记为$\lambda_{n} = {\rm i} \omega_{n}$), 将$\lambda_{n} = {\rm i} \omega_{n}$代入特征方程(2.4)中, 则有

$$$\left\{\begin{array}{l}{C_{n} \sin \omega_{n} \tau+D \omega_{n} \cos \omega_{n} \tau = A_{n} \omega_{n}}, \\ {C_{n} \cos \omega_{n} \tau-D \omega_{n} \sin \omega_{n} \tau = \omega_{n}^{2}-B_{n}}.\end{array}\right.$$$

$\begin{eqnarray} &&\cos \omega_{n} \tau = \frac{A_{n} D_{n} \omega_{n}^{2}+C_{n} \omega_{n}^{2}-B_{n} C_{n}}{C_{n}^{2}+D^{2} \omega_{n}^{2}} : = M_{n}(\omega), \end{eqnarray}$

$\begin{eqnarray} && \sin \omega_{n} \tau = \frac{C_{n} A_{n} \omega_{n}-D \omega_{n}\left(\omega_{n}^{2}-B_{n}\right)}{C_{n}^{2}+D^{2} \omega_{n}^{2}} : = N_{n}(\omega). \end{eqnarray}$

$$$\omega_{n}^{4}+\left(A_{n}^{2}-D^{2}-2 B_{n}\right) \omega_{n}^{2}+B_{n}^{2}-C_{n}^{2} = 0.$$$

$$$z_{n}^{2}+\left(A_{n}^{2}-D^{2}-2 B_{n}\right) z_{n}+B_{n}^{2}-C_{n}^{2} = 0.$$$

$$$z^{\pm}_{n} = \frac{-\left(A_{n}^{2}-D^{2}-2 B_{n}\right) \pm \sqrt{\left(A_{n}^{2}-D^{2}-2 B_{n}\right)^{2}-4\left(B_{n}^{2}-C_{n}^{2}\right)}}{2}.$$$

$z_{n} = \omega_{n}^{2}$, 经计算, 可得(2.21)式存在正实根的条件如下.

$$$B_{n}-C_{n} = \rho d^{2}_{2}\left(\frac{n^{2}}{l^{2}}\right)^{2}+\left( (D+1) \rho -A \right)d_{2} \left(\frac{n^{2}}{l^{2}}\right)-(B C+(1+D) A).$$$

$n \in\left[0, n_{0}\right), B_{n}^{2}-C_{n}^{2}<0$；当$n \in\left(n_{0}, \infty\right), B_{n}^{2}-C_{n}^{2}>0$.

$$$\tau_{n, j} = \left\{\begin{array}{ll} \frac{\arccos M_{n}\left(\omega^*_{n}\right)+2 j \pi}{\omega^*_{n}}, & 当 N_{n}\left(\omega^*_{n}\right) \geq 0, \\ \frac{2 (j+1) \pi-\arccos M_{n}\left(\omega^*_{n}\right)}{\omega^*_{n}}, & 当 N_{n}\left(\omega^*_{n}\right) < 0. \end{array}\right.$$$

$\tau$充分接近$\tau_{n, j}$时, 特征方程(2.4)的根$\lambda_{n}(\tau) = \alpha_{n}(\tau)+i \omega_{n}(\tau)$满足$\alpha_{n}\left(\tau_{n, j}\right) = 0$, $\omega_{n}\left(\tau_{n, j}\right) = \omega^*_{n}$. 下面定义

$$$\tau_{0}^{*} = \min \tau_{n, 0}, n\in N, 0 \leq n<n^{0}.$$$

$\rm(1) $$0 \leq \tau<\tau_{0}^{*} ( \tau_{0}^{*} 如(2.25)式所示)时, 则模型(1.4)正稳态解 E_{+} 是渐近稳定的. 当 \tau>\tau_{0}^{*} 时, 模型(1.4)正稳态解 E_{+} 是不稳定的. \rm(2)$$ \tau = \tau_{n, j}^{+}$时, 则模型(1.4)在正稳态解$E_{+}$经历波数为$n$的Hopf分支, 出现周期振动.

3 数值模拟

$$$\alpha = 0.75, \eta = 0.05, \gamma = 0.05, m = 1.5, \beta_{s} = 10^{(-10)}, \beta_{u} = 0.1.$$$

$$$\alpha = 1.75, \eta = 0.05, \gamma = 0.05, m = 1.5, \beta_{s} = 10^{(-10)}, \beta_{u} = 0.1.$$$

图 1

$$$\rho_{2}\left(d_{2}\right) = 0.1198 /\left(d_{2}+0.1296\right),$$$

$$$\rho\left(d_{2}, n\right) = \frac{0.1198\left(d_{2} \frac{n^{2}}{9}-2.4418\right)}{\left(d_{2} \frac{n^{2}}{9}\left(d_{2} \frac{n^{2}}{9}+0.1296\right)\right)}.$$$

参考文献 原文顺序 文献年度倒序 文中引用次数倒序 被引期刊影响因子

Naheed A , Singh M , Lucy D .

Numerical study of SARS epidemic model with the inclusion of diffusion in the system

Appl Math Comput, 2014, 229, 480- 498

Beauchemin C A A , Handel A .

A review of mathematical models of influenza a infections within a host or cell culture: lessons learned and challenges ahead

BMC Public Health, 2011, 11 (Suppl 1): S7

Han X L , Wang W G , Mo J Q .

Bionomics dynamic system for epidemic virus transmission

Acta Math Sci, 2019, 39A (1): 200- 208

Tang S Y , Tang B , Bragazzi N L , et al.

Analysis of COVID-19 epidemic traced data and stochastic discrete transmission dynamic model

Sci Sin Math, 2020, 50 (8): 1071- 1086

Cheng X X , Rao Y Q , Huang G .

COVID-19 transmission model in an enclosed space: a case study of Japan Diamond Princess Cruises

Acta Math Sci, 2020, 40A (2): 540- 544

Gilchrist M A , Sasaki A .

Modeling host-parasite coevolution: a nested approach based on mechanistic models

J Theor Biol, 2002, 218 (3): 289- 308

Mohtashemi M , Levins R .

Transient dynamics and early diagnostics in infectious disease

J Math Biol, 2001, 43 (5): 446- 470

Pugliese A , Gandolfi A .

A simple model of pathogen-immune dynamics including specific and non-specific immunity

Math Biosci, 2008, 214 (1): 73- 80

Wang W , Ma W .

Hepatitis C virus infection is blocked by HMGB1: a new nonlocal and time-delayed reaction-diffusion model

Appl Math Comput, 2018, 320, 633- 653

Stancevic O , Angstmann C N , Murray J M , et al.

Turing patterns from dynamics of early HIV infection

B Math Biol, 2013, 75 (5): 774- 795

Lee M R , Huang Y T , Lee P I , et al.

Healthcare-associated bacteraemia caused by Leuconostoc species at a university hospital in Taiwan between 1995 and 2008

J Hosp Infect, 2011, 78 (1): 45- 49

Cui Q M , Yuan C Y , Li C L , et al.

Control of white - spot disease and scuticociliatida disease of some main cultured sea fishes

Reservoir Fisheries, 2007, 27 (6): 85- 87

Sun W S . Medical Immunology. Bei Jing: Higher Education Press, 2010

Han X L , Jin Z .

A dynamic model of hepatitis B virus with delayed immune response

J North University of China, 2011, 32 (1): 197- 208

Bai Z , Peng R , Zhao X Q .

A reaction-diffusion malaria model with seasonality and incubation period

J Math Biol, 2018, 77 (1): 201- 228

Zhu D D , Ren J L , Zhu H P .

Spatial-temporal basic reproduction number and dynamics for a dengue disease diffusion model for a dengue disease diffusion model

Math Meth Appl Sci, 2018, 41, 5388- 5403

Yamazaki K .

Threshold dynamics of reaction-diffusion partial differential equations model of Ebola virus disease

Int J Biomath, 2018, 11 (8): 1850108

Diggles B K , Lester R J G .

Influence of temperature and host species on the development of cryptocaryon irritans

J Parasitol, 1996, 82 (1): 45- 51

Wang K , Wang W , Pang H , et al.

Complex dynamic behavior in a viral model with delayed immune response

Physica D, 2007, 226 (2): 197- 208

Xie Q , Huang D , Zhang S , et al.

Analysis of a viral infection model with delayed immune response

Appl Math Model, 2010, 34 (9): 2388- 2395

Niu B , Guo Y X , Du Y F .

Hopf bifurcation induced by delay effect in a diffusive tumor-immune system

Int J Bifurcat Chaos, 2018, 28 (11): 1850136

Canabarro A A , Gléria I M , Lyra M L .

Periodic solutions and chaos in a nonlinear model for the delayed cellular immune response

Physica A, 2004, 342 (1/2): 234- 241

Wu J H . Theory and Applications of Partial Functional Differential Equations. New York: Springer-Verlag, 1996

Jiang W , Wang H , Cao X .

Turing instability and Turing-Hopf bifurcation in diffusive Schnakenberg systems with gene expression time delay

J Dynam Differential Equations, 2019, 31 (4): 2223- 2247

Wang W , Liu Q X , Jin Z .

Spatiotemporal complexity of a ratio-dependent predator-prey system

Phys Rev E, 2007, 75 (5): 051913

Baurmann M , Gross T , Feudel U .

Instabilities in spatially extended predator-prey systems: spatio-temporal patterns in the neighborhood of Turing-Hopf bifurcations

J Theor Biol, 2007, 245 (2): 220- 229

Song Y , Jiang H , Liu Q X , et al.

Spatiotemporal dynamics of the diffusive Mussel-Algae model near Turing-Hopf bifurcation

SIAM J on Appl Dyn Syst, 2017, 16 (4): 2030- 2062

Garvie M R .

Finite-difference schemes for reaction-diffusion equations modeling predator-prey interactions in MATLAB

B Math Biol, 2007, 69 (3): 931- 956

/

 〈 〉