Processing math: 56%
Research article Special Issues

TKRD: Trusted kernel rootkit detection for cybersecurity of VMs based on machine learning and memory forensic analysis

  • The promotion of cloud computing makes the virtual machine (VM) increasingly a target of malware attacks in cybersecurity such as those by kernel rootkits. Memory forensic, which observes the malicious tracks from the memory aspect, is a useful way for malware detection. In this paper, we propose a novel TKRD method to automatically detect kernel rootkits in VMs from private cloud, by combining VM memory forensic analysis with bio-inspired machine learning technology. Malicious features are extracted from the memory dumps of the VM through memory forensic analysis method. Based on these features, various machine learning classifiers are trained including Decision tree, Rule based classifiers, Bayesian and Support vector machines (SVM). The experiment results show that the Random Forest classifier has the best performance which can effectively detect unknown kernel rootkits with an Accuracy of 0.986 and an AUC value (the area under the receiver operating characteristic curve) of 0.998.

    Citation: Xiao Wang, Jianbiao Zhang, Ai Zhang, Jinchang Ren. TKRD: Trusted kernel rootkit detection for cybersecurity of VMs based on machine learning and memory forensic analysis[J]. Mathematical Biosciences and Engineering, 2019, 16(4): 2650-2667. doi: 10.3934/mbe.2019132

    Related Papers:

    [1] Qi Cui, Ruohan Meng, Zhili Zhou, Xingming Sun, Kaiwen Zhu . An anti-forensic scheme on computer graphic images and natural images using generative adversarial networks. Mathematical Biosciences and Engineering, 2019, 16(5): 4923-4935. doi: 10.3934/mbe.2019248
    [2] Zairong Wang, Xuan Tang, Haohuai Liu, Lingxi Peng . Artificial immune intelligence-inspired dynamic real-time computer forensics model. Mathematical Biosciences and Engineering, 2020, 17(6): 7221-7233. doi: 10.3934/mbe.2020370
    [3] Liang Yu, Zhengkuan Zhang, Yangbing Lai, Yang Zhao, Fu Mo . Edge computing-based intelligent monitoring system for manhole cover. Mathematical Biosciences and Engineering, 2023, 20(10): 18792-18819. doi: 10.3934/mbe.2023833
    [4] Giuseppe Ciaburro . Machine fault detection methods based on machine learning algorithms: A review. Mathematical Biosciences and Engineering, 2022, 19(11): 11453-11490. doi: 10.3934/mbe.2022534
    [5] Gayathri Vivekanandhan, Mahtab Mehrabbeik, Karthikeyan Rajagopal, Sajad Jafari, Stephen G. Lomber, Yaser Merrikhi . Applying machine learning techniques to detect the deployment of spatial working memory from the spiking activity of MT neurons. Mathematical Biosciences and Engineering, 2023, 20(2): 3216-3236. doi: 10.3934/mbe.2023151
    [6] Sadia Anjum, Lal Hussain, Mushtaq Ali, Adeel Ahmed Abbasi, Tim Q. Duong . Automated multi-class brain tumor types detection by extracting RICA based features and employing machine learning techniques. Mathematical Biosciences and Engineering, 2021, 18(3): 2882-2908. doi: 10.3934/mbe.2021146
    [7] Liang Tian, Fengjun Shang, Chenquan Gan . Optimal control analysis of malware propagation in cloud environments. Mathematical Biosciences and Engineering, 2023, 20(8): 14502-14517. doi: 10.3934/mbe.2023649
    [8] Zhongwei Li, Wenqi Jiang, Xiaosheng Liu, Kai Tan, Xianji Jin, Ming Yang . GAN model using field fuzz mutation for in-vehicle CAN bus intrusion detection. Mathematical Biosciences and Engineering, 2022, 19(7): 6996-7018. doi: 10.3934/mbe.2022330
    [9] Keyue Yan, Tengyue Li, João Alexandre Lobo Marques, Juntao Gao, Simon James Fong . A review on multimodal machine learning in medical diagnostics. Mathematical Biosciences and Engineering, 2023, 20(5): 8708-8726. doi: 10.3934/mbe.2023382
    [10] Yufeng Qian . Exploration of machine algorithms based on deep learning model and feature extraction. Mathematical Biosciences and Engineering, 2021, 18(6): 7602-7618. doi: 10.3934/mbe.2021376
  • The promotion of cloud computing makes the virtual machine (VM) increasingly a target of malware attacks in cybersecurity such as those by kernel rootkits. Memory forensic, which observes the malicious tracks from the memory aspect, is a useful way for malware detection. In this paper, we propose a novel TKRD method to automatically detect kernel rootkits in VMs from private cloud, by combining VM memory forensic analysis with bio-inspired machine learning technology. Malicious features are extracted from the memory dumps of the VM through memory forensic analysis method. Based on these features, various machine learning classifiers are trained including Decision tree, Rule based classifiers, Bayesian and Support vector machines (SVM). The experiment results show that the Random Forest classifier has the best performance which can effectively detect unknown kernel rootkits with an Accuracy of 0.986 and an AUC value (the area under the receiver operating characteristic curve) of 0.998.


    According to clinical data, drug treatment against persistent human infections sometimes fails to consistently eradicate the infection from the host [1,2,3], such as human immunodeficiency virus (HIV), hepatitis B virus (HBV), and hepatitis C virus (HCV). For all these viruses, drug treatment is not effective, and for HIV, lifelong therapy is commonly required to control viral replication [4,5]. The main reason is that these viruses are able to suppress and impair immune responses, resulting in persistent infection [6,7]. An alternative strategy is to use drug therapies to boost virus-specific immune response and then induce sustained viral suppression in the absence of drugs. Antiviral therapy aimed at boosting specific immune responses has attracted more and more attentions recently. In 2003, Komarova et al. [6] used mathematical models to study immune response dynamics during therapy in the context of immunosuppressive infections. They showed that a single phase of antiviral drug therapy is also possible to establish sustained immunity. They assumed that the degree of immune expansion depends on virus load and the response inhibits virus growth. This assumption is natural for any branch of the adaptive immune system, such as CD4 T cells, CD8 T cells, or antibodies. Denote by y(t) and z(t), respectively, the virus population and immune cell population at time t. Komarova et al. [6] proposed the following system of ordinary differential equations:

    y(t)=ry(t)(1y(t)K)py(t)z(t)by(t),z(t)=cy(t)z(t)1+dy(t)qy(t)z(t)μz(t). (1.1)

    The virus population grows at a rate described by the logistic function. The viral replication rate is a linear decreasing function of viral loads with a maximum value of r and it vanishes at a viral load K. The treatment is modeled as a reduction in r. Immune cells kill virus at a rate pyz. b and μ are the natural decay rates of viral particles and immune cells, respectively. Immune expansion is modeled by the growth function cyz1+dy. When the virus load is low, the level of immune response is simply proportional to both the viral load and the immune cells. The immune response saturates when the virus load is sufficiently high. Immune cells are inhibited by the virus at a rate qyz. Komarova et al. presented a simple relationship between the timing of therapy and efficacy of the drugs. In the presence of strong viral suppression, they showed that therapy should be stopped relatively early because a longer duration of treatment may lead to a failure. On the other hand, in the presence of weaker viral suppression, stopping treatment too early is detrimental, and therapy has to be continued beyond a time threshold. Essentially, model (1.1) has two stable equilibria: a virus dominant equilibrium with no sustained immunity and an immune control equilibrium with sustained immunity. This bistability allows a solution from the basin of the attraction of the virus dominant equilibrium to be lifted to that of the immune control equilibrium via a single phase of therapy.

    Model (1.1) assumes that the immune response is instantaneous. Nevertheless, it has been established that the immune response process involves a sequence of events such as antigenic activation, selection, and proliferation of the immune cells [8]. Moreover, oscillatory viral loads and immune cells were observed from clinical data [9], while the oscillatory behavior is not exhibited in model (1.1). A natural question is what is the cause of the oscillations. In 2014, Shu et al. [10] incorporated the time lag during the immune response process into model (1.1). By studying the dynamics between an immunosuppressive infection and antiviral immune response, they demonstrated that the time lag is the main factor causing oscillations. Here, τ is denoted as the time lag for the immune system to trigger a sequence of events. The activation rate of immune cells at time t depends on the virus load and the number of immune cells at time tτ. The above assumptions lead to the following system:

    y(t)=ry(t)(1y(t)K)py(t)z(t)by(t),z(t)=cy(tτ)z(tτ)1+dy(tτ)qy(t)z(t)μz(t). (1.2)

    It was shown that the delayed antiviral immune response may induce sustained periodic oscillations, transient oscillations and even sustained aperiodic oscillations (chaos), and the model can admit two stable equilibria and can also allow a stable equilibrium to coexist with a stable periodic solution.

    In the aforementioned work, the mortality of the population during the time lag has been ignored. Consideration of the survival probability during the time lag requires an additional delay-dependent multiplier in the nonlinear delayed feedback term. The death of the free viral particles depends on the mature immune cells instead of all the newly generated immune cells. We denote z(t) as the population of mature immune cells at time t, and consider the following delay differential equation with a delay-dependent coefficient:

    y(t)=f(y(t))py(t)z(t)by(t),z(t)=eδτh(y(tτ))z(tτ)qy(t)z(t)μz(t). (1.3)

    Here, f(y) describes the virus population growth function; δ>0 is the death rate of the immature immune cells during the time lag; eδτ represents the survival probability of the immature immune cells surviving τ time units before becoming mature; h(y)z denotes the immune cell growth function; the other parameters have the same meaning as in model (1.2). In addition to more general virus and immune cell growth functions, model (1.3) differs from model (1.2) in the sense that it introduces an additional delay-dependent survival probability eδτ. As we shall see later, this quantity will make a significant difference in the bifurcation analysis. For instance, model (1.3) does not possess any positive equilibrium for sufficiently large τ and its global Hopf branches always have bounded τ components. However, for model (1.2), the existence condition of positive equilibrium is independent of τ and the global Hopf branches always have unbounded τ components. Furthermore, there exists only one stability switch of positive equilibrium for model (1.2), but the stability of positive equilibrium for our model (1.3) may switch multiple times.

    Actually, the model (1.3) can be derived from a stage structured population model for u(t,a) as below

    y(t)=f(y(t))py(t)z(t)by(t),tu(t,a)+au(t,a)=(μ(a)+q(a)y(t))u(t,a),

    with the stage-specific decay rate of immune cells, μ(a), inhibition rate of immune cells, q(a),

    μ(a)={μ,aτ,δ,a<τ,    q(a)={q,aτ,0,a<τ,

    where u(t,a) is the population of immune cells at age a and time t. Note that the population of mature immune cells is z(t)=τu(t,a)da. We integrate along characteristic lines to find

    z(t)=u(t,τ)u(t,)qy(t)z(t)μz(t)=u(tτ,0)eδτu(t,)qy(t)z(t)μz(t).

    Note that u(t,)=0 and u(tτ,0)=h(y(tτ))z(tτ). The equation for z(t) follows.

    Throughout this paper, we assume that the virus population growth function f(y) and the immune expansion rate h(y) satisfy the following conditions:

    (H1) f(y) is continuously differentiable, f(0)=0 and f(0)>b; there exists ˉy>0 such that f(ˉy)=0, and (yˉy)f(y)<0 for positive yˉy; f(y)<0 for y[0,ˉy].

    (H2) h(y) is continuously differentiable, h(0)=0 and h(0)>q; h(y)>0 and h(y)<0 for all y0.

    The assumptions (H1)-(H2) are biological conditions for infection dynamics. For instance, f(0)=0 means no new infection in the absence of virus; f(0)>b indicates that the intrinsic growth rate is greater than the decay rate. Similarly, h(0)=0 implies that immune cells cannot reproduce without the presence of virus; h(0)>q guarantees that the initial expansion rate is greater than the inhibition rate.

    The rest of this paper is organized as follows. We first present some preliminary results concerning the well-posedness of model (1.3), and describe global dynamics of the ordinary differential equation system (2.4) in Section 2. Local and global stability analysis of equilibria as well as local and global Hopf bifurcation analysis are given in Section 3. In Section 4, numerical simulations based on the bifurcation analysis are reported and discussed. Finally, we give a brief summary in Section 5.

    To establish the well-posedness of model (1.3), we choose the phase space C×C, where C is the Banach space of continuous functions on [τ,0] defined by C:={ϕ:[τ,0]R is continuous}, and the norm is defined by ϕ=supτθ0ϕ(θ). The nonnegative cone of C is denoted by C+. As usual, ϕtC is defined by ϕt(θ)=ϕ(t+θ) for θ[τ,0]. For biological applications, the initial condition of (1.3) is given as

    (y0,z0)X:=C+×C+.

    The existence and uniqueness of the solution of model (1.3) follows from the standard theory of functional differential equations [11]. Using the same manner as in proving [12,Proposition 2.1], we can easily show that the solution of (1.3) with initial conditions in X is nonnegative, which implies that X is positively invariant under the solution map of (1.3).

    Next, we shall prove that all solutions of (1.3) are ultimately bounded. It follows from the first equation in (1.3) that y(t)f(y), which yields lim supty(t)ˉy. By the nonnegativity of the solution of (1.3) and (H1)-(H2), we obtain

    (y(t)+ph(0)z(t+τ))Mf+(eδτh(y)h(0)y1)py(t)z(t)by(t)μph(0)z(t+τ)Mfγ(y(t)+ph(0)z(t+τ)),

    where Mf=maxy0f(y) and γ=min{b,μ}. Thus,

    lim supt(y(t)+ph(0)z(t+τ))Mfγ,

    which implies that y(t) and z(t) are ultimately bounded in X. To summarize, we obtain the following proposition.

    Proposition 2.1. The solutions of model (1.3) with the initial conditions in X are nonnegative, and the region

    Γ={(ϕ,ψ)X:ϕˉy, ϕ(τ)+ph(0)ψ(0)Mfγ}

    is positively invariant and absorbing in X; namely, all solutions ultimately enter Γ.

    Clearly, model (1.3) always has a trivial equilibrium E0=(0,0). Since f(0)>b, there also exists a boundary equilibrium Eb=(˜y,0), where ˜y is the unique positive root of f(y)=by. We call Eb the virus dominant equilibrium (VDE). Assume that E=(y,z) is a positive equilibrium with y>0 and z>0, then it must satisfy the following equilibrium equations

    eδτh(y)=qy+μ,  z=f(y)bypy. (2.1)

    Obviously, z>0 if and only if y<˜y. Hence, E exists if and only if y is a positive root of

    h(y)y+μ/q=qeδτ (2.2)

    with y<˜y. It is easily to calculate that sign(h(y)y+μ/q)=sign g(y), where g(y)=h(y)(y+μq)h(y). Note that g(0)>0,g()<0 and g(y)=h(y)(y+μ/q)<0. There exists a unique yc>0 such that g(yc)=0, and (yyc)g(y)<0 for yyc. Therefore, h(y)y+μ/q is strictly increasing on [0,yc), and strictly decreasing on (yc,). We denote

    s:=h(yc)yc+μ/q=supy0h(y)y+μ/q. (2.3)

    If s<qeδτ, (2.2) has no positive root. If sqeδτ, (2.2) has exactly two positive roots (counting multiplicity) y1 and y2 such that y1ycy2. These two roots coincide (i.e., y1=y2=yc) if and only if s=qeδτ. Throughout this paper, we make the following assumption to ensure the existence of positive roots for (2.2).

    (H3) s>q and ττmax, where τmax=1δlnsq.

    If (H3) is violated, there does not exist any positive equilibrium for our model and the boundary equilibrium is locally asymptotically stable. Summarizing the above analysis, we obtain the following result.

    Proposition 2.2. Assume that (H1)-(H3) are satisfied.

    (i) If y1˜y holds, then there are two equilibria: E0=(0,0) and Eb=(˜y,0).

    (ii) If y1<˜yy2 or  y1=y2<˜y holds, then there are three equilibria: E0, Eb and E1=(y1,z1).

    (iii) If y1<y2<˜y holds, then there are four equilibria: E0, Eb, E1 and E2=(y2,z2).

    When τ=0, model (1.3) reduces to an ODE system which generalizes (1.1)

    y(t)=f(y(t))py(t)z(t)by(t),z(t)=h(y(t))z(t)qy(t)z(t)μz(t). (2.4)

    We next provide a complete description about its global dynamics.

    Lemma 2.3. System (2.4) has no closed orbits.

    Proof. Using an argument similar to Proposition 2.1, we can show that

    Γ0={(ϕ,ψ)R2+:ϕˉy, ϕ+ph(0)ψMfγ}

    is positively invariant and absorbing in R2+. Thus, it suffices to show that System (2.4) has no closed orbits in Γ0.

    Let (P(y,z),Q(y,z)) denote the vector field of (2.4), then we have

    (BP)y+(BQ)z=1z(t)(f(y)y),  where  B(y,z)=1yz.

    Note that (H1) implied that (f(y)y)<0 for all yˉy. Thus, (BP)y+(BQ)z<0 in Γ0. The nonexistence of closed orbits for (2.4) follows from the classical Dulac-Bendixson criterion. This ends the proof.

    We linearize system (2.4) at E0=(0,0) and calculate that the two eigenvalues are f(0)b>0 and μ<0. Thus, E0 is a saddle point. For the linearized system at Eb=(˜y,0), the two eigenvalues are

    λ1=f(˜y)b<0  and  λ2=h(˜y)q˜yμ.

    From the proof of Proposition 2.2, we know h(˜y)q˜yμ<0 if ˜y<y1 or ˜y>y2, and h(˜y)q˜yμ>0 provided that y1<˜y<y2. This shows that Eb is stable if y1>˜y or y2<˜y, and Eb is a saddle point if y1<˜y<y2. For the linearized system at E1=(y1,z1), the associated characteristic equation is

    λ2+(f(y1)+pz1+b)λ+py1z1(h(y1)q)=0.

    It follows from (2.1) and (H1) that f(y1)+pz1+b=f(y1)+f(y1)/y1>0. This, together with h(y1)q>0, implies that all eigenvalues must have negative real parts, and hence E1 is asymptotically stable. At E2=(y2,z2), the two eigenvalues λ1 and λ2 satisfy

    λ1+λ2=f(y2)f(y2)y2<0,  λ1λ2=py2z2(h(y2)q)<0.

    This implies that one eigenvalue must be positive and the other negative. Thus, E2 is a saddle point whenever it exists. Summarizing the above analysis, we obtain the following global stability results on model (2.4).

    Theorem 2.4. Assume that (H1)-(H2) are satisfied and s>q.

    (i) If y1>˜y holds, then E0 is a saddle point and Eb attracts all solutions with initial conditions in {(y(t),z(t))R2+:y(0)>0,z(0)0}.

    (ii) If y1<˜y<y2 holds, then both E0 and Eb are saddle points and E1 attracts all solutions with initial conditions in {(y(t),z(t))R2+:y(0)>0,z(0)>0}.

    (iii) If y1<y2<˜y holds, then both E0 and E2 are saddle points, and for any initial condition (y(0),z(0)) with y(0)>0, the solution of (2.4) approaches either Eb or E1, the stable manifold of E2 separate the basins of attraction of Eb and E1.

    In this section, we study the dynamics of model (1.3) with τ>0.

    By linearizing (1.3) at E0=(0,0), we find two eigenvalues f(0)b>0 and μ<0. Hence, E0 is a saddle point. The characteristic equation associated with the linearization of model (1.3) at Eb is

    (λf(˜y)+b)(λ+q˜y+μeδτh(˜y)eλτ)=0.

    One eigenvalue is λ1=f(˜y)b<0. Hence Eb is locally asymptotically stable if all zeros of

    Δb(λ)=λ+q˜y+μeδτh(˜y)eλτ

    have negative real parts. By [13,Lemma 6], there exists at least one positive eigenvalue if q˜y+μ<eδτh(˜y), or equivalently, y1<˜y<y2. On the other hand, if q˜y+μ>eδτh(˜y), that is, either ˜y<y1 or ˜y>y2, then all eigenvalues have negative real parts. For the critical case q˜y+μ=eδτh(˜y); namely, either ˜y=y1 or ˜y=y2, one eigenvalue is 0, and all other eigenvalues have negative real parts. By using the method of normal forms, we obtain that Eb is locally asymptotically stable. To summarize, we have the following result on the local stability of E0 and Eb.

    Theorem 3.1. Consider model (1.3) under the assumptions (H1)-(H3). E0 is a saddle point. If either ˜yy1 or ˜yy2, then Eb is locally asymptotically stable; while Eb is unstable if y1<˜y<y2.

    To establish the global stability of Eb, we first show that the infection is persistent.

    Lemma 3.2. lim infty(t)>0 for any solution of (1.3) with the initial condition in X1={(ϕ,ψ)X:ϕ(0)>0,ψ(0)0}.

    Proof. We claim that lim supty(t)>0. If not, then y(t)0 as t. It then follows from (1.3) that z(t)0 as t. This contradicts to the fact that E0 is unstable. Thus, y(t) is weakly persistent. By using [14,Theorem 3.4.6] and Proposition 2.1, there exists a global attractor for the solution semiflow of (1.3). This, together with [15,Theorem 2.3], implies that y(t) is actually strongly persistent; that is, lim infty(t)>0.

    Our next result is concerned with global stability of Eb.

    Theorem 3.3. Assume that (H1)-(H3) hold. Then Eb of model (1.3) is globally asymptotically stable in the region X1 provided that ττb, where τb:=1δlnh(0)˜yq˜y+μ. Moreover, the condition ττb implies that ˜yy1.

    Proof. By Proposition 2.1 and Theorem 3.1, it suffices to show that Eb is globally attractive in X1Γ, which is a positively invariant set of (1.3). Let (y(t),z(t)) be a solution of (1.3) with the initial condition in X1, it follows from Lemma 3.2 that y(t)>0 for t>0. Motivated by the earlier work in [16], we construct a Lyapunov functional L:X1ΓR by

    L(yt,zt)=yt(0)˜yln(yt(0))+czt(0)+ceδτ0τh(yt(θ))zt(θ)dθ,

    where c is a constant to be determined later. Calculating the time derivative of L along the solution, and using b=f(˜y)/˜y, h(y)h(0)y for all y0, we obtain

    dLdt=f(y)y(y˜y)b(y˜y)+(p˜ycμ)z+ceδτh(y)z(cq+p)yz,(f(y)yf(˜y)˜y)(y˜y)+(p˜ycμ)z+(ceδτh(0)cqp)yz.

    Assumption (H1) implies that (f(y)yf(˜y)˜y)(y˜y)0 for all yˉy. Since h(0)eδτq>0, the condition ττb is the same as p˜yμpeδτh(0)q. We may choose a constant c[p˜yμ,peδτh(0)q] such that dL/dt0. Moreover, dL/dt=0 if and only if y(t)=˜y and z(t)=0. Thus the maximal compact invariant set in {dL/dt=0} is the singleton {Eb}. By the LaSalle invariance principle [11], Eb is globally attractive in X1Γ. Since X1Γ is absorbing in X1, we conclude that Eb is globally attractive in X1. By Theorem 3.2, Eb is globally asymptotically stable in X1 if ττb.

    If ττb, then eδτh(0)˜yq˜yμ0=eδτh(y1)qy1μeδτh(0)y1qy1μ. Since eδτh(0)>q, we obtain ˜yy1.

    Here, we remark that if y1=y2<˜y, then two positive equilibria E1 and E2 coincide, and Eb is locally asymptotically stable.

    In this subsection, we investigate the stability of E1 and identify parameter regions in which the time delay can destabilize E1, lead to Hopf bifurcation and induce sustained oscillations.

    The characteristic equation associated with the linearization of system (1.3) at equilibrium Ei (i=1,2) is

    Δi(λ)=λ2+a1,iλ+a0,i+(b1,iλ+b0,i)eλτ=0, (3.1)

    where

    a0,i=(f(yi)yif(yi))h(yi)eδτpqyizi,  a1,i=f(yi)yif(yi)+h(yi)eδτ>0,b0,i=(f(yi)f(yi)yi)h(yi)eδτ+pyizih(yi)eδτ,  b1,i=h(yi)eδτ<0.

    We first investigate the stability of E1, regarding τ as the bifurcation parameter. In the proof of Theorem 2.4, we have demonstrated that, when τ=0, all eigenvalues of (3.1) with i=1 lie to the left of the imaginary axis and E1 is locally asymptotically stable. As τ increases, E1 may lose its stability only when some eigenvalues cross the imaginary axis to the right. In view of

    a0,1+b0,1=py1z1(h(y1)eδτq)>0,

    0 is not an eigenvalue for any τ0. For simplicity, we drop the subscript 1 in the following arguments. Substituting λ=iω(ω>0) into (3.1) and separating the real and imaginary parts, we obtain

    ω2a0=b0cos(ωτ)+b1ωsin(ωτ),a1ω=b1ωcos(ωτ)b0sin(ωτ).

    Squaring and adding both the above equations lead to

    G(ω,τ):=ω4+c1(τ)ω2+c0(τ)=0, (3.2)

    where

    c1(τ)=(f(y1)y1f(y1))2+2pqy1z1>0,  c0(τ)=a20b20.

    Then iω(ω>0) is an imaginary root of (3.1) if and only if G(ω,τ)=0. Since c1(τ)>0 for all 0ττmax. we know that G(ω,τ)=0 has exactly one positive root if and only if c0(τ)<0, or equivalently, a0<b0.

    If a0<b0, the implicit function theorem implies that there exists a unique C1 function ω=ω(τ) such that G(ω(τ),τ)=0 for 0<ττmax. For iω(τ) to be a root of (3.1), ω(τ) needs to satisfy the system

    sin(ω(τ)τ)=b1ω(τ)3+(a1b0a0b1)ω(τ)b21ω(τ)2+b20:=g1(τ),cos(ω(τ)τ)=(b0a1b1)ω(τ)2a0b0b21ω(τ)2+b20:=g2(τ). (3.3)

    Set

    I={τ:0ττmax  satisfies a0<b0}. (3.4)

    For τI, let θ(τ) be the unique solution of sinθ=g1 and cosθ=g2 in (0,2π]; that is,

    θ(τ)={arccosg2(τ),   if  ω(τ)2a0a1b0b1,2πarccosg2(τ),   if  ω(τ)2>a0a1b0b1.

    Following Beretta and Kuang [17], we define

    Sn(τ)=τθ(τ)+2nπω(τ)   for  τI  with  nN. (3.5)

    One can check that ±iω(τ) are a pair of purely imaginary eigenvalues of (3.1) if and only if τ is a zero of function Sn(τ) for some nN. From [17,Theorem 2.2], we have

    Sign(dReλ(τ)dτ|τ=τ)=Sign(Gω(ω(τ),τ))SignSn(τ).

    Note that Gω(ω(τ),τ)>0, thus we have the following result concerning the transversality condition.

    Lemma 3.4. Assume that Sn(τ) has a positive root τI for some nN, then a pair of simple purely imaginary roots ±iω(τ) of (3.1) exist at τ=τ, and

    Sign(dReλ(τ)dτ|τ=τ)=SignSn(τ).

    Moreover, this pair of simple purely imaginary roots ±iω(τ) cross the imaginary axis from left to right at τ=τ if Sn(τ)>0, and from right to left if Sn(τ)<0.

    It is easily seen that Sn(τ)>Sn+1(τ) for all τI,nN. When τ=0, the asymptotic stability of E1 implies that S0(0)<0. Denote

    ˆτ=supI:=sup{τ:0ττmax  satisfies a0<b0}. (3.6)

    If ˆτ<τmax, then c0(ˆτ)=0, and thus ω(τ)0 as τˆτ. This, together with (3.3), yields limτˆτsinθ(τ)=0 and limτˆτcosθ(τ)=1. Therefore, limτˆτθ(τ)=π, limτˆτSn(τ)=.

    If supτIS0(τ)<0, Sn(τ) has no zeros in I for all nN, which excludes the existence of purely imaginary eigenvalues and thus implies that E1 is locally asymptotically stable for all τ[0,τmax].

    If supτIS0(τ)=0, S0(τ) has a double zero in I, denoted by τ, and S0(τ)=0. This, together with Lemma 3.4, implies that the transversality condition at τ is not satisfied and all eigenvalues remain to the left of the pure imaginary axis. Thus, E1 is still locally asymptotically stable for τ[0,τmax].

    If supτIS0(τ)>0, it then follows from S0(0)<0 and limτˆτSn(τ)= that S0(τ) has at least two zeros in I. For simplicity, we assume that

    (H4) supτIS0(τ)>0 and all zeros of Sn(τ) with odd multiplicity are simple.

    Under assumption (H4), each Sn(τ) has even number of simple zeros. Now, we collect all simple zeros of Sn(τ) with n0 and list them in increasing order: 0<τ0<τ1<<τ2K1<ˆτ (KN+). For each integer 0iK1, we apply Lemma 3.4 to obtain Sm(τi)>0 and Sm(τ2Ki1)<0 for some mN. Therefore, the pair of simple conjugate purely imaginary eigenvalues ±iω(τi) crosses the imaginary axis from left to right, and the pair of simple conjugate purely imaginary eigenvalues ±iω(τ2Ki1) crosses the imaginary axis from right to left. Thus, system (1.3) undergoes a Hopf bifurcation at E1 when τ=τj (0j2K1). Moreover, E1 is asymptotically stable for τ[0,τ0)(τ2K1,τmax], and unstable for τ(τ0,τ2K1).

    For each k=0,,2K1, there exists nk such that τk is a simple zero of Snk(τ). Let Tk be the period of periodic solution bifurcated at τk. Applying the Hopf bifurcation theorem for delay differential equations [11,18], we have

    Tk=2πω(τk)=2πτkθ(τk)+2nkπ.

    This, together with θ(τk)(0,2π], implies that Tk>τk if nk=0, and {\tau_k\over n_k+1}\le T_k < {\tau_k\over n_k} if n_k > 0 . To summarize, we have the following results on stability of E_1^* and Hopf bifurcation.

    Theorem 3.5. Assume that ( \textbf{H}_{1} )-( \textbf{H}_{3} ) hold, y_1^* < \tilde y and \hat\tau < \tau_{max} . Let I, S_n(\tau), \tau_{max}, \hat\tau be defined in (3.4), (3.5), ( \textbf{H}_{3} ) and (3.6), respectively.

    (i) If either I = \emptyset or \sup\limits_{\tau\in I} S_0(\tau)\le 0 , then E_1^* is locally asymptotically stable for all \tau\in[0, \tau_{max}] .

    (ii) If ( \textbf{H}_{4} ) holds, then there exist exactly 2K local Hopf bifurcation values, namely, 0 < \tau_0 < \tau_1 < \cdots < \tau_{2K-1} < \hat\tau such that model (3.1) undergoes a Hopf bifurcation at E_1^* when \tau = \tau_j for 0\le j\le 2K-1 . E_1^* is locally asymptotically stable for \tau\in[0, \tau_0)\cup(\tau_{2K-1}, \tau_{max}] , and unstable for \tau\in(\tau_0, \tau_{2K-1}) . Moreover, the period T_k of periodic solution bifurcated at \tau_k satisfies T_k > \tau_k if \tau_k is a simple zero of S_0(\tau) , and {\tau_k\over n_k+1}\le T_k < {\tau_k\over n_k} if \tau_k is a simple zero of S_{n_k}(\tau) with n_k > 0 .

    Theorem 3.6. Assume that ( \textbf{H}_{1} )-( \textbf{H}_{3} ) are satisfied. If y_1^* < y_2^* < \tilde{y} holds, then E_2^* is unstable for all \tau\geq 0 . Moreover, if a_{0, 2}\le b_{0, 2} , then the characteristic equation (3.1) with i = 2 has no purely imaginary eigenvalues.

    Proof. Note that in (3.1) with i = 2 , \Delta_2(0) = p y_{2}^{*}z_{2}^{*}(h'(y_{2}^{*})e^{-\delta\tau}-q) < 0 . For any \tau\geq 0 , we know that \Delta_2(\infty) > 0 , thus as long as E_2^* exists, the associated characteristic equation must have at least one positive real root and hence E_2^* is always unstable for \tau\geq 0 .

    Using a similar argument in the study of root distribution for (3.1) at E_1^* , we note that c_{1, 2}(\tau) > 0, c_{0, 2}(\tau) = a_{0, 2}^2-b_{0, 2}^2 and a_{0, 2}+b_{0, 2} = p y_{2}^{*}z_{2}^{*}(h'(y_{2}^{*})e^{-\delta\tau}-q) < 0 . Thus, G(\omega, \tau) in (3.2) has no positive zeros if a_{0, 2}\le b_{0, 2} . Therefore, (3.1) with i = 2 has no purely imaginary eigenvalues if a_{0, 2}\le b_{0, 2} .

    Since E_2^* is unstable (and thus biologically irrelevant) whenever it exists, we name E_1^* , instead of E_2^* , the immune control equilibrium (ICE).

    Note that Theorem 3.5 states that if y_1^* < \tilde y and ( \textbf{H}_{4} ) holds, then periodic solutions can bifurcate from E_1^* when \tau is near the local Hopf bifurcation values \tau_k, k = 0, 1, 2, \ldots, 2K-1. For integer n\ge0 , we denote

    \begin{equation} J: = \{\tau_0, \tau_1,\cdots,\tau_{2K-1}\}, \ \ J_n: = \{\tau\in J: S_n(\tau) = 0\}, \ \ J^0 = J\setminus J_0. \end{equation} (3.7)

    According to Theorem 3.5, J_n has even number of bifurcation values, and the period of any periodic solution bifurcated near a bifurcation point in J^0 is bounded from both below and above, while the periodic solution bifurcated near a bifurcation point in J_0 has a lower bound for its period. As we shall see later in the numerical simulation, it seems impossible to find an upper bound for such period. In what follows, we will restrict our investigation on the set J^0 and assume that J^0\neq\emptyset . Especially, we will discuss the global continuation of periodic solutions bifurcated from the point (E_1^*, \tau_k) with \tau_k\in J^0 as the bifurcation parameter \tau varies. We shall use the global Hopf bifurcation theorem for delay differential equations [19] and show that model (1.3) admits periodic solutions globally for all \tau\in(\underline\tau, \bar\tau) , where \underline\tau = \min J^0 and \bar\tau = \max J^0 .

    Let x(t) = (y(\tau t), z(\tau t))^T , model (1.3) can be rewritten as a general functional differential equation:

    \begin{equation} x'(t) = F(x_{t},\tau,T), \; \; (t,\tau,T)\in\mathbb{R}_{+}\times I\times\mathbb{R}_{+}, \end{equation} (3.8)

    where x_{t}(\theta) = x(t+\theta) for \theta\in[-1, 0] , and x_{t}\in \mathcal{X}: = C([-1, 0], \mathbb{R}_{+}^2) , and

    \begin{equation} F(x_t,\tau,T) = \left( \begin{array}{ll} \tau f(x_{1t}(0))-p\tau x_{1t}(0)x_{2t}(0)-b\tau x_{1t}(0)\\ \tau e^{-\delta\tau} h(x_{1t}(-1))x_{2t}(-1)-q\tau x_{1t}(0)x_{2t}(0)-\mu\tau x_{2t}(0) \end{array} \right) \end{equation} (3.9)

    with x_t = (x_{1t}, x_{2t})\in\mathcal{X} . Identifying the subspace of \mathcal{X} consisting of all constant functions with \mathbb{R}_{+}^2 , we obtain a restricted function

    \begin{equation*} \widetilde{F}(x,\tau,T): = F\mid_{\mathbb{R}_{+}^2\times I\times\mathbb{R}_{+}} = \left( \begin{array}{ll} \tau f(x_1)-p\tau x_1 x_2-b\tau x_1\\ \tau e^{-\delta\tau} h(x_1)x_2-q\tau x_1 x_2-\mu\tau x_2 \end{array} \right). \end{equation*}

    Obviously, \widetilde{F} is twice continuously differentiable, i.e., the condition (\textbf{A1}) in [19] holds. We denote the set of stationary solutions of system (3.8) by {\bf{E}}(F) = \{(\widetilde{x}, \widetilde{\tau}, \widetilde{T})\in\mathbb{R}_{+}^2\times I\times\mathbb{R}_{+}: \; \widetilde{F}(\widetilde{x}, \widetilde{\tau}, \widetilde{T}) = 0 \}. It follows from Proposition 2.2 that (ⅰ) if y_1^* < \tilde{y} < y_2^* holds,

    {\bf{E}}(F) = \{(E_0,\tau,T), (E_b,\tau,T), (E_1^*,\tau,T); \; (\tau,T)\in I\times\mathbb{R}_{+}\};

    (ⅱ) if y_1^* < y_2^* < \tilde{y} holds,

    {\bf{E}}(F) = \{(E_0,\tau,T), (E_b,\tau,T), (E_1^*,\tau,T), (E_2^*,\tau,T); \; (\tau,T)\in I\times\mathbb{R}_{+}\}.

    For any stationary solution (\widetilde{x}, \tau, T) , the characteristic function is

    \Delta_{(\widetilde{x},\tau,T)}(\lambda) = \lambda \text{Id}-DF(\widetilde{x},\tau,T)(e^{\lambda\cdot}\text{Id}).

    Note that if either y_1^* < \tilde{y} < y_2^* or y_1^* < y_2^* < \tilde{y} holds, 0 is not a eigenvalue of any stationary solution of (3.8) and hence the condition (A2) in [19] holds. It can be checked easily from (3.9) that the smoothness condition (A3) in [19] is satisfied.

    Theorem 3.5 implies that if y_1^* < \tilde y and ( \textbf{H}_{4} ) holds, then for each integer 0\le k\le 2K-1 the stationary solution (E_1^*, \tau_k, 2\pi/(\omega_k\tau_k)) is an isolated center of (3.8), where \omega_k = \omega(\tau_k) is the unique positive zero of G(\omega, \tau) in (3.2), and there is only one purely imaginary characteristic value of the form im (2\pi/\widetilde{T}) with m = 1 and \widetilde{T} = 2\pi/(\omega_k\tau_k) . Moreover, if follows from Lemma 3.4 that the crossing number \zeta_{1}(E_1^*, \tau_k, 2\pi/(\omega_k\tau_k)) at each of these centers is

    \begin{equation} \zeta_{1}(E_1^*, \tau_k,\frac{2\pi}{\omega_k\tau_k}) = -\text{Sign}\left({d Re\lambda(\tau)\over d\tau}\big|_{\tau = \tau_k}\right) = \begin{cases} -1,&\ 0\le k\le K-1,\\ 1,&\ K\le k\le 2K-1. \end{cases} \end{equation} (3.10)

    Thus the condition (A4) in [19] holds. We then define a closed subset \Sigma(F) of \mathcal{X}\times I \times \mathbb{R}_+ by

    \Sigma(F) = \mathcal{C}l\{(x,\tau,T)\in \mathcal{X}\times I \times \mathbb{R}_+: x \mbox{ is a nontrivial T-periodic soltuion of (3.8)}\},

    and for each integer 0\le k\le 2K-1 , let {\bf{C}}(E_1^*, \tau_k, 2\pi/(\omega_k\tau_k)) denote the connected component of (E_1^*, \tau_k, 2\pi/(\omega_k\tau_k)) in \Sigma(F) . By Theorem 3.5, {\bf{C}}(E_1^*, \tau_k, 2\pi/(\omega_k\tau_k)) is a nonempty subset of \Sigma(F) .

    To find the interval of \tau in which periodic solutions exist, we shall further investigate the properties of periodic solutions of (3.8).

    Lemma 3.7. All nonconstant and nonnegative periodic solutions of (3.8) are uniformly bounded. Actually, we have 0 < y(t), z(t)\le M for any t\in\mathbb{R} , where M = \max\{\bar y, h'(0)M_f/(p\gamma)\} .

    Proof. Since y(t) is nonconstant and nonnegative, there exists t_0 > 0 such that y(t_0) > 0 . Integrating the equation for y'(t) gives

    e^{\int_{t_0}^t[pz(s)+b]ds}y(t) = y(t_0)+\int_{t_0}^t e^{\int_{t_0}^r[pz(s)+b]ds}f(y(r))dr \gt 0,\; t \gt t_0.

    Hence, y(t) > 0 for all t > t_0 . Since y is periodic, we have y(t) > 0 for all t > 0 . Similarly, if z(t_0) > 0 for some t_0 > 0 , we integrate the equation for z'(t) to obtain

    e^{\int_{t_0}^t[qy(s)+\mu]ds}z(t) = z(t_0)+\int_{t_0}^t e^{\int_{t_0}^r[qy(s)+\mu]ds}e^{-\delta\tau}h(y(r-\tau))z(r-\tau)dr \gt 0

    for all t > t_0 . This together with periodicity of z implies that z(t) > 0 for all t > 0 .

    By Proposition 2.1, we have \limsup\limits_{t\rightarrow\infty}y(t)\le\bar y . We claim y(t)\le M for all t\ge 0 . Otherwise, if y(t_1) > M for some t_1\ge0 , then

    \lim\limits_{n\rightarrow\infty}y(t_1+nT) = y(t_1) \gt M,

    where T is the period of the nonconstant and nonnegative periodic solutions. This contradicts the fact that y(t) is eventually bounded above by M as t\rightarrow\infty . Thus, M is the upper bound of y(t) . Similarly, from Proposition 2.1, we have \limsup\limits_{t\rightarrow\infty}z(t)\le h'(0)M_f/(p\gamma)\le M . This implies that z(t) is uniformly bounded by M .

    Lemma 3.8. System (3.8) has no nontrivial periodic solution of period 1.

    Proof. Note that any nontrivial 1 -periodic solution x(t) of system (3.8) is also a nontrivial periodic solution of model (2.4), which does not admit any nontrivial periodic solution via Lemma 2.3.

    We are now in the position to state our result concerning the properties of the global Hopf branches.

    Theorem 3.9. Assume that ( \textbf{H}_{1} )-( \textbf{H}_{4} ), either y_1^* < \tilde{y} < y_2^* or y_1^* < y_2^* < \tilde{y} holds, \hat\tau < \tau_{max}, J^0\neq\emptyset and a_{0, 2}\le b_{0, 2} . Then for system (3.8), we have the following results:

    (i) The connected component {\bf{C}}(E_1^*, \tau_k, 2\pi/(\omega_k\tau_k)) is bounded for all \tau_k\in J^0 .

    (ii) Let n\ge1 . If J_n has 2k_n bifurcation values, ordered as \tau_{n, 1} < \cdots < \tau_{n, 2k_n} , then for each i = 1, \cdots, k_n , there exists j > k_n such that \tau_{n, i} and \tau_{n, j} are connected by a global Hopf branch. Similarly, for each j = k_n+1, \cdots, 2k_n , there exists i\le k_n such that \tau_{n, i} and \tau_{n, j} are connected by a global Hopf branch. Especially, if k_n = 1 , then the two bifurcation values of J_n are connected.

    (iii) For each \tau\in(\min J_n, \max J_n) with integer n\ge1 , system (3.8) has at least one periodic solution with period in (1/(n+1), 1/n) .

    Proof. Lemma 3.7 implies that the projection of {\bf{C}}(E_1^*, \tau_k, 2\pi/(\omega_k\tau_k)) onto \mathcal{X} is bounded. Note that system (3.8) has no periodic solutions of period 1 , and thus no periodic solutions of period 1/n_{k} or 1/(n_{k}+1) for any positive integer n_{k} . It follows from Lemma 3.8 that the period T of a periodic solution on the connected component {\bf{C}}(E_1^*, \tau_k, 2\pi/(\omega_k\tau_k)) satisfies

    {1\over n_{k}+1} \lt T \lt {1\over n_{k}}

    with \tau_{k}\in J_{n_k} and n_k\ge1 . Hence, the projection of {\bf{C}}(E_1^*, \tau_k, 2\pi/(\omega_k\tau_k)) onto the T -space is bounded for \tau_k\in J^0 . Note that \tau\in I and I is a bounded interval. Therefore, {\bf{C}}(E_1^*, \tau_k, 2\pi/(\omega_k\tau_k)) is bounded in \mathcal{X}\times I\times \mathbb{R}_+ for any \tau_k\in J^0 . This proves (ⅰ).

    From the proof of Theorems 3.1 and 3.6, we know that the stationary solutions (E_0, \tau, T) , (E_b, \tau, T) and (E_2^*, \tau, T) cannot be a center for any \tau and T if either y_1^* < \tilde{y} < y_2^* or y_1^* < y_2^* < \tilde{y} holds and a_{0, 2}\le b_{0, 2} . It then follows from the global bifurcation theorem [19,Theorem 3.3] that \mathcal{E}: = {\bf{C}}(E_1^*, \tau_k, 2\pi/(\omega_k\tau_k))\cap {\bf{E}}(F) is finite and

    \sum\limits_{(\widetilde{x},\tau,T)\in \mathcal{E}}\zeta_1(\widetilde{x},\tau,T) = 0.

    If J_n has 2k_n bifurcation values, ordered as \tau_{n, 1} < \cdots < \tau_{n, 2k_n} , then by (3.10), \zeta_1(\widetilde{x}, \tau_{n, i}, T) = -1 for each i = 1, \cdots, k_n , and \zeta_1(\widetilde{x}, \tau_{n, j}, T) = 1 for each j = k_n+1, \cdots, 2k_n . This together with the above equation implies (ⅱ).

    Note that \tau_{n, 1} = \min J_n and \tau_{n, 2k_n} = \max J_n . The \tau projection of the global bifurcation branch bifurcated from \tau_{n, 1} includes (\tau_{n, 1}, \tau_{n, k_n+1}) , while the \tau projection of the global bifurcation branch bifurcated from \tau_{n, 2k_n} includes (\tau_{n, k_n}, \tau_{n, 2k_n}) . Thus, for any \tau\in(\tau_{n, 1}, \tau_{n, 2k}) , system (3.8) has at least one periodic solution. According to the proof of (ⅰ), the period of this periodic solution lies in (1/(n+1), 1/n) . Thus, (ⅲ) is proved.

    In this section, we present some numerical simulations to demonstrate our theoretical results. We will use the Matlab package DDE-BIFTOOL developed by Engelborghs et al. [20,21] to sketch the global Hopf branches {\bf{C}}(E_1^*, \tau_k, 2\pi/(\omega_k\tau_k)) for 0\le k\le 2K-1 . Following the study of ODE model in [6], we choose

    f(y) = ry(1-{y\over K}), \ \ h(y) = {cy\over 1+dy},

    and set the parameter values as follows:

    \begin{alignat} {3} r& = 6 \; day^{-1}, & \quad K& = 3 \; virus\; mm^{-3}, &\quad p& = 1 \; mm^{3}\ cells^{-1}\ day^{-1}, \\ b& = 3 \; day^{-1}, & \quad c& = 4 \; mm^3\ virus^{-1}day^{-1}, &\quad d& = 0.5 \; mm^3\ virus^{-1}, \\ \delta & = 0.03 \; day^{-1}, & \quad q& = 1 \; mm^3\ virus^{-1}day^{-1}, &\quad \mu& = 0.5 \; day^{-1}. \end{alignat} (4.1)

    Here, the units for viral loads y and immune cell density z are virus per \text{mm}^3 and cells per \text{mm}^3 , respectively. The delay \tau has the unit in days. It is easy to calculate

    \hat\tau = 17.1939 \lt \tau_{max} = 19.1788 \lt \tau_b = 36.6204.

    From Proposition 2.2, there are exactly two equilibria E_{0} = (0, 0) and E_{b} = (1.5, 0) when \tau > \tau_{max} , and, if \tau\leq\tau_{max} , there exists a third equilibrium E_1^* = (y_1^*, z_1^*) , where y_1^*, z_1^* depend on \tau . This is the unique positive equilibrium when \tau\in[0, \tau_e] and there exists another positive equilibrium E_2^* = (y_2^*, z_2^*) when \tau\in(\tau_e, \tau_{max}) ; see Figure 1. Here, \tau_e is the critical value when y_2^* = \tilde y ; that is,

    \tau_e = {1\over\delta}\ln{c\tilde y\over(q\tilde y+\mu)(1+d\tilde y)} = 17.9666.
    Figure 1.  The graphs of y_1^*(\tau) , y_2^*(\tau) and \tilde y(\tau) , which characterize the existence of positive equilibria.

    At \tau = \tau_{max} , the two positive equilibria coincide: y_1^* = y_2^* .

    Theorem 3.1 implies that E_0 is a saddle point, and E_b is locally asymptotically stable for \tau\ge\tau_e , and unstable for \tau < \tau_e , Moreover, by Theorem 3.3, E_b is globally asymptotically stable in the region X_1 provided that \tau\geq\tau_b . It can be verified that the conditions of case (ⅱ) in Theorem 3.5 are satisfied only if 0\leq\tau < \hat\tau . By Theorem 3.5, we know that there are exactly 4 local Hopf bifurcation values with K = 2 , namely,

    \tau_0\approx0.2994 \lt \tau_1\approx8.9274 \lt \tau_2\approx13.9287 \lt \tau_3\approx17.0322,

    as shown in Figure 2. Correspondingly,

    \omega_0\approx0.9537 \gt \omega_1\approx0.7540 \gt \omega_2\approx0.4994 \gt \omega_3\approx0.1069.
    Figure 2.  The graphs of S_{n}(\tau) (n = 0, 1, 2) . This gives the solution \tau_{j}, j = 0, 1, 2, 3 .

    Clearly, E_1^* is stable for \tau\in[0, \tau_0)\cup(\tau_3, \tau_{max}] and unstable for \tau\in(\tau_0, \tau_3) , which implies that \tau_0 and \tau_3 are stability switches. Remark that \tau_3 < \tau_e .

    In Figure 3, we plot the four global Hopf branches {\bf{C}}(E_1^*, \tau_k, 2\pi/(\omega_k\tau_k)) with 0\le k\le 3 . It is observed that the branches {\bf{C}}(E_1^*, \tau_1, 2\pi/(\omega_1\tau_1)) and {\bf{C}}(E_1^*, \tau_2, 2\pi/(\omega_2\tau_2)) are bounded and connected, which agrees with Theorem 3.9. The periodic solutions on these global Hopf branches are plotted on the phase plane of y and z ; see Figures 46. As we know, the stability of the periodic solution is completely determined by the associated principal Floquet multiplier: if the principal Floquet multiplier is larger than 1 , then the corresponding periodic solution is unstable, otherwise, the bifurcated periodic solution is stable [11].

    Figure 3.  All global Hopf branches of model (1.3) with parameter values given in (4.1).
    Figure 4.  The periodic solutions on the global Hopf branch {\bf{C}}(E_1^*, \tau_0, 2\pi/(\omega_0\tau_0)) .
    Figure 5.  The periodic solutions on the global Hopf branch {\bf{C}}(E_1^*,\tau_1, 2\pi/(\omega_1\tau_1)) = {\bf{C}}(E_1^*,\tau_2, 2\pi/(\omega_2\tau_2)) .
    Figure 6.  The periodic solutions on the global Hopf branch {\bf{C}}(E_1^*, \tau_3, 2\pi/(\omega_3\tau_3)) .

    By using DDE-BIFTOOL, we can further calculate the associated principal Floquet multipliers; see Figure 7. Note that the periodic solution on the first branch {\bf{C}}(E_1^*, \tau_0, 2\pi/(\omega_0\tau_0)) is stable for small \tau , and becomes unstable as \tau increases; the periodic solution on the second branch {\bf{C}}(E_1^*, \tau_1, 2\pi/(\omega_1\tau_1)) –or equivalently, the third branch {\bf{C}}(E_1^*, \tau_2, 2\pi/(\omega_2\tau_2)) –is always unstable; the periodic solution on the fourth branch {\bf{C}}(E_1^*, \tau_3, 2\pi/(\omega_3\tau_3)) is stable for \tau near \tau_3 , and becomes unstable as \tau decreases. Moreover, for any \tau\in(\tau_1, \tau_2) , model (1.3) has at least one periodic solution on the second (or third) global Hopf branch with period in the interval (\tau/2, \tau) , as shown in Figure 8.

    Figure 7.  The principal Floquet multipliers of periodic solutions on all Hopf branches of model (1.3).
    Figure 8.  The periods of periodic solutions on the global Hopf branches of model (1.3).

    It is interesting to note from Figure 3 that the first and the fourth Hopf branches do not connect. This phenomenon is quite different from that in the existing literature. We have proved that both the periodic solutions and delays on any global Hopf branch are uniformly bounded. According to the global Hopf bifurcation theory [19], either the first and the fourth Hopf branches connect, or these two branches have unbounded periods. From Figure 8, we observe that the periods of periodic oscillations on the fourth branch {\bf{C}}(E_1^*, \tau_3, 2\pi/(\omega_3\tau_3)) seem to be unbounded. To further explore this interesting phenomenon, we numerically sketch in Figure 9 a bifurcation diagram for model (1.3). It is suggested that model (1.3) can have chaotic solutions and thus allow aperiodic oscillations when \tau lies in the interval where the first and fourth Hopf branches disappear. Figure 9 also confirms that \tau_0 and \tau_3 are stability switches: (ⅰ) when \tau crosses \tau_0 from left to right, the positive equilibrium E_1^* becomes unstable, while a stable periodic solution bifurcates; and (ⅱ) when \tau crosses \tau_3 from left to right, E_1^* regains its stability while there exists a stable periodic solution when \tau is close to \tau_3 from the left. It should be mentioned that the bifurcation diagram has a jump at \tau_{max} , which is due to the fact that there exists no positive equilibrium when \tau > \tau_{max} but the positive equilibrium at \tau = \tau_{max} is nonzero.

    Figure 9.  Bifurcation diagram of (1.3) with \tau as the bifurcation parameter.

    Considering the fact the immune response process involves a sequence of events such as antigenic activation, selection, and proliferation of the immune cells and the mortality of immune cells during the time lag, we derive an immunosuppressive infection model from a stage structured population model. We find that the time lag \tau plays a key role in the infection and immune response dynamics. As \tau increases, the model dynamics shifts through four possible outcomes: (ⅰ) virus persists without immune response when \tau > \tau_{max} ; (ⅱ) the outcome is initial condition dependent and delay dependent, either the immunity can be completely destroyed, or the immune mediated control can be established when \tau\in(\tau_e, \tau_{max}] ; (ⅲ) immune response controls the virus when \tau\in[0, \tau_0)\cup(\tau_{2K-1}, \tau_e) ; and (ⅳ) sustained oscillation appears or chaotic dynamics occurs when \tau\in(\tau_0, \tau_{2K-1}) . Our theoretical results coincide with the phenomenon of oscillatory viral loads and immune cells observed in clinical data [9].

    From numerical simulation, we also observe an unbounded global Hopf bifurcation branch with bounded periodic solutions, bounded delays, but unbounded periods. To the best of our knowledge, this phenomenon is new, because, in the existing literature, either the global Hopf branch is bounded or it is unbounded with unbounded delays. A detailed theoretical study of this interesting numerical phenomenon should enrich the global Hopf bifurcation theory for delay differential equations, and we leave it as a future work.

    We are very grateful to the anonymous referees for their valuable comments and suggestions, which help to improve the presentation of this manuscript. H. Shu was partially supported by the National Natural Science Foundation of China (No.11971285, No.11601392), and the Fundamental Research Funds for the Central Universities.

    The authors declare there is no conflict of interest.



    [1] Y. Ding, H. M. Wang, P. C. Shi, et al., Trusted cloud service, Chin. J. Comput., 38 (2015), 133–149.
    [2] M. Ali, S. U. Khan and A. V. Vasilakos, Security in cloud computing: Opportunities and challenges, Inform. Sciences., 305 (2015), 357–383.
    [3] Y. Q. Zhang, X. F. Wang, X. F. Liu, et al., Survey on cloud computing security, J. Software, 27 (2016), 1328−1348.
    [4] J. Wilhelm and T. C. Chiueh, A forced sampled execution approach to kernel rootkit identification, In: International Workshop on Recent Advances in Intrusion Detection; 2007 Sept 5–7; Gold Goast, Australia. Berlin: Springer; 2007: 219–235.
    [5] N. Zhang, R. Zhang, K. Sun, et al., Memory Forensic Challenges Under Misused Architectural Features, IEEE T. Inf. Foren. Sec., 13 (2018), 2345–2358.
    [6] A. Cohen, N. Nissim, L. Rokach, et al., SFEM: Structural feature extraction methodology for the detection of malicious office documents using machine learning methods, Expert Syst. Appl., 63 (2016), 324–343.
    [7] N. Nissim, R. Moskovitch, O. BarAd, et al., ALDROID:Efficient update of Android anti-virus software using designated active learning methods, Knowl. Inf. Syst., 49 (2016), 795–833.
    [8] N. Nissim, A. Cohen, C. Glezer, et al., Detection of malicious PDF files and directions for enhancements: A state-of-the art survey, Comput. Secur., 48 (2015), 246–266.
    [9] G. Hoglund and J. Butler, Rootkits: subverting the Windows kernel, Addison-Wesley Professional, New Jersey, 2006.
    [10] A. Case and G. G. Richard III, Advancing Mac OS X rootkit detection, Digit. Invest., 14 (2015), S25–S33.
    [11] H. Yang, J. Zhuge, H. Liu, et al., A tool for volatile memory acquisition from Android devices, In: IFIP International Conference on Digital Forensics; 2016 Jan 4-6; New Delhi, India. Cham: Springer; 2016: 365–378.
    [12] A. Kumara and C. D. Jaidhar, Execution time measurement of virtual machine volatile artifacts analyzers, In: 2015 IEEE 21st International Conference on Parallel and Distributed Systems (ICPADS);2015 Dec14-17; Melbourne, VIC, Australia. IEEE; 2015: 314–319.
    [13] Q. Hua and Y. Zhang, Detecting Malware and Rootkit via Memory Forensics, In:2015 International Conference on Computer Science and Mechanical Automation (CSMA); 2015 Oct 23–25; Hangzhou, China. IEEE; 2015: 92–96.
    [14] C. W. Tien, J. W. Liao, S. C. Chang, et al., Memory forensics using virtual machine introspection for Malware analysis, In:2017 IEEE Conference on Dependable and Secure Computing; 2017 Aug 7-10; Taipei, Taiwan. IEEE; 2017: 518–519.
    [15] A. Cohen and N. Nissim, Trusted detection of ransomware in a private cloud using machine learning methods leveraging meta-features from volatile memory, Expert Syst. Appl., 102 (2018), 158–178.
    [16] N. Nissim, Y. Lapidot, A. Cohen, et al., Trusted system-calls analysis methodology aimed at detection of compromised virtual machines using sequential mining, Knowl.-Based Syst., 153 (2018), 147–175.
    [17] A. Kumara and C. D. Jaidhar, Automated multi-level malware detection system based on reconstructed semantic view of executables using machine learning techniques at VMM, Future Gener. Comp. Sy., 79 (2018), 431–446.
    [18] H. Upadhyay, H. A. Gohel, A. Pons, et al., Windows Virtualization Architecture For Cyber Threats Detection. In:2018 1st International Conference on Data Intelligence and Security (ICDIS). 2018 Apr 8-10; South Padre Island, TX, USA.IEEE; 2018: 119–122.
    [19] R. Mosli, R. Li, B. Yuan, et al., Automated malware detection using artifacts in forensic memory images. In:2016 IEEE Symposium on Technologies for Homeland Security (HST). 2016 May 10-11; Waltham, MA, USA. IEEE; 2016: 1–6.
    [20] M. A. Kumara and C. D. Jaidhar, Leveraging virtual machine introspection with memory forensics to detect and characterize unknown malware using machine learning techniques at hypervisor, Digit. Invest., 23 (2017), 99–123.
    [21] J. Bai and J. Wang, Improving malware detection using multi view ensemble learning, Secur. Commun. Netw., 9 (2016), 4227–4241.
    [22] OpenStack. Available from: https://docs.openstack.org/rocky/.
    [23] Volatility. Available from: https://www.volatilityfoundation.org/.
    [24] M. H. Ligh, A. Case, J. Levy, et al., The art of memory forensics: detecting malware and threats in windows, linux, and Mac memory, John Wiley & Sons, New Jersey, 2014.
    [25] Malshare. Available from: http://www.malshare.com
    [26] T. Kim, B. Kang, M. Rho, et al., A Multimodal Deep Learning Method for Android Malware Detection Using Various Features. IEEE T. Inform. Foren. Sec., 14 (2019), 773–788.
    [27] Virustotal. Available from: https://www.virustotal.com/
    [28] M. Hall, E. Frank, G. Holmes, et al., The WEKA data mining software: an update. ACM SIGKDD explorations newsletter, 11 (2009): 10–18.
    [29] Z. Wang, J. Ren, D. Zhang, et al., A deep-learning based feature hybrid framework for spatiotemporal saliency detection inside videos, Neurocomputing, 287 (2018), 68–83.
    [30] J. Zabalza, J. Ren, J. Zheng, et al., Novel segmented stacked autoencoder for effective dimensionality reduction and feature extraction in hyperspectral imaging, Neurocomputing, 185 (2016), 1–10.
    [31] S. Md Noor, J. Ren, S. Marshall, et al., Hyperspectral Image Enhancement and Mixture Deep-Learning Classification of Corneal Epithelium Injuries, Sensors, 17 (2017), 2644.
    [32] J. Ren, D. Wang and J Jiang, Effective recognition of MCCs in mammograms using an improved neural classifier, Eng. Appl. Artif. Intel., 24 (2011), 638–645.
  • This article has been cited by:

    1. Wonjun Lee, Mohammad Nadim, 2020, Kernel-Level Rootkits Features to Train Learning Models Against Namespace Attacks on Containers, 978-1-7281-6550-9, 50, 10.1109/CSCloud-EdgeCom49738.2020.00018
    2. Rongfu Zhou, Jun Lin, Lan Liu, Min Ye, Shunhe Wei, 2020, Chapter 46, 978-3-030-39430-1, 479, 10.1007/978-3-030-39431-8_46
    3. Suresh Kumar S, SudalaiMuthu T, 2022, Kernel Rootkit Secret Detection in Cloud Computing, 978-1-6654-7655-3, 276, 10.1109/ICCST55948.2022.10040354
    4. Maryam Mohammadzad, Jaber Karimpour, Using rootkits hiding techniques to conceal honeypot functionality, 2023, 10848045, 103606, 10.1016/j.jnca.2023.103606
    5. Rahul N. Vaza, Ramesh Prajapati, Dushyantsinh Rathod, Dineshkumar Vaghela, Developing a novel methodology for virtual machine introspection to classify unknown malware functions, 2022, 15, 1936-6442, 793, 10.1007/s12083-021-01281-5
    6. Mohammad Nadim, David Akopian, Wonjun Lee, 2021, A Review on Learning-based Detection Approaches of the Kernel-level Rootkit, 978-1-6654-2714-2, 1, 10.1109/ICEET53442.2021.9659710
    7. Tomer Panker, Nir Nissim, Leveraging malicious behavior traces from volatile memory using machine learning methods for trusted unknown malware detection in Linux cloud environments, 2021, 226, 09507051, 107095, 10.1016/j.knosys.2021.107095
    8. Luxin Zheng, Jian Zhang, Qi Jiang, A New Malware Detection Method Based on VMCADR in Cloud Environments, 2022, 2022, 1939-0122, 1, 10.1155/2022/4208066
    9. Jean-Paul A. Yaacoub, Hassan N. Noura, Ola Salman, Ali Chehab, Advanced digital forensics and anti-digital forensics for IoT systems: Techniques, limitations and recommendations, 2022, 19, 25426605, 100544, 10.1016/j.iot.2022.100544
    10. Tom Landman, Nir Nissim, Deep-Hook: A trusted deep learning-based framework for unknown malware detection and classification in Linux cloud environments, 2021, 144, 08936080, 648, 10.1016/j.neunet.2021.09.019
    11. Faouzi Kamoun, Farkhund Iqbal, Mohamed Amir Esseghir, Thar Baker, 2020, AI and machine learning: A mixed blessing for cybersecurity, 978-1-7281-5628-6, 1, 10.1109/ISNCC49221.2020.9297323
    12. Asad Arfeen, Muhammad Asim Khan, Obad Zafar, Usama Ahsan, Process based volatile memory forensics for ransomware detection, 2022, 34, 1532-0626, 10.1002/cpe.6672
    13. Amir Djenna, Ahmed Bouridane, Saddaf Rubab, Ibrahim Moussa Marou, Artificial Intelligence-Based Malware Detection, Analysis, and Mitigation, 2023, 15, 2073-8994, 677, 10.3390/sym15030677
    14. Hamad Naeem, Shi Dong, Olorunjube James Falana, Farhan Ullah, Development of a Deep Stacked Ensemble With Process Based Volatile Memory Forensics for Platform Independent Malware Detection and Classification, 2023, 09574174, 119952, 10.1016/j.eswa.2023.119952
    15. Mathieu Drolet, Vincent Roberge, Enterprise Malware Detection using Digital Forensic Artifacts and Machine Learning, 2024, 12, 2415-1521, 336, 10.37394/232018.2024.12.33
    16. S Suresh Kumar, T SudalaiMuthu, 2023, Advance Kernel Rootkit Detection: Survey, 979-8-3503-9725-3, 944, 10.1109/ICICCS56967.2023.10142360
    17. 2023, 9781119898870, 1, 10.1002/9781119898900.ch1
    18. Suresh Kumar S, Stephen S, Suhainul Rumysia M, 2024, Rootkit Detection Using Hybrid Machine Learning Models and Deep Learning Model: Implementation, 979-8-3503-8943-2, 1, 10.1109/ACCAI61061.2024.10602165
    19. S Suresh Kumar, S Stephen, M Suhainul Rumysia, 2024, Rootkit Detection using Deep Learning: A Comprehensive Survey, 979-8-3503-5306-8, 365, 10.1109/ICCSP60870.2024.10543963
    20. Tomer Panker, Aviad Cohen, Tom Landman, Chen Bery, Nir Nissim, MinCloud: Trusted and transferable MinHash-based framework for unknown malware detection for Linux cloud environments, 2024, 87, 22142126, 103907, 10.1016/j.jisa.2024.103907
    21. Basirah Noor, Sana Qadir, Machine Learning and Deep Learning Based Model for the Detection of Rootkits Using Memory Analysis, 2023, 13, 2076-3417, 10730, 10.3390/app131910730
    22. R. M. Muthulakshmi, T. P. Anithaashri, C. Nataraj, V. S. N. Talasila, 2024, 3161, 0094-243X, 020339, 10.1063/5.0229420
    23. Nour Moustafa, Nickolaos Koroniotis, Marwa Keshk, Albert Y. Zomaya, Zahir Tari, Explainable Intrusion Detection for Cyber Defences in the Internet of Things: Opportunities and Solutions, 2023, 25, 1553-877X, 1775, 10.1109/COMST.2023.3280465
    24. Eszter Kail, Anna Banati, Rita Fleiner, Amirhosein Mosavi, Csaba Mako, 2023, Machine Learning Methods for Cybersecurity: Review and Bibliometric Analysis, 979-8-3503-4336-6, 000527, 10.1109/SISY60376.2023.10417894
    25. Gaurav Gogia, Parag Rughani, 2024, Chapter 23, 978-981-97-2838-1, 333, 10.1007/978-981-97-2839-8_23
    26. Wenjie Liu, Liming Wang, 2023, A Novel Malware Classification Method Based on Memory Image Representation, 979-8-3503-0048-2, 1404, 10.1109/ISCC58397.2023.10217992
    27. Teh Beng Yen, Joey Li, Shih-Wei Li, 2024, SECvma: Virtualization-based Linux Kernel Protection for Arm, 979-8-3315-2088-5, 579, 10.1109/ACSAC63791.2024.00056
  • Reader Comments
  • © 2019 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(8179) PDF downloads(1152) Cited by(27)

Figures and Tables

Figures(8)  /  Tables(7)

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog