Innovative approaches to fractional modeling: Aboodh transform for the Keller-Segel equation

1.   Introduction

1.1.   Motivation and literature review

One of the most important health issues in the world is the Hepatitis B virus (HBV). All over the world, about $500$ million chronic infections are caused by HBV or HCV ^[1]. The term "hepatitis" describes liver swelling. Alcohol, specific drugs, chemicals, or a viral infection are the main factors that contribute to liver inflammation. The use of contaminated razors and needles, blood transfusion, and saliva exchange can also cause infection in a susceptible person. The two main viral infections that affect liver cells are HCV and HBV. Although HBV and HCV have similar names, they are genetically and clinically distinct viruses. Both DNA viruses known as HBV and single strand RNA viruses known as HCV damage the liver cells. HBV can cause a short-term infection called an acute infection, when the virus is eradicated from the body by the immune system. Cirrhosis and hepatocellular carcinoma are two liver diseases that are ultimately caused by a chronic infection. Those who just have the acute disease still suffer severe symptoms, such as jaundice, high weakness, nausea, vomiting, and stomach pain, for up to a year. It is noted here that $1\%$ to $5\%$ of infected people are chronically infected; however, this percentage is significantly higher in newborns and young children. The virus causes liver cancer in around $25\%$ of chronic carriers. Chronic HBV infections frequently happen in our early years, and the virus survives in the body because significant immune responses are destroyed. Both acute and chronic infections caused by HBV have a high death rate. The immune system works to protect the body against threats from outside particles. By eliminating virus-infected cells, the cytotoxic T lymphocytes (CTLs) play a crucial role in antiviral defense. Cytotoxic T lymphocytes (CTLs) are a crucial component of the immune system responsible for recognizing and eliminating virus-infected cells and cancerous cells. They are a type of T lymphocyte, a subset of white blood cells that play a central role in adaptive immunity. The term "cytotoxic" refers to their ability to induce cell death, particularly in target cells marked for destruction. CTLs are equipped with a specialized receptor known as the T-cell receptor (TCR), which allows them to recognize specific antigens presented on the surface of infected or abnormal cells. When a CTL encounters a target cell displaying the matching antigen, it releases cytotoxic molecules such as perforin and granzymes, which penetrate the target cell's membrane and induce apoptosis, or programmed cell death. This process effectively eliminates the infected or abnormal cell, preventing the spread of pathogens and the development of cancer. Additionally, CTLs play a vital role in immune surveillance, continuously patrolling the body to identify and eliminate threats to maintain overall health and immunity. Moreover, CTLs and B-cells are main parts of the immune response. CTLs target and kill virus-infected cells, whereas B-cells create antibodies that neutralize the viruses, referred to as antibodies for immune response. A typical immune response to a virus must include antibodies, cytokines, natural killer cells, and T-cells ^[2,3]. Further, it is important here to mention that within-host dynamics between free-viruses, susceptible/infected cells, and immune response refers to the interactions and changes that occur within an infected organism at the cellular and molecular level. Here is a breakdown:

(i) Free-viruses: These are viruses that are circulating freely within the host organism. They can infect susceptible cells and replicate, leading to the spread of infection.

(ii) Susceptible/infected cells: Cells within the host organism that are vulnerable to viral infection are referred to as susceptible cells. When a virus enters a susceptible cell and begins to replicate, that cell becomes infected. The dynamics between susceptible and infected cells involve the process of viral replication, cell damage, and potential cell death.

(iii) Immune response: When the host organism detects the presence of viruses or infected cells, the immune system is activated. This response involves various components, including immune cells such as T- and B-cells, as well as signaling molecules like cytokines. The immune response aims to eliminate the virus and infected cells, thus controlling the infection.

Understanding the within-host dynamics between these components is crucial for developing strategies to combat viral infections, such as vaccines or antiviral therapies. It involves studying how viruses interact with host cells, evade immune detection, and how the immune system responds to limit viral spread and clear the infection. On the other hand, in recent years, many scientists have proposed and studied the dynamics of infectious disease models. For instance, Qesmi et al. ^[4] investigated backward bifurcation analysis of continuous-time Hepatitis B and C virus mathematical models. Li et al.^[5] explored the dynamics of a Hepatitis B deterministic model. Chen and Xu ^[6] explored bifurcation analysis of a viral infection model. For more intersecting results in this direction, we refer the reader to the work done by eminent researchers ^{[7,8,9,10,11,12,13]}. See also ^[14] for the analysis through the basic reproduction number of a general discrete-time model.

1.2.   Main contribution

Our main findings in this paper include:

● Model formulation of a discrete time HBV model.

● Examination of the local behavior at equilibrium states with a basic reproduction number (BRN) of the HBV model.

● A study of the rate of convergence of the discrete HBV model.

● Bifurcation analysis at equilibrium states.

● Examination of the chaos of the state feedback method.

● Numerical validation of theoretical results.

1.3.   Paper layout

The paper is structured as follows: In Section 2, we give the model formulation of a discrete time HBV model. In Section 3, we investigate the equilibrium states, linearized form, and BRN for a discrete HBV model. Whereas in Section 4, we examine local stability analysis at the equilibrium states and bifurcation sets. In Section 5, detailed bifurcation analysis at equilibrium states is investigated. In Sections 6 and 7, we study convergence rate and chaos control, respectively, whereas theoretical findings are numerically verified in Section 8. The paper's conclusion with future work is presented in Section 9.

