Research article

On the zero forcing number and propagation time of oriented graphs

  • Zero forcing is a process of coloring in a graph in time steps known as propagation time. These graph-theoretic parameters have diverse applications in computer science, electrical engineering and mathematics itself. The problem of evaluating these parameters for a network is known to be NP-hard. Therefore, it is interesting to study these parameters for special families of networks. Perila et al. (2017) studied properties of these parameters for some basic oriented graph families such as cycles, stars and caterpillar networks. In this paper, we extend their study to more non-trivial structures such as oriented wheel graphs, fan graphs, friendship graphs, helm graphs and generalized comb graphs. We also investigate the change in propagation time when the orientation of one edge is flipped.

    Citation: Sakander Hayat, Hafiz Muhammad Afzal Siddiqui, Muhammad Imran, Hafiz Muhammad Ikhlaq, Jinde Cao. On the zero forcing number and propagation time of oriented graphs[J]. AIMS Mathematics, 2021, 6(2): 1833-1850. doi: 10.3934/math.2021111

    Related Papers:

    [1] Gaixiang Cai, Fengru Xiao, Guidong Yu . The identification numbers of lollipop graphs. AIMS Mathematics, 2025, 10(4): 7813-7827. doi: 10.3934/math.2025358
    [2] Ali N. A. Koam, Ali Ahmad, Azeem Haider, Moin A. Ansari . Computation of eccentric topological indices of zero-divisor graphs based on their edges. AIMS Mathematics, 2022, 7(7): 11509-11518. doi: 10.3934/math.2022641
    [3] Bilal A. Rather, M. Aijaz, Fawad Ali, Nabil Mlaiki, Asad Ullah . On distance signless Laplacian eigenvalues of zero divisor graph of commutative rings. AIMS Mathematics, 2022, 7(7): 12635-12649. doi: 10.3934/math.2022699
    [4] Mohammad Hamidi, Irina Cristea . Hyperideal-based zero-divisor graph of the general hyperring Zn. AIMS Mathematics, 2024, 9(6): 15891-15910. doi: 10.3934/math.2024768
    [5] Meijin Luo, Qiutao Qin . Exponents of a class of special three-colored primitive digraphs with n vertices in graph theory. AIMS Mathematics, 2025, 10(4): 9415-9434. doi: 10.3934/math.2025435
    [6] Yunfeng Tang, Huixin Yin, Miaomiao Han . Star edge coloring of K2,t-free planar graphs. AIMS Mathematics, 2023, 8(6): 13154-13161. doi: 10.3934/math.2023664
    [7] Hongyu Chen, Li Zhang . A smaller upper bound for the list injective chromatic number of planar graphs. AIMS Mathematics, 2025, 10(1): 289-310. doi: 10.3934/math.2025014
    [8] Jin Cai, Shuangliang Tian, Lizhen Peng . On star and acyclic coloring of generalized lexicographic product of graphs. AIMS Mathematics, 2022, 7(8): 14270-14281. doi: 10.3934/math.2022786
    [9] Xiaoxue Hu, Jiangxu Kong . An improved upper bound for the dynamic list coloring of 1-planar graphs. AIMS Mathematics, 2022, 7(5): 7337-7348. doi: 10.3934/math.2022409
    [10] Jianxin Luo, Jiangxu Kong . An upper bound for the semistrong chromatic index of Halin graphs. AIMS Mathematics, 2025, 10(7): 15811-15820. doi: 10.3934/math.2025708
  • Zero forcing is a process of coloring in a graph in time steps known as propagation time. These graph-theoretic parameters have diverse applications in computer science, electrical engineering and mathematics itself. The problem of evaluating these parameters for a network is known to be NP-hard. Therefore, it is interesting to study these parameters for special families of networks. Perila et al. (2017) studied properties of these parameters for some basic oriented graph families such as cycles, stars and caterpillar networks. In this paper, we extend their study to more non-trivial structures such as oriented wheel graphs, fan graphs, friendship graphs, helm graphs and generalized comb graphs. We also investigate the change in propagation time when the orientation of one edge is flipped.


    Some human pathogens, such as hepatitis B virus (HBV), hepatitis C virus (HCV) and human immunodeficiency virus (HIV), possess the ability to suppress the immune response and damage the immune system, ultimately leading to the development of certain chronic diseases and even tumors [1]. In recent years, researchers have extensively studied numerous models pertaining to these infections [2,3,4]. In recent years, there has been notable progress in treating these infections. For example, the latest direct-acting antiviral drugs (DAA) for HCV can cure the majority of cases quickly, with cure rates reaching up to 95%[5,6]. For HBV, the optimal approach involves long-term medication to manage the condition, reduce viral transmission and stabilize host immune system[7,8]. Similarly, antiretroviral therapy (ART) is commonly used to suppress HIV replication[9]. Notably, a small group of elite controllers can naturally manage HIV without treatment[10]. However, challenges persist, including drug resistance and the risk of viral replication or recurrence[5,6,7,8,9]. In the past few years, immune checkpoint therapy, a method that shows significant efficacy in completely eliminating tumor cells for certain types of tumors, has been discovered, which activates the host's own immune system to inhibit tumor growth and spread [11,12,13]. Enhancing antiviral treatment by strengthening specific immune responses has become a hot topic among researchers.

    In 2003, Komarova et al. [1] first established an immunosuppressive infection model to investigate the relationship between drug efficacy and treatment duration, the model is given as follows:

    {dydt=ry(1yK)aypyz,dzdt=cyz1+ηybzqyz:=zf(y), (1.1)

    where y, z denote the population size of viruses and immune cells at time t, respectively. Viral reproduction follows logistic growth, r represents the viral replication rate at a viral load K. It is hypothesized that the virus population undergoes decay at a rate a, and immune cells eliminate the viruses at a rate pyz. Furthermore, immune cells are suppressed by the viruses at a rate qyz and die at a rate b. In the context of antigenic stimulation, the proliferation of immune cells relies on both the viruses and immune cells, and it can be quantified by the expansion function

    Φ(y,z)=cyz1+ηy,

    where η represents the inhibitory effects of the viruses on the proliferation of immune cells.

    In 2014, Shu et al. [14] introduced a time delay in the immune expansion function of the model described in [1] and studied the existence of Hopf bifurcation and periodic solutions. In 2017, Tian et al. [15] incorporated diffusion into the model presented in [14] and studied the direction and stability of Hopf bifurcations. However, viruses require a certain amount of time to effectively stimulate immune cells. We refer to this time period as the virus stimulation delay, denoted as τ. The immune cells at time t are essentially activated by the viruses at time tτ. Therefore, in 2021, Li et al. [16] studied the function as follows:

    Φ(yτ,z)=cy(tτ)z1+ηy(tτ).

    However, it has been noted by Wodarz (see pp 31, Chapter 2 in [17]) that the growth rate of immune cells, denoted as zf(y), in all of these models was assumed to be directly proportional to z, but nonlinearity can better capture the realistic dynamics. Thus, a more realistic model, proposed by De Boer and Perelson [18,19] in 1995–1997, has been suggested to provide a better representation of the actual dynamics. In this model, the immune expansion function is described by the well-known Beddington-DeAngelis functional response [20,21], the form of which is as follows:

    Φ(y,z)=cyz1+ηy+sz.

    It was initially used in predator-prey models [22,23,24]. The parameter s in the function represents the strength of competition between immune responses. This immune expansion function provides a more biologically realistic description of the interactions and competitive relationships between virus populations and immune cell populations within a host [25]. Thus the Beddington-DeAngelis functional response has also been widely discussed in immunosuppressive infection models [26,27,28].

    In 2023, Chen et al. [29] introduced the Beddington-DeAngelis functional response into the immunosuppressive infection model established by Komarova et al. [1], the model is as follows:

    {dydt=ry(1yK)aypyz,dzdt=cy(tτ)z1+ηy(tτ)+szbzqyz. (1.2)

    In [29], Chen et al. discussed the local and global stability of equilibrium points in the absence of delay, as well as the stability and Hopf bifurcation of model (1.2) when delay is present. They designed antiviral drug treatment strategies from the perspectives of adaptivity and feasibility of immune competition strength and viral suppression intensity.

    In practical problems, the spread of viruses within a host is possible. Thus we introduce the diffusion into the model proposed by Chen et al. in [29]. As far as we know, this task has not been attempted yet. Therefore, we study the reaction-diffusion system under the Neumann boundary condition as follows:

    {y(x,t)t=d1Δy(x,t)+ry(x,t)(1y(x,t)K)ay(x,t)py(x,t)z(x,t),z(x,t)t=d2Δz(x,t)+cy(x,tτ)z(x,t)1+ηy(x,tτ)+sz(x,t)bz(x,t)qy(x,t)z(x,t),νy(x,t)=νz(x,t)=0,xΩ,t>0,y(x,t)=ψ1(x,t),z(x,t)=ψ2(x,t),xΩ,t[τ,0], (1.3)

    where Ω=(0,lπ) with l>0, y=y(x,t) and z=z(x,t) denoting the population densities of viruses and immune cells at location x and time t, respectively. The parameters s and τ stand for the competition intensity between immune responses and the viral stimulation delay, respectively. The diffusion rates of viruses and immune cells are d1 and d2, respectively. The same parameters as in model (1.1) have the same biological meanings. Note that all parameters in (1.3) remain nonnegative.

    The remaining sections of the paper are organized as follows: Section 2 discusses the local asymptotic stability of the positive constant steady state solution and the existence of Hopf bifurcation. In Section 3, we derive the normal form of the Hopf bifurcation by using the multiple time scales (MTS) method. Section 4 presents numerical simulations to illustrate the theoretical analysis, and the conclusion is given in Section 5.

    In this section, we concentrate on studying the local asymptotic stability of the positive constant steady state solution and the existence of Hopf bifurcation for system (1.3). From a biological perspective, we only focus on the positive constant steady state solution (y,z). It is evident that (y,z) satisfies the following equations:

    ry(1yK)aypyz=0,cyz1+ηy+szbzqyz=0. (2.1)

    Then we can write y in terms of z as

    y=K(1a+pzr), (2.2)

    substituting Eq (2.2) into Eq (2.1), we obtain

    Az2+Bz+C=0,

    where

    A=qηK2p2+Kpqrs,B=2qη(ra)pK2+[r(bηc+q)p(ra)rqs]Kbr2s,C=qη(ra)2K2+r(bηc+q)(ar)Kbr2. (2.3)

    Denote

    Δ1=B24AC. (2.4)

    Lemma 1. For the positive constant steady state solution of system (1.3), we have the following results.

    (i) If Δ1>0, AC>0,AB<0,anda+pz<r holds, then system (1.3) has two positive constant steady state solutions, E1=(y1,z1) and E2=(y2,z2). Here z1=B+Δ12A, z2=BΔ12A, with y1=K(1a+pz1r) and y2=K(1a+pz2r);

    (ii) If Δ1>0, AC<0,anda+pz<r holds, then system (1.3) has one positive constant steady state solution E3=(y3,z3). Here z3=B+Δ12A, with y3=K(1a+pz3r);

    (iii) If Δ1=0, AB<0,anda+pz<r holds, then system (1.3) has one positive constant steady state solution E4=(y4,z4). Here z4=B2A, with y4=K(1a+pz4r).

    Without loss of generality, we denote E1, E2, E3, and E4 as E. Therefore, when the parameters of system (1.3) satisfy any one of the conditions (ⅰ)–(ⅱ) of Lemma 1, the positive constant steady state solution E=(y,z) exists.

    Now, we analyze the stability of E=(y,z). Set U(x,t)=(y(x,t),z(x,t))T, then the linearization of system (1.3) at E can be written as

    U(x,t)t=DΔU(x,t)+AU(x,t)+BU(x,tτ), (2.5)

    where

    D=(d100d2),A=(ryKpyqzzs(b+qy)2cy),B=(00z(b+qy)[cη(b+qy)]cy0).

    Therefore, the characteristic equation of the linearization of system (1.3) at E can be obtained as

    λ2+Tnλ+Nn+C0(eλτb0q)=0,nN0={0,1,2,3,...}, (2.6)

    with

    Tn=n2l2d1+n2l2d2+ryK+zs(b+qy)2cy,Nn=n4l4d1d2+n2l2d1zs(b+qy)2cy+n2l2d2ryK+rzs(b+qy)2Kc,C0=pyz,b0=(b+qy)[cη(b+qy)]cy.

    Lemma 2. When τ=0, the positive constant steady state solution E of system (1.3) is locally asymptotically stable if b0q>0.

    Proof. When τ=0, characteristic Eq (2.6) becomes

    λ2+Tnλ+Nn+C0(b0q)=0,nN0. (2.7)

    It is clear that Tn>0, Nn>0, thus if b0q>0 holds, the real parts of the roots of Eq (2.6) are negative. Thus, E is locally asymptotically stable.

    This completes the proof.

    Suppose that λ=iω(ω>0) is a root of Eq (2.6), then we have

    {C0b0cos(ωτ)=ω2+C0qNn,C0b0sin(ωτ)=Tnω. (2.8)

    Square and add above equations, then denote g=ω2. We obtain

    φ(g)=g2+(2C0q2Nn+T2n)g+(C0qNn)2C20b20=0. (2.9)

    Next, we will discuss the existence of ω. For the sake of convenience in our discussion, we present the following assumed conditions:

    (H1)Δ2=0,2C0q2Nn+T2n<0,(H2)Δ2>0,(C0qNn)2C20b20<0,(H3)Δ2>0,2C0q2Nn+T2n<0,(C0qNn)2C20b20>0,(H4)Δ2<0,(H5)Δ2>0,2C0q2Nn+T2n>0,(C0qNn)2C20b20>0,

    where Δ2=(2C0q2Nn+T2n)24[(C0qNn)2C20b20], and Nn,Tn,b0,C0 are given in Eq (2.6).

    Lemma 3. When τ>0, concerning the cases of positive roots of Eq (2.9), we have the following results:

    (i) Equation (2.9) has a unique positive root if (H1) or (H2) is satisfied, then the positive constant steady state solution of system (1.3) will not exhibit stability switches,

    (ii) Equation (2.9) has two positive roots if (H3) is satisfied, then the positive constant steady state solution of system (1.3) will exhibit stability switches,

    (iii) Equation (2.9) has no positive root if (H4) or (H5) is satisfied.

    Proof. If the condition (H1) is satisfied, we have g1=g2, g1+g2=(2C0q2Nn+T2n)>0. If the condition (H2) is satisfied, we have g1g2=(C0qNn)2C20b20<0. Both of them imply that there is a unique positive root g satisfying Eq (2.9). This proves (ⅰ).

    Similarly, we conclude (ⅱ)–(ⅲ).

    Considering the scenario with only a pair of purely imaginary roots, the unique positive ω that satisfies Eq (2.9) can be described as follows:

    ω+n=(2C0q2Nn+T2n)+Δ22,nN0. (2.10)

    Denote

    Rn=sin(ω+nτ)=Tnω+nC0b0,Qn=cos(ω+nτ)=(ω+n)2+C0qNnC0b0. (2.11)

    Now, define a set

    I1={nN0|such that (H1) or (H2) holds}. (2.12)

    From Eq (2.11), we obtain

    τ+n,j={1ω+n(arccos(Qn)+2jπ),forRn0,1ω+n(2πarccos(Qn)+2jπ),forRn<0.nI1,jN0. (2.13)

    From Eq (2.13), we find that {τ+n,j}j=0 is increasing about j for some fixed nI1. Thus, for fixed n, one has τ+n,0=minjN0{τ+n,j}. To investigate the stability of E, we define the smallest critical value of time delay

    τc=min{τ+n,0|nI1}. (2.14)

    Theorem 1. If b0q>0 and one of the two conditions (H1) or (H2) is satisfied, then the positive constant steady state solution E of system (1.3) is locally asymptotically stable when τ[0,τc), and unstable when τ(τc,). Moreover, system (1.3) undergoes the Hopf bifurcation at E when τ=τc, where τc is given in Eq (2.14).

    Proof. Differentiating the both sides of the characteristic Eq (2.6) with respect to time delay τ gives rise to

    (dλdτ)1=2λ+Tnλ[C0q(λ2+Tnλ+Nn)]τλ.

    As a result, one has

    (Redλdτ)1|τ=τ+n,j=Δ2T2n(ω+n,j)2+[NnC0q(ω+n,j)2]2>0,nI1,jN0. (2.15)

    This means that the transversality condition is satisfied.

    This completes the proof.

    We have obtained that under certain conditions, the characteristic Eq (2.6) of system (1.3) has a unique pair of purely imaginary roots λ=±iω+n, then system (1.3) undergoes the Hopf bifurcation at E when τ=τc, with τc as given at (2.14). In this section, we will utilize the method of MTS to derive the normal form of the Hopf bifurcation.

    Denote that

    f(y(x,tτ),z(x,t))=cy(x,tτ)z(x,t)1+ηy(x,tτ)+sz(x,t) (3.1)

    and

    L=1+ηy+sz,α=1+sz,β=1+ηy. (3.2)

    Perform a Taylor expansion of f(y(x,tτ),z(x,t)) at E, then we obtain

    f(y(x,tτ),z(x,t))=cyzL+czαL2(y(x,tτ)y)+cyβL2(z(x,t)z)cηzαL3(y(x,tτ)y)2+c(2ηyzs+ηy+α)L3(y(x,tτ)y)(z(x,t)z)csyβL3(z(x,t)z)2+cη2zαL4(y(x,tτ)y)3+cη[3szαL41+2szL3](y(x,tτ)y)2(z(x,t)z)+cs[1zs(1+2ηy)+η2y2]L4(y(x,tτ)y)(z(x,t)z)2+cs2yβL4(z(x,t)z)3+. (3.3)

    Now, we set ˆy(x,t)=y(x,τt)y, ˆz(x,t)=z(x,τt)z, and still denote ˆy(x,t) and ˆz(x,t) by y(x,t) and z(x,t), respectively. Then system (1.3) takes the form

    {y(x,t)t=τ[d1Δy(x,t)+r(y(x,t)+y)(1y(x,t)+yK)a(y(x,t)+y)p(y(x,t)+y)(z(x,t)+z)],z(x,t)t=τ[d2Δz(x,t)+czαL2y(x,t1)+cyβL2z(x,t)cηzαL3y2(x,t1)+c(2ηyzs+ηy+α)L3y(x,t1)z(x,t)csyβL3z2(x,t)+cη2zαL4y3(x,t1)+cη[3szαL41+2szL3]y2(x,t1)z(x,t)+cs[1zs(1+2ηy)+η2y2]L4y(x,t1)z2(x,t)+cs2yβL4z3(x,t)bz(x,t)q[y(x,t)z(x,t)+y(x,t)z+z(x,t)y]],νy(x,t)=νz(x,t)=0,xΩ,t>0,y(x,t)=ψ1(x,t),z(x,t)=ψ2(x,t),xΩ=(0,lπ),t[τ,0]. (3.4)

    Let h=(h11,h12)T be the eigenvector of the characteristic matrix of the linearization of system Eq (3.4) at E corresponding to the eigenvalue iωτ and h=d(h21,h22)T be the eigenvector of the adjoint matrix corresponding to the eigenvalue iωτ. They satisfy the inner product

    h,h=ˉhTh=1.

    By calculating, we have

    h11=1,h12=1py(iω+n2l2d1+rKy),h21=1py[iω+n2l2d2+zs(b+qy)2cy],h22=1,d=(h21ˉh11+h22ˉh12)1. (3.5)

    Consider the delay τ as a bifurcation parameter, let τ=τc+εμ, where τc is the critical point of the Hopf bifurcation, given as Eq (2.14). μ is a small perturbation parameter, and ε is the characteristic scale used for non-dimensionalization. The solution of Eq (3.4) is assumed to take the form:

    U(x,t)=U(x,T0,T1,T2,)=+k=1εkUk(x,T0,T1,T2,), (3.6)

    where

    U(x,T0,T1,T2,)=(y(x,T0,T1,T2,),z(x,T0,T1,T2,))T,Uk(x,T0,T1,T2,)=(yk(x,T0,T1,T2,),zk(x,T0,T1,T2,))T. (3.7)

    The derivative with respect to t is now transformed into

    t=T0+εT1+ε2T2+=D0+εD1+ε2D2+, (3.8)

    with the differential operator Di=Ti,iN0. Let

    yj=yj(x,T0,T1,T2,),yj,1=yj(x,T01,T1,T2,),zj=zj(x,T0,T1,T2,),zj,1=zj(x,T01,T1,T2,), (3.9)

    where jN={1,2,3,...}. According to Eqs (3.7) and (3.8), we get

    Ut=εD0U1+ε2D1U1+ε3D2U1+ε2D0U2+ε3D1U2+ε3D0U3+. (3.10)

    To deal with the delayed terms, we expand y(x,t1) at y(x,T01,T1,T2,), then we have

    y(x,t1)=εy11+ε2y21+ε3y31ε2D1y11ε3D1y21ε3D2y11+, (3.11)

    where yj1=yj(x,T01,T1,T2,),jN.

    Substituting Eqs (3.6)–(3.11) into Eq (3.4) will lead to some perturbation equations. By collecting the coefficients with respect to different powers of ε, we can obtain a set of ordered linear differential equations. First, for the ε-order terms, we have

    {D0y1τc[ry1ay1+d1Δy12ryy1Kpy1zpyz1]=0,D0z1τc[d2Δz1+czαy11L2+cyβz1L2bz1q(yz1+y1z)]=0, (3.12)

    where L,α,β are given by Eq (3.2). Thus the solution of Eq (3.12) can be expressed in the form as follows:

    y1=G(T1,T2,)eiωτcT0h11cos(nx)+c.c.,z1=G(T1,T2,)eiωτcT0h12cos(nx)+c.c., (3.13)

    where h11,h12 are given by Eq (3.5), and c.c. represents the complex conjugate of the previous term. Next, for the ε2-order terms, we obtain

    {D0y2τc[ry2ay2+d1Δy22ryy2Kpy2zpyz2]=D1y1τcry21Kμay1+μry1+μd1Δy1μpy1zμpyz12μryy1Kτcpy1z1,D0z2τc[d2Δz2+czαy21L2+cyβz2L2bz2q(yz2+y2z)]=τc[czαD1y11L2cηzαy211L3+c(2ηyzs+ηy+α)y11z1L3cysβz21L3qy1z1]+μ[d2Δz1+czαy11L2+cyβz1L2bz1q(yz1+y1z)]D1z1, (3.14)

    where L,α,β are given by Eq (3.2). We substitute solution Eq (3.13) into the right side of Eq (3.14), and denote the coefficient vector of eiωτcT0 as m1, which satisfies the solvability condition

    h,(m1,cos(nx))=0.

    Then one has

    GT1=MμG, (3.15)

    where

    M=ˉh21M1+ˉh22M2h11ˉh21+h12ˉh22+τczcαL2eiωτch11ˉh22, (3.16)

    with

    M1=ah11+rh11n2d1h11pzh11pyh122ryKh11,M2=n2d2h12+czαL2eiωτch11+cyβL2h12bh12qyh12qzh11.

    Suppose that the solution of Eq (3.14) is

    {y2=+k=0(η0,kGˉG+η1,ke2iωτcT0G2+ˉη1,ke2iωτcT0ˉG2)cos(kx),z2=+k=0(ζ0,kGˉG+ζ1,ke2iωτcT0G2+ˉζ1,ke2iωτcT0ˉG2)cos(kx). (3.17)

    Let Ck=cos(nx)cos(nx),cos(kx)=lπ0cos(nx)cos(nx)cos(kx)dx. Besides, we have

    lπ0cos(nx)cos(kx)dx={0kn,lπ2k=n0,lπk=n=0.

    We substitute solution Eqs (3.13) and (3.17) into Eq (3.14), then make the inner product with cos(kx) on both sides of the obtained equation, where kN0. By adding the obtained results, we can compare the coefficients in front of terms G2 and GˉG, then we get

    η1,k=Ck(W(J1+2iω)J2py)lπ0cos(kx)cos(kx)dx[(2iω+J1)(l1+2iω)pyczαL2e2iωτcqz],η0,k=Ck(SJ1l2py)lπ0cos(kx)cos(kx)dx[J1l1pyczαL2qz],ζ1,k=Cklπ0cos(kx)cos(kx)dxJ2pyl1+2iωpyη1,k,ζ0,k=Cklπ0cos(kx)cos(kx)dxl2pyl1pyη0,k, (3.18)

    where

    W=cηzαL3e2iωτch211+c(2ηyzs+ηy+α)L3eiωτccysβL3h212qh11h12,S=cηzαL32h11ˉh11+c(2ηyzs+ηy+α)L3(eiωτch11ˉh12+eiωτcˉh11h12)cysβL32h12ˉh12q(h11ˉh12+h12ˉh11),J1=d2k2cyβL2+b+qy,J2=rKh211ph11h12,l1=ar+2ryK+pz+d1k2,l2=rK2h11ˉh11p(h11ˉh12+h12ˉh11),L=1+ηy+sz,α=1+sz,β=1+ηy.

    Specifically, if n=0, the condition for Ck0 is k=0. Thus we consider n=k=0, then we have Ck=lπ, lπ0cos(nx)cos(kx)dx=lπ. Thus we can simply get the value of η1,0, η0,0, ζ1,0, ζ0,0. For the ε3-order terms, we get

    {D0y3+τc(ar+2ryK+pz)y3d1τcΔy3+τcpyz3=μd1Δy2μ(ar+2ryK+pz)y2D2y1D1y2μpy1z1μpyz2μry21K2τcry2y1Kτcpy1z2τcpy2z1,D0z3τcd2Δz3τcczαL2y31(τccyβL2bqy)z3=D2z1D1z2+τc[czαD1y21czαD2y11L22cηzαD1y211+2cηzαy11y21L3+c(2ηyzs+ηy+α)y11z2L3+c(2ηyzs+ηy+α)(D1y11z1+y21z1)L32cysβz1z2L3+c(1+z(2ηy1)s)sy11z21L4+cys2βz31L4+czη2αy311L4+cη[3szαL42sz+1L3]y211z1qy1z2qy2z1]+μ[d2Δz2+czαD1y11+czαy21L2+(cyβL2bqy)z2cηzαy211L3+c(2ηyzs+ηy+α)y11z1L3cyβsz21L3qy1z1qy2z], (3.19)

    where L,α,β are given by Eq (3.2). We substitute solution Eqs (3.13) and (3.17) into Eq (3.19), then denote the coefficient vector of eiωτcT0 as m2, which satisfies the solvability condition

    h,(m2,cos(nx))=0.

    Since μ is the disturbance parameter and μ2 has very little influence on the result, the items containing μ2G can be ignored. Then we obtain

    GT2=χG2ˉG,

    where

    χ=RCk+QV, (3.20)

    with

    R=ˉh21[2τcrK(h11+k=0η0,k+ˉh11+k=0η1,k)τcp(h11+k=0ζ0,k+ˉh11+k=0ζ1,k+h12+k=0η0,k+ˉh12+k=0η1,k)]+ˉh22τc[2cηzαL3(eiωτch11+k=0η0,k+eiωτcˉh11+k=0η1,k)+c(2ηyzs+ηy+α)L3(eiωτch11+k=0ζ0,k+eiωτcˉh11+k=0ζ1,k+h12+k=0η0,k+ˉh12+k=0η1,k)2cysβL3(h12+k=0ζ0,k+ˉh12+k=0ζ1,k)],Q=ˉh22τc[c(1+z(2ηy1)s)sL4(2eiωτch11h12ˉh12+eiωτcˉh11h212)+cyβs2L4(2h212ˉh12)+czη2αL4(2eiωτch211ˉh11)+cη[3szαL42sz+1L3](e2iωτch211ˉh12+2h11ˉh11h12)q(h11+k=0ζ0,k+ˉh11+k=0ζ1,k+h12+k=0η0,k+ˉh12+k=0η1,k)]lπ0cos4(nx)dx,V=[h11ˉh21+h12ˉh22+τcczαL2eiωτch11ˉh22]lπ0cos(nx)cos(nx)dx,L=1+ηy+sz,α=1+sz,β=1+ηy.

    Therefore, the normal form of the Hopf bifurcation for system Eq (1.3) reduced on the center manifold is

    ˙G=MμG+χG2ˉG, (3.21)

    where M and χ are given by Eqs (3.16) and (3.20), respectively. Now introduce G=reiθ and put it into the normal form Eq (3.21), then one has

    {˙r=Re(M)μr+Re(χ)r3,˙θ=Im(M)μ+Im(χ)r2. (3.22)

    For the stability of the periodic solutions of system (1.3), we have the following result.

    Theorem 2. When Re(M)μRe(χ)<0, system Eq (1.3) exists as periodic solutions near the positive constant steady state solution E. Moreover,

    (i) If Re(M)μ<0, the bifurcating periodic solutions reduced on the center manifold are unstable, and when μ>0 (μ<0), the direction of bifurcation is forward (backward);

    (ii) If Re(M)μ>0, the bifurcating periodic solutions reduced on the center manifold are stable, and when μ>0 (μ<0), the direction of bifurcation is forward (backward).

    In this section, we will perform numerical simulations for different values of τ. Based on the work of Komarova et al. [1], we fix

    a=3.0day1,p=1.0mm3cells1day1,q=1.0mm3virus1day1,r=6.0day1.

    For the remaining parameters, we consider the actual meaning of the parameters and choose

    d1=0.4mm3day1,K=10virusmm3,b=2day1,d2=0.1mm3day1,c=6.5mm3virus1day1,η=0.35mm1virus,s=0.02mm3immunecells1,l=2.

    It is easy to check that the condition (ⅱ) of Lemma 1 is satisfied, which implies that the system (1.3) has only one positive constant steady state solution,

    E=(y,z)=(0.4600,2.7239).

    In addition, we can also find that this set of parameters satisfies the condition (H2) of Lemma 3. Furthermore, by Eqs (2.10), (2.12) and (2.14), we obtain

    ω+0=2.1289,I1={n=0,1,2,3},τc=τ+0,0=0.0667.

    Therefore, by Theorem 2.1, E is locally asymptotically stable when τ[0,τc)=[0,0.0667), and unstable when τ(0.0667,+). Take τ=0.055[0,0.0667) as an example, then E is locally asymptotically stable as is illustrated by the computer simulations (see Figure 1).

    Figure 1.  Choose the initial function as (0.4600+0.0001cosx,2.7239+0.002cosx), then the constant steady state solution E=(0.4600,2.7239) of system (1.3) is locally asymptotically stable when τ=0.055<τc.

    Biologically, when the virus stimulation delay does not exceed the threshold, immune cells can control the virus in the long-term at low levels. Although the virus is not completely eliminated, it does not cause serious disease deterioration, this is the ideal scenario that we hope for.

    According to Eqs (3.16) and (3.20), we can get

    M=0.387410.6020i,χ=0.42210.2652i,

    which implies that Re(M)>0, and Re(χ)<0. From Theorem 2, when μ>0, system (1.3) has stable bifurcating periodic solutions. In addition, the direction of the Hopf bifurcation is forward. Take τ=0.067(0.0667,+) as an example, then system (1.3) generates a spatially homogeneous stable periodic solution (see Figure 2, where Figure 2(a), (c) depict simulations of the virus and immune cell dynamics, respectively, while Figure 2(b), (d) represent their respective local zoom-in plots).

    Figure 2.  Choose the initial function as (0.4600+0.0001cosx, 2.7239+0.002cosx), then system (1.3) has a spatially homogeneous stable periodic solution when τ=0.067>τc.

    Biologically, when the virus stimulation delay slightly exceeds the threshold τc, the densities of immune cells and virus exhibits periodic oscillations. As the number of immune cells increases, they effectively control the virus and reduce the viral load. After a while, an immunoevasion mechanism for the virus may emerge and cause the number of immune cells to drop, allowing the virus to remultiply. This scenario is not ideal, but it is still manageable.

    Once τ significantly exceeds the threshold, take τ=0.076(0.0667,+) as an example, the virus proliferates rapidly (see Figure 3). This is the worst-case scenario, as the virus could spread rapidly and pose significant risks to the host.

    Figure 3.  Choose the initial function as (0.4600+0.0001cosx, 2.7239+0.002cosx), then the constant steady state solution E of system (1.3) is unstable when τ=0.076>τc.

    In this paper, we incorporated diffusion into an immunosuppressive infection model with virus stimulation delay and Beddington-DeAngelis functional response. We studied the stability of the constant steady state solution of the model and provided conditions for the occurrence of Hopf bifurcation. Moreover, we analyzed the stability and direction of the Hopf bifurcation, and further derived the normal form. By analyzing the model, we obtained a threshold τc for the virus stimulation delay τ and discovered that when τ was less than this threshold, the constant steady state solution was locally asymptotically stable. However, once τ exceeded the threshold, the constant steady state solution lost stability, and the model exhibited stable periodic solutions. This implied that when the virus could quickly stimulate immune cells, the host possessed strong immunity, which could inhibit virus growth and spread. Conversely, rapid virus proliferation posed significant risks to the host. The discovery of the threshold could help biologists predict, judge, and control the growth of virus cells and select effective treatment strategies.

    The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

    This study was funded by Fundamental Research Funds for the Central Universities of China (Grant No.2572022DJ07).

    The authors declare there is no conflicts of interest.



    [1] A. Aazami, Hardness results and approximation algorithms for some problems on graphs, Thesis (Ph.D.), University of Waterloo, 2008.
    [2] F. Barioli, W. Barrett, S. Butler, S. M. Cioabă, D. Cvetković, S. M. Fallat, et al., Zero forcing sets and the minimum rank of graphs, Linear Algebra Appl., 428 (2007), 1628-1648.
    [3] F. Barioli, W. Barrett, S. Fallat, H. T. Hall, L. Hogben, B. Shader, et al., Parameters related to tree-width, zero forcing, and maximum nullity of a graph, J. Graph Theory, 72 (2013), 146177.
    [4] F. Barioli, S. M. Fallat, H. T. Hall, D. Hershkowitz, L. Hogben, H. van der Holst, et al., On the minimum rank of not necessarily symmetric matrices: a preliminary study, Electron. J. Linear Algebra, 18 (2009), 126-145.
    [5] K. F. Benson, D. Ferrero, M. Flagg, V. Furst, L. Hogben, V. Vasilevska, et al., Power domination and zero forcing, arXiv preprint arXiv: 1510.02421, 2015.
    [6] A. Berliner, C. Bozeman, S. Butler, M. Catral, L. Hogben, B. Kroschel, et al., Zero forcing propagation time on oriented graphs, Discrete Appl. Math., 224 (2017), 45-59. doi: 10.1016/j.dam.2017.02.017
    [7] A. Berliner, C. Brown, J. Carlson, N. Cox, L. Hogben, J. Hu, et al., Path cover number, maximum nullity, and zero forcing number of oriented graphs and other simple digraphs, Involve, 8 (2015), 147-167. doi: 10.2140/involve.2015.8.147
    [8] A. Berliner, M. Catral, L. Hogben, M. Huynh, K. Lied, M. Young. Mini- mum rank, maximum nullity and zero forcing number of simple digraphs, Electron. J. Linear Algebra, 26 (2013), 762- 780.
    [9] D. Burgarth, D. D'Alessandro, L. Hogben, S. Severini, M. Young, Zero forcing, linear and quantum controllability for systems evolving on networks, IEEE Trans. Autom. Control, 58 (2013), 2349- 2354.
    [10] D. Burgarth, V. Giovannetti, Full control by locally induced relaxation, Phys. Rev. Lett., 99 (2007), 100501. doi: 10.1103/PhysRevLett.99.100501
    [11] K. B. Chilakamarri, N. Dean, C. X. Kang, E. Yi, Iteration index of a zero forcing set in a graph, Bull. Inst. Combin. Appl., 64 (2012), 57-72.
    [12] S. Fallat, L. Hogben, Minimum Rank, maximum nullity, and zero forcing number of graphs. In: Handbook of Linear Algebra, 2nd edition, L. Hogben editor, CRC Press, Boca Raton, 2014.
    [13] S. Fallat, K. Meagher, B. Yang, On the complexity of the positive semidefinite zero forcing number, Linear Algebra Appl., 491 (2016), 101-122. doi: 10.1016/j.laa.2015.03.011
    [14] L. Hogben, Minimum rank problems, Linear Algebra Appl., 432 (2010), 1961-1974.
    [15] L. Hogben, M. Huynh, N. Kingsley, S. Meyer, S. Walker, M. Young, Propagation time for zero forcing on a graph, Discrete Appl. Math., 160 (2012), 1994-2005. doi: 10.1016/j.dam.2012.04.003
    [16] I. J. Kim, B. P. Barthel, Y. Park, J. R. Tait, J. L. Dobmeier, S. Kim, D. Shin, Network analysis for active and passive propagation models, Networks, 63 (2014), 160-169. doi: 10.1002/net.21532
    [17] M. Perila, K. Rutschke, B. Kroschel, Zero forcing and propagation time of oriented graphs, Submitted.
  • This article has been cited by:

    1. K. Ishwarya, K. G. Rani, K. Appathurai, R. Surendran, R. Selvanarayanan, C. Y. Lau, 2024, 3161, 0094-243X, 020301, 10.1063/5.0229271
    2. Ishwarya K, Saranya G, Appathurai K, Surendran R, 2023, Exploring the Significance and Paradoxical Nature of the Zero Value in Mathematics Using Artificial Intelligence, 979-8-3503-0088-8, 861, 10.1109/ICOSEC58147.2023.10276215
    3. Shiying Wang, Hongmei Li, Lina Zhao, The m-Component Connectivity of Leaf-Sort Graphs, 2024, 12, 2227-7390, 404, 10.3390/math12030404
    4. Krittawit Limkul, Sayan Panma, On the Independence Number of Cayley Digraphs of Clifford Semigroups, 2023, 11, 2227-7390, 3445, 10.3390/math11163445
    5. S. T. Vikram, S. Balaji, An Analysis of the Factors Influencing the Strong Chromatic Index of Graphs Derived by Inflating a Few Common Classes of Graphs, 2023, 15, 2073-8994, 1301, 10.3390/sym15071301
    6. J. Anitha, Indra Rajasingh, Hossein Rashmanlou, Propagating Sets in Sierpiński Fractal Graphs, 2025, 59, 0988-3754, 1, 10.1051/ita/2024016
  • Reader Comments
  • © 2021 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(3750) PDF downloads(281) Cited by(6)

Figures and Tables

Figures(22)

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog