Fractional differential equation modeling of the HBV infection with time delay and logistic proliferation

In this article, a fractional-order differential equation model of HBV infection was proposed with a Caputo derivative, delayed immune response, and logistic proliferation. Initially, infection-free and infection equilibriums and the basic reproduction number were computed. Thereafter, the stability of the two equilibriums was analyzed based on the fractional Routh–Hurwitz stability criterion, and the results indicated that the stability will change if the time delay or fractional order changes. In addition, the sensitivity of the basic reproduction number was analyzed to find out the most sensitive parameter. Lastly, the theoretical analysis was verified by numerical simulations. The results showed that the time delay of immune response and fractional order can significantly affect the dynamic behavior in the HBV infection process. Therefore, it is necessary to consider time delay and fractional order in modeling HBV infection and studying its dynamics.

Introduction

Hepatitis B virus (HBV) can attack the liver and cause both acute and chronic diseases and further lead to fibrosis, cirrhosis, or even cancer. It is estimated that 296 million people have chronic hepatitis B, and 1.5 million new infections are reported each year; 820 000 people died of hepatitis B infections in 2019 (1). Therefore, HBV has become a major public health problem affecting human health (2).

Mathematical modeling and analysis of infectious viruses help understand the infection mechanism and realize the disease progression (3–5). Furthermore, mathematical modeling can also provide new insights to find the key factors to treat infectious diseases (6). In 1996, the basic ordinary differential equation (ODE) model of HBV infection was established with uninfected cells, infected cells, and free viruses (7). This is an early mathematical model for studying the spread of viruses. As research progresses, the mathematical modeling of virus transmission has become more and more complicated. For instance, Peter et al. (3) established a deterministic ODE model with six compartments to study the transmission dynamics of measles and obtained the best fit using available data, which could help health workers in decision-making and policymakers to frame policies to eradicate the spread of measles in Nigeria. Mayowa et al. (5) divided the population into six classes and formulated a six-compartmental deterministic model to investigate the effect of vaccination on the dynamics of tuberculosis in a given population. All the aforementioned mathematical models are based on ordinary differential equations with bilinear incidence rate.

Subsequently, a large number of dynamic models were proposed to describe and analyze virus infection according to different biological mechanisms (8–13); for example, because hepatocytes have the ability to regenerate, the models are constrained by the number of healthy and infected hepatocytes. Li et al. (10) developed a logistic growth model of HBV. Moreover, in order to characterize the time of a body's immune response after the virus infection of target cells, time delay has been considered. Therefore, Zhang et al. (13) proposed a susceptible-vaccinated-exposed-infectious-removed (SVEIR) epidemic model with two time delays and constructed a Lyapunov function to discuss the asymptotic stability of the positive equilibrium point. Babasola et al. (14) modeled the spread of COVID-19 with a convex incidence rate incorporated with a time delay and proved that delay can destabilize the system and lead to periodic oscillation.

In recent years, a fractional derivative for describing memory, history, and heredity effects in modeling physical, chemical, financial, and biological systems has received increasing attention (15–28). For example, Diethelm (29) used a fractional-order model to simulate the dynamics of a dengue fever outbreak. The results showed that the simulation accuracy of the fractional-order model is much higher than that of the integer-order derivative. Gilberto et al. (30) proposed a fractional-order model to research the dynamics of influenza A (H1N1), and the results showed that the fractional-order model was in good agreement with real data. Similarly, Ogunrinde et al. (27) divided the population into five classes and proposed a fractional-order differential equation model to study COVID-19. The basic reproduction number was calculated by the spectral radius method, and the stability analysis of the model was carried out by constructing the Lyapunov function. Finally, the parameters were estimated by collected data, and the model can offer guidance to policymakers.

In addition to the mathematical modeling of fractional differential equations for the aforementioned infectious diseases, there are also many studies that use the fractional-order model to characterize the process of HBV infections (31–33). For example, Simelane and Dlamini (33) established a fractional-order HBV model with a saturated incidence rate by using the Caputo fractional derivatives. Then, the basic reproduction number was calculated, and the stability of the equilibriums was discussed. The simulation results demonstrated that the fractional-order model is more appropriate for modeling HBV transmission dynamics than the integer-order model. The time of HBV entry into the healthy liver cells and the production of new virus particles should be taken into account; therefore, Gao et al. (32) established a three-dimensional delayed fractional-order HBV model, which included healthy hepatocytes, infected hepatocytes, and free viruses, as follows:

$\begin{array}{ll}\{\begin{array}{l}{{\displaystyle {}_0{}^{CF}D}}_t^{{\sigma }_1}x(t)={{\displaystyle {\lambda }_1}}^{{\sigma }_1}-{{\displaystyle {\mu }_1}}^{{\sigma }_1}x(t)-{{\displaystyle {\beta }_1}}^{{\sigma }_1}x(t)v(t)+{{\displaystyle {\delta }_1}}^{{\sigma }_1}y(t),\hfill \\ {{\displaystyle {}_0{}^{CF}D}}_t^{{\sigma }_1}y(t)={{\displaystyle {\beta }_1}}^{{\sigma }_1}x(t)v(t)-({{\displaystyle {\alpha }_1}}^{{\sigma }_1}+{{\displaystyle {\delta }_1}}^{{\sigma }_1})y(t),\hfill \\ {{\displaystyle {}_0{}^{CF}D}}_t^{{\sigma }_2}v(t)={{\displaystyle c_1}}^{{\sigma }_2}y(t-\tau )e^{-\rho \tau }-{{\displaystyle {\gamma }_1}}^{{\sigma }_2}v(t).\hfill \end{array}& (1)\end{array}$

This model has not considered the cytotoxic T lymphocyte (CTL) and alanine aminotransferase (ALT) levels, which reflect the extent of liver damage. Therefore, the items of CTL and ALT will be considered in our established model.

However, until now, no study has been designed to analyze the dynamics of HBV involving logistic proliferation, time delay, and items of CTL and ALT by fractional-order differential equations. Motivated by the aforementioned discussion, we proposed a fractional-order differential equation model with time delay and logistic proliferation in order to better understand the transmission mechanism of HBV in the human body.

The remaining part of this article is organized as follows: Section Mathematical model deals with the formulation of the model. Section Equilibriums and the basic reproduction number discusses the infection-free and infection equilibriums and the basic reproduction number. Section Equilibriums and the basic reproduction number discusses the stability analysis of the two equilibriums and analyzes the sensitivity of the basic reproduction number. Section Numerical simulation gives an account of the numerical simulations of equilibriums and the Hopf bifurcation. Finally, Section Conclusion and discussion comprises the conclusion and discussion.

Mathematical model

Therefore, based on our work (34), we proposed a fractional-order differential equation model with time delay and logistic proliferation as follows:

$\begin{array}{ll}\{\begin{array}{l}\frac{d^{\alpha }x(t)}{dt}={\xi }^{\alpha }+r^{\alpha }x(t)\left(1-\frac{x(t)+y(t)}{T_{{\max }^{\alpha }}}\right)-d^{\alpha }x(t)-b^{\alpha }x(t)v(t),\hfill \\ \frac{d^{\alpha }y(t)}{dt}=b^{\alpha }x(t)v(t)-a^{\alpha }y(t)-k_1^{\alpha }y(t-\tau )z(t-\tau ),\hfill \\ \frac{d^{\alpha }v(t)}{dt}=k^{\alpha }y(t)-{\epsilon }^{\alpha }v(t)-k_2^{\alpha }y(t-\tau )z(t-\tau ),\hfill \\ \frac{d^{\alpha }z(t)}{dt}=k_3^{\alpha }y(t-\tau )z(t-\tau )-k_4^{\alpha }z(t),\hfill \\ \frac{d^{\alpha }w(t)}{dt}=k_5^{\alpha }+k_6^{\alpha }y(t)z(t)-k_7^{\alpha }w(t).\hfill \end{array}& (2)\end{array}$

where the variables x, y, v, z, and w represent the uninfected cells, infected cells, viruses, CTL level, and ALT level, respectively; ξ and r are the production rate and proliferation rate of uninfected cells, respectively; T_max is the maximum hepatocyte count in the liver; d is the death rate of uninfected cells; b is the infection rate of uninfected cells to become infected cells; a is the death rate of infected cells; k₁ represents the cure rate of infected cells by CTL; k and ε are the production rate and death rate of free viruses, respectively; k₂ represents the clearance rate of free viruses by CTL; k₃ and k₄ are the production rate and death rate of CTL, respectively; k₅ is the natural production rate of ALT; k₆ is the production rate of ALT from infected cells; k₇ is the death rate; τ is time delay with the order of α (0 < α ≤ 1); and $\frac{d^{\alpha }x(t)}{dt}$, $\frac{d^{\alpha }y(t)}{dt}$, $\frac{d^{\alpha }v(t)}{dt}$, $\frac{d^{\alpha }z(t)}{dt}$, and $\frac{d^{\alpha }w(t)}{dt}$ denote the Caputo fractional derivatives. Hence, $\frac{d^{\alpha }x(t)}{dt}$ is defined as follows:

$\begin{array}{ll}\frac{d^{\alpha }{x}_i}{dt}=I^{n-\alpha }\frac{d^{n}x}{d{t}^n}=\frac{1}{\Gamma (n-\alpha )}{\displaystyle {\int }_0^t}{(t-s)}^{(n-\alpha -1)}{x}^{(n)}(s)ds& (3)\end{array}$

where n − 1 < α < n, n ∈ ℕ and Γ(▪) is the gamma function. When 0 < α < 1,

$\begin{array}{ll}\frac{d^{\alpha }x}{dt}=\frac{1}{\Gamma (1-\alpha )}{\displaystyle {\int }_0^t}\frac{x^{\prime }(s)}{{(t-s)}^{\alpha }}ds& (4)\end{array}$

Based on the aforementioned model, the equilibriums and stability analysis are discussed in Section Equilibriums and the basic reproduction number.

Equilibriums and the basic reproduction number

In the following paragraphs, the equilibriums and the basic reproduction number are discussed.

Equilibriums

The method to compute the equilibrium is to set $\frac{d^{\alpha }x(t)}{dt}=0$, $\frac{d^{\alpha }y(t)}{dt}=0$, $\frac{d^{\alpha }v(t)}{dt}=0$, $\frac{d^{\alpha }z(t)}{dt}=0$, and $\frac{d^{\alpha }w(t)}{dt}=0$. Hence, we get the following equations:

$\begin{array}{ll}\{\begin{array}{l}{\xi }^{\alpha }-d^{\alpha }x(t)+r^{\alpha }x(t)\left(1-\frac{x(t)+y(t)}{T_{{\max }^{\alpha }}}\right)-b^{\alpha }x(t)v(t)=0,\hfill \\ b^{\alpha }x(t)v(t)-a^{\alpha }y(t)-k_1^{\alpha }y(t-\tau )z(t-\tau )=0,\hfill \\ k^{\alpha }y(t)-{\epsilon }^{\alpha }v(t)-k_2^{\alpha }y(t-\tau )z(t-\tau )=0,\hfill \\ k_3^{\alpha }y(t-\tau )z(t-\tau )-k_4^{\alpha }z(t)=0,\hfill \\ k_5^{\alpha }+k_6^{\alpha }y(t)z(t)-k_7^{\alpha }w(t)=0.\hfill \end{array}& (5)\end{array}$

The infection-free equilibrium E₀ denotes x ≠ 0, w ≠ 0, y = v = z = 0; thus, the infection-free equilibrium is as follows:

$\begin{array}{l}{E}_0=\left(x_0,y_0,v_0,z_0,w_0\right)=\\ \left(\frac{T_{max}^{\alpha }}{2r^{\alpha }}\left[-(d^{\alpha }-r^{\alpha })+\sqrt{{(d^{\alpha }-r^{\alpha })}^2+\frac{4{\xi }^{\alpha }{r}^{\alpha }}{T_{max}^{\alpha }}}\right],0,0,0,\frac{k_5^{\alpha }}{k_7^{\alpha }}\right)\end{array}$

Similarly, the infection equilibrium E₁, which denotes x ≠ 0, y ≠ 0, v ≠ 0, z ≠ 0, w ≠ 0, was computed by the following equations:

$\begin{array}{ll}\{\begin{array}{l}{\xi }^{\alpha }-d^{\alpha }{x}^{*}+r^{\alpha }{x}^{*}\left(1-\frac{x^{*}+y^{*}}{T_{{\max }^{\alpha }}}\right)-b^{\alpha }xv=0,\hfill \\ b^{\alpha }{x}^{*}{v}^{*}-a^{\alpha }{y}^{*}-k_1^{\alpha }{y}^{*}{z}^{*}=0,\hfill \\ k^{\alpha }{y}^{*}-{\epsilon }^{\alpha }{v}^{*}-k_2^{\alpha }{y}^{*}{z}^{*}=0,\hfill \\ k_3^{\alpha }{y}^{*}{z}^{*}-k_4^{\alpha }{z}^{*}=0,\hfill \\ k_5^{\alpha }+k_6^{\alpha }{y}^{*}{z}^{*}-k_7^{\alpha }{w}^{*}=0.\hfill \end{array}& (6)\end{array}$

The previous equations were solved, and the infection equilibrium was obtained as follows:

$x^{*}=-\frac{B}{3A}+\sqrt[3]{-\frac{q}{2}+\sqrt{\frac{q^2}{4}+\frac{p^3}{27}}}+\sqrt[3]{-\frac{q}{2}-\sqrt{\frac{q^2}{4}+\frac{p^3}{27}}},y^{*}\mathrm{\text{=}}\frac{k_4^{\alpha }}{k_3^{\alpha }},v^{*}\mathrm{\text{=}}\frac{a^{\alpha }{k}_2^{\alpha }{k}_4^{\alpha }+k^{\alpha }{k}_1^{\alpha }{k}_4^{\alpha }}{b^{\alpha }{k}_2^{\alpha }{k}_3^{\alpha }{x}^{*}+{\epsilon }^{\alpha }{k}_1^{\alpha }{k}_3^{\alpha }},z^{*}=\frac{k^{\alpha }}{k_2^{\alpha }}-\frac{{\epsilon }^{\alpha }(a^{\alpha }{k}_2^{\alpha }-k^{\alpha }{k}_1^{\alpha })}{k_2^{\alpha }(b^{\alpha }{k}_2^{\alpha }{x}^{*}+{\epsilon }^{\alpha }{k}_1^{\alpha })},w^{*}=\frac{k_5^{\alpha }}{k_7^{\alpha }}+\frac{k^{\alpha }{k}_4^{\alpha }{k}_6^{\alpha }}{k_2^{\alpha }{k}_3^{\alpha }{k}_7^{\alpha }}-\frac{{\epsilon }^{\alpha }{k}_4^{\alpha }{k}_6^{\alpha }(a^{\alpha }{k}_2^{\alpha }-k^{\alpha }{k}_1^{\alpha })}{k_2^{\alpha }{k}_3^{\alpha }{k}_7^{\alpha }(b^{\alpha }{k}_2^{\alpha }{x}^{*}+{\epsilon }^{\alpha }{k}_1^{\alpha })}$.

where

$\begin{array}{l}\begin{array}{c}A=\frac{b^{\alpha }{r}^{\alpha }{k}_2^{\alpha }}{T_{max}^{\alpha }},B=-\left[b^{\alpha }{k}_2^{\alpha }(r^{\alpha }-d^{\alpha }-\frac{r^{\alpha }{k}_4^{\alpha }}{T_{max}^{\alpha }{k}_3^{\alpha }})-\frac{{\epsilon }^{\alpha }{r}^{\alpha }{k}_1^{\alpha }}{T_{max}^{\alpha }}\right],\hfill \\ C=-\left[b^{\alpha }{\xi }^{\alpha }{k}_2^{\alpha }+{\epsilon }^{\alpha }{k}_1^{\alpha }(r^{\alpha }-d^{\alpha }-\frac{r^{\alpha }{k}_4^{\alpha }}{T_{max}^{\alpha }{k}_3^{\alpha }})-\frac{b^{\alpha }{k}_4^{\alpha }(a^{\alpha }{k}_2^{\alpha }+k^{\alpha }{k}_1^{\alpha })}{k_3^{\alpha }}\right],\\ D=-{\epsilon }^{\alpha }{\xi }^{\alpha }{k}_1^{\alpha },p=\frac{3AC-B^2}{3A^2},q=\frac{27A^{2}D-9ABC+2B^3}{27A^3}.\hfill \end{array}\hfill \end{array}$

Thus, the infection equilibrium is as follows:

$\begin{array}{l}{E}_1=(x^{*},y^{*},v^{*},z^{*},w^{*})\\ =(x^{*},\frac{k_4^{\alpha }}{k_3^{\alpha }},\frac{a^{\alpha }{k}_2^{\alpha }{k}_4^{\alpha }+k^{\alpha }{k}_1^{\alpha }{k}_4^{\alpha }}{b^{\alpha }{k}_2^{\alpha }{k}_3^{\alpha }{x}^{*}+{\epsilon }^{\alpha }{k}_1^{\alpha }{k}_3^{\alpha }},\frac{k^{\alpha }}{k_2^{\alpha }}-\frac{{\epsilon }^{\alpha }(a^{\alpha }{k}_2^{\alpha }-k^{\alpha }{k}_1^{\alpha })}{k_2^{\alpha }(b^{\alpha }{k}_2^{\alpha }{x}^{*}+{\epsilon }^{\alpha }{k}_1^{\alpha })},\frac{k_5^{\alpha }}{k_7^{\alpha }}\\ +\frac{k^{\alpha }{k}_4^{\alpha }{k}_6^{\alpha }}{k_2^{\alpha }{k}_3^{\alpha }{k}_7^{\alpha }}-\frac{{\epsilon }^{\alpha }{k}_4^{\alpha }{k}_6^{\alpha }(a^{\alpha }{k}_2^{\alpha }-k^{\alpha }{k}_1^{\alpha })}{k_2^{\alpha }{k}_3^{\alpha }{k}_7^{\alpha }(b^{\alpha }{k}_2^{\alpha }{x}^{*}+{\epsilon }^{\alpha }{k}_1^{\alpha })})\end{array}$

Basic reproduction number

The basic reproduction number can be calculated by the method of integral operator spectral radius given as follows:

$\begin{array}{l}{R}_0=\rho (F{V}^{-1})\end{array}$

Thus, the basic reproduction number of E₀ is as follows:

$\begin{array}{l}{R}_0=\frac{b^{\alpha }{k}^{\alpha }{x}_0}{a^{\alpha }{\epsilon }^{\alpha }},\end{array}$

where

$F=\left[\begin{array}{ccc}{b}^{\alpha }v& 0& b^{\alpha }x\\ 0& k^{\alpha }y& 0\end{array}\right],V=\left[\begin{array}{ccc}0& a^{\alpha }+k_1^{\alpha }z{e}^{-\lambda \tau }& 0\\ 0& k_2^{\alpha }z{e}^{-\lambda \tau }& {\epsilon }^{\alpha }\end{array}\right].$

Similarly, the basic reproduction number of E₁ is as follows:

$\begin{array}{l}{R}_1=\frac{b^{\alpha }{k}^{\alpha }{x}^{*}}{{\epsilon }^{\alpha }\left(a^{\alpha }+\frac{k^{\alpha }{k}_1^{\alpha }}{k_2^{\alpha }}-\frac{{\epsilon }^{\alpha }{k}_1^{\alpha }(a^{\alpha }{k}_2^{\alpha }-k^{\alpha }{k}_1^{\alpha })}{k_2^{\alpha }(b^{\alpha }{k}_2^{\alpha }{x}^{*}+{\epsilon }^{\alpha }{k}_1^{\alpha })}\right)}\end{array}$

Stability and sensitivity analyses

The local asymptotic stability of E₀ and E₁ is discussed in this part.

First, the Jacobi matrix was computed as follows:

$\begin{array}{ll}Jac=\left[\begin{array}{ccccc}-d^{\alpha }+r^{\alpha }-\frac{2r^{\alpha }x+y}{T_{max}^{\alpha }}& -\frac{r^{\alpha }x}{T_{max}^{\alpha }}& -b^{\alpha }x& 0& 0\\ b^{\alpha }v& -a^{\alpha }-k_1^{\alpha }z{e}^{-S^{\alpha }\tau }& b^{\alpha }x& -k_1^{\alpha }y{e}^{-S^{\alpha }\tau }& 0\\ 0& k^{\alpha }-k_2^{\alpha }z{e}^{-S^{\alpha }\tau }& -{\epsilon }^{\alpha }& -k_2^{\alpha }y{e}^{-S^{\alpha }\tau }& 0\\ 0& k_3^{\alpha }z{e}^{-S^{\alpha }\tau }& 0& k_3^{\alpha }y{e}^{-S^{\alpha }\tau }-k_4^{\alpha }& 0\\ 0& k_6^{\alpha }z& 0& k_6^{\alpha }y& -k_7^{\alpha }\end{array}\right]& (7)\end{array}$

Based on the previous Jacobi matrix, we got the characteristic determinant:

$\begin{array}{l}|S^{\alpha }I-Jac|\\ =\left|\begin{array}{ccccc}{S}^{\alpha }+d^{\alpha }-r^{\alpha }+\frac{2r^{\alpha }x+y}{T_{max}^{\alpha }}& \frac{r^{\alpha }x}{T_{max}^{\alpha }}& b^{\alpha }x& 0& 0\\ -b^{\alpha }v& S^{\alpha }+a^{\alpha }+k_1^{\alpha }z{e}^{-S^{\alpha }\tau }& -b^{\alpha }x& k_1^{\alpha }y{e}^{-S^{\alpha }\tau }& 0\\ 0& -k^{\alpha }+k_2^{\alpha }z{e}^{-S^{\alpha }\tau }& S^{\alpha }+{\epsilon }^{\alpha }& k_2^{\alpha }y{e}^{-S^{\alpha }\tau }& 0\\ 0& -k_3^{\alpha }z{e}^{-S^{\alpha }\tau }& 0& S^{\alpha }+k_4^{\alpha }-k_3^{\alpha }y{e}^{-S^{\alpha }\tau }& 0\\ 0& -k_6^{\alpha }z& 0& -k_6^{\alpha }y& S^{\alpha }+k_7^{\alpha }\end{array}\right|\end{array}$

Let S^α = λ, then the simplified characteristic determinant is as follows:

$\begin{array}{l}|\lambda I-Jac|=\left|\begin{array}{ccccc}\lambda +d^{\alpha }-r^{\alpha }+\frac{2r^{\alpha }x+y}{T_{max}}& \frac{r^{\alpha }x}{T_{max}}& b^{\alpha }x& 0& 0\\ -b^{\alpha }v& \lambda +a^{\alpha }+k_1^{\alpha }z{e}^{-\lambda \tau }& -b^{\alpha }x& k_1^{\alpha }y{e}^{-\lambda \tau }& 0\\ 0& -k^{\alpha }+k_2^{\alpha }z{e}^{-\lambda \tau }& \lambda +{\epsilon }^{\alpha }& k_2^{\alpha }y{e}^{-\lambda \tau }& 0\\ 0& -k_3^{\alpha }z{e}^{-\lambda \tau }& 0& \lambda +k_4^{\alpha }-k_3^{\alpha }y{e}^{-\lambda \tau }& 0\\ 0& -k_6^{\alpha }z& 0& -k_6^{\alpha }y& \lambda +k_7^{\alpha }\end{array}\right|\end{array}$

Local asymptotic stability of the infection-free equilibrium

The characteristic determinant at the infection-free equilibrium (E₀) is as follows:

$\begin{array}{l}\left|\lambda I-Jac\right|=\left|\begin{array}{ccccc}\lambda +d^{\alpha }-r^{\alpha }+\frac{2r^{\alpha }{x}_0}{T_{{\max }^{\alpha }}}& \frac{r^{\alpha }{x}_0}{T_{{\max }^{\alpha }}}& b^{\alpha }x& 0& 0\\ 0& \lambda +a^{\alpha }& -b^{\alpha }{x}_0& 0& 0\\ 0& -k^{\alpha }& \lambda +{\epsilon }^{\alpha }& 0& 0\\ 0& 0& 0& \lambda +k_4^{\alpha }& 0\\ 0& 0& 0& 0& \lambda +k_7^{\alpha }\end{array}\right|\\ =(\lambda +d^{\alpha }-r^{\alpha }+\frac{2r^{\alpha }{x}_0}{T_{{\max }^{\alpha }}})(\lambda +k_4^{\alpha })(\lambda +k_7^{\alpha })[(\lambda +a^{\alpha })(\lambda +{\epsilon }^{\alpha })\\ -b^{\alpha }{k}^{\alpha }{x}_0]\end{array}$

When |λI − Jac| = 0, the eigenvalues are ${\lambda }_1=-d^{\alpha }+r^{\alpha }-\frac{2r^{\alpha }{x}_0}{T_{max}^{\alpha }}$, ${\lambda }_2=k_4^{\alpha }$, ${\lambda }_3=k_7^{\alpha }$, ${\lambda }_4=\frac{-(a^{\alpha }+{\epsilon }^{\alpha })+\sqrt{{(a^{\alpha }+{\epsilon }^{\alpha })}^2-4(a^{\alpha }{\epsilon }^{\alpha }-b^{\alpha }{k}^{\alpha }{x}_0)}}{2}$, and ${\lambda }_5=\frac{-(a^{\alpha }+{\epsilon }^{\alpha })-\sqrt{{(a^{\alpha }+{\epsilon }^{\alpha })}^2-4(a^{\alpha }{\epsilon }^{\alpha }-b^{\alpha }{k}^{\alpha }{x}_0)}}{2}$.

Since d > r and $R_0=\frac{b^{\alpha }{k}^{\alpha }{x}_0}{a^{\alpha }{\epsilon }^{\alpha }}<1$, we have λ_1,2,3,4,5 < 0. Thus, $|arg(S_{1,2,3,4,5})|>\frac{\alpha \pi }{2}$.

Thus, we get the conclusion that when $R_0=\frac{b^{\alpha }{k}^{\alpha }{x}_0}{a^{\alpha }{\epsilon }^{\alpha }}<1$, E₀ is locally asymptotically stable.

Local asymptotic stability of the infection equilibrium

The characteristic determinant at the infection equilibrium (E₁) is as follows:

$\begin{array}{l}\left|\begin{array}{ccccc}\lambda +d^{\alpha }-r^{\alpha }+\frac{2r^{\alpha }x+y}{T_{max}^{\alpha }}& \frac{r^{\alpha }x}{T_{max}^{\alpha }}& b^{\alpha }x& 0& 0\\ -b^{\alpha }v& \lambda +a^{\alpha }+k_1^{\alpha }z{e}^{-\lambda \tau }& -b^{\alpha }x& k_1^{\alpha }y{e}^{-\lambda \tau }& 0\\ 0& -k^{\alpha }+k_2^{\alpha }z{e}^{-\lambda \tau }& \lambda +{\epsilon }^{\alpha }& k_2^{\alpha }y{e}^{-\lambda \tau }& 0\\ 0& -k_3^{\alpha }z{e}^{-\lambda \tau }& 0& \lambda +k_4^{\alpha }-k_3^{\alpha }y{e}^{-\lambda \tau }& 0\\ 0& -k_6^{\alpha }z& 0& -k_6^{\alpha }y& \lambda +k_7^{\alpha }\end{array}\right|\end{array}$

$\begin{array}{l}=\left(\lambda +k_7^{\alpha }\right)\left\{\begin{array}{c}\left(\lambda +A_0\right)\left[\begin{array}{c}\left({\lambda }^2+\left(a^{\alpha }+{\epsilon }^{\alpha }\right)\lambda +a^{\alpha }{\epsilon }^{\alpha }+k_1^{\alpha }z\lambda e^{-\lambda \tau }+{\epsilon }^{\alpha }{k}_1^{\alpha }z{e}^{-\lambda \tau }\right)\left(\lambda +k_4^{\alpha }-k_3^{\alpha }y{e}^{-\lambda \tau }\right)\\ +b^{\alpha }x\left(-k^{\alpha }+k_2^{\alpha }z{e}^{-\lambda \tau }\right)\left(\lambda +k_4^{\alpha }-k_3^{\alpha }y{e}^{-\lambda \tau }\right)+k_1^{\alpha }{k}_3^{\alpha }yz\left(\lambda +{\epsilon }^{\alpha }\right)e^{-2\lambda \tau }\\ +b^{\alpha }{k}_2^{\alpha }xy{e}^{-\lambda \tau }\end{array}\right]\\ +b^{\alpha }v\left[\frac{r\alpha x}{T_{max}^{\alpha }}\left(\lambda +{\epsilon }^{\alpha }\right)\left(\lambda +k_4^{\alpha }-k_3^{\alpha }y{e}^{-\lambda \tau }\right)-b^{\alpha }x\left(\begin{array}{c}\left(\lambda +k_4^{\alpha }-k_3^{\alpha }y{e}^{-\lambda \tau }\right)\left(k_2^{\alpha }z{e}^{-\lambda \tau }-k^{\alpha }\right)\\ +k_2^{\alpha }{k}_3^{\alpha }yz{e}^{-2\lambda \tau }\end{array}\right)\right]\end{array}\right\}\end{array}$

where $A_0=d^{\alpha }-r^{\alpha }+\frac{2r^{\alpha }x+y}{T_{max}^{\alpha }}$.

For convenience, we made the following simplifications:

$\begin{array}{l}\begin{array}{c}\left({\lambda }^2+\left(a^{\alpha }+{\epsilon }^{\alpha }\right)\lambda +a^{\alpha }{\epsilon }^{\alpha }+k_1^{\alpha }z\lambda e^{-\lambda \tau }+{\epsilon }^{\alpha }{k}_1^{\alpha }z{e}^{-\lambda \tau }\right)\left(\lambda +k_4^{\alpha }-k_3^{\alpha }y{e}^{-\lambda \tau }\right)\hfill \\ ={\lambda }^3+A_1{\lambda }^2+A_2\lambda +A_3+A_4{\lambda }^2{e}^{-\lambda \tau }+A_5\lambda e^{-\lambda \tau }+A_6{e}^{-\lambda \tau }+A_7\lambda e^{-2\lambda \tau }+A_8{e}^{-2\lambda \tau }\end{array}\hfill \end{array}$

where

$\begin{array}{l}\begin{array}{c}{A}_1=a^{\alpha }+{\epsilon }^{\alpha }+k_4^{\alpha },A_2=a^{\alpha }{\epsilon }^{\alpha }+\left(a^{\alpha }+{\epsilon }^{\alpha }\right)k_4^{\alpha },A_3=a^{\alpha }{\epsilon }^{\alpha }{k}_4^{\alpha },A_4=\left(k_1^{\alpha }z-k_3^{\alpha }y\right),\hfill \\ A_5={\epsilon }^{\alpha }{k}_1^{\alpha }z+k_1^{\alpha }{k}_4^{\alpha }z-\left(a^{\alpha }+{\epsilon }^{\alpha }\right)k_3^{\alpha }y,A_6={\epsilon }^{\alpha }{k}_1^{\alpha }{k}_4^{\alpha }z-a^{\alpha }{\epsilon }^{\alpha }{k}_3^{\alpha }y,A_7=-k_1^{\alpha }{k}_3^{\alpha }yz,\hfill \\ A_8=-{\epsilon }^{\alpha }{k}_1^{\alpha }{k}_3^{\alpha }yz.\hfill \end{array}\hfill \end{array}$

$\begin{array}{l}\begin{array}{c}{b}^{\alpha }x\left(\lambda +k_4^{\alpha }-k_3^{\alpha }y{e}^{-\lambda \tau }\right)\left(-k^{\alpha }+k_2^{\alpha }z{e}^{-\lambda \tau }\right)+k_1^{\alpha }{k}_3^{\alpha }yz{\left(\lambda +\epsilon \right)}^{\alpha }{e}^{-2\lambda \tau }+b^{\alpha }{k}_2^{\alpha }xy{e}^{-\lambda \tau }\hfill \\ =A_9\lambda +A_{10}+A_{11}\lambda e^{-\lambda \tau }+A_{12}{e}^{-\lambda \tau }+A_{13}\lambda e^{-2\lambda \tau }+A_{14}{e}^{-2\lambda \tau }\hfill \end{array}\hfill \end{array}$

where

$\begin{array}{l}\begin{array}{c}{A}_9=-b^{\alpha }{k}^{\alpha }x,A_{10}=-b^{\alpha }{k}^{\alpha }{k}_4^{\alpha }x,A_{11}=b^{\alpha }{k}_2^{\alpha }xz,A_{12}=b^{\alpha }{k}^{\alpha }{k}_3^{\alpha }xy+b^{\alpha }{k}_2^{\alpha }{k}_4^{\alpha }xz+b^{\alpha }{k}_2^{\alpha }xy,\hfill \\ A_{13}=k_1^{\alpha }{k}_3^{\alpha }yz,A_{14}={\epsilon }^{\alpha }{k}_1^{\alpha }{k}_3^{\alpha }yz-b{k}_2^{\alpha }{k}_3^{\alpha }xyz\hfill \end{array}\hfill \end{array}$

$\begin{array}{l}\begin{array}{c}{b}^{\alpha }v\left[\frac{r^{\alpha }x}{T_{max}^{\alpha }}\left(\lambda +{\epsilon }^{\alpha }\right)\left(\lambda +k_4^{\alpha }-k_3^{\alpha }y{e}^{-\lambda \tau }\right)-b^{\alpha }x\left(\begin{array}{c}\left(\lambda +k_4^{\alpha }-k_3^{\alpha }y{e}^{-\lambda \tau }\right)\left(k_2^{\alpha }z{e}^{-\lambda \tau }-k^{\alpha }\right)\hfill \\ +k_2^{\alpha }{k}_3^{\alpha }yz{e}^{-2\lambda \tau }\hfill \end{array}\right)\right]\hfill \\ =A_{15}{\lambda }^2+A_{16}\lambda +A_{17}+A_{18}\lambda e^{-\lambda \tau }+A_{19}{e}^{-\lambda \tau }+A_{20}{e}^{-2\lambda \tau }\hfill \end{array}\hfill \end{array}$

where

$A_{15}=\frac{b^{\alpha }{r}^{\alpha }xv}{T_{max}^{\alpha }},A_{16}=b^{\alpha }v\left(\frac{r^{\alpha }{k}_4^{\alpha }x}{T_{max}^{\alpha }}+\frac{{\epsilon }^{\alpha }{r}^{\alpha }x}{T_{max}}+b^{\alpha }{k}^{\alpha }x\right),A_{17}=b^{\alpha }v\left(\frac{{\epsilon }^{\alpha }{r}^{\alpha }{k}_4^{\alpha }x}{T_{max}^{\alpha }}+b^{\alpha }{k}^{\alpha }{k}_4^{\alpha }x\right),$

$\begin{array}{c}{A}_{18}=-b^{\alpha }v\left(\frac{r^{\alpha }{k}_3^{\alpha }xy}{T_{max}^{\alpha }}+b^{\alpha }{k}_2^{\alpha }xz\right),A_{19}=-b^{\alpha }v\left(\frac{{\epsilon }^{\alpha }{r}^{\alpha }{k}_3^{\alpha }xy}{T_{max}^{\alpha }}+b^{\alpha }{k}_2^{\alpha }{k}_4^{\alpha }xz+b^{\alpha }{k}^{\alpha }{k}_3^{\alpha }xy\right),\hfill \\ A_{20}=-b^{2\alpha }xv\left(k_2^{\alpha }{k}_3^{\alpha }yz-k_2^{\alpha }{k}_3^{\alpha }yz\right).\hfill \end{array}$

Thus, the characteristic determinant becomes as follows:

$\begin{array}{l}\left(\lambda +k_7\right)\left\{\begin{array}{c}\left(\lambda +A_0\right)\left[\begin{array}{c}{\lambda }^3+A_1{\lambda }^2+A_2\lambda +A_3+A_4{\lambda }^2{e}^{-\lambda \tau }+A_5\lambda e^{-\lambda \tau }+A_6{e}^{-\lambda \tau }+A_7\lambda e^{-2\lambda \tau }\hfill \\ +A_8{e}^{-2\lambda \tau }+A_9\lambda +A_{10}+A_{11}\lambda e^{-\lambda \tau }+A_{12}{e}^{-\lambda \tau }+A_{13}\lambda e^{-2\lambda \tau }+A_{14}{e}^{-2\lambda \tau }\hfill \end{array}\right]\hfill \\ +A_{15}{\lambda }^2+A_{16}\lambda +A_{17}+A_{18}\lambda e^{-\lambda \tau }+A_{19}{e}^{-\lambda \tau }+A_{20}{e}^{-2\lambda \tau }\hfill \end{array}\right\}\hfill \end{array}$

$\begin{array}{l}=\left(\lambda +k_7^{\alpha }\right)\left\{\left[\begin{array}{c}{\lambda }^4+B_1{\lambda }^3+B_2{\lambda }^2+B_3\lambda +B_4+B_5{\lambda }^3{e}^{-\lambda \tau }+B_6{\lambda }^2{e}^{-\lambda \tau }+B_7\lambda e^{-\lambda \tau }+B_8{e}^{-\lambda \tau }\hfill \\ +B_9{\lambda }^2{e}^{-2\lambda \tau }+B_{10}\lambda e^{-2\lambda \tau }+B_{11}{e}^{-2\lambda \tau }\hfill \end{array}\right]\right\}\hfill \end{array}$

where

$\begin{array}{l}\begin{array}{c}{B}_1=A_1+A_0,B_2=A_2+A_9+A_0{A}_1+A_{15},B_3=A_3+A_{10}+A_0{A}_2+A_0{A}_9+A_{16},\hfill \\ B_4=A_0{A}_3+A_0{A}_{10}+A_{17},B_5=A_4,B_6=A_5+A_{11}+A_0{A}_4,B_7=A_6+A_{12}+A_0{A}_5+A_0{A}_{11}+A_{18},\hfill \\ B_8=A_0{A}_6+A_0{A}_{12}+A_{19},B_9=A_7+A_{13},B_{10}=A_8+A_{14}+A_0{A}_7+A_0{A}_{13},B_{11}=A_0{A}_8+A_0{A}_{14}+A_{20}.\hfill \end{array}\hfill \end{array}$

For further simplification, we derived the following assignment:

$\begin{array}{l}\begin{array}{c}{C}_1=B_1\text{+}{k}_7,C_2=B_2+k_7{B}_1,C_3=B_3+k_7{B}_2,C_4=B_4+k_7{B}_3,C_5=k_7,C_6=B_5,C_7=B_6+k_7{B}_5,\hfill \\ C_8=B_7+k_7{B}_6,C_9=B_8+k_7{B}_7,C_{10}=k_7{B}_8,C_{11}=B_9,C_{12}=B_{10}+k_7{B}_9,C_{13}=B_{11}+k_7{B}_{10},C_{14}=k_7{B}_{11}.\hfill \end{array}\hfill \end{array}$

The characteristic determinant is as follows:

$\begin{array}{ll}\begin{array}{c}H\left(\lambda ;\tau \right)\text{=}{\lambda }^5+C_1{\lambda }^4+C_2{\lambda }^3+C_3{\lambda }^2+C_4\lambda +C_5+\left(C_6{\lambda }^4+C_7{\lambda }^3+C_8{\lambda }^2+C_9\lambda +C_{10}\right)e^{-\lambda \tau }\hfill \\ +\left(C_{11}{\lambda }^3+C_{12}{\lambda }^2+C_{13}\lambda +C_{14}\right)e^{-2\lambda \tau }\text{=0}\hfill \end{array}\hfill & \left(8\right)\end{array}$

When τ = 0, the previous equation becomes as follows:

$\begin{array}{ll}{\lambda }^5+D_1{\lambda }^4+D_2{\lambda }^3+D_3{\lambda }^2+D_4\lambda +D_5=0\hfill & \left(9\right)\end{array}$

where

$\begin{array}{l}{D}_1=C_1+C_6,D_2=C_2+C_7+C_{11},D_3=C_3+C_8+C_{12},D_4=C_4+C_9+C_{13},D_5=C_5+C_{10}+C_{14}.\hfill \end{array}$

Based on equation (9), we get the following lemma by applying the Routh–Hurwitz criterion.

Lemma If equation (9) satisfies Δ₁ ≡ D₁ > 0, ${\Delta }_2\equiv \left|\begin{array}{cc}{D}_1\hfill & 1\\ D_3\hfill & D_2\end{array}\right|>0$, and ${\Delta }_3\equiv \left|\begin{array}{ccc}{D}_1\hfill & 1& 0\\ D_3\hfill & D_2& D_1\\ D_5\hfill & D_4& D_3\end{array}\right|>0$, E₁ is locally asymptotically stable when τ = 0.

Proof. The detailed proof can be referred to Peter et al. (26), Ogunrinde et al. (27).

The aforementioned lemma indicated that when τ = 0, all roots of H(λ; τ) are to the left of the imaginary axis, and some roots may cross to the right from the imaginary axis as τ increases. Thus, E₁ is unstable because of its positive real parts.

Then, the stability of system (2) was investigated when τ > 0.

Both sides of equation (8) were multiplied by e^λτ:

$\begin{array}{ll}\begin{array}{c}\left(C_6{\lambda }^4+C_7{\lambda }^3+C_8{\lambda }^2+C_9\lambda +C_{10}\right)+\left({\lambda }^5+C_1{\lambda }^4+C_2{\lambda }^3+C_3{\lambda }^2+C_4\lambda +C_5\right)e^{\lambda \tau }\hfill \\ +\left(C_{11}{\lambda }^3+C_{12}{\lambda }^2+C_{13}\lambda +C_{14}\right)e^{-\lambda \tau }\text{=0}\hfill \end{array}\hfill & \left(10\right)\end{array}$

Suppose the aforementioned equation has a purely imaginary root λ = iω (ω > 0), then we have e^iω = cosω+isinω, e^−iω = cosω − isinω. Substituting λ = iω into equation (10), we have

$\begin{array}{ll}{C}_6{\lambda }^4+C_7{\lambda }^3+C_8{\lambda }^2+C_9\lambda +C_{10}=C_6{\omega }^4-C_7{\omega }^{3}i-C_8{\omega }^2+C_9\omega i+C_{10}\hfill & \left(11\right)\end{array}$

$\begin{array}{ll}\begin{array}{l}({\lambda }^5+C_1{\lambda }^4+C_2{\lambda }^3+C_3{\lambda }^2+C_4\lambda +C_5)e^{\lambda \tau }=(C_1{\omega }^4-C_3{\omega }^2+C_5)\cos \omega \tau +({\omega }^5-\hfill \\ C_2{\omega }^3+C_4\omega )\cos \omega \tau i+(-{\omega }^5+C_2{\omega }^3-C_4\omega )\sin \omega \tau +(C_1{\omega }^4-C_3{\omega }^2+C_5)\sin \omega \tau i\hfill \end{array}\hfill & \left(12\right)\end{array}$

$\begin{array}{ll}\begin{array}{c}\left(C_{11}{\lambda }^3+C_{12}{\lambda }^2+C_{13}\lambda +C_{14}\right)e^{-\lambda \tau }=\left(-C_{12}{\omega }^2+C_{14}\right)cos\omega \tau +\left(-C_{11}{\omega }^3+C_{13}\omega \right)\hfill \\ \times cos\omega \tau i+\left(-C_{11}{\omega }^3+C_{13}\omega \right)sin\omega \tau +\left(C_{12}{\omega }^2-C_{14}\right)sin\omega \tau i\hfill \end{array}\hfill & \left(13\right)\end{array}$

Therefore, equation (10) becomes as follows:

$\begin{array}{l}\begin{array}{c}\left(C_1{\omega }^4-C_3{\omega }^2-C_{12}{\omega }^2+C_5+C_{14}\right)cos\omega \tau +\left({\omega }^5-C_2{\omega }^3-C_{11}{\omega }^3+C_4\omega +C_{13}\omega \right)cos\omega \tau i\hfill \\ +\left(-{\omega }^5+C_2{\omega }^3-C_{11}{\omega }^3-C_4\omega +C_{13}\omega \right)sin\omega \tau +\left(C_1{\omega }^4-C_3{\omega }^2+C_{12}{\omega }^2+C_5-C_{14}\right)sin\omega \tau i\hfill \\ +C_6{\omega }^4-C_7{\omega }^{3}i-C_8{\omega }^2+C_9\omega i+C_{10}=0\hfill \end{array}\hfill \end{array}$

For convenience, we assumed the following:

$\begin{array}{l}\begin{array}{c}{a}_1=C_1{\omega }^4-C_3{\omega }^2-C_{12}{\omega }^2+C_5+C_{14},a_2={\omega }^5-C_2{\omega }^3-C_{11}{\omega }^3+C_4\omega +C_{13}\omega ,\hfill \\ a_3=-{\omega }^5+C_2{\omega }^3-C_{11}{\omega }^3-C_4\omega +C_{13}\omega ,a_4=C_1{\omega }^4-C_3{\omega }^2+C_{12}{\omega }^2+C_5-C_{14}.\hfill \end{array}\hfill \end{array}$

Thus, we get

$\begin{array}{ll}\begin{array}{c}{a}_{1}cos\omega \tau +a_{2}cos\omega \tau i+a_{3}sin\omega \tau +a_{4}sin\omega \tau i+C_6{\omega }^4-C_7{\omega }^{3}i-C_8{\omega }^2\hfill \\ +C_9\omega i+C_{10}=0\hfill \end{array}\hfill & \left(14\right)\end{array}$

The real part after separating the real and imaginary parts is as follows:

$\begin{array}{ll}{a}_{1}cos\omega \tau +a_{3}sin\omega \tau =-C_6{\omega }^4+C_8{\omega }^2-C_{10}=D_1\hfill & \left(15\right)\end{array}$

and the imaginary part is as follows:

$\begin{array}{ll}{a}_{2}cos\omega \tau +a_{4}sin\omega \tau =C_7{\omega }^3-C_9\omega =D_2\hfill & \left(16\right)\end{array}$

It follows from the real part and imaginary part that

$\begin{array}{ll}cos\omega \tau =\frac{a_4{D}_1-a_3{D}_2}{a_1{a}_4-a_2{a}_3};sin\omega \tau =\frac{a_1{D}_2-a_2{D}_1}{a_1{a}_4-a_2{a}_3}\hfill & \left(17\right)\end{array}$

Suppose equation (10) has ñ(1 ≤ ñ ≤ 5) positive real roots, denoted by x_n(1 ≤ n ≤ ñ).

Let $\sqrt{x_n}=\omega$, we get

$\begin{array}{l}\begin{array}{c}cos\left(\sqrt{x_n}\tau \right)=Q_n\hfill \\ =\frac{\left(C_1{\omega }^4-C_3{\omega }^2+C_{12}{\omega }^2+C_5-C_{14}\right)\left(-C_6{\omega }^4+C_8{\omega }^2-C_{10}\right)-\left(-{\omega }^5+C_2{\omega }^3-C_{11}{\omega }^3-C_4\omega +C_{13}\omega \right)\left(C_7{\omega }^3-C_9\omega \right)}{\left(C_1{\omega }^4-C_3{\omega }^2-C_{12}{\omega }^2+C_5+C_{14}\right)\left(C_1{\omega }^4-C_3{\omega }^2+C_{12}{\omega }^2+C_5-C_{14}\right)-\left({\omega }^5-C_2{\omega }^3-C_{11}{\omega }^3+C_4\omega +C_{13}\omega \right)\left(-{\omega }^5+C_2{\omega }^3-C_{11}{\omega }^3-C_4\omega +C_{13}\omega \right)}\hfill \end{array}\hfill \end{array}$

$\begin{array}{l}\begin{array}{c}sin\left(\sqrt{x_n}\tau \right)=P_n\hfill \\ =\frac{\left(C_1{\omega }^4-C_3{\omega }^2-C_{12}{\omega }^2+C_5+C_{14}\right)\left(C_7{\omega }^3-C_9\omega \right)-\left({\omega }^5-C_2{\omega }^3-C_{11}{\omega }^3+C_4\omega +C_{13}\omega \right)\left(-C_6{\omega }^4+C_8{\omega }^2-C_{10}\right)}{\left(C_1{\omega }^4-C_3{\omega }^2-C_{12}{\omega }^2+C_5+C_{14}\right)\left(C_1{\omega }^4-C_3{\omega }^2+C_{12}{\omega }^2+C_5-C_{14}\right)-\left({\omega }^5-C_2{\omega }^3-C_{11}{\omega }^3+C_4\omega +C_{13}\omega \right)\left(-{\omega }^5+C_2{\omega }^3-C_{11}{\omega }^3-C_4\omega +C_{13}\omega \right)}\hfill \end{array}\hfill \end{array}$

Let

$\begin{array}{l}{\tau }_n^{\left(j\right)}=\{\begin{array}{ll}\frac{1}{\sqrt{x_n}}\left[\arccos \left(Q_n\right)+2j\pi \right],\hfill & if{P}_n\ge 0\hfill \\ \frac{1}{\sqrt{x_n}}\left[2\pi -\arccos \left(Q_n\right)+2j\pi \right],\hfill & if{P}_n<0\hfill \end{array}\hfill \end{array}$

Here, the positive integer n satisfies 1 ≤ n ≤ ñ, j = 0, 1, 2, ...

Thus, from the aforementioned equation, we know that the characteristic equation has a pair of purely imaginary roots $\pm i\sqrt{x_n}$. We define ${\lambda }_n^{(j)}(\tau )={\alpha }_n^{(j)}(\tau )+i{\omega }_n^{(j)}(\tau )$ as the root of equation (10) near ${\tau }_n^{(j)}$ for every 1 ≤ n ≤ ñ and j, satisfying ${\alpha }_n^{(j)}({\tau }_n^{(j)})=0$ and ${\omega }_n^{(j)}({\tau }_n^{(j)})=\sqrt{x_n}$. In summary, we arrived at the following theorem:

Theorem 1 When $\tau \in [0,{\tau }_{n_0}^{(0)})$ and there are positive real roots in equation (10), infection equilibrium E₁ is locally asymptotically stable, where

$\begin{array}{l}{\tau }_{n_0}^{\left(0\right)}=min\left\{{\displaystyle {\tau }_n^{\left(j\right)}}|1\le n\le ñ,j=0,1,2,...\right\}.\hfill \end{array}$

Proof. When $\tau \in [0,{\tau }_{n_0}^{(0)})$ and equation (10) have no positive real roots, where ${\tau }_{n_0}^{(0)}=min\left\{{\tau }_n^{(j)}|1\le n\le ñ,j=0,1,2,...\right\}$, all the roots have strictly negative real parts. Thus, E₁ is locally asymptotically stable for $\tau \in [0,{\tau }_{n_0}^{(0)})$.

Sensitivity of the basic reproduction number

In this part, the sensitivity index of the basic reproduction number is explored in order to find out the most sensitive parameter that can significantly affect the basic reproduction number and give proper treatment strategies (3).

The sensitivity index can be computed by using the following equation:

$\begin{array}{ll}{K}_q^{R_0}=\frac{\partial R_0}{\partial q}\times \frac{q}{R_0}\hfill & \left(18\right)\end{array}$

The basic reproduction number of E₀ is as follows:

$\begin{array}{l}{R}_0=\frac{b^{\alpha }{k}^{\alpha }{x}_0}{a^{\alpha }{\epsilon }^{\alpha }},where{x}_0=\frac{T_{max}^{\alpha }}{2r^{\alpha }}\left[-\left(d^{\alpha }-r^{\alpha }\right)+\sqrt{{\left(d^{\alpha }-r^{\alpha }\right)}^2+\frac{4{\xi }^{\alpha }{r}^{\alpha }}{T_{max}^{\alpha }}}\right].\hfill \end{array}$

The results of sensitivity indexes (Table 1) demonstrated that the infection rate of uninfected cells to become infected cells (b), production rate of free viruses (k), maximum hepatocyte counts in the liver (T_max), and production rate of uninfected cells (ξ) have the highest positive index. Therefore, decreasing the infection rate, the production rate of free viruses, and the production rate of uninfected cells can help treat patients with hepatitis B. On the contrary, the death rate of infected cells (a), the death rate of free viruses (ε), and the death rate of uninfected cells (d) have the highest negative index. This also suggests that increasing the death rate of infected cells, the death rate of free viruses, and the death rate of uninfected cells can also keep R₀ < 1 and help the treatment of patients with hepatitis B.

TABLE 1

Table 1. Sensitivity indexes of R₀ to model parameters.

Numerical simulation

In this section, a simulation is carried out to prove the accuracy of the aforementioned theoretical analysis.

Algorithm

Before the simulation, first, we provide the algorithm to solve the fractional-order differential equation (35, 36):

$\begin{array}{ll}\{\begin{array}{l}x(t_k)=\left[{\xi }^{\alpha }-d^{\alpha }x(t_{k-1})+r^{\alpha }x(t_{k-1})\left(1-\frac{x(t_{k-1})+y(t_{k-1})}{T_{{\max }^{\alpha }}}\right)-b^{\alpha }x(t_{k-1})v(t_{k-1})\right]h^{q_1}-{\displaystyle \sum _{j=v}^k{c}_j^{(q_1)}x(t_{k-j})},\hfill \\ y(t_k)=\left[b^{\alpha }x(t_{k-1})v(t_{k-1})-ay(t_{k-1})-k_1^{\alpha }y(t_{k-m-1})z(t_{k-m-1})\right]h^{q_1}-{\displaystyle \sum _{j=v}^k{c}_j^{(q_1)}y(t_{k-j})},\hfill \\ v(t_k)=\left[k^{\alpha }y(t_{k-1})-{\epsilon }^{\alpha }v(t_{k-1})-k_2^{\alpha }y(t_{k-m-1})z(t_{k-m-1})\right]h^{q_1}-{\displaystyle \sum _{j=v}^k{c}_j^{(q_1)}v(t_{k-j})},\hfill \\ z(t_k)=\left[k_3^{\alpha }y(t_{k-m-1})z(t_{k-m-1})-k_4^{\alpha }z(t_{k-1})\right]h^{q_1}-{\displaystyle \sum _{j=v}^k{c}_j^{(q_1)}z(t_{k-j})},\hfill \\ w(t_k)=\left[k_5^{\alpha }+k_6^{\alpha }y(t_{k-1})z(t_{k-1})-k_7^{\alpha }w(t_{k-1})\right]h^{q_1}-{\displaystyle \sum _{j=v}^k{c}_j^{(q_1)}w(t_{k-j})},\hfill \end{array}\hfill & (19)\end{array}$

where T_sim is time length, k = 1, 2, 3, ..., N, N = [T_sim/h], m = [τ/h], and x(0) = x₀, v(0) = v₀, w(0) = w₀, y(t) = y₀, z(t) = z₀, t ∈ [−τ, 0] are the initial conditions. $c_0^{(q)}=1,c_j^{(q)}=\left(1-\frac{1+q}{j}\right)c_{j-1}^{(q)}$.

Simulation of asymptotically stable infection-free equilibrium

First, we simulate the case of infection-free. The parameters are shown in Table 2.

TABLE 2

Table 2. Description of parameters and values when R₀ < 1.

The time length is 400, and the initial conditions are x(0) = 1, v(0) = 1, w(0) = 1, y(t) = 1, z(t) = 1, t ∈ [−τ, 0]. The order α = 0.9 and the time delay τ = 0.7. Therefore, we have E₀ = (x₀, y₀, v₀, z₀, w₀) = (2.0289, 0, 0, 0, 1.4372), and the basic reproduction number R₀ = 0.7546.

The behaviors of the uninfected cells (x), infected cells (y), free viruses (v), CTLs (z), and ALT (w) are shown in Figure 1. In Figure 1, all individuals converge to the infection-free equilibrium E₀, and the basic reproduction number R₀ is 0.7546, which is smaller than 1. This coincides with our theoretical analysis, which showed that when R₀ < 1, the infection-free equilibrium E₀ is asymptotically stable.

FIGURE 1

Figure 1. Dynamic change trend of uninfected cells (x), infected cells (y), free viruses (v), CTLs (z), and ALT (w) with R₀ < 1 and τ = 0.7.

Simulation of asymptotically stable infection equilibrium

The theoretical analysis of the infection equilibrium is verified in this section. Similarly, the parameters are shown in Table 3.

TABLE 3

Table 3. Description of parameters and values when R₁ > 1.

The initial conditions are the same as in the previous section. The time length is 400. The order α = 0.9, and the time delay τ = 1.2. Therefore, the infection equilibrium is as follows: E₁ = (x₁, y₁, v₁, z₁, w₁) = (1.2762, 0.3341, 0.5338, 1.0787, 1.4802), and the basic reproduction number is R₁ = 2.0132.

Figure 2 is the behavior of the uninfected cells (x), infected cells (y), free viruses (v), CTLs (z), and ALT (w) with R₁ > 1 and τ = 1.2. From Figure 2, we know that although all the individuals oscillate at the beginning, they converge to infection equilibrium E₁shortly. Figure 3 shows the phase portraits of the uninfected cell–infected cell–free virus space; the arrow indicates the direction of convergence of the phase portraits, and it converges to the infection equilibrium E₁ (red dot).

FIGURE 2

Figure 2. Dynamic change trend of uninfected cells (x), infected cells (y), free viruses (v), CTLs (z), ALT (w), and phase portraits of the xyz-space with R₁ > 1 and τ = 1.2.

FIGURE 3

Figure 3. Phase portraits of the xyz-space with R₁ > 1 and τ = 1.2.

Simulation of the Hopf bifurcation of the infection equilibrium

In this subsection, the Hopf bifurcation of the infection equilibrium is simulated. All parameters are the same as those in Section Simulation of asymptotically stable infection equilibrium, except τ = 3.2.

Figure 4 shows that when τ = 3.2, the uninfected cells (x), infected cells (y), free viruses (v), CTLs (z), and ALT (w) oscillate periodically around the infection equilibrium E₁. Figure 5 shows the phase portraits of the uninfected cell–infected cell–free virus space, and when τ = 3.2, the phase portraits are a stable limit cycle which is around the infection equilibrium E₁. The bifurcation diagram (Figure 6) shows that the stability of infection equilibrium E₁ changes at τ = 1.2.

FIGURE 4

Figure 4. Dynamic change trend of uninfected cells (x), infected cells (y), free viruses (v), CTLs (z), ALT (w), and phase portraits of the xyz-space with τ = 3.2.

FIGURE 5

Figure 5. Phase portraits of the xyz-space with τ = 3.2.

FIGURE 6

Figure 6. One parameter bifurcation diagram with respect to τ.

Simulation of phase portraits with different orders

In this section, the phase portraits with different orders are studied by using the method of numerical simulation. The initial order is α = 0.75, and with step of 0.05, the order increases to 0.99. The phase portraits also used the uninfected cell–infected cell–free virus space. As shown in Figure 7, when τ = 1.2 and the order (α) increases from 0.75 to 0.99, the volume of the phase portraits becomes bigger and the phase portraits become more complicated. Furthermore, the numerical simulations indicated that when the order increases from 0.75 to 0.95, the uninfected cell–infected cell–free virus space converges to the infection equilibrium E₁. However, when α = 0.99, the phase portrait is a stable limit cycle, which is around the infection equilibrium E₁. This indicated that the order can significantly affect the stability of the system.

FIGURE 7

Figure 7. Phase portraits of the xyz-space with τ = 1.2 and α = 0.75, α = 0.80, α = 0.85, α = 0.90, and α = 0.99.

Conclusion and discussion

In this study, a fractional differential model of HBV infection with time delay and logistic proliferation was proposed in order to better understand the infection mechanism and realize the infection progression. First, the infection-free equilibrium, infection equilibrium, and the basic reproduction number were computed. In epidemiology, R₀ is considered the most important parameter, which provides an insight into how the disease spreads and helps us understand how to control the disease. Therefore, we proved that if the basic reproduction number $R_0=\frac{b^{\alpha }{k}^{\alpha }{x}_0}{a^{\alpha }{\epsilon }^{\alpha }}<1$, the infection-free equilibrium (E₀) is locally asymptotically stable, which indicated that if the basic reproduction number R₀ < 1 can be controlled in patients, hepatitis B will disappear. Similarly, the stability analysis of the infection-free equilibrium (E₁) was discussed. In addition, the Hopf bifurcation of the infection equilibrium was studied at the theoretical level. Furthermore, sensitivity was analyzed to screen out the parameters that can significantly affect the basic reproduction number in our model. The results indicated that decreasing the infection rate (b), production rate of free viruses (k) and production rate of uninfected cells (ξ) can significantly decrease the basic reproduction number (R₀). Similarly, increasing the death rate of infected cells (a), the death rate of free viruses (ε) and he death rate of uninfected cells (d) can also decrease the basic reproduction number (R₀). Therefore, in order to keep R₀ < 1, the patient can decrease parameters b, k, and ξ or increase a, ε, and d to achieve the purpose of treatment.

In order to verify the accuracy of the aforementioned theoretical analysis, the numerical simulations were carried out. The simulation results showed that when R₀ < 1 and τ < 1.2, the infection-free equilibrium E₀ is asymptotically stable, which indicates that the disease will disappear. When R₁ > 1 and τ < 1.2, the infection equilibrium E₁ is asymptotically stable, which indicates that the disease could be mitigated and will lead to a lower infectious class over a period. However, with the increase in τ, the uninfected cells, infected cells, free viruses, CTL levels, and ALT levels oscillate periodically around the infection equilibrium E₁, and the phase portrait is a stable limit cycle, which around the infection equilibrium E₁ indicate that the disease would be out of control. Furthermore, the simulations also indicated that the order can significantly affect the stability of the system. For example, if the order is in the range of 0.75–0.95, the phase portraits converge to the infection equilibrium E₁, and when α = 0.99, the phase portrait is a stable limit cycle.

Therefore, time delay and fractional order are necessary factors that should be considered in modeling HBV infection and for researching dynamic characteristics. Although the process of HBV infection is more complicated than is established in this study, we believe that the model and analysis can play an important role in improving the HBV treatment regimen.

Data availability statement

The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding authors.

Author contributions

DS: conceptualization, methodology, software, and writing—reviewing and editing. JL: software, methodology, and writing. XS: supervision, writing—review editing. GP: writing—review editing. All authors contributed to the article and approved the submitted version.

Funding

The work is supported by The National Natural Science Foundation of China (62103287 and 62003071), Basic Research General Project of Shenzhen (JCYJ20210324103209026), PhD Basic Research Initiation Project (RCBS20200714114856171), and Clinical Research Project of Shenzhen Second People's Hospital (20213357007).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

1. WHO. Hepatitis B. (2021). Available online at: https://www.who.int/news-room/fact-sheets/detail/hepatitis-b (accessed June 24, 2022).

Google Scholar

2. Li M, Zu J. The review of differential equation models of HBV infection dynamics. J Virol Methods. (2019) 266:103–13. doi: 10.1016/j.jviromet.2019.01.014

PubMed Abstract | CrossRef Full Text | Google Scholar

3. James Peter O, Ojo MM, Viriyapong R, Abiodun Oguntolu F. Mathematical model of measles transmission dynamics using real data from Nigeria. J Diff Equ Appl. (2022) 28:753–70. doi: 10.1080/10236198.2022.2079411

CrossRef Full Text | Google Scholar

4. Peter OJ, Qureshi S, Yusuf A, Al-Shomrani M, Idowu AA. A new mathematical model of COVID-19 using real data from Pakistan. Result Phys. (2021) 24:104098. doi: 10.1016/j.rinp.2021.104098

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Ojo MM, Peter OJ, Goufo EFD, Panigoro HS, Oguntolu FA. Mathematical model for control of tuberculosis epidemiology. J Appl Math Comput. (2022). doi: 10.1007/s12190-022-01734-x

CrossRef Full Text | Google Scholar

6. Elaiw AM, AlShamrani NH. Global stability of humoral immunity virus dynamics models with nonlinear infection rate and removal. Nonlinear Anal Real World Appl. (2015) 26:161–90. doi: 10.1016/j.nonrwa.2015.05.007

CrossRef Full Text | Google Scholar

7. Nowak MA, BS, Hill AM. Viral dynamics in hepatitis B virus infection. Proc Natl Acad Sci U S A. (1996) 93:4398–402. doi: 10.1073/pnas.93.9.4398

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Su Y, Sun D, Zhao L. Global analysis of a humoral and cellular immunity virus dynamics model with the Beddington–DeAngelis incidence rate. Math Methods Appl Sci. (2015) 38:2984–93. doi: 10.1002/mma.3274

CrossRef Full Text | Google Scholar

9. Manna K, Chakrabarty SP. Chronic hepatitis B infection and HBV DNA-containing capsids: modeling and analysis. Commun Nonlinear Sci Numer Simul. (2015) 22:383–95. doi: 10.1016/j.cnsns.2014.08.036

CrossRef Full Text | Google Scholar

10. Li J, Wang K, Yang Y. Dynamical behaviors of an HBV infection model with logistic hepatocyte growth. Math Comput Model. (2011) 54:704–11. doi: 10.1016/j.mcm.2011.03.013

CrossRef Full Text | Google Scholar

11. Prakash M, Balasubramaniam P. Bifurcation analysis of macrophages infection model with delayed immune response. Commun Nonlinear Sci Numer Simul. (2016) 35:1–16. doi: 10.1016/j.cnsns.2015.10.012

CrossRef Full Text | Google Scholar

12. Sun C, Cao Z, Lin Y. Analysis of stability and Hopf bifurcation for a viral infectious model with delay. Chaos, Soliton, and Fractals. (2007) 33:234–45. doi: 10.1016/j.chaos.2005.12.029

CrossRef Full Text | Google Scholar

13. Zhang Z, Kundu S, Tripathi JP, Bugalia S. Stability and Hopf bifurcation analysis of an SVEIR epidemic model with vaccination and multiple time delays. Chaos, Solitons, and Fractals. (2020) 131:109483. doi: 10.1016/j.chaos.2019.109483

CrossRef Full Text | Google Scholar

14. Babasola O, Oshinubi K, Peter OJ, Onwuegbuche FC, Khan I. Time-delayed modelling of the COVID-19 dynamics with a convex incidence rate. Res Sq. (2022). doi: 10.21203/rs.3.rs-1814397/v1

CrossRef Full Text | Google Scholar

15. Chinnathambi R, Rihan FA. Stability of fractional-order prey–predator system with time-delay and Monod–Haldane functional response. Nonlinear Dyn. (2018) 92:1637–48. doi: 10.1007/s11071-018-4151-z

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Sun D, Lu L, Liu F, Duan L, Wang D, Xiong J. Analysis of an improved fractional-order model of boundary formation in the Drosophila large intestine dependent on Delta–Notch pathway. Adv Differ Equ. (2020) 2020:377. doi: 10.1186/s13662-020-02836-1

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Mondal S, Lahiri A, Bairagi N. Analysis of a fractional order eco-epidemiological model with prey infection and type 2 functional response. Math Methods Appl Sci. (2017) 40:6776–89. doi: 10.1002/mma.4490

CrossRef Full Text | Google Scholar

18. Almeida R. Analysis of a fractional SEIR model with treatment. Appl Math Lett. (2018) 84:56–62. doi: 10.1016/j.aml.2018.04.015

CrossRef Full Text | Google Scholar

19. George Maria Selvam A, Abraham Vianny D. Bifurcation and dynamical behaviour of a fractional order Lorenz-Chen-Lu like chaotic system with discretization. J Phys Conf Ser. (2019) 1377:012002. doi: 10.1088/1742-6596/1377/1/012002

CrossRef Full Text | Google Scholar

20. Wang X, Wang Z, Huang X, Li Y. Dynamic analysis of a delayed fractional-order SIR model with saturated incidence and treatment functions. Int J Bifurc Chaos. (2019) 28:1850180. doi: 10.1142/S0218127418501808

CrossRef Full Text | Google Scholar

21. Li H, Huang C, Li T. Dynamic complexity of a fractional-order predator–prey system with double delays. Phys A Stat Mech Appl. (2019) 526:120852. doi: 10.1016/j.physa.2019.04.088

CrossRef Full Text | Google Scholar

22. Balci E, Öztürk I, Kartal S. Dynamical behaviour of fractional order tumor model with Caputo and conformable fractional derivative. Chaos, Solitons, and Fractals. (2019) 123:43–51. doi: 10.1016/j.chaos.2019.03.032

CrossRef Full Text | Google Scholar

23. Ye Y, Kou C. Stability analysis for a fractional HIV-1 model with time delay. In 4th International Conference on Biomedical Engineering and Informatics. IEEE, New York, NY (2011).

Google Scholar

24. Nwajeri UK, Fadugba D, Ohaeri EO, Oshinubi K, Ogunrinde RR, Ogunrinde PB. Co-Dynamic Model of Drug Trafficking and Money Laundering Coupled with Atangana-Baleanu Derivative. (2022). doi: 10.2139/ssrn.4165832

CrossRef Full Text | Google Scholar

25. Peter OJ. Transmission Dynamics of Fractional Order Brucellosis Model Using Caputo–Fabrizio Operator. Int. J. Diff. Equ. (2020) 2020:2791380. doi: 10.1155/2020/2791380

CrossRef Full Text | Google Scholar

26. Peter OJ, Shaikh AS, Ibrahim MO, Nisar KS, Baleanu D, Khan I, et al. Analysis and dynamics of fractional order mathematical model of COVID-19 in Nigeria using Atangana–Baleanu operator. Comput Mater Contin. (2021) 66:1823–48. doi: 10.32604/cmc.2020.012314

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Ogunrinde RB, Nwajeri UK, Fadugba SE, Ogunrinde RR, Oshinubi KI. Dynamic model of COVID-19 and citizens reaction using fractional derivative. Alex Eng J. (2021) 60:2001–12. doi: 10.1016/j.aej.2020.09.016

CrossRef Full Text | Google Scholar

28. Peter OJ, Yusuf A, Oshinubi K, Oguntolu FA, Lawal JO, Abioye AI, et al. Fractional order of pneumococcal pneumonia infection model with Caputo Fabrizio operator. Result Phys. (2021) 29:104581. doi: 10.1016/j.rinp.2021.104581

CrossRef Full Text | Google Scholar

29. Diethelm K, A. fractional calculus based model for the simulation of an outbreak of dengue fever. Nonlinear Dyn. (2012) 71:613–9. doi: 10.1007/s11071-012-0475-2

CrossRef Full Text | Google Scholar

30. González-Parra G, Arenas AJ, Chen-Charpentier BM. A fractional order epidemic model for the simulation of outbreaks of influenza A(H1N1). Math Methods Appl Sci. (2014) 37:2218–26. doi: 10.1002/mma.2968

CrossRef Full Text | Google Scholar

31. Danane J, Allali K, Hammouch Z. Mathematical analysis of a fractional differential model of HBV infection with antibody immune response. Chaos, Solitons, and Fractals. (2020) 136:109787. doi: 10.1016/j.chaos.2020.109787

CrossRef Full Text | Google Scholar

32. Gao F, Li X, Li W, Zhou X. Stability analysis of a fractional-order novel hepatitis B virus model with immune delay based on Caputo–Fabrizio derivative. Chaos, Solitons, and Fractals. (2021) 142:110436. doi: 10.1016/j.chaos.2020.110436

CrossRef Full Text | Google Scholar

33. Simelane SM, Dlamini PG, A. fractional order differential equation model for Hepatitis B virus with saturated incidence. Result Phys. (2021) 24:104114. doi: 10.1016/j.rinp.2021.104114

CrossRef Full Text | Google Scholar

34. Su Y, Sun D. Optimal control of anti-HBV treatment based on combination of traditional Chinese medicine and western medicine. Biomed Signal Process Control. (2015) 15:41–8. doi: 10.1016/j.bspc.2014.09.007

CrossRef Full Text | Google Scholar

35. Adebisi AF, Ojurongbe TA, Okunlola KA, Peter OJ. Application of Chebyshev polynomial basis function on the solution of Volterra integro-differential equations using Galerkin method. Math Comput Sci. (2021) 2:41–51. doi: 10.30511/mcs.2021.540133.1047

CrossRef Full Text | Google Scholar

36. Christie I, Adebayo T, Peter O, Folasade A. Numerical solution of two-dimensional Fredholm integro-differential equations by Chebyshev integral operational matrix method. J Appl Math Comput Mech. (2022) 21:29–40. doi: 10.17512/jamcm.2022.1.03

CrossRef Full Text | Google Scholar

37. Sun D, Liu F. Analysis of a new delayed HBV model with exposed state and immune response to infected cells and viruses. Biomed Res Int. (2017) 2017:7805675. doi: 10.1155/2017/7805675

PubMed Abstract | CrossRef Full Text | Google Scholar