2.   Model formulation

Motivated from the aforementioned studies, in this paper, we explore dynamical characteristics of a discrete Hepatitis B virus mathematical model. For this, we first reformulate the mathematical formulation of a continuous-time HBV model based on Figure 1, where $z$, $v$, $s,$ and $I$, respectively, denote the size of the CTL response or the quantity of virus-specific CTLs, free virus particles, uninfected cells, and infected cells. Furthermore, it is worthwhile to mention that the investigation consists of the simple dynamics of virus-host cell interaction and the kinetic interaction of CTLs with infected cells, as well as the impact of immune responses on viral burden and antigenic diversity. Viral organisms rely on their host cell for both survival and reproduction. Based on the amount of virus present, the tissues affected, and the length of the infection, the antiviral immune response's beneficial and harmful effects must be balanced. The virus is the immune system's reaction to the infection and may directly destroy the host cell. At the rate $\beta s v$, the cells $s$ are infected by the $v$, where the constant rate $\beta$ represents the efficiency of this process, including the rate of successful infection and probability. The cells $I$ produce $v$ at the rate $kI$, whereas the cells $I$ die at a rate $\alpha I$. $v$ is eradicated from the system at the rate $uv$, whereas $\alpha$, $\delta$, and $u$ denote death rates of cells $I$, $s$, and free virus $v$, respectively. Further, $\frac{1}{u}$ and $\frac{1}{\alpha}$ are the average life times of $v$ and $I$, respectively. $\vartheta$ and $\delta s$ are the constant rates at which cells $s$ are produced and then die, respectively, while $\frac{k}{\alpha}$ is the amount of virus produced from the $I$. The CTLs proliferation rate in response to antigens is represented by $cIz$. CTLs decrease at a rate of $bz$ in the absence of stimulus. CTLs kill $I$ at a rate of $pIz$. The CTLs responsiveness is indicated by the parameter $c$. The parameter $p$ determines how quickly CTLs eliminate $I$. Based on these presumptions, the continuous-time HBV mathematical model designated by a system of differential equations takes the following form ^[15]:

In the context of populations with non-overlapping generations, discrete-time models governed by maps are preferable to continuous ones. Discrete models can also produce effective computational results for numerical simulations. For instance, after discretization by the Euler-forward formula, HBV model (2.1) takes the following form:

After simplification, the required discrete HBV model (2.2) takes the form:

where $h > 0$ denotes the integral step size such that the quantities $1-h \delta$, $1-\alpha h$, $1- h u$, and $1- h \mu$ are positive and this statement underscores a crucial aspect of discretizing continuous models: the choice of step size $h$ and its impact on the validity of numerical solutions. Discretization involves approximating continuous systems into a finite set of points or intervals, a fundamental process in numerical simulations and computational methods. However, the discretization scheme must be carefully chosen to ensure the stability and accuracy of the numerical solution. The necessity for a small step size $h > 0$ arises from the preservation of certain properties of the continuous model during discretization. Specifically, the quantities $1-h \delta$, $1-\alpha h$, $1- h u$, and $1- h \mu$, which likely represent rates or constants in the model equations, must remain positive. This condition ensures that the numerical solution accurately reflects the behavior of the system without introducing instability or negative values. Choosing an appropriate step size is crucial because it directly impacts the fidelity of the numerical solution. A small enough step size ensures that the discretization accurately captures the dynamics of the system over time. If the step size is too large, it may lead to numerical instability or inaccuracies in the solution, compromising the reliability of the results. Moreover, the statement highlights the non-uniqueness of discretization schemes. Different discretization methods can yield varying results and behaviors of the numerical solution. Some schemes may be more suitable than others for preserving stability and positivity, depending on the specific characteristics of the model equations and the system being studied. In summary, ensuring that the quantities $1-h \delta$, $1-\alpha h$, $1- h u$, and $1- h \mu$ remain positive for small enough step sizes $h > 0$ is essential for maintaining the stability and validity of numerical solutions obtained through the discretization of continuous models. The choice of the discretization scheme and step size play a crucial role in accurately representing the dynamics of the system and obtaining reliable results in numerical simulations and computational analysis.

3.   Linearized form, equilibrium states, and BRN

We will explore equilibrium states of a discrete HBV model (2.3), linearized form, and BRN in this section.

Lemma 3.1. For the HBV model's equilibrium states, the following statements hold:

(i) For all $\vartheta, \delta, \beta, \alpha, p, k, u, c, \mu, h$, discrete HBV model (2.3) has disease-free equilibrium state (DFES) $\Phi_1 = \left(\frac{\vartheta}{\delta}, 0, 0, 0\right)$;

(ii) Discrete HBV model (2.3) has boundary equilibrium state (ES) $\Phi_2 = \left(\frac{u \alpha}{k \beta}, \frac{\vartheta \beta k-\delta u \alpha}{\alpha k \beta}, \frac{\vartheta \beta k-\delta u \alpha}{\alpha u \beta}, 0 \right)$ if $\beta > \frac{\alpha u \delta}{\vartheta k}$. Moreover, biologically, ES represents the state where the pathogens are present while CTLs are absents;

(iii) If $\beta > \frac{c \alpha u \delta}{c \vartheta k-k \mu \alpha}$ then the discrete HBV model (2.3) has an epidemic equilibrium state (EES) $\Phi_3 = \left(\frac{c u \vartheta}{c \delta u+ \beta k \mu}, \frac{\mu}{c}, \frac{k \mu}{c u}, \frac{c k \beta \vartheta-c \delta u \alpha-k \alpha \beta \mu}{p (c\delta u+k \beta \mu)} \right)$.

Proof. If HBV model (2.3) has an equilibrium state $\Phi = (s, I, v, z)$, then

It is noted that algebraic system (3.1) is satisfied obviously if $\Phi = (s, I, v, z) = \left(\frac{\vartheta}{\delta}, 0, 0, 0\right)$, and so the HBV model (2.3) has a disease-free equilibrium state (DFES) $\Phi_1 = \left(\frac{\vartheta}{\delta}, 0, 0, 0\right)$ for all model parameters $\vartheta$, $\delta$, $\beta$, $\alpha$, $p$, $k$, $u$, $c$, $\mu$, and $h$. Now, in order to find equilibrium state $\Phi_2$, we will solve the following system simultaneously:

From the last equation of (3.2), we get

On utilizing (3.3) into the second equation of system (3.2), we get

Now, on solving the third equation of (3.2) and (3.4) simultaneously, one gets

Using (3.5) into the first equation of system (3.2), we get

From the third equation of system (3.2) and (3.6), we get

Equations (3.3) and (3.5)–(3.7) imply that the HBV model (2.3) has ES $\Phi_2 = \left(\frac{u \alpha}{k \beta}, \frac{\vartheta \beta k-\delta u \alpha}{\alpha k \beta}, \frac{\vartheta \beta k-\delta u \alpha}{\alpha u \beta}, 0\right)$ if $\beta > \frac{\alpha u \delta}{\vartheta k}$. More importantly, it is noted here that if $\beta > \frac{\alpha u \delta}{\vartheta k}$, that is, $\frac{\beta\vartheta k}{\alpha u \delta} > 1$, then the HBV model (2.3) has ES $\Phi_2 = \left(\frac{u \alpha}{k \beta}, \frac{\vartheta \beta k-\delta u \alpha}{\alpha k \beta}, \frac{\vartheta \beta k-\delta u \alpha}{\alpha u \beta}, 0\right)$, and so we define $\mathfrak{R}_0 = :\frac{\beta\vartheta k}{\alpha u \delta} > 1$ as a BRN of the virus, which indicates the typical number of newly infected cells produced from one infected cell at the start of the infectious process ^[16]. On the other hand, we can say that the HBV model (2.3) has ES $\Phi_2 = \left(\frac{u \alpha}{k \beta}, \frac{\vartheta \beta k-\delta u \alpha}{\alpha k \beta}, \frac{\vartheta \beta k-\delta u \alpha}{\alpha u \beta}, 0\right)$ if $\mathfrak{R}_0 > 1$. Finally, for EES in $\mathbb{R}_+^4$, one needs to solve following algebraic system (3.2). From the last equation of (3.2), we get

By utilizing (3.8) in the third equation of system (3.2), we get

Using (3.8) and (3.9) in the first equation of system (3.2), we get

Using (3.8)–(3.10) in the second equation of system (3.2), we get

Equations (3.8)–(3.11) imply that if $\beta > \frac{\alpha c u \delta}{c \vartheta k-k \alpha \mu}$, then the HBV model (2.3) has EES $\Phi_3 = \left(\frac{c u \vartheta}{c \delta u+\beta k \mu}, \frac{\mu}{c}, \frac{k \mu}{c u}, \frac{c k \beta \vartheta-c \delta u\alpha-k \alpha \beta \mu}{p(c \delta u+k \beta \mu)}\right)$. □ Now a linearized form for the HBV model (2.3) at equilibrium state $\Phi$ is formulated. So, under the map $(e, f, g, h)\mapsto (s_{t+1}, I_{t+1}, v_{t+1}, z_{t+1})$, one has the following variational matrix $J|_{\Phi}$, which is evaluated at equilibrium state $\Phi$:

where

Remark 3.1. It is important here to explain that the continuous-time model, which is depicted in (2.1), can be modified by reconsidering the following key aspects and assumptions:

$(1)$ The proportionality of $k$ to the burst size: Parameter $k$ is proportional to the burst size, which is the number of virion released from an infected cell upon lysis. This indicates that $k$ represents the rate of virion production per infected cell.

$(2)$ Adsorption rate $\beta$: The parameter $\beta$ denotes the rate at which viruses attach to and infect susceptible cells.

$(3)$ Free viruses' lifetime and loss: In the original model, the free viruses $v$ do not account for losses due to adsorption by uninfected cells. In biological terms, when a virus particle adsorbs to a susceptible cell and successfully infects it, the virus is effectively removed from the pool of free viruses. To incorporate this loss, the equation for free viruses should include a term that represents the reduction of free viruses due to adsorption by the susceptible cells. The modified equation would be:

where $\beta sv$ represents the rate at which viruses are lost due to successful infection of susceptible cells.

Based on the key aspects and assumptions, which are explained above, the modified version of the continuous-time model (2.1) now becomes the following:

Additionally, basic reproduction number, $\mathfrak{R}_0$, in the context of within-host viral dynamics, is commonly interpreted as the average number of new infections caused by a single infected cell in a fully susceptible population. However, given that free viruses are the primary infecting agents, a more appropriate interpretation of $\mathfrak{R}_0$ would be "the expected number of virion produced by each virus during its lifetime". This considers the entire cycle from infection to the release of new virion and their subsequent interactions with susceptible cells. In conclusion, incorporating the loss of viruses due to adsorption by susceptible cells provides a more accurate representation of the within-host viral dynamics. This adjustment ensures a more biologically realistic and mathematically consistent framework for understanding the dynamics of viral infections within a host.

4.   Local stability analysis at equilibrium states and bifurcation sets

Here, we will study the local behavior along with the identification of bifurcation sets of the HBV model (2.3) at the equilibrium states using existing theory ^{[17,18,19,20,21,22,23]}. For DFES, (3.12) becomes

with

Now, at DFES, we give the local behavior as follows:

Theorem 4.1. DFES of the HBV model (2.3) is

(i) a locally asymptotically stable if

(ii) an unstable if

(iii) a saddle if

or

or

or

or

or

(iv) non-hyperbolic if

or

or

or

where in original parameters (4.13) and (4.14), respectively, it will be shown in the following form:

and

Proof. DFES of the HBV model (2.3) is asymptotically stable if $|\lambda_1| = |1-\delta h| < 1$, $|\lambda_2| = |1-h \mu| < 1$, and $|\lambda_{3, 4}| = \left|\frac{1}{2}\left(2-\alpha h-h u\pm \frac{h \sqrt{\delta \left(u-\alpha\right)^2+4 k \beta \vartheta}}{\sqrt{\delta}}\right)\right| < 1$, that is, $0 < \delta < \frac{2}{h}$, $0 < \mu < \frac{2}{h}$, and $\frac{4-2\alpha h-2h u+\alpha u h^2}{\alpha u h^2} < \mathfrak{R}_0 < 1$. Therefore, DFES is asymptotically stable if $0 < \delta < \frac{2}{h}$, $0 < \mu < \frac{2}{h}$, and $\frac{4-2\alpha h-2h u+\alpha u h^2}{\alpha u h^2} < \mathfrak{R}_0 < 1$. A similar calculation shows that DFES of the HBV model (2.3) is unstable if $\delta > \frac{2}{h}$, $\mu > \frac{2}{h}$, and $1 < \mathfrak{R}_0 < \frac{4-2\alpha h-2h u+\beta u h^2}{\alpha u h^2}$, a saddle if $0 < \delta < \frac{2}{h}$, $\mu > \frac{2}{h}$, and $1 < \mathfrak{R}_0 < \frac{4-2\alpha h-2h u+\alpha u h^2}{\alpha u h^2}$ or $\delta > \frac{2}{h}$, $0 < \mu < \frac{2}{h}$, and $1 < \mathfrak{R}_0 < \frac{4-2\alpha h-2h u+\alpha u h^2}{\alpha u h^2}$ or $\delta > \frac{2}{h}$, $0 < \mu < \frac{2}{h}$, and $\frac{4-2\alpha h-2h u+\alpha u h^2}{\alpha u h^2} < \mathfrak{R}_0 < 1$ or $0 < \delta < \frac{2}{h}$, $\mu > \frac{2}{h}$, and $\frac{4-2\alpha h-2h u+\alpha u h^2}{\alpha u h^2} < \mathfrak{R}_0 < 1$ or $\delta > \frac{2}{h}$, $\mu > \frac{2}{h}$, and $\frac{4-2\alpha h-2h u+\alpha u h^2}{\alpha u h^2} < \mathfrak{R}_0 < 1$ or $0 < \delta < \frac{2}{h}$, $0 < \mu < \frac{2}{h}$, and $1 < \mathfrak{R}_0 < \frac{4-2\alpha h-2h u+\alpha u h^2}{\alpha u h^2}$, and finally, non-hyperbolic if $\delta = \frac{2}{h}$ or $\delta = \frac{2}{h}$ or $\mathfrak{R}_0 = 1$ or $\mathfrak{R}_0 = \frac{4-2\alpha h-2h u+\alpha u h^2}{\alpha u h^2}$. □

Hereafter, based on Theorem 4.1, we will discuss the existence of possible bifurcation sets. It is noted that if parameters $\chi = (\vartheta, \delta, \beta, \alpha, p, k, u, c, \mu, h)$ pass certain regions, then a number of possible bifurcation sets may occur, which are summarized as a following theorem:

Theorem 4.2. We have the following bifurcation sets if parameter $\chi = (\vartheta, \delta, \beta, \alpha, p, k, u, c, \mu, h)$ passes certain regions:

