In this paper, we investigate a delayed HIV-1 infection model with immune response. Though a logistic growth is incorporated in the growth of the target cells, our focus is on the effect of delays on the infection dynamics. We first study the existence of steady states, which depends on the basic reproduction number R0. When R0≤1, there is only the infection-free steady state, which is globally asymptotically stable if R0<1. When R0>1, besides the unstable infection-free steady state, there is a unique infected steady state. We then study the local stability of the infected steady state and local Hopf bifurcation at it. The theoretical analysis indicates that the dynamics scenario is complicated. For example, there can be three sequences of critical values, stability switches and double Hopf bifurcation can occur. Some of the theoretical results are also illustrated with numerical simulations.
1.
Introduction
AIDS (acquired immunodeficiency syndrome) is a syndrome caused by HIV (human immunodeficiency virus). HIV is a lentivirus (a subetaoup of retrovirus). It infects vital cells in the human immune system, such as helper T cells (specifically CD4+ T cells), macrophages, and dendritic cells [1]. When the number of CD4+ T cells declines below a critical level, cell-mediated immunity is lost, and the body becomes progressively more susceptible to opportunistic infections, leading to the development of AIDS. Mathematical modeling has contributed a lot to the understanding of HIV infection (see, for example, the review by Perelson and Ribeiro [2] for within-host models).
In the simplest and earliest models of HIV infection, only the key players were taken into account. These models include uninfected target cells (T), productively infected cells (T∗), and free viruses (V). One such model is described by the following system of ordinary differential equations,
For more detail, we refer the readers to Ribeiro and Perelson [3]. Inspired by this model, researchers have proposed many other HIV models by considering, for example, different uninfected target cell growth and incidence, latently infected CD4+ T cells, treatment, drug resistance, and immune response (to name a few, see [4,5,6,7,8,9,10,11,12,13,14,15,16]).
Time delay is commonly observed in many biological processes. For HIV infection, on the one hand, it roughly takes about 1 day for a newly infected cell to become productive and then to be able to produce new virus particles [17]. On the other hand, during CTL response, effector CTLs need time to recognize infected cells and destroy them. Herz et al. [18] were the first to introduce an intracellular delay to describe the time between the initial viral entry into a target cell and subsequent viral production. They obtained the effect of the delay on viral load change. Since then, delayed HIV models have attracted the attention of many researchers. See, for example, [19,20,21,22,23,24,25,26] and the references therein.
In this paper, motivated by the studies in [7,10,27], we propose and study the following delayed HIV model,
Here T(t), T∗(t), V(t), and E(t) represent the densities of uninfected CD4+ T-cells, productively infected CD4+ T-cells, free viruses, and immune effectors at time t, respectively. As in [7,28], k1=ke−ατ1, where α∈[d,δ] is the death rate of infected cells before becoming productive. τ1 denotes the time delay between viral entry and viral production while τ2 stands for the time needed for the CTLs immune response to emerge to control viral replication. The interpretations of the parameters are summarized in Table 1, where their units and ranges will be given in Section 4. The logistic growth in target cells and natural growth of immune effectors combined is a new feature of Model (1.1). Our main focus is on the effects of delays, especially τ2, on the dynamics of (1.1).
The rest of the paper is organized as follows. In Section 2, we present some preliminary results of (1.1a), which include the positivity and boundedness of solution, the existence of steady states. Then we analyze the stability of steady states and possible Hopf bifurcation in Section 3. We conclude the paper with some numerical simulations to illustrate the main theoretical results.
2.
Preliminaries
The suitable phase space for (1.1) is C=C1×C2×C1×R, where Ci=C([−τi,0],R) is the Banach space of all continuous functions from [−τi,0] to R equipped with the supremum norm, i=1, 2. The norm on C is the usual product norm. The nonnegative cone of Ci is C+i=C([−τi,0],R+). To be biologically meaningful, in the sequel, the initial conditions of (1.1) will be always from C+=C+1×C+2×C+1×R+.
For each Φ=(ϕ1,ϕ2,ϕ3,ϕ4)∈C+, by the standard theory of functional differential equations [29], Model (1.1) has a unique and global solution through it. For such a solution, we first claim that T(t)>0 for t>0. In fact, it is clear that there exists t0>0 such that T(t)>0 for t∈(0,t0). Suppose to the contrary that there exists t1>t0 such that T(t)>0 for t∈(0,t1) and T(t1)=0. Then by (1.1), dT(t1)dt=s>0 and hence there exists ε∈(0,t1) such that T(t)<0 for t∈(t1−ε,t1), a contradiction. This proves the claim. Next, with step-by-step method we show that T∗(t)≥0 and V(t)≥0 for t≥0. Note that, for t≥0, by (1.1b) and (1.1c), we have
and
respectively. It follows from (2.1) that T∗(t)≥0 for t∈[0,τ1]. This, combined with (2.2), gives V(t)≥0 for t∈[0,τ1], which together with (2.1) yields T∗(t)≥0 for t∈[0,2τ1]. In turn from (2.2) we have V(t)≥0 for t∈[0,2τ1]. Continuing this way gives the desired result. Finally, from (1.1d), we get
for t≥0 and hence E(t)≥0 for t≥0. Therefore, the solution of (1.1a) with initial condition in C+ is nonnegative.
Next, we consider the boundedness of solutions. Firstly, we obtain from (1.1a) that
for t≥0. It follows that
where
is the unique positive zero of s−dT+rT(1−TTmax). Moreover, if T(0)≤T0 then T(t)≤T0 for t≥0. Secondly, consider the Lyapunov functional
The derivative of L1 along solutions of (1.1) is
where d1=min{δ,d} and M0=rTmax+4s4 (>0). Then lim supt→∞L1(t)≤M0d1. In particular, lim supt→∞T∗(t)≤k1M0kd1. Finally, this combined with (1.1c) and (1.1d) immediately gives
Lastly, we study the lower boundedness of T. For any ε>0, there exists t2>0 such that V(t)≤Nδk1M0ckd1+ε for t≥t2. This, together with (1.1a), gives us
Hence, as ε is arbitrary, we get
To summarize, we have shown the following result.
Proposition 1. The solutions of (1.1) with initial conditions in C+ are nonnegative and bounded. Moreover, the region
is a positively invariant and attracting subset of (1.1) in C+.
In the remaining of this section, we consider the steady states of (1.1). Note that a steady state is a solution of the following system of algebraic equations,
It follows from (2.3c) that V=NδT∗c. Sustituting it into (2.3b) gives
Then T∗=0 or T=c(δ+dxE)Nδk1. When T∗=0, we get the infection-free steady state P0=(T0,0,0,E0), where E0=λEdE. Now assume T=c(δ+dxE)Nδk1. Combining it with E=λE+pT∗dE obtained from (2.3d), we can get after a little computation that
Then
Substituting it into (2.3a) yields
where
Note that G always has a positive zero and it only has one positive zero. However, for infected steady states, we have T>c(δ+dxλEdE)Nδk1 from (2.4), or equivalently, G(c(δ+dxλEdE)Nδk1)>0 or c(δ+dxλEdE)Nδk1<T0. Thus there is an infected steady state if and only if R0>1, where
In summary we have obtained the following result.
Theorem 2.1. (i) If R0≤1 then (1.1) only has the infection-free steady state P0.
(ii) If R0>1 then, besides P0, (1.1) also has a unique infected steady state P1=(T1,T∗1,V1,E1), where
Note that, in epidemiology, R0 is called the basic reproduction number, whose expression can also be derived by the procedure in [30].
3.
Stability and bifurcation
3.1. Global stability of P0
We start with the local stability of the infection-free steady state P0.
Theorem 3.1. (i) If R0<1, then the infection-free steady state P0 of (1.1) is locally asymptotically stable.
(ii) If R0>1, then P0 is unstable.
Proof. The characteristic equation at P0 is
Clearly, −dE and −(d−r+2rT0Tmax)=−sT0−rT0Tmax are eigenvalues and both are negative. The other eigenvalues are roots of the following transcendental equation,
Noting
we can rewrite (3.1) as
or equivalently
(ⅰ) Assume R0<1. We claim that all roots of (3.3) have negative real parts. Otherwise, (3.3) has a root λ=σ+ωi with σ≥0 and σ2+ω2>0 since 0 is not a root by R0>1. Taking moduli of both sides of (3.3) gives
This is impossible as the right side of the above is >1 and R0<1. This proves the claim and hence P0 is locally asymptotically stable if R0<1.
(ⅱ) Assume R0>1. In this case, (3.2) has a positive root. In fact, this follows from the Intermediate Value Theorem and
Therefore, P0 is unstable if R0>1. This completes the proof.
In fact, the local stability of P0 implies its global stability.
Theorem 3.2. If R0<1, then the infection-free steady state P0 of (1.1) is globally asymptotically stable.
Proof. Define the Lyapunov functional
Then the time derivative of W0 along solutions of (1.1) is
Moreover, dW0dt=0 if and only if T∗(t)=0 and V(t)(T(t)−T0)=0. Then one can see that the largest invariant subset of {dW0dt=0} is {P0}. By the Lyapunov-LaSalle invariance principle (see [29,Theorem 5.3.1] or [31,Theorem 3.4.7]) and Theorem 3.1, we conclude that if R0<1 then P0 is globally asymptotically stable.
3.2. Stability of P1 and bifurcation analysis
Recall that P1 exists only when R0>1, which implies that necessarily Nδk1T0c(δ+dxλEdE)>1 as k1=ke−ατ1. The purpose of this paper is to consider the effects of delays on the dynamics. As a result, in the sequel of this section, we always assume that Nδk1T0c(δ+dxλEdE)>1 and denote
Then R0>1 is equivalent to τ1<ˆτ1.
The characteristic equation at P1 is
In the following, we follow the arguments in [23] to first show that P1 is locally stable for τ1∈[0,ˆτ1) and τ2=0. Then for given τ1∈[0,ˆτ1), we discuss the possible bifurcations.
Theorem 3.3. Suppose τ1∈[0,ˆτ1), τ2=0, and 0≤r<d1−T1Tmax. Then the infected steady state P1 is locally asymptotically stable.
Proof. When τ2=0, the characteristic equation (3.4) reduces to
We will prove that all roots of (3.5) have negative real parts in three steps.
Firstly, we show that (3.5) has no roots on the imaginary axis with contradictory arguments. Let λ=iω0 with ω0≥0 be a root of (3.5). Then
Note that the modulus of the right hand side of (3.6) is
However, since
and
it follows that the modulus of the left hand side of (3.6) is strictly larger than that of its right hand side, a contradiction. Thus we have verified that (3.5) has no roots on the imaginary axis.
Secondly, we show that (3.5) has no nonnegative real roots. Again, by contradiction, assume that (3.5) has a root λ0≥0 and we know that e−λ0τ1∈(e−λ0ˆτ1,1]. Noting
we get from (3.5) that
But
as Nδk1T1=c(δ+dxE1), which contradicts with (3.7). This proved that (3.5) has no nonnegative real root.
Finally, we claim that there exists η0>0 such that all roots of (3.5) have negative real parts when 1<R0<1+η0. If this is not true, then there exists a sequence of values for the parameters where R0 (>1) →1 such that for each set of values for the parameters there exists a pair of conjugate roots for (3.5) with positive real parts (which follows from the results just proved above). Note that roots of (3.5) having nonnegative real parts are uniformly bounded. Without loss of generality, we suppose that the sequence of the conjugate roots converges to α0±iβ0, otherwise just consider a subsequence. Then α0≥0 and α0±iβ0 are roots of the characteristic equation of the infection-free steady state P0 when R0=1. However, when R0=1, this characteristic equation has no roots with nonnegative real parts except the simple root 0. Then α0=β0=0, which implies that 0 is a root of at least multiplicity 2 of this characteristic equation, a contradiction. This proves the claim.
Now the proof is done by noting the fact that the roots of (3.5) depend continuously on the parameters.
Theorem 3.4. If τ1=τ2=0 and 0≤r<d1−T1Tmax holds, then the infected steady state P1 is global asymptotically stable.
Proof. We define a Lyapunov functional
Then the time derivative of W1(t) along solutions of system (1.1) is
It is clear that dW1dt=0 if and only if (T,T∗,V,E)=P1. By the Laypunov asymptotic stability theorem, we conclude that P1 of system (1.1) is globally asymptotically stable.
In the following, we study the effects of delays by fixing τ1=τ01∈[0,ˆτ1) and changing τ2. Then the characteristic equation (3.4) can be rewritten as
where
Since
we know that λ=0 is not a root of (3.8). Therefore, for stability changes of P1 to occur, we first look for τ2 where (3.8) has a pair of conjugate roots λ=±iω(τ01,τ2) with ω(τ01,τ2)>0. Substitute λ=iω(τ01,τ2) into (3.8) and then separate the real and imaginary parts to obtain
where
It follows that
Using sin2ωτ2+cos2ωτ2=1, we see that ω(τ01,τ2) satisfies F(ω,τ01)=0, where
with
Therefore, ω(τ01,τ2) is independent of τ2. Denote
which is a finite set. If Iτ01=∅ then P1 is stable for τ1=τ01 and τ2≥0. Now, suppose Iτ01≠∅. For example, this is true if A0+B0<C0 since
and A0+B0+C0>0.
Assume Iτ01={ω1(τ01),…,ωj(τ01)(τ01)}. For j∈N(τ01)={1,…,j(τ01)}, choose the unique angle θj(τ01)∈[0,2π) such that
Now, define
Then the characteristic equation (3.8) at τ2=τn2j has a pair of conjugate eigenvalues λ=±iωj(τ01) for j∈N(τ01) and n∈N. The following result comes from [32,Theorem 2.2].
Theorem 3.5. Let τ01∈[0,ˆτ1). Then the following two statements are true.
(i) If Iτ01=∅, then P1 is locallay asymptotically stable for τ1=τ01 and τ2≥0.
(ii) If Iτ01≠∅, then a pair of simple conjugate pure imaginary roots λ(τn2j(τ01))=±iωj(τ01) of (3.8) exists at τ2=τn2j(τ01) for j∈N(τ01) and n∈N, which crosses the imaginary axis from left to right if δ(τn2j(τ01))>0 and crosses the imaginary axis from right to left if δ(τn2j(τ01))<0, where
By Theorem 3.5, for any τ1∈[0,^τ1), there exists a τ∗2(τ1)∈(0,∞] such that P1 is locally asymptotically stable for τ<τ∗2(τ1).
In general, it is hard to determine whether Iτ1 is empty or not. Moreover, if Iτ1≠∅ and has more than one element, then Theorem 3.5 indicates that there may be stability switches for P1. To get a clear picture of it, we consider the case where τ01=0. Moreover, F(ω,0) in (3.10) reduces to h(ω2), where
with
In this case, δ(τn2j(0)) in Theorem 3.5 is sign(h′(ω2j(0))). As a result, for Hopf bifurcation to occur, we only focus on the situations where h(z) defined by (3.12) has simple positive real zeros.
Though a simple calculation gives b3=d2E+(c+δ+dxE1)2+(sT1+rT1Tmax)2>0, we cannot easily get the signs of the other coefficients. By Descartes' rule of sign, the polynomial h(z) has at most three positive real zeros. In fact, Yan and Li [33] have obtained the conditions on the existence of at least one positive real zero for h(z). To cite the result, let
where ρ is one of the cubic roots of the complex number −q2+√Δ and ϵ=−1+√3i2. Note that when Δ>0, z∗1 is the only real zero of h′(z); when Δ=0, z1=−b34−23√q2, z2=z3=−b34+3√q2 are the three real zeros of h′(z); when Δ<0, −b34+2Re{ρ}, −b34+2Re{ρϵ}, and −b34+2Re{ρˉϵ} are the three real zeros of h′(z) and we arrange them as ˆz1<ˆz2<z∗3.
Proposition 3.1 ([33], Lemma 2.1]).
(i) If b0<0, then h(z) has at least one positive real zero. (ii) If b0≥0, then h(z) has no positive real zero if one of the following conditions holds.
(ii−1) Δ>0 and z∗1<0;
(ii−2) Δ=0 and z∗2<0;
(ii−3) Δ<0 and z∗3<0.
(iii) If b0≥0, then h(z) has at least one positive real zero if one of the following conditions holds.
(iii−1) Δ>0, z∗1>0 and h(z∗1)<0;
(iii−2) Δ=0, z∗2>0 and h(z∗2)<0;
(iii−3) Δ<0, z∗3>0 and h(z∗3)<0.
When τ1=0, by Proposition 3.1 we can have the following result.
Theorem 3.6. Assume τ1=0 and one of the conditions in statement (ii) of Proposition 3.1 holds. Then the infected steady state P1 is locally asymptotically stable for all τ2≥0.
In the following result, we characterize the situations where h(z) has simple positive real zeros, which is not difficult to see by considering the possible graphs of h(z) and h′(z). Recall that h(z) can have at most three positive real zeros.
Proposition 3.2. For the polynomial h(z) defined by (3.12), the following results hold.
(i) h(z) has one simple positive zero and no other positive zeros if and only if (H1): one of the following conditions hold.
(i−1) Δ≥0 and b0<0.
(i−2) Δ>0, b0=0, and z∗1>0.
(i−3) Δ=0, b0=0 and (z∗2=z1>0 or z∗2=z2>z1>0).
(i−4) Δ<0, b0<0 and ˆz2≤0.
(i−5) Δ<0, b0<0, ˆz2>0 and h(ˆz2)<0.
(i−6) Δ<0, b0<0, ˆz2>0, h(ˆz2)>0 and h(z∗3)>0.
(i−7) Δ<0, b0=0 and ˆz2<0<z∗3.
(i−8) Δ<0, b0=0, ˆz1>0, h(ˆz2)>0 and h(z∗3)>0.
(i−9) Δ<0, b0=0, ˆz1>0 and h(ˆz2)<0.
(ii) h(z) has two simple positive zeros and no other positive zeros if and only if (H2) : one of the following conditions holds.
(ii−1) Δ>0, b0>0, z∗1>0 and h(z∗1)<0.
(ii−2) Δ=0, b0>0, z∗2=z1>0 and h(z∗2)<0.
(ii−3) Δ=0, b0>0, z∗2=z2>z1>0 and h(z1)<0.
(ii−4) Δ<0, b0=0, ˆz1≤0<ˆz2 and h(z∗3)<0.
(ii−5) Δ<0, b0>0, ˆz2≤0<z∗3 and h(z∗3)<0.
(ii−6) Δ<0, b0>0, ˆz1≤0<ˆz2 and h(z∗3)<0.
(ii−7) Δ<0, b0>0, ˆz1>0, h(ˆz1)>0 and h(z∗3)<0.
(ii−8) Δ<0, b0>0, ˆz1>0, h(ˆz1)<0 and h(ˆz2)<0.
(ii−9) Δ<0, b0>0, ˆz1>0, h(ˆz1)<0, h(ˆz2)>0 and h(z∗3)>0.
(iii) h(z) has three simple positive zeros and no other positive zeros if and only if (H3): one of the following conditions holds.
(iii−1) Δ<0, b0<0, ˆz2>0, h(ˆz2)>0, and h(z∗3)<0.
(iii−2) Δ<0, b0=0, ˆz1>0, h(ˆz2)>0 and h(z∗3)<0.
If (H1) holds, let ˜z>0 be the unique simple positive zero of h(z) and denote ˜ω=√˜z. Solving (3.11) to obtain the unique ˜θ∈[0,2π). Define
As h′(˜z)>0, we have δ(τn2)=1 and hence the following result holds.
Theorem 3.7. Assume that τ1=0 and assumption (H1) holds. Then there exists a sequence 0<τ02<τ12<τ22<⋯ such that P1 is locally asymptotically stable for τ2∈[0,τ02) and unstable for τ>τ02, and system (1.1) undergoes a Hopf bifurcation at P1 when τ2=τn2 for n∈N.
Now, assume (H2) holds. Let ˜z2<˜z1 be the only positive real zeros of h(z), which are also simple. Similarly as for the case of (H1), we can get two increasing positive sequences {τn21} and {τn22}, associated with ˜z1 and ˜z2, respectively. Since h′(˜z1)>0>h′(˜z2), we easily see that τ021<τ022. Since ˜ω1=√˜z1>√˜z2=˜ω2, we have 2π˜ω1<2π˜ω2. Thus we define
Such k exists due to τl+121−τl22→−∞ as l→∞. Then the first few Hopf bifurcation values can be ordered as
Theorem 3.8. Assume τ1=0 and (H2) holds. Given n∈N and j∈{1,2}, system (1.1) undergoes Hopf bifurcation at τ2=τn2j if τn2j≠τl2(3−j) for all l∈N. Furthermore, the stability of P1 switches off (namely, it becomes unstable) when τ2 crosses τ021, …, τk21 and switches on (namely, it becomes stable) when τ2 crosses τ022, …, τk−122. In other words, P1 is stable when τ2∈[0,τ021)∪(τ022,τ121)∪⋯∪(τk−122,τk21) and unstable when τ2∈(τ021,τ022)∪⋯∪(τk−121,τk−122)∪(τk21,∞).
We mention that we can study the global continuation of Hopf bifurcation in Theorem 3.8 as in Li and Shu [34]. We believe that the Hopf bifurcation branches are bounded and each joins a pair of τn21 and τn22 for n∈N. As a result, for τ2∈(τk+121,τk22), there will be two stable periodic orbits.
Also in Theorem 3.8, we exclude the situation where τn21=τl22 for some n, l∈N. If this happens, then n>l since τ021<τ022 and 2π˜ω1<2π˜ω2. In this critical situation, as τ2 crosses this common critical value, two pairs of purely imaginary eigenvalues ±i˜ω1 and ±i˜ω2 appear and all other eigenvalues have nonzero real parts. Therefore, a double Hopf bifurcation occurs.
Theorem 3.9. Assume τ1=0 and (H2) holds. If there exist integers n>l≥0 such that τn21=τl22=τ20, then (1.1) undergoes a double Hopf bifurcation at P1 when τ2=τ20.
When (H3) holds, we can similarly get three sequences of critical values for τ2. Similar results as those in Theorem 3.8 and Theorem 3.9 can be obtained. Moreover, the global Hopf bifurcation associated with the third sequence is unbounded (one can refer to Li et al. [35] for a similar discussion).
4.
Numerical simulations
In this paper, we rigorously analyzed an HIV-infection model with CTL-immune response and two time delays. The model incorporates a logistic growth term for the target cell growth and a natural resource for the immune effectors. The basic reproduction number R0 played an important in the infection dynamics. If R0<1 then the infection-free steady state is globally asymptotically stable. Note that R0 explicitly depends on τ1. It follows that if τ1 is large enough then the virus will be cleared. We emphasize that this has not been noted in most existing study. Of course, in the real situation, this delay between viral entry and subsequent viral production usually is not very big. This leads to the complicated dynamics when R0>1. In this case, the unique infected steady state could be stable or unstable, depending on the parameter values. In particular, we focused on the effects of time delays. Theoretical results indicate that there can be Hopf bifurcation, double Hopf bifurcation and stability switches.
We conclude this paper with some numerical simulations to illustrate the above mentioned main results. The ranges of the parameters except the delays are summarized in Table 2.
For simplicity, we use the same initial condition (T0,T∗0,V0,E0)=(100,0,10−2,0.6) in all simulations.
First, we take s=10, d=0.01, r=0.03, Tmax=1500, α=0.02, δ=0.3, dx=0.01, N=21, c=3, λE=1, p=0.3, dE=0.1, τ1=1.2 and k=2.4×10−5. Then R0=0.1680<1. By Theorem 3.2, the infection-free steady state P0=(1366,0,0,10) is globally asymptotically stable (see Figure 1 with τ2=4).
Next, we only change k to k=2.4×10−3 and keep the others as above. In this case, we have ˆτ1=142.2801. Take τ1=1.2∈[0,ˆτ1). Then R0=16.8038>1. Through numerical calculations, we get the first few critical values τ021=2.9940 and τ121=32.9360 associated with ω1=0.2098 and τ022=28.8271 associated with ω2=0.1066. By Theorem 3.5, the infected steady state P1=(171.7615,14.8383,31.1604,54.5149) is locally asymptotically stable for τ2<τ021 (see Figure 2 with τ2=2).
In fact, numerical simulations indicates that there is Hopf bifurcation for τ2∈(τ021,τ022) (see Figure 3 with τ2=5) and there is a stability switch at τ2=τ022, that is, P1 is locally asymptotically stable for τ2∈(τ022,τ121) (see Figure 4 with τ2=32).
In the following we illustrate this more clearly with the special case where τ1=0. We distinguish three cases.
Case 1: (ⅱ) of Proposition 3.1 holds. We take s=5, d=0.2, r=0.03, Tmax=1500, α=0.2, δ=0.3, dx=0.01, N=2800, c=15, λE=1, p=0.3, dE=0.1, τ1=0 and k=2.4×10−3. Then R0=9.8484>1 and system (1.1) has the unique infected steady state P1=(4.5277,6.9510,389.2547,30.8529). In this case, Δ=4.1507×106>0, b0=0.6078>0, and z∗1=−1.9553×10−2<0. This means that (ⅱ-1) of Proposition 3.1 holds. It follows from Theorem 3.6 that the infected steady state P1 is locally asymptotically stable for all τ2≥0 (see Figure 5 with τ2=1).
Case 2: Assumption (H1) holds. We choose the parameter values s=10, d=0.01, r=0.25, Tmax=1500, α=0.02, δ=0.3, dx=0.01, N=21, c=3, λE=1, p=0.3, dE=0.1 and k=2.4×10−4. Then R0=1.8655>1 and system (1.1) has the unique infected steady state P1=(1448.2000,10.9961,23.0918,42.9882). Note that Δ=−0.2247<0, b0=−6.0220×10−4<0, and ˆz2=−0.0442<0, namely, assumption (H1) (ⅰ-4) holds. We can get ω=0.1519 and τj2=4.3175+2jπω for j∈N. It follows from Theorem 3.7 that the infected steady state P1 is locally asymptotically stable for τ2∈[0,τ02) (see Figure 6 with τ2=4<τ02) and unstable for τ>τ02. Moreover, system (1.1) undergoes Hopf bifurcation at τ2=τj2 for j∈N. Figure 7 supports this with τ2=5>τ02. Figure 8 provides the bifurcation diagram.
Case 3: Assumption (H2) holds. This time we replace k with k=2.4×10−3, and r with r=0.03, and keep the other parameter values as in Case 2. It follows that R0=17.2119>1 and system (1.1) has the unique infected steady state P1=(168.9145,15.0443,31.5930,55.1329). Moreover, Δ=−0.4441<0, b0=3.0321×10−4>0, ˆz1=−11.1936<0<ˆz2=0.0018>0, and h(z∗3)=−9.3650×10−4<0. Therefore, assumption (H2) (ⅱ-6) holds. In this case, we have ˜ω1=0.2845, ˜ω2=0.1401, τj21=2.0033+2jπ˜ω1, and τj22=22.3237+2jπ˜ω2 for j∈N. Then the first few Hopf bifurcation values are ordered as τ021<τ022<τ121<τ221<τ122<⋯. By Theorem 3.8, the infected steady state P1 is locally asymptotically stable for τ2∈[0,τ021)∪(τ022,τ121) (see Figure 9 and Figure 10 with τ2=1.5<τ022 and τ2=24∈(τ022,τ121), respectively) and is unstable for τ2∈(τ021,τ022)∪(τ121,∞). Thus there are stability switches for P1. Moreover, there are supercritical Hopf bifurcation at τ2=τj21 and subcritical Hopf bifurcation at τ2=τj22 (see Figures 11 and 12). The numerical simulations also strongly indicate that the global Hopf bifurcation branches are bounded and each branch connects a pair τj21 and τj22, which we will confirm in a future work.
Acknowledgment
The authors are grateful to the anonymous reviewers for their constructive comments. JW is supported by the Nanhu Scholar Program for Young Scholars of Xinyang Normal University. CQ is supported by Scientific Research Foundation of Graduate School of Xinyang Normal University (No. 2018KYJJ12). XW is supported by the NSFC (No. 11771374), the Nanhu Scholar Program for Young Scholars of Xinyang Normal University, the Program for Science and Technology Innovation Talents in Universities of Henan Province (17HASTIT011). YC is supported by NSERC of Canada.
Conflict of interest
All authors declare no conflicts of interest in this paper.