
Atopic dermatitis (AD, eczema) is an inflammatory skin condition whose histopathology involves remodeling. Few preclinical AD studies are performed using male mice. The histopathological mechanisms underlying AD development were investigated here in male mice at a pre-lesional stage using a human AD-like mouse model. Hypodermal cellular infiltration without thickening of skin layers was observed after one epicutaneous exposure to antigen ovalbumin (OVA), compared to controls. In contrast to our previous report using female mice, OVA treatment did not activate skin mast cells (MC) or elevate sphingosine-1-phosphate (S1P) levels while increasing systemic but not local levels of CCL2, CCL3 and CCL5 chemokines. In contrast to the pathogenic AD mechanisms we recently uncovered in female, S1P-mediated skin MC activation with subsequent local chemokine production is not observed in male mice, supporting sex differences in pre-lesional stages of AD. We are proposing that differential involvement of the MC/S1P axis in early pathogenic skin changes contributes to the well documented yet still incompletely understood sex-dimorphic susceptibility to AD in humans.
Citation: Ross M. Tanis, Piper A. Wedman-Robida, Alena P. Chumanevich, John W. Fuseler, Carole A. Oskeritzian. The mast cell/S1P axis is not linked to pre-lesional male skin remodeling in a mouse model of eczema[J]. AIMS Allergy and Immunology, 2021, 5(3): 160-174. doi: 10.3934/Allergy.2021012
[1] | Honghua Bin, Daifeng Duan, Junjie Wei . Bifurcation analysis of a reaction-diffusion-advection predator-prey system with delay. Mathematical Biosciences and Engineering, 2023, 20(7): 12194-12210. doi: 10.3934/mbe.2023543 |
[2] | Yu Yang, Gang Huang, Yueping Dong . Stability and Hopf bifurcation of an HIV infection model with two time delays. Mathematical Biosciences and Engineering, 2023, 20(2): 1938-1959. doi: 10.3934/mbe.2023089 |
[3] | Dan Liu, Shigui Ruan, Deming Zhu . Stable periodic oscillations in a two-stage cancer model of tumor and immune system interactions. Mathematical Biosciences and Engineering, 2012, 9(2): 347-368. doi: 10.3934/mbe.2012.9.347 |
[4] | Hongying Shu, Wanxiao Xu, Zenghui Hao . Global dynamics of an immunosuppressive infection model with stage structure. Mathematical Biosciences and Engineering, 2020, 17(3): 2082-2102. doi: 10.3934/mbe.2020111 |
[5] | Fang Liu, Yanfei Du . Spatiotemporal dynamics of a diffusive predator-prey model with delay and Allee effect in predator. Mathematical Biosciences and Engineering, 2023, 20(11): 19372-19400. doi: 10.3934/mbe.2023857 |
[6] | Jinhu Xu, Yicang Zhou . Bifurcation analysis of HIV-1 infection model with cell-to-cell transmission and immune response delay. Mathematical Biosciences and Engineering, 2016, 13(2): 343-367. doi: 10.3934/mbe.2015006 |
[7] | Anuj Kumar, Yasuhiro Takeuchi, Prashant K Srivastava . Stability switches, periodic oscillations and global stability in an infectious disease model with multiple time delays. Mathematical Biosciences and Engineering, 2023, 20(6): 11000-11032. doi: 10.3934/mbe.2023487 |
[8] | Huan Kong, Guohong Zhang, Kaifa Wang . Stability and Hopf bifurcation in a virus model with self-proliferation and delayed activation of immune cells. Mathematical Biosciences and Engineering, 2020, 17(5): 4384-4405. doi: 10.3934/mbe.2020242 |
[9] | Meiling Zhu, Huijun Xu . Dynamics of a delayed reaction-diffusion predator-prey model with the effect of the toxins. Mathematical Biosciences and Engineering, 2023, 20(4): 6894-6911. doi: 10.3934/mbe.2023297 |
[10] | Huan Dai, Yuying Liu, Junjie Wei . Stability analysis and Hopf bifurcation in a diffusive epidemic model with two delays. Mathematical Biosciences and Engineering, 2020, 17(4): 4127-4146. doi: 10.3934/mbe.2020229 |
Atopic dermatitis (AD, eczema) is an inflammatory skin condition whose histopathology involves remodeling. Few preclinical AD studies are performed using male mice. The histopathological mechanisms underlying AD development were investigated here in male mice at a pre-lesional stage using a human AD-like mouse model. Hypodermal cellular infiltration without thickening of skin layers was observed after one epicutaneous exposure to antigen ovalbumin (OVA), compared to controls. In contrast to our previous report using female mice, OVA treatment did not activate skin mast cells (MC) or elevate sphingosine-1-phosphate (S1P) levels while increasing systemic but not local levels of CCL2, CCL3 and CCL5 chemokines. In contrast to the pathogenic AD mechanisms we recently uncovered in female, S1P-mediated skin MC activation with subsequent local chemokine production is not observed in male mice, supporting sex differences in pre-lesional stages of AD. We are proposing that differential involvement of the MC/S1P axis in early pathogenic skin changes contributes to the well documented yet still incompletely understood sex-dimorphic susceptibility to AD in humans.
atopic dermatitis;
ovalbumin;
mast cells;
sphingosine-1-phosphate;
blood vessels
Cancer, as the second-most common fatal disease, remains a serious threat to human health, and millions of people suffer from cancer around the world. Cancer cells proliferate and spread rapidly, which causes many difficulties in treatment. The function of the immune system is to prevent tumor cells from invading the body, and the body restricts the growth of tumor cells. In view of molecular cell biology, the immune system has a natural defense for the body, since it can recognize and destroy tumor cells [1,2]. It is obvious that immune system treatment is effective for the treatment of cancer and has more advantages than traditional treatment methods, such as chemotherapy, operation, and so on. Just owing to the character of the immune system, immunotherapy arises gradually, and immunotherapy has become the primary treatment instead of other methods in recent decades. Due to the complexity associated with the function of the immune system in the process of tumor proliferation, pathologists and clinicians think that treating malignant tumors is still very hard, such that human beings cannot overcome it even today.
It is necessary to construct a mathematical model for understanding the mechanism of the immune system, numerous linear and nonlinear mathematical models endowed with functional response have been considered by mathematicians, which play a crucial role in predicting or controlling tumor growth. It is also necessary to study the mechanisms underlying the interaction between the tumor-immune system to understand immunotherapy. Further development of immunotherapy depends on a better comprehension of the interaction between the tumor system and immune system. Tumor cells can escape the monitoring of the immune system to proliferate infinitely, which is called immune suppression. Many researchers are devoted to removing immune suppression so that cancer cells can be controlled by the immune system. Mathematical models, for example, ordinary differential equation (ODE), partial differential equation (PDE) and partial functional differential equation (PFDE), are extensively used to investigate the dynamics of the immune system. Many researches study how the immune system influences tumor cells by modeling and analyzing mathematical models. Experts attempt to solve the problem based on theoretical analysis and experimental data. There is no doubt that constructing mathematical models provides us with a new view to recognize the tumor-immune system.
In recent years, immunotherapy models have been extensively utilized to predict the behaviors of tumor cells. Ref. [3] constructed a tumor-immune system with two distributed delays to describe the interactions between tumor cells and immune cells. Ref. [4] proposed a tumor-immune model with time delay, and explained that the time delay of the immune reaction could change the stability of equilibrium. Ref. [5] analyzed a two-dimensional tumor-immune system with two delays. By using bifurcation theory, the stability of equilibrium and the existence of Hopf bifurcation were given when the time delays were regarded as the bifurcation parameters. Relevant scholars proposed all kinds of differential equations [6-8]. Refs. [9,10] showed that different sizes of delays could cause a change in equilibrium stability. Refs. [11-13] studied the various dynamic properties associated with the changing time delay between tumor cells and immune cells. Ref. [14] studied a tumor-immune system in view of optimal control, and showed the method associated with minimizing the cost of immune-chemotherapy and reducing the load of tumor cells. Ref. [15] investigated a tumour-immune system with delays and fractional-order, and provided the necessary and sufficient conditions for stability of the steady states and Hopf bifurcating periodic solutions. Ref. [16] presented a delay differential equation to analyze the dynamics between effector cells and tumor cells. The existence and stability of possible steady states and the local stability of Hopf bifurcating periodic solutions are discussed. Das et al. studied a series of studies associated with tumor-immune interaction. For example, Ref. [17] studied the deterministic and stochastic dynamics of tumor-immune interactions. Refs. [18,19] considered the dynamics of time-delayed tumor-immune systems. Refs. [20,21] compared the deterministic and noise-induced chaotic dynamics, and described the growth and proliferation process of tumor cells. The analysis showed that the onset of chaos in the system can be predicted. Refs.[22,23] established the optimal treatment strategies that maximize the number of immune-effector cells, minimized the number of cancer cells, and detrimental effects caused by the amount of drugs.
The motivation of this paper is as follows. First, Ref. [24] showed that there existed a time delay during the immune system recognition of tumor cells. Ref. [25] pointed out that delay played an important role in the interaction between the immune system and tumor cells. The different delays could cause different phenomena of the tumor-immune system, which is useful for understanding the tumor-immune system and providing treatment for tumors. Second, Ref. [26] indicated that tumor cells could spread in the body and revealed some factors causing the diffusion of tumor cells. Therefore, it is necessary and reasonable to introduce delay and diffusion into the tumor-immune system. Third, Ref. [27] considered some results on the nonlinear dynamics of delayed differential equation models describing the interaction between tumor cells and effector cells of the immune system, however, it did not consider the influence of diffusion of tumor cells and immune cells. For the immune reaction of organisms, it must be enough time to identify nonself cells, and the growth of tumor cells also needs some time. Based on the model of Ref. [27], we introduce two delays and diffusion, and propose a tumor-immune system with two delays and diffusion.
The paper is organized as follows. In Section 2, we establish a partial functional differential equation (PFDE) associated with tumor-immune system, and study the dynamics of this system. In Section 3, the stability of nonnegative equilibrium and the existence of Hopf bifurcation are presented. The normal form of Hopf bifurcation is deduced in Section 4. Numerical simulations to test our findings are shown in Section 5. In Section 6, we give a conclusion for our work.
The essential mechanisms of interaction between tumor cells and immune cells is shown as follows (see Figure 1).
Immune cells assault tumor cells, and the multiplication of immune cells is stimulated by the appearance of tumor cells. The body can not recognize tumor cells immediately, and there exists a time delay τ1 for responding to the entrance of nonself cells, and the response time may be short, however, it cannot be ignored. Meanwhile, there also exists a time delay τ2 when tumor cells are recognized by immune cells. Furthermore, tumor cells and immune cells diffuse in the body under conditions of limited resources and space, and the growth and interaction between tumor cells and immune cells not only depend on time, but are also affected and restricted by the living spatial environment. Hence, it is also indispensable to investigate the complex interaction caused by spatial factors. Considering the above factors, we construct the following model to analyze the dynamics of the tumor-immune system:
{∂u(x,t)∂t=d1Δu(x,t)+ru(x,t)(1−u(x,t−τ1)α)−mu(x,t)v(x,t),t>0,x∈Ω,∂v(x,t)∂t=d2Δv(x,t)+βu(x,t−τ2)v(x,t−τ2)−dv(x,t)−pu(x,t)v(x,t),t>0,x∈Ω,∂u(x,t)∂x=∂v(x,t)∂x=0,t≥0,x∈∂Ω,u(x,t)=u0(x,t)≥0,v(x,t)=v0(x,t)≥0,(t,x)∈[−τ,0]×ˉΩ. | (2.1) |
where Ω=[0,π] is a bounded open domain in RN(N≥1) with a smooth boundary ∂Ω, and Δ denotes the Laplacian operator in RN(N≥1). u and v are the numbers of tumor cells and immune cells, respectively. The description of the parameters is given in Table 1, and these parameters are all positive. We also assume that u0,v0∈C=C([−τ,0],X) and X is defined by X={u,v∈W2,2(Ω):∂u(x,t)∂x=∂v(x,t)∂x=0,x∈∂Ω}. There exists a time delay for producing a proper immune response in host cells, and considering the time delay of tumor growth τ1, the growth of tumor cells follows ru(1−u(x,t−τ1)α). βu(x,t−τ2)v(x,t−τ2) means that immune cells are activated due to the stimulus of tumor cells.
Parameter | description |
d1 | Diffusion rate of tumor cells |
d2 | Diffusion rate of immune cells |
r | Growth rate of malignant cells |
α | Environmental carrying capacity |
m | Death rate of tumor cells |
p | Death rate of immune cells |
d | Intrinsic death rate of the immune cells |
β | Activating rate of the immune cells |
τ1 | The time delay of tumor cells proliferation |
τ2 | The time of immune cells recognizing tumor cells |
System (2.1) has two boundary equilibria E0=(0,0),E1=(α,0), and it has a unique positive constant steady equilibrium E∗=(u∗,v∗)=(dβ−p,rm(1−dα(β−p))) when d<α(β−p) and β>p. For convenience, we show this condition associated with the existence of positive equilibrium E∗ by using the following assumption.
(H1):d<α(β−p),β>p. |
Denote U(x,t) = (u(x,t),v(x,t))T, then linearized system for Eq. (2.1) can be rewritten as a differential equation at equilibrium E=(u0,v0) (here (u0,v0)=(0,0),(α,0),(u∗,v∗)):
∂U(x,t)∂t=DΔU(x,t)+AU(x,t)+BU(x,t−τ1)+CU(x,t−τ2) |
where
D=(d100d2),A=(r−ru0α−mv0−mu0−pv0−d−pu0),B=(−ru0α000),C=(00βv0βu0). |
Hence, the characteristic equation of (2.1) at (u0,v0) is given as follows:
λ2+Anλ+(Bnλ+Cn)e−λτ1+(Dnλ+En)e−λτ2+Fne−λ(τ1+τ2)+Gn=0,n=0,1,2,⋯, | (3.1) |
where
An=d+d1n2+d2n2+pu0−r+ru0α+mv0,Bn=ru0α,Cn=ru0α(d2n2+d+pu0),Dn=−βu0,En=mu20β−βu0(d1n2−r+ru0α+mv0),Fn=−ru20βα,Gn=(d1n2−r+ru0α+mv0)(d2n2+d+pu0)−mpu0v0. |
When τ1=0, τ2=0, Eq. (3.1) for equilibrium E0=(0,0) becomes
λ2+(d1n2+d2n2+d−r)λ+(d1n2−r)(d2n2+d)=0,n=0,1,2,⋯. | (3.2) |
For Eq. (3.1) with n=0, the product of two eigenvalues −rd<0, obviously, equilibrium E0=(0,0) is always unstable.
When τ1=0, τ2=0, Eq. (3.1) for equilibrium E1=(α,0) becomes
λ2+(d1n2+d2n2+d+pα+r−βα)λ+(d1n2+r)(d2n2+d+pα−βα)+mα2β=0. | (3.3) |
When (H1) does not hold,
d1n2+d2n2+d+pα+r−βα≥d+pα+r−βα>d+pα−βα>0,(d1n2+r)(d2n2+d+pα−βα)+mα2β>0, |
thus, equilibrium E1=(α,0) is always locally asymptotically stable when (H1) does not hold.
When (H1) holds, there also exists positive equilibrium E2. The stability analysis for E1=(α,0) is similar to E∗=(u∗,v∗), and we only show the general stability results for E=(u0,v0) (here (u0,v0)=(α,0),(u∗,v∗)). Thus, we show following assumption.
(H2):d+pu0+2ru0α−r+mv0−βu0>0. |
When τ1=0, τ2=0, Eq. (3.1) becomes
λ2+(An+Bn+Dn)λ+Cn+En+Fn+Gn=0,n=0,1,2,⋯, |
note that under (H2), An+Bn+Dn≥A0+B0+D0=d+pu0+2ru0α−r+mv0−βu0>0. Denote
F(n)=Cn+En+Fn+Gn=ξ2n4+ξ1n2+ξ0, | (3.4) |
where
ξ2=d1d2,ξ1=d2(2ru0α−r+mv0)−βu0d1+d1(d+pu0),ξ0=ru0α(d+pu0)+mu20β−ru20βα−mpu0v0−(ru0α−r+mv0)βu0. |
Thus, we obtain following theorem.
Theorem 3.1. For system (2.1) with τ1=0, τ2=0, the stability results for equilibria are given as follows.
(1) Equilibrium E0=(0,0) is always unstable.
(2) Equilibrium E1=(α,0) is always locally asymptotically stable when (H1) does not hold.
(3) When (H1) and (H2) hold, the stability results for equilibrium E1 (or E∗) are given as follows:
(3-1) If ξ0<0, equilibrium E1 (or E∗) is unstable.
(3-2) If ξ0>0 and −ξ12ξ2<0, equilibrium E1 (or E∗) is locally asymptotically stable.
(3-3) If ξ0>0, −ξ12ξ2>0, F(⌊−ξ12ξ2⌋)>0 and F(⌊−ξ12ξ2⌋+1)>0, equilibrium E1 (or E∗) is locally asymptotically stable.
(3-4) If ξ0>0, −ξ12ξ2>0, and F(⌊−ξ12ξ2⌋)<0 or F(⌊−ξ12ξ2⌋+1)<0, equilibrium E1 (or E∗) is unstable.
The stability analysis for E1=(α,0) is similar to E∗=(u∗,v∗), thus, we only show the general stability results for E=(u0,v0) (here (u0,v0)=(α,0),(u∗,v∗)).
When τ2=0, τ1≠0, the characteristic equation of system (2.1) is given as follows:
λ2+L1λ+(L3λ+L4)e−λτ1+L2=0,n=0,1,2,⋯, | (3.5) |
where L1=An+Dn, L2=En+Gn, L3=Bn,L4=Cn+Fn, with
An=d+d1n2+d2n2+pu0−r+ru0α+mv0,Bn=ru0α,Cn=ru0α(d2n2+d+pu0),Dn=−βu0,En=mu20β−βu0(d1n2−r+ru0α+mv0),Fn=−ru20βα,Gn=(d1n2−r+ru0α+mv0)(d2n2+d+pu0)−mpu0v0. |
We might suppose that ±iω(ω>0) are a pair of purely imaginary roots of Eq. (3.5). Substituting them into Eq. (3.5) and separating the real and imaginary parts, we obtain
{ω2−L2=L3ωsin(ωτ1)+L4cos(ωτ1),−L1ω=L3ωcos(ωτ1)−L4sin(ωτ1), |
thus, R:=sin(ωτ1)=(ω2−L2)L3ω+L1L4ωL23ω2+L24, Q:=cos(ωτ1)=(ω2−L2)L4−L1L3ω2L23ω2+L24, which implies that
ω4−(L23+2L2−L21)ω2+L22−L24=0. |
Let z=ω2, then the above equation can be rewritten in the following form
h(z)=z2−(L23+2L2−L21)z+L22−L24. | (3.6) |
Under assumption
(H3):L23−L21+2L2<0,L22−L24>0, |
Eq. (3.6) has no positive roots.
Under assumption
(H4):L22−L24<0, |
Eq. (3.6) has one unique positive root ω0.
Under assumption
(H5):L23−L21+2L2>0,L22−L24>0,(L23−L21+2L2)2−4(L22−L24)2>0, |
Eq. (3.6) has two positive roots ω1,2, and the critical time delay can be determined as,
τ(j)1l={1ωl[arccos(Q)+2jπ],R≥0,1ωl[2π−arccos(Q)+2jπ],R<0, | (3.7) |
where R:=sin(ωτ1)=(ω2−L2)L3ω+L1L4ωL23ω2+L24, Q:=cos(ωτ1)=(ω2−L2)L4−L1L3ω2L23ω2+L24, l=0,1,2;j=0,1,2,⋯. Then the following transversality conditions yield,
Re[d(λ)dτ1]|−1τ1=τ(j)1l=Re((2λ+L1)eλτ1λ(L3λ+L4)+L3λ(L3λ+L4))=zΔh′(z), |
where Δ=L22z+L23z2,l=0,1,2;j=0,1,2,⋯. Obviously, under (H1)and(H4), h(z)=0 has one positive real root z0 and h′(z0)>0; under (H1)and(H5), h(z)=0 has two positive real roots z1 and z2. Suppose z1<z2, then h′(z1)<0 and h′(z2)>0.
Theorem 3.2. For system (2.1) with τ2=0,τ1>0, the following conclusions hold when assumptions (H1),(H2) and the stability conditions of Theorem 3.1 (3-2) or (3-3) hold,
(1) If (H3) is satisfied, then equilibrium E1 (or E∗) is locally asymptotically stable for all τ1>0.
(2) If (H4) is satisfied, then equilibrium equilibrium E1 (or E∗) is locally asymptotically stable for all 0<τ1<τ(0)10, and unstable for τ1>τ(0)10, furthermore, system (2.1) undergoes a Hopf bifurcation at equilibrium E1 (or E∗) when τ1=τ(j)10,j=0,1,2,⋯.
(3) If (H5) is satisfied, then there is a positive integer m such that the equilibrium E1 (or E∗) is locally asymptotically stable when τ1∈[0,τ(0)12)⋃∪(m−1)k=0(τ(k)11,τ(k+1)12), and unstable τ1∈∪(m−1)k=0(τ(k)12,τ(k)11)⋃(τ(m)12,+∞), furthermore, system (2.1) undergoes a Hopf bifurcation at equilibrium E1 (or E∗) when τ1=τ(j)1l,l=1,2;j=0,1,2,⋯.
The first stable interval mentioned in Theorem 3.2, on which the equilibrium E1 (or E∗) is locally asymptotically stable, is denoted by stable region I. That is, assumptions (H1),(H2) and the stability conditions of Theorem 3.1 (3-2) or (3-3) hold, and if (H3) also holds, I={τ1|τ1>0}. If (H4) also holds, I={τ1|0<τ1<τ(0)10}. If (H5) holds, I={τ1|0<τ1<τ(0)12}. Next, we choose τ1=τ∗1∈I. Regarding τ2 as a parameter, the characteristic equation of (2.1) at (u0,v0) (here (u0,v0)=(α,0),(u∗,v∗)) is rewritten as follows:
λ2+(An+Bne−λτ∗1)λ+(Dnλ+En)e−λτ2+Fne−λτ∗1+Cne−λτ∗1+Gn=0, |
where
An=d+d1n2+d2n2+pu0−r+ru0α+mv0,Bn=ru0α,Cn=ru0α(d2n2+d+pu0),Dn=−βu0,En=mu20β−βu0(d1n2−r+ru0α+mv0),Fn=−ru20βα,Gn=(d1n2−r+ru0α+mv0)(d2n2+d+pu0)−mpu0v0. |
Letting λ=iω(τ2)(ω>0) be the root of the above equation, then separate the real and imaginary parts for above equation,
ω2−Gn−Cncos(ωτ∗1)−Bnωsin(ωτ∗1)=Encos(ωτ2)+Dnωsin(ωτ2)+Fncos(ω(τ∗1+τ2)),Anω+Cnsin(ωτ∗1)+Bnωcos(ωτ∗1)=−Dnωcos(ωτ2)+Ensin(ωτ2)+Fnsin(ω(τ∗1+τ2)), |
which leads to,
F(ω)=ω4+(B2n−2Gn+A2n−D2n)ω2+2(BnGnω−AnCnω−Bnω3+DnFnω)sin(ωτ∗1)+G2n+C2n−E2n−F2n+2(CnGn+AnBnω2−Cnω2−EnFn)cos(ωτ∗1)=0. |
Suppose
(H6):(Gn+Cn)2−(En+Fn)2<0, |
then, we get F(0)=G2n+C2n−E2n−F2n+2CnGn−2EnFn<0,F(+∞)>0. Hence, F(ω)=0 has definite positive roots ω2i,i=0,1,2. For every fixed ωi, there is a sequence of τ(j)2i,j=0,1,2,⋯, defined by
τ(j)2i={1ω2iarccos[Q+2jπ],P≥0,1ω2iarccos[Q+2jπ],P<0, | (3.8) |
where j=0,1,2,⋯,
Q=cos(ωτ2)=f(En+Fncos(ωτ∗1))−g(Dnω−Fnsin(ωτ∗1))(En+Fncos(ωτ∗1))2+(Dnω−Fnsin(ωτ∗1))2,P=sin(ωτ2)=g(En+Fncos(ωτ∗1))+f(Dnω−Fnsin(ωτ∗1))(En+Fncos(ωτ∗1))2+(Dnω−Fnsin(ωτ∗1))2, |
with
f=ω2−Bnωsin(ωτ∗1)−Cncos(ωτ∗1)−Gn,g=Anω+Bncos(ωτ∗1)ω+Cnsin(ωτ∗1). |
Let τ20=minτ(j)2i,i=0,1,2,j=0,1,2,⋯, when τ2=τ20, Eq. (3.1) has a pair of purely imaginary roots ±iω20 for τ∗1∈I. Assume
(H7):Re[dλdτ2]τ2=τ(j)20≠0, |
then we have,
Theorem 3.3. Assuming that (H1),(H6)and(H7) are satisfied, and one of the assumptions (H3),(H4) and (H5) holds, then for τ∗1∈I,system (2.1) undergoes a Hopf bifurcation at E=(u0,v0) when τ2=τ(j)20, here (u0,v0)=(α,0) or (u∗,v∗) with u∗=dβ−p,v∗=rm(1−dα(β−p)), and the equilibrium E is locally asymptotically stable when τ2∈[0,τ20).
In this section, we derive the normal formal of Hopf bifurcation at E∗=(u∗,v∗) when τ2=0,τ1≠0 for the system (2.1). We denote the critical value τ=τ2=τc. Define ¯u(x,t)=u(x,τt)−u∗,¯v(x,t)=v(x,τt)−v∗, and drop the bar for convenience, the system (2.1) can be rewritten as,
{∂u(x,t)∂t=τ[d1△u(x,t)−ru(x,t−1)α(u∗+u(x,t))−rαu∗u(x,t)+ru(x,t)−m(u(x,t)v(x,t)+u∗v(x,t)+v∗u(x,t))],t>0,x∈Ω,∂v(x,t)∂t=τ[d2△v(x,t)+(β−p)(u(x,t)v(x,t)+u(x,t)v∗+u∗v(x,t))−dv(x,t)],t>0,x∈Ω,∂u(x,t)∂x=∂v(x,t)∂x=0,t≥0,x∈∂Ω,u(x,t)=u0(x,t)≥0,v(x,t)=v0(x,t)≥0,(t,x)∈[−τ,0]×ˉΩ. | (4.1) |
Let h=(h11,h12)T be the eigenvector of the linear operator corresponding to the eigenvalue iωτ, and let h∗=(h21,h22)T be the normalized eigenvector of the adjoint operator of the linear operator corresponding to the eigenvalues −iωτ satisfying the inner product <¯h∗,h>=¯h∗T⋅h=1. By a simple calculation, we get
h=(h11,h12)T=(1,(β−p)v∗iω+d2n2+d+pu∗−βu∗)T,h∗=(h21,h22)T=l(iω−d2n2−d−pu∗+βu∗mu∗,1)T,l=[iω−d2n2−d−pu∗+βu∗mu∗+(β−p)v∗−iω+d2n2+d+pu∗−βu∗]−1. | (4.2) |
Suppose the solution of Eq. (4.1) is
U(x,t)=U(x,T0,T1,T2,⋯)=+∞∑k=1εkUk(x,T0,T1,T2⋯), | (4.3) |
where
U(x,T0,T1,T2,⋯)=(u(x,T0,T1,T2,⋯),v(x,T0,T1,T2,⋯))T, |
Uk(x,T0,T1,T2,⋯)=(uk(x,T0,T1,T2,⋯),vk(x,T0,T1,T2,⋯))T. |
The derivation with respect to t is
ddt=∂∂T0+ε∂∂T1+ε2∂∂T2+⋯=D0+εD1+ε2D2+⋯, |
where the differential operator Di=∂∂Ti,i=0,1,2,⋯.
Denote
uj=uj(x,T0,T1,T2,⋯),uj,1=uj(x,T0−1,T1,T2,⋯), |
vj=vj(x,T0,T1,T2,⋯),vj,1=vj(x,T0−1,T1,T2,⋯),j=1,2,3,⋯. |
We obtain
dU(x,t)dt=εD0U1+ε2D0U2+ε2D1U1+ε3D0U3+ε3D1U2+ε3D2U1+⋯,ΔU(x,t)=εΔU1(x,t)+ε2ΔU2(x,t)+ε3ΔU3(x,t)+⋯. | (4.4) |
We take perturbations as τ=τc+εμ, to deal with the delayed terms, we expend u(x,t−1),v(x,t−1) at u(x,T0−1,T1,T2,⋯) and v(x,T0−1,T1,T2,⋯), respectively, that is,
u(x,t−1)=εu1,1+ε2u2,1−ε2D1u1,1+ε3u3,1−ε3D1u2,1−ε3D2u1,1+⋯,v(x,t−1)=εv1,1+ε2v2,1−ε2D1v1,1+ε3v3,1−ε3D1v2,1−ε3D2v1,1+⋯, | (4.5) |
where uj,1=uj(x,T0−1,T1,T2,⋯),vj,1=vj(x,T0−1,T1,T2,⋯),j=1,2,3,⋯.
Substituting Eqs. (4.3)-(4.5) into Eq. (4.1), for the ε-order terms, we obtain
{D0u1−τcd1Δu1+ru∗ατcu1,1+ru∗ατcu1+rτcu1+mτc(u∗v1+v∗u1)=0,D0v1−τcd2Δv1−τc(β−p)(v∗u1+u∗v1)+dτcv1=0. | (4.6) |
Since ±iωτ are the eigenvalues of the linear part of Eq. (4.1), the solution of Eq. (4.6) can be expressed in the following form
U1(x,T0,T1,T2,⋯)=G(T1,T2,⋯)eiωτcT0hcos(nx)+c.c., | (4.7) |
where c.c. means the complex conjugate of the preceding terms, and h is given in Eq. (4.2).
For the ε2-order terms, we obtain
{D0u2−τcd1Δu2+ru∗ατcu2,1+ru∗ατcu2+rτcu2+mτc(u∗v2+v∗u2)=−D1u1+μd1Δu1+rα(τcu∗D1u1,1−τcu1u1,1−μu∗u1,1−μu∗u1)−rμu1−m(τcu1v1+μu∗v1+μv∗u1),D0v2−τcd2Δv2−τc(β−p)(v∗u2+u∗v2)+dτcv2=−D1v1+μd2Δv1+(β−p)μ(u∗v1+u1v∗)+τc(β−p)u1v1−dμv1. | (4.8) |
Substituting Eq. (4.7) into the right side of Eq. (4.8), we obtain the coefficient vector of term eiωτcT0, denoted as m1, by solvability conditions, that is, <h∗,m1>=0, we obtain
∂G∂T1=MμG, | (4.9) |
where
M=α¯h21h11(v∗−n2d1−u∗e−iωτc−u∗−r)−mαu∗¯h21h12−α¯h22[(β−p)(u∗h12+v∗h11)−h12(d2n2+d)]α(h11¯h21+h12¯h22)−ru∗τch11¯h21e−iωτc. |
Suppose the solution of Eq. (4.8) is
{u2=+∞∑k=0(η0kG¯G+η1kG2e2iωτcT0+¯η1k¯G2e−2iωτcT0)cos(kx),v2=+∞∑k=0(ζ0kG¯G+ζ1kG2e2iωτcT0+¯ζ1k¯G2e−2iωτcT0)cos(kx). | (4.10) |
Denoted ck=<cos(nx)cos(nx),cos(kx)>=∫π0cos(nx)cos(nx)cos(kx)dx, thus
η1k=−ckF1,k[(re−iωτch211+αmh11h12)(d2k2+2iω−βu∗+pu∗+d)+mα(β−p)h11h12],ζ1k=ckF1,k[(2iαω+d1k2α+ru∗e−2iωτc+ru∗+rα+mv∗)(β−p)h11h12−(β−p)v∗(re−iωτch211+mαh11h12)],η0k=−ckF2,k{[1+u∗(β−p)]αm(h11¯h12+h12¯h11)+(d2k2−βu∗+pu∗+d)rh11¯h11(eiωτc+eiωτc)},ζ0k=−ckF2,k{(β−p)v∗rh11¯h12(eiωτc+eiωτc)+(h11¯h12+h12¯h11)[mα−(β−p)(d1k2α+ru∗+rα+v∗mα)]}, | (4.11) |
with
F1,k=[(2iαω+d1k2α+ru∗e−2iωτc+ru∗+rα+mαv∗)(2iω+d2k2−βu∗+pu∗+d)+α(β−p)v∗m]×∫π0cos(kx)cos(kx)dx,F2,k=(∫π0cos(kx)cos(kx)dx)[mu∗(β−p)v∗α+(d2k2−βu∗+pu∗+d)(d1k2α+ru∗+rα+v∗mα)]. |
For the ε3-order terms, we obtain
{D0u3−τcd1Δu3+ru∗ατcu3,1+ru∗ατcu3+rτcu3+mτc(u∗v3+v∗u3)=−D1u2−D2u1+μd1Δu2+rατc(u∗D1u2,1+u∗D2u1,1−u1u2,1+u1D1u1,1−u2u1,1)−rαμ(u∗u2,1−u∗D1u1,1+u1u1,1)−rαμu∗u2−rμu2−mτc(u2v1+u1v2)−mμ(u1v1+u∗v2+v∗u2),D0v3−τcd2Δv3−τc(β−p)(v∗u3+u∗v3)+dτcv3=−D1v2−D2v1+μd2Δv2+(β−p)τc(u2v1+u1v2)+μ(β−p)(u1v1+v∗u2+u∗v2)−dμv2. | (4.12) |
Substituting solutions Eq. (4.7) and Eq. (4.10) into the right side of Eq. (4.12), we obtain the coefficient vector of term eiωτcT0, denoted as m2, by solvability conditions, let <h∗,m2>=0, we obtain
∂G∂T2=χG2¯G, | (4.13) |
where
χ=1F3[−rτc¯h21∑k≥0(ckη0kh11+¯h11η0ke−2iωτc+η0kh11e−iωτc+η1k¯h11eiωτc)+¯h22(β−p)τcα∑k≥0ck(η0kh12+η1k¯h12+ζ0kh11+¯h11ζ1k)], |
with
F3=[¯h21(αh11−ru∗τch11e−iωτc)+α¯h22h12]∫π0cos(nx)cos(nx)dx. |
According to the above analysis, the normal form of Hopf bifurcation for system (2.1) reduced on the center manifold is
∂G∂T=ε∂G∂T1+ε2∂G∂T2+⋯, | (4.14) |
making G→G/ε, thus, Eq. (4.14) becomes:
˙G=MμG+χG2¯G, | (4.15) |
where M and χ are given by Eq. (4.9) and Eq. (4.13), respectively.
Let G=reiθ and substitute it into Eq. (4.15), and we obtain the Hopf bifurcation normal form in polar coordinates:
{⋅r=Re(M)μr+Re(χ)r3,⋅θ=Im(M)μ+Im(χ)r2. | (4.16) |
According to the bifurcation normal form Eq. (4.16) in polar coordinates, we have following theorem.
Theorem 4.1. For system (4.16), if Re(M)μRe(χ)<0 holds, system (2.1) exists periodic solutions near equilibrium E=(u0,v0), here (u0,v0)=(α,0) or (u∗,v∗) with u∗=dβ−p,v∗=rm(1−dα(β−p)).
(1) If Re(M)μ<0, the bifurcating periodic solutions are unstable.
(2) If Re(M)μ>0, the bifurcating periodic solutions are locally asymptotically stable.
According to Ref. [12] we choose the following parameter (see Table 2) for numerical simulation. Obviously, (H1) holds, and system (2.1) has two boundary equilibria E0=(0,0) and E1=(α,0)=(5×106,0), and one nontrivial equilibrium E∗=(u∗,v∗)=(2.5948×105,7.8643×105). Actually, we only care about the stability of the nontrivial equilibrium E∗. In this section, we plot numerical simulations by using MATLAB software.
Parameter | Value |
d1 | 0.02day−1 |
d2 | 0.2day−1 |
r | 0.18day−1 |
α | 5.0×106viablecells−1 |
m | 1.101×10−7day−1viablecells−1 |
p | 3.422×10−10day−1viablecells−1 |
d | 0.00152day−1 |
β | 6.2×10−9day−1viablecells−1 |
When τ1=0,τ2=0, according to Theorem 3.1, E0=(0,0) is always unstable. This means that if system (2.1) without delays, the immune system cannot restrain the growth of tumor cells, and tumor cells proliferate quickly in the body. For equilibrium E∗, assumption (H2) and the conditions in Theorem 3.1 (3-2) hold, thus equilibrium E∗=(u∗,v∗) is locally asymptotically stable. We choose the initial value (2.5947×106,7.8642×105), E∗=(u∗,v∗) is locally asymptotically stable (see Figure 2). Although tumor cells and immune cells can coexist at this moment, the immune system can suppress the growth of tumor cells effectively.
By a simple calculation, (H5) holds, from Eq. (3.7), we obtain τ(0)12=14.8403, we choose τ1=3∈(0,τ(0)12) and choose the initial value (2.5947×106,7.8642×105), according to Theorem 3.2, E∗=(u∗,v∗) is locally asymptotically stable (see Figure 3). From Figure 3, we conclude that if tumor cells need a long time to proliferate, the tumor cells cannot diffuse, and it is easy to cure the tumor at this time. Thus, the diffusion of tumor cells can be controlled by making the size of the delay small. We reduce the diffusion of tumor cells by controlling the time of the proliferation of tumor cells.
We choose τ1=15 and the initial value (2.5947×106,7.8642×105), from (4.15), we obtain Re(M)>0,Re(χ)<0, thus, according to Theorem 4.1, system (2.1) exhibits stable periodic solutions near E∗, and the direction of Hopf bifurcation is forward. This means that, when τ is close to τc, the number of tumor cells can be controlled, and the number of tumor cells and immune cells varies periodically (see Figure 4).
When τ2≠0,τ1≠0, the stable interval mentioned in Theorem 3.2 I={τ1|0<τ1<τ(0)12=14.8403}, on which the equilibrium E∗ is locally asymptotically stable for τ2=0, and (H2) holds, from (3.8), we obtain τ(0)21=4.3239. We choose τ1=1∈I,τ2=1, according to Theorem 3.3, we choose the initial value (2.5947×106,7.8642×105), E∗=(u∗,v∗) is locally asymptotically stable (see Figure 5). The time of tumor growth is fixed within a certain range, and the tumor cells cannot escape the control of the immune system.
Remark: Comparing the above figures, we have the following results. When the time delay is less than the critical value, the equilibrium is locally asymptotically stable, that is, the immune system can restrain the growth of tumor cells effectively at this moment, and tumor cells may not proliferate insanely in the body. Tumor cells originate from the malignant growth of normal cells. It is shown that the immune system can restrain tumor growth, although immune cells and tumor cells could coexist, and the immune system has higher efficiency (see Figures 1, 2 and 5). As the time delay of immune cells identifying tumor cells increased, immune system is able to identify itself from nonself substances, and immunotherapy can work effectively. Then, the numbers of tumor cells and immune cells demonstrate stable periodic variation (see Figure 3). When the time delay of immune cells identifying tumor cells is larger, the numbers of tumor cells are larger, that is, the immune system loses function, and the body may be damaged. On the other hand, the growing speed of tumor cells depends on the mood and diet of the patient, and it is very important to maintain an optimistic mood and a regular diet.
We constructed a tumor-immune system with two delays and diffusion, and investigated how the delays affect the dynamics of the system. We analyzed the existence and stability of equilibria, and studied the dynamic properties of Hopf bifurcation. The nontrivial equilibrium is locally asymptotically stable under suitable parameters for the system with two delays, and the system occurred stable periodical solutions when the delay of tumor cell proliferation passed through the critical value.
The authors declare that they have no competing interests.
This study was funded by Fundamental Research Funds for the Central Universities of China (Grant No. 2572021DJ01) and the Heilongjiang Provincial Natural Science Foundation of China (Grant No. LH2019A001).
The datasets used or analyzed during the current study are available from the corresponding author on reasonable request.
The idea of this research was introduced by Yuting Ding. All authors contributed to the main results and numerical simulations.
In this section, we show the code of MATLAB software for simulating partial differential equations (that is, Figure 2).
functionpdex43m=0;x=0:0.1:1∗pi;t=0:0.25:500;sol=pdepe(m,@pdex44pde,@pdex44ic,@pdex44bc,x,t);u1=sol(:,:,1);u2=sol(:,:,2);figuresurf(x,t,u1);title(′u(x,t)′);xlabel(′x′);ylabel(′t′)figuresurf(x,t,u2);title(′v(x,t)′);xlabel(′x′);ylabel(′t′)function[c,f,s]=pdex44pde(x,t,u,DuDx)d1=0.02;d2=0.2;r=0.18;alpha=5.0∗106;m=1.101∗10−7;beta=6.2∗10−9;d=0.0152;s=3.422∗10−10;c=[1;1];f=[d1;d2].∗DuDx;s=[r∗u(1)∗(1−u(1)/alpha)−m∗u(1)∗u(2);beta∗u(1)∗u(2)−d∗u(2)−s∗u(1)∗u(2)];functionu0=pdex44ic(x);u0=[300000;800000];function[pl,ql,pr,qr]=pdex44bc(xl,ul,xr,ur,t)pl=[0;0];ql=[1;1];pr=[0;0];qr=[1;1]; |
In this section, we show the code of MATLAB software for simulating partial functional differential equations (that is, Figures 3-5).
functiontumor−delayd1=0.02;d2=0.2;r=0.18;alpha=5.0∗106;m=1.101∗10−7;beta=6.2∗10−9;d=0.0152;s=3.422∗10−10;q=3;l=1∗pi;N=10;h=q/N;K=6000;T=h∗K;p=K+N+1;M=100;dx=l/M;u=zeros(M+1,p);v=zeros(M+1,p);fori=1:M+1forj=1:N+1u(i,j)=(2594700);v(i,j)=(786420);endendf=zeros(M+1,p);g=zeros(M+1,p);forj=N+1:p−1fori=2:Mf(i,j)=r∗u(i,j)∗(1−u(i,j−N)/alpha)−m∗u(i,j)∗v(i,j);g(i,j)=beta∗u(i,j)∗v(i,j)−d∗v(i,j)−s∗u(i,j)∗v(i,j);u(i,j+1)=u(i,j)+h∗d1∗(u(i+1,j)−2∗u(i,j)+u(i−1,j))/dx/dx+h∗f(i,j);v(i,j+1)=v(i,j)+h∗d2∗(v(i+1,j)−2∗v(i,j)+v(i−1,j))/dx/dx+h∗g(i,j);endu(1,j+1)=u(2,j+1);u(M+1,j+1)=u(M,j+1);v(1,j+1)=v(2,j+1);v(M+1,j+1)=v(M,j+1);endz1=zeros(p,M+1);z2=zeros(p,M+1);fori=1:M+1forj=1:pz1(j,i)=u(i,j);z2(j,i)=v(i,j);endendx=0:dx:1∗pi;t=−q:h:T;[xx,tt]=meshgrid(x,t);figure;mesh(xx,tt,z1);title(′u(x,t)′);xlabel(′x′);ylabel(′t′);figure;mesh(xx,tt,z2);title(′v(x,t)′);xlabel(′x′);ylabel(′t′); |
1. | Larysa Dzyubak, Oleksandr Dzyubak, Jan Awrejcewicz, Nonlinear multiscale diffusion cancer invasion model with memory of states, 2023, 168, 09600779, 113091, 10.1016/j.chaos.2022.113091 | |
2. | Yanchuang Hou, Chunyue Wei, Yuting Ding, Dynamic analysis of reaction-diffusion dual carbon model considering economic development in China, 2023, 31, 2688-1594, 2483, 10.3934/era.2023126 | |
3. | Zhaoxuan Huo, Jicai Huang, Yang Kuang, Shigui Ruan, Yuyue Zhang, Oscillations in a tumor–immune system interaction model with immune response delay, 2024, 1477-8599, 10.1093/imammb/dqae016 | |
4. | Weili Liu, Hongpeng Zhang, Weipeng Zhang, Xuenan Sun, DYNAMICAL BEHAVIORS OF A TUMOR-IMMUNE-VITAMIN MODEL WITH RANDOM PERTURBATION, 2023, 13, 2156-907X, 2739, 10.11948/20220558 |
Parameter | description |
d1 | Diffusion rate of tumor cells |
d2 | Diffusion rate of immune cells |
r | Growth rate of malignant cells |
α | Environmental carrying capacity |
m | Death rate of tumor cells |
p | Death rate of immune cells |
d | Intrinsic death rate of the immune cells |
β | Activating rate of the immune cells |
τ1 | The time delay of tumor cells proliferation |
τ2 | The time of immune cells recognizing tumor cells |
Parameter | Value |
d1 | 0.02day−1 |
d2 | 0.2day−1 |
r | 0.18day−1 |
α | 5.0×106viablecells−1 |
m | 1.101×10−7day−1viablecells−1 |
p | 3.422×10−10day−1viablecells−1 |
d | 0.00152day−1 |
β | 6.2×10−9day−1viablecells−1 |
Parameter | description |
d1 | Diffusion rate of tumor cells |
d2 | Diffusion rate of immune cells |
r | Growth rate of malignant cells |
α | Environmental carrying capacity |
m | Death rate of tumor cells |
p | Death rate of immune cells |
d | Intrinsic death rate of the immune cells |
β | Activating rate of the immune cells |
τ1 | The time delay of tumor cells proliferation |
τ2 | The time of immune cells recognizing tumor cells |
Parameter | Value |
d1 | 0.02day−1 |
d2 | 0.2day−1 |
r | 0.18day−1 |
α | 5.0×106viablecells−1 |
m | 1.101×10−7day−1viablecells−1 |
p | 3.422×10−10day−1viablecells−1 |
d | 0.00152day−1 |
β | 6.2×10−9day−1viablecells−1 |