(i) If (4.11) holds, then $\lambda_1|_{(4.11)} = -1$ but $\lambda_{3, 4}|_{(4.11)} = \frac{2-h \alpha-h u\pm\sqrt{h^2((u-\alpha)^2+2 h k \beta \vartheta)}}{2}\ne 1\ \mathrm {or }\ -1$, as well as $\lambda_2|_{(4.11)} = 1-h \mu \ne 1\ \mathrm {or }\ -1$. So, for DFES, one has the following flip bifurcation set associated with Theorem 4.1:

(ii) If (4.12) holds, then $\lambda_2|_{(4.12)} = -1$ but $\lambda_{3, 4}|_{(4.12)} = \frac{1}{2}\left(2-\alpha h-h u\pm \frac{h \sqrt{\delta \left(u-\alpha\right)^2+4 k \beta \vartheta}}{\sqrt{\delta}}\right)\ne 1\ \mathrm {or }\ -1$, as well as $\lambda_{1}|_{(4.12)} = 1-h d \ne 1\ \mathrm {or }\ -1$. So, for DFES, one has the following flip bifurcation set associated with Theorem 4.1:

(iii) If (4.15) holds, then $\lambda_3|_{(4.15)} = 1$ but $\lambda_{1, 4}|_{(4.15)} = 1-h \left(\frac{k \beta\vartheta}{\alpha u}\right), 1-h\left(\alpha+u\right)\ne 1\ \mathrm {or }\ -1$, as well as $\lambda_2|_{(4.15)} = 1-h \mu \ne 1\ \mathrm {or }\ -1$. So, for DFES, one has the following fold bifurcation set associated with Theorem 4.1:

(iv) If (4.16) holds, then $\lambda_4|_{(4.16)} = -1$ but $\lambda_{1, 3}|_{(4.16)} = 1-h \left(\frac{k\beta\vartheta h^2}{4-2h u-2\alpha h+\alpha h^2 u}\right), 3-h \left(u+\alpha\right)\ne 1\ \mathrm {or }\ -1$, as well as $\lambda_2|_{(4.16)} = 1-h \mu \ne 1\ \mathrm {or }\ -1$. So, for DFES, one has the following flip bifurcation set associated with Theorem 4.1:

Hereafter, we will study the local behavior at ES and EES of the HBV model (2.3) by Theorem 1.4 of ^[18]. At ES, (3.12) gives

with

where

Theorem 4.3. If

then the ES of the HBV model (2.3) is a sink where $\mathcal{S}_i(i = 1, 2, 3, 4)$ are depicted in (4.23).

Finally, at EES, (3.12) gives

with

where

Theorem 4.4. If

then the EES of the HBV model (2.3) is a sink where $\mathcal{U}_i(i = 1, 2, 3, 4)$ are depicted in (4.27).

It is noted that by explicit criterion, we explore the existence of bifurcation at the ES and EES of HBV model (2.3).

5.   Bifurcation analysis

We will examine bifurcations at DFES, ES, and EES of the HBV model (2.3) in this section by bifurcation theory ^{[23,24,25,26,27,28,29]}.

Theorem 5.1. At DFES, the HBV model (2.3) does not undergo flip bifurcation if (4.17) holds.

Proof. For DFES, HBV model (2.3) is invariant subject to $I = v = z = 0$. Therefore, the HBV model (2.3) restricted to $I = v = z = 0$ becomes

From (5.1), we define

Since $\delta = \delta^* = \frac{2}{h}, \ s = s^* = \frac{\vartheta}{\delta}$, and therefore, from (5.2), one has that $f_s|_{\delta = \delta^* = \frac{2}{h}, \ s = s^* = \frac{\vartheta}{\delta}} = -1$ and $f_\delta|_{\delta = \delta^* = \frac{2}{h}, \ s = s^* = \frac{\vartheta}{\delta}} = -\frac{h^2\vartheta}{2}\neq0$ hold but $f_{ss}|_{\delta = \delta^* = \frac{2}{h}, \ s = s^* = \frac{\vartheta}{\delta}} = 0$ violates the non-degenerate condition for the occurrence of flip bifurcation at DFES of HBV model (2.3). □

Theorem 5.2. At DFES, HBV model (2.3) does not undergo flip, flip, and fold bifurcations if (4.18)–(4.20) hold, respectively.

Proof. This is the same as the proof of Theorem 5.1. □

Hereafter, bifurcation analysis at ES and EES of the HBV model (2.3) is investigated by implementing explicit criterion ^{[30,31,32,33]}.

Theorem 5.3. At ES, the discrete HBV model (2.3) undergoes flip bifurcation if

where $h_0$ is a real roots of $1 - \mathcal{S}_1(h) +\mathcal{S}_2(h) -\mathcal{S}_3(h)+ \mathcal{S}_4(h) = 0$.

Proof. For $n = 4$, the explicit criterion yields

and

□

Theorem 5.4. At EES, the discrete HBV model (2.3) undergoes N-S bifurcation if

where $h_0$ is a real root of $1 - \mathcal{U}_4(h) -\mathcal{U}_2(h)-\mathcal{U}^2_3(h)+ \mathcal{U}^3_4(h)-\mathcal{U}^2_4(h)\left(1+\mathcal{U}_2(h) \right)+2 \mathcal{U}_2(h)\mathcal{U}_4(h)-\mathcal{U}^2_1(h)\mathcal{U}_4(h)+\mathcal{U}_1(h)\mathcal{U}_3(h)\left(1+\mathcal{U}_4(h)\right) = 0$.

Proof. For $n = 4$, the explicit criterion yields

and

□

Theorem 5.5. At EES, the discrete HBV model (2.3) undergoes flip bifurcation if

where $h_0$ is a real root of $1 - \mathcal{U}_1(h) +\mathcal{U}_2(h) -\mathcal{U}_3(h)+ \mathcal{U}_4(h) = 0$.

Proof. This is the same as the proof of Theorem 5.3. □

Remark 5.1. The existence of Neimark-Sacker and flip bifurcations in the context of the epidemic equilibrium state within a discrete Hepatitis B virus (HBV) model carries significant biological and epidemiological implications:

$(1)$ A Neimark-Sacker bifurcation typically results in the emergence of stable periodic orbits from an epidemic equilibrium state. In the context of an HBV model, the presence of a Neimark-Sacker bifurcation could indicate the onset of sustained oscillations in the prevalence of the virus within the population. These oscillations may signify periodic fluctuations in the number of infected individuals over time. Biologically, this could reflect cyclical patterns in disease transmission, possibly influenced by seasonal factors or periodic changes in host immunity. Epidemiologically, the appearance of stable periodic orbits could have implications for the design and timing of control measures, as interventions may need to be adapted to account for the cyclical nature of Hepatitis B virus transmission.

$(2)$ A flip bifurcation typically occurs when a stable equilibrium point undergoes a sudden change in stability, leading to the emergence of a stable limit cycle or chaotic behavior. In the context of a Hepatitis B virus model, a flip bifurcation at the epidemic equilibrium state could indicate a critical transition in disease dynamics. This transition may manifest as a sudden shift from stable disease control to sustained endemicity or periodic outbreaks. Biologically, the flip bifurcation could signify a threshold effect, where small changes in key parameters lead to qualitative changes in disease behavior. Epidemiologically, this transition could have profound implications for public health strategies, necessitating a reassessment of control measures and surveillance efforts to effectively manage the evolving dynamics of Hepatitis B virus transmission.

Overall, the presence of Neimark-Sacker and flip bifurcations in a discrete Hepatitis B virus model underscores the complex and nonlinear nature of Hepatitis B virus transmission dynamics. Understanding these bifurcations provides valuable insights into the potential for oscillatory behavior, critical transitions, and the effectiveness of control strategies, ultimately informing decision-making processes aimed at mitigating the burden of Hepatitis B virus infection on public health.

6.   Convergence rate of HBV model (2.3)

We will study the convergence rate of the HBV model (2.3) in this section as follows:

Theorem 6.1. If the positive solution of (2.3) is $\{(s_t, I_t, v_t, z_t)\}_{t = 0}^\infty$ such that $\lim\limits_{t\to\infty}\{(s_t, I_t, v_t, z_t)\} = \Phi$ then

satisfying

Proof. If the HBV model (2.3) has a positive solution $\{(s_t, I_t, v_t, z_t)\}_{t = 0}^\infty$ such that $\lim\limits_{t\to\infty}\{(s_t, I_t, v_t, z_t)\} = \Phi$ then

on setting

From (6.3) and (6.4), one has

where

From (6.6), one has

that is

where $\eta_{11}, \eta_{13}, \eta_{21}, \eta_{22}, \eta_{23}, \eta_{24}, \eta_{32}, \eta_{33}, \eta_{42}, \eta_{44} \rightarrow0$ as $t\rightarrow \infty$. From existing literature (see ^[34]), we have an error system

where $A = J|_{\Phi}$ and $B_t = \left(\begin{array}{cccc} \eta_{11} & 0 & \eta_{13} & 0\\ \eta_{21} & \eta_{22} & \eta_{23} & \eta_{24}\\ 0 & \eta_{32} & \eta_{33} & 0\\ 0 & \eta_{42} & 0 & \eta_{44} \end{array} \right)$. So, the error system becomes

which is the same as a linearized system of discrete HBV model (2.3) at $\Phi$. In particular, (6.10) gives

and

which are similar to the linearized system of the discrete HBV model (2.3) obtained at DFES, ES, and EES, respectively. □

7.   Chaos control

In this section, the chaos control is studied by the feedback control strategy for the discrete HBV model (2.3) about EES by existing theory ^[35,36]. By utilizing feedback control strategy, the discrete HBV model (2.3) takes the form

where the control parameter is $\varrho$. The $J|_{\mathrm{EES}}$ of the discrete HBV model (7.1) is

with

where

So the dynamics of (7.1) can be concluded as follows:

Theorem 7.1. EES is a sink if

where $\mathcal{U}_1^{*}, \mathcal{U}_2^{*}$, $\mathcal{U}_3^{*}$, and $\mathcal{U}_4^{*}$ are depicted in (7.4).

8.   Numerical simulations

The following cases are to be presented in order to verify the correctness of the theoretical results.

Case I. In this case, it is examined that at $h = 0.3186837362435554$, the discrete HBV model (2.3) undergoes a flip bifurcation if $\vartheta = 13.9, \ \beta = 2.654, \ \alpha = 3.21, \ k = 1.976, \ u = 3, \ \delta = 1.004, \ c = 1.4, \ \mu = 1.7432768, \ p = 1.9$, and $h\in[0.01, 0.3]$ with $(s_0, I_0, v_0, z_0) = (0.0147, 0.169, 0.2, 0.0)$. If $\vartheta = 13.9, \ \beta = 2.654, \ \alpha = 3.21, \ k = 1.976, \ u = 3, \ \delta = 1.004, \ c = 1.4, \ \mu = 1.7432768, \ p = 1.9$, and $h = 0.3186837362435554$, then from (4.22), we get

with $\lambda_{1} = -1$ but $\lambda_{2, 3, 4} = 0.49978202821265805$, $0.8068705153090874$, $0.8250500673962229\not = 1\ \mathrm{or} \ -1$, and for the occurrence of flip bifurcation eigenvalues criterion holds. Therefore, the discrete HBV model (2.3) may undergo flip bifurcation. But following simulation ensures that flip bifurcation must exist, and so, if $\vartheta = 13.9, \ \beta = 2.654, \ \alpha = 3.21, \ k = 1.976, \ u = 3, \ \delta = 1.004, \ c = 1.4, \ \mu = 1.7432768, \ p = 1.9$, and $h = 0.3186837362435554$, then from (5.3), one gets

Thus, from (8.2), all conditions of Theorem 5.3 hold, and hence, flip bifurcation takes place at ES. Finally, MLE along with flip bifurcation diagrams are shown in Figure 2.

Case II. If $h = 0.0734578$, $\vartheta = 4.3$, $\delta = 2.004$, $\alpha = 0.41$, $\beta = 2.5$, $p = 1.9$, $k = 5.5$, $u = 4.6$, $c = 8.4$, and $\mu = 4.432768$, then from (4.28), one has $\left|\mathcal{U}_1 + \mathcal{U}_3\right| = 4.954082970547482 < 1+\mathcal{U}_2+\mathcal{U}_4$ $= 4.960842450652595$, $\left|\mathcal{U}_3-\mathcal{U}_1\right| = 1.3166773911382283 < 2(1-\mathcal{U}_4) = 1.3591781755420462$, $\mathcal{U}_3-3 \mathcal{U}_4 = 2.679198801736687 < 3, \ \mathrm{and} \ \mathcal{U}_2+{\mathcal{U}_4}^2(1+\mathcal{U}_2)+{\mathcal{U}_3}^2+\mathcal{U}_4(1+{\mathcal{U}_1}^2) = 10.894757977942465 < 1+2\mathcal{U}_2\mathcal{U}_4+ \mathcal{U}_1\mathcal{U}_3(1+\mathcal{U}_4)+{\mathcal{U}_4}^3 = 10.895174109435445$, and therefore, $\mathrm{EES} = \Bigr(1.200649318366021, 0.5277104761904762, 0.6309581780538303, 1.6731039047520355\Bigl)$ of the discrete HBV model (2.3) is a sink (see Figure 3).

Case III. Now, we will examine that at $h = 0.00032327767915996635$, the discrete HBV model (2.3) undergoes a N-S bifurcation if $\vartheta = 4.3, \ \beta = 2.5, \ \alpha = 0.41, \ k = 5.5, \ u = 4.6, \ \delta = 3.004, \ c = 8.4, \ \mu = 4.432768, \ p = 1.9$, and $h\in[0.01, 0.2]$ with $(s_0, I_0, v_0, z_0) = \left(0.147, 3.169, 1.2, 1.3\right)$. If $\vartheta = 4.3, \ \beta = 2.5, \ \alpha = 0.41, \ k = 5.5, \ u = 4.6, \ \delta = 3.004, \ c = 8.4, \ \mu = 4.432768, \ p = 1.9$, and $h = 0.00032327767915996635$, then from (4.26), one has

with $\lambda_{1, 2} = 0.9998769965433948\pm0.000806187461780547\iota$, $\lambda_3 = 0.9978504688163795$, and $\lambda_4 = 0.9990487427589927$, where $|\lambda_{1, 2}| = 1$. Therefore, the eigenvalues criterion holds for the occurrence of said bifurcation, and so, the HBV model (2.3) may undergo a N-S bifurcation. But the rest of the simulation shows that N-S bifurcation must exist, that is, if $\vartheta = 4.3, \ \beta = 2.5, \ \alpha = 0.41, \ k = 5.5, \ u = 4.6, \ \delta = 3.004, \ c = 8.4, \ \mu = 4.432768, \ p = 1.9$, and $h = 0.00027921473998766453$, then from (5.6), we get

Moreover, $\cos\frac{2\pi}{l} = 0.987$ implies $l = \pm4.33246986002$. From (8.4), it is concluded that all conditions of Theorem 5.4 hold, and hence, at EES, the HBV model (2.3) undergoes N-S bifurcation, where the MLE along with N-S bifurcation diagrams are plotted in Figure 4.

Case IV. If $\vartheta = 12.3, \ \beta = 2.45, \ \alpha = 4.1, \ k = 4.256, \ u = 1, \ \delta = 1.004, \ c = 3.84$, $\mu = 1.7432312322768, \ p = 1.9$, and $h = 0.08978168084826603$, then from (4.26), one gets:

with $\lambda_{1} = -1$ but $\lambda_{2, 3, 4} = 0.6929389869902541$, $0.7851203029606955$, $0.9101064264695002\not = 1\ \mathrm{or} \ -1$. This implies that eigenvalues criterion holds for the appearance of flip bifurcation if $\vartheta = 12.3, \ \beta = 2.45, \ \alpha = 4.1, \ k = 4.256, \ u = 1, \ \delta = 1.004, \ c = 3.84, \ \mu = 1.7432312322768, \ p = 1.9$, and $h = 0.08978168084826603$. Furthermore, if $\vartheta = 12.3, \ \beta = 2.45, \ \alpha = 4.1, \ k = 4.256, \ u = 1, \ \delta = 1.004, \ c = 3.84, \ \mu = 1.7432312322768, \ p = 1.9$, and $h = 0.08978168084826603$, then from (5.9), one has

From (8.6), all parametric conditions of Theorem 5.5 hold, and at EES, the HBV model (2.3) undergoes flip bifurcation (see Figure 5).

Case V. If $h = 0.07834578$, $\vartheta = 4.3$, $\delta = 3.004$, $\alpha = 0.41$, $\beta = 2.5$, $p = 1.9$, $k = 5.5$, $u = 4.6$, $c = 8.4$, $\mu = 4.432768$, and $\varrho = 0.007$, then from (7.5), one gets $\left|\mathcal{U}_1^{*} + \mathcal{U}_3^{*}\right| = 4.810518929365905 < 1+\mathcal{U}_2^{*}+\mathcal{U}_4^{*} = 4.8185666034767305$, $\left|\mathcal{U}_3^{*}-\mathcal{U}_1^{*}\right| = 1.3672303900825242 < 2(1-\mathcal{U}_4^{*}) = 1.4041181034286887$, $\mathcal{U}_2^{*}-3 \mathcal{U}_4^{*} = 2.626802810334108 < 3$, and finally, $\mathcal{U}_2^{*}+{\mathcal{U}_4^{*}}^2(1+\mathcal{U}_2^{*})+{\mathcal{U}_3^{*}}^2+ \mathcal{U}_4^{*}(1+{\mathcal{U}_1^{*}}^2) = 10.026614433151806 < 1+2\mathcal{U}_2^{*}\mathcal{U}_4^{*}+\mathcal{U}_1^{*} \mathcal{U}_3^{*}(1+\mathcal{U}_4^{*})+{\mathcal{U}_4^{*}}^3 = 10.026701399999961$, which implies that $\mathrm{EES} = \Bigl(0.516585, 2.7202, 1.10213, 8.60536\Bigr)$ of the control HBV model (7.1) is a sink. The plots for control HBV model (7.1) with $(s_0, I_0, v_0, z_0) = (3.06276, 0.9821, 0.9, 0.1)$ are shown in Figure 6 which demonstrates that $\mathrm{EES} = \Bigl(0.9385786604748523, 0.5277104761904762, 0.6309581780538303, 1.2608 073891909866\Bigr)$ is a sink.

9.   Conclusions

This work is about the local dynamics at equilibrium states, and we have examined the basic reproduction number, rate of convergence, and bifurcation analysis of discrete HBV model (2.3). We examined that for all of the model's parameters, discrete HBV model (2.3) had a disease-free equilibrium state. Furthermore, under certain parametric restriction(s), discrete HBV model (2.3) had a boundary equilibrium state, that is, the state where the pathogens are present while CTLs are absent, and an epidemic equilibrium state. Next, by liner stability theory, we examined the local dynamical properties at equilibrium states of discrete HBV model (2.3). We have also studied the convergence rate for discrete HBV model (2.3). Further, in order to understand the dynamics of discrete HBV model (2.3) deeply, we have studied the possible bifurcation scenarios. We have proved that with a disease-free equilibrium state, there exists no flip and fold bifurcations, but discrete HBV model (2.3) undergoes flip bifurcations about the boundary equilibrium state, and Neimark-Sacker and flip bifurcations about the epidemic equilibrium state. We have studied flip bifurcations about the boundary equilibrium state, and Neimark-Sacker and flip bifurcations about the epidemic equilibrium state of discrete HBV model (2.3) by utilizing explicit criterion. The rate of convergence and chaos in discrete HBV model (2.3) were also examined. Finally, our main findings are illustrated numerically.

Boundedness, persistence, and global dynamics at equilibrium states of a discrete HBV model (2.3) are our next aim to study.

Author contributions

Abdul Qadeer Khan: Conceptualization, formal analysis, investigation, methodology, resources, software, supervision, validation, visualization, writing original draft, writing review & editing; Fakhra Bibi: Formal analysis, investigation, software, writing original draft, writing review & editing; Saud Fahad Aldosary: Conceptualization, funding acquisition, resources, writing original draft, writing review & editing. All authors have read and approved the final version of the manuscript for publication.

Use of AI tools declaration

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

Acknowledgements

This study is supported via funding from Prince Sattam bin Abdulaziz University, project number (PSAU/2024/R/1445).

Conflict of interest

The authors declare that they have no conflicts of interest.