Neimark-Sacker bifurcation, chaos, and local stability of a discrete Hepatitis C virus model

1.   Introduction

Chronic Hepatitis C virus (HCV) infection is a significant global health concern, affecting an estimated 71 million individuals worldwide. This bloodborne infection, primarily targeting the liver, led to approximately 399,000 deaths in 2016 due to severe complications such as liver cirrhosis and hepatocellular carcinoma. A large proportion (60%–80%) of those infected with HCV develop chronic infection, with 15%–30% of these individuals at risk of progressing to liver cirrhosis over time. The primary modes of transmission include unscreened blood transfusions, injection drug use, and the use of improperly sterilized medical equipment, making it a preventable yet persistent health issue in many parts of the world. Despite the recent advancements in antiviral treatments, particularly with the introduction of direct-acting antivirals (DAAs), which have significantly improved cure rates, interferon-alpha IFN-$\alpha$ therapy has long played a central role in managing HCV infection. Traditionally administered either alone or in combination with ribavirin, IFN-$\alpha$ has been instrumental in achieving sustained virologic response (SVR) in chronic HCV patients, thereby preventing disease progression. However, the precise mechanisms through which IFN-$\alpha$ exerts its therapeutic effects, especially in terms of HCV dynamics within the host, remain not fully understood. HCV infection triggers a complex interplay of biological processes within the liver. Once the virus enters hepatocytes, it initiates viral replication, leading to an immune response from the host. At the same time, the liver's natural regenerative capacity attempts to maintain homeostasis by proliferating healthy hepatocytes. This dynamic interaction between viral replication, immune activation, and hepatocellular proliferation presents a challenging scenario for effective treatment. The availability of uninfected hepatocytes, the rate of hepatocyte proliferation, and the overall homeostatic balance of the liver become critical factors in determining the progression of the infection and the host's ability to control it. To better understand these intricate processes, mathematical modeling has become a crucial tool in studying HCV dynamics. Mathematical models allow researchers to simulate the spread of the virus within the liver, analyze how the infection interacts with treatment protocols like IFN-$\alpha$ and DAAs, and explore the role of immune response and hepatocyte regeneration. These models not only provide insights into the transmission dynamics of HCV but also help to refine treatment strategies by identifying key factors influencing viral clearance or persistence. Through such modeling, it becomes possible to investigate the intricate balance between viral proliferation, host immune defense, and the liver's regenerative capabilities, offering valuable perspectives on both disease progression and therapeutic intervention. To illustrate, Pan and Chakrabarty ^[1] explored the Hopf bifurcation of the HCV model

where $T$, $V$, $I$, and $Z$ respectively denote the densities of uninfected hepatocytes, virions, actively infected hepatocytes, and humoral immunity. Mojaver and Kheiri ^[2] explored the behavior of a class of HCV infection models with application of optimal control. Gümüş and Türk ^[3] examined the dynamical analysis of a HBV model. Kong et al. ^[4] examined the Hopf bifurcation of the HCV model

where $T$ and $I$ respectively denote uninfected and infected cells, while $V$ is free virus. Furthermore, cells $T$ are generated, die, and are infected at rates $s$, $d_1 T$, and $\beta T V$, cells $I$ emerge and die at rates $\beta T V$ and $b I$, and finally, free virus are released and decay at rates $k I$ and $u V$. Jiang et al. ^[5] examined the bifurcation, chaos, and circuit realisation of a four-dimensional memristor system. Chen et al. ^[6] explored the bifurcation, chaos and fixed-time synchronization of memristor cellular neural networks. On the other hand, many authors have performed theoretical and mathematical analysis of the mathematical models ^[7,8,9]. Furthermore, different mathematicians ^[10,11,12] have first analyzed the continuous-time HCV model

and then showed that it takes the following form after decomposing the viral dynamics to two distinct time scales

It is worthwhile to mention here that discrete-time models for Hepatitis C virus (HCV) dynamics are often preferred over continuous-time models due to several practical advantages. One of the primary reasons is their compatibility with the nature of clinical data collection. In practice, data such as viral load measurements, liver enzyme levels, and patient responses to treatments are typically recorded at discrete time intervals, such as daily, weekly, or monthly. Discrete-time models naturally align with these sampling intervals, allowing for direct integration of empirical data without the need for interpolation or extrapolation, which can introduce errors and complicate model validation. Furthermore, discrete-time models tend to be computationally more efficient. Solving continuous-time models involves complex numerical methods and often requires fine time discretization to ensure accuracy, which increases computational costs. In contrast, discrete-time models can achieve accurate results with simpler algorithms and larger time steps, making them more practical for large-scale simulations or real-time applications. This efficiency is particularly beneficial when modeling the long-term progression of HCV or evaluating various treatment strategies over extended periods. The flexibility of discrete-time models is another significant advantage. They can more easily incorporate complex biological processes and treatment protocols. For instance, antiviral drug effects, immune responses, and patient adherence to therapy, which often change at specific time points (e.g., weekly dosing schedules), can be modeled more intuitively in a discrete framework. This capability is crucial for accurately reflecting the real-world dynamics of HCV and its treatment, where interventions are frequently time-dependent. From an implementation perspective, discrete-time models are generally easier to understand and use. Researchers and clinicians can work with models that update state variables at each time step based on straightforward recurrence relations. This simplicity facilitates broader adoption of the models in both clinical and research settings, promoting better communication and collaboration among multidisciplinary teams. The ease of implementation also lowers the barrier for developing and testing new hypotheses or treatment strategies, accelerating progress in HCV research. Lastly, discrete-time models offer greater robustness to uncertainties and variabilities inherent to biological systems and clinical data. They can handle irregular sampling intervals and missing data more effectively, which is critical in real-world scenarios where ideal, continuous data collection is rarely feasible. This robustness ensures that the models remain reliable and useful even when the data is imperfect, enhancing their applicability and trustworthiness in guiding clinical decisions and policy-making. In summary, while continuous-time models have their theoretical advantages, the practical benefits of discrete-time models such as data compatibility, computational efficiency, flexibility, ease of implementation, and robustness to uncertainties make them a preferable choice for modeling HCV dynamics in many contexts. Due to the aforementioned facts, our aim in this paper is to explore the dynamics of a discrete-time HCV model instead of continuous-time HCV model, which is depicted in (1.4), where the discrete version is obtained by a non-standard finite difference scheme. The dynamics of a HCV infection model can vary depending on whether it is a discrete-time model or a continuous-time model. Below, we shall highlight the key differences in their dynamical behaviors, focusing on stability, bifurcations, and chaotic behaviors. In a continuous-time model, the system evolves continuously over time, with changes happening smoothly at every instant; in a discrete-time model, the system evolves at specific time intervals, and changes happen at discrete points. Discrete systems are prone to phenomena like Neimark-Sacker and period-doubling bifurcations, which can lead to quasi-periodic behavior and chaos more readily than in continuous models. Numerically, continuous models are often simulated using methods like Euler or Runge-Kutta, requiring fine time steps to approximate smooth dynamics. Discrete models, however, are simpler to implement and can reveal complex behaviors with fewer time steps. The contrast between the smooth progression of continuous models and the stepwise, often abrupt behavior of discrete models highlights how the mathematical formulation influences the biological insights derived from HCV infection dynamics. In order to derive the discrete version, we replace $\dot{u}$ by $\frac{u_{t+1}-u_t}{h}$, $\dot{v}$ by $\frac{v_{t+1}-v_t}{h}$, $u^2$ by $u_{t}u_{t+1}$, $v^2$ by $v_{t}v_{t+1}$, $u$ by $u_t$, and $v$ by $v_t$, and so continuous HCV (1.4) becomes

with $h$ being the integral step size. Furthermore, after simplification, (1.5) takes the form

More precisely, our main findings in this paper include:

● Existence of fixed points of a discrete model (1.6).

● Dynamics at fixed points of a discrete model (1.6).

● Identification of bifurcation sets for fixed points.

● Bifurcation analysis of a discrete model (1.6).

● Study of chaos.

● Numerical validation of theoretical results.

The rest of the paper is organized as follows: Existence of fixed points are given in Section 2. Section 3 regards the local behavior at fixed points of the model (1.6). Bifurcation analysis is given in Section 4, whereas in Section 5 chaos control is studied. The numerical simulation of obtained results is given in Section 6, while the conclusions are given in Section 7.

2.   Existence of fixed points

In this section, we will study the existence of fixed points of the discrete model (1.6). If $(\hat{u}, \hat{v})$ is a fixed point of model (1.6), then

Since system (2.1) is obviously satisfied if $(\hat{u}, \hat{v}) = (0, 0), \left(1-g, 0\right), \left(0, 1-\zeta\right)$. Therefore, $\forall$ $b$, $\alpha$, $g$, $\zeta$, $h$, discrete HCV model (1.6) has liver-free fixed point (LFFP), disease-free fixed point (DFFP) $\left(1-g, 0\right)$ if $0 < g < 1$, and total infection fixed point (TIFP) $\left(0, 1-\zeta\right)$ if $0 < \zeta < 1$. In order to find the partial-infection fixed point (PIFP), from (2.1), one has

From the first equation of (2.2), one gets

Using (2.3) into the second equation of (2.2), we get

Finally, using (2.4) into (2.3), one gets

Equations (2.4) and (2.5) imply that if $\frac{\zeta-g}{(1-\alpha)(1-g)} < b < \frac{\zeta-g}{(1-\alpha)(1-\zeta)}$, then HCV model (1.6) has PIFP $\left(\frac{\zeta-g-(1-\alpha)(1-\zeta)b}{(1-\alpha)^2 b^2}, \frac{g-\zeta+(1-\alpha)(1-g)b}{(1-\alpha)^2 b^2}\right)$.

3.   Local stability

In this section, the local stability at fixed points of the HCV model (1.6) is analyzed using linear stability theory ^[13,14]. Specifically, we explore local stability by evaluating the eigenvalues of the Jacobian matrix at fixed points. This approach allows us to assess whether perturbations around the fixed points are stable or unstable. For this, we first compute the Jacobian matrix at fixed points, and then eigenvalues will be determined in order to apply linear stability theory ^[15,16]. So, Jacobian $J|_{(\hat{u}, \hat{v})}$ evaluated at $(\hat{u}, \hat{v})$ under the map $\left(g_1, g_2\right)\mapsto \left(u_{t+1}, v_{t+1}\right)$ is

where

Now, at LFFP, (3.1) becomes

with

From (3.4), it is easy to see that LFFP is a saddle if $\frac{2}{\zeta-1} < h < \frac{2}{g-1}$ or $\frac{2}{g-1} < h < \frac{2}{\zeta-1}$, a source if $0 < h < \min\left\{\frac{2}{g-1}, \frac{2}{\zeta-1}\right\}$, a sink if $h > \max\left\{\frac{2}{g-1}, \frac{2}{\zeta-1}\right\}$, and finally, non-hyperbolic if

or

At DFFP, (3.1) becomes

with

From (3.8), one can also obtain that DFFP is a saddle if $\frac{2}{g-1} < h < \frac{2}{-g+\zeta-(1-g)(1-\alpha)b}$ or $\frac{2}{-g+\zeta-(1-g)(1-\alpha)b} < h < \frac{2}{g-1}$, a source if $0 < h < \min\left\{\frac{2}{g-1}, \frac{2}{-g+\zeta-(1-g)(1-\alpha)b}\right\}$, a sink if $h > \max\left\{\frac{2}{g-1}, \frac{2}{-g+\zeta-(1-g)(1-\alpha)b}\right\}$, and non-hyperbolic if

or

Furthermore, at TIFP, (3.1) becomes

with

from which it is easy to prove that TIFP is a saddle if $\frac{2}{\zeta-1} < h < \frac{2}{-\zeta+g+(1-\alpha)(1-\zeta)b}$ or $\frac{2}{-\zeta+g+(1-\alpha)(1-\zeta)b} < h < \frac{2}{\zeta-1}$, a source if $0 < h < \min\left\{\frac{2}{\zeta-1}, \frac{2}{-\zeta+g+(1-\alpha)(1-\zeta)b}\right\}$, a sink if $h > \max\left\{\frac{2}{\zeta-1}, \frac{2}{-\zeta+g+(1-\alpha)(1-\zeta)b}\right\}$, and non-hyperbolic if

or

Finally, at PIFP, (3.1) becomes

with

where

The $J|_{\mathrm{PIFP}}$ has characteristics roots

where

If $\Delta < 0$, then PIFP is an unstable focus if $h > \frac{(1-\alpha)^3(\zeta-g)b^3}{(2-(1-\alpha)^2b^2)(g+b(1-g)(1-\alpha)-\zeta)(g+b(1-\zeta)(1-\alpha)-\zeta)}$, a stable focus if $0 < h < \frac{(1-\alpha)^3(\zeta-g)b^3}{(2-(1-\alpha)^2b^2)(g+b(1-g)(1-\alpha)-\zeta)(g+b(1-\zeta)(1-\alpha)-\zeta)}$, and non-hyperbolic if

In addition, if $\Delta > 0$, then PIFP is an unstable node if

a stable node if

and non-hyperbolic if

or

Example 1. If $b = 12.3, g = 0.0006, \zeta = 4.5$, then at PIFP, the dynamics of the discrete HCV model (1.6) are given in Figure 1.

Example 2. If $b = 0.5, g = 0.09, \zeta = 0.55$, then at PIFP, the dynamics of the discrete HCV model (1.6) are given in Figure 2.

4.   Bifurcation analysis

In this section, detailed bifurcation analysis is given based on bifurcation theory ^[17,18,19]. Before detailed bifurcation analysis, we first give the existence of possible bifurcation sets for fixed points as follows:

(a) If (3.5) holds, then roots of $J|_{\mathrm{LFFP}}$ satisfying $\lambda_1|_{(3.5)} = -1$ but $\lambda_2|_{(3.5)} = \frac{g-1+2(1-\zeta)}{g-1}\not = 1 \ \mathrm{or}\ -1$. So, period-doubling bifurcation may take place if $\left(b, \alpha, g, \zeta, h \right) = \mathfrak{L}$ crosses through the curve

(b) if (3.6) holds, then roots of $J|_{\mathrm{LFFP}}$ satisfying $\lambda_1|_{(3.6)} = \frac{1-2g+\zeta}{\zeta-1}\not = 1 \ \mathrm{or}\ -1$ but $\lambda_2|_{(3.6)} = -1$. So, period-doubling bifurcation may take place if $\mathfrak{L}$ crosses through the curve

(c) if (3.9) holds, then roots of $J|_{\mathrm{DFFP}}$ satisfying $\lambda_1|_{(3.9)} = -1$ but $\lambda_2|_{(3.9)} = 1+\frac{2(g-\zeta+b(1-\alpha)(1-g))}{g-1} \not = 1 \ \mathrm{or}\ -1$. So, period-doubling bifurcation may take place if $\mathfrak{L}$ crosses through the curve

(d) if (3.10) holds, then roots of $J|_{\mathrm{DFFP}}$ satisfying $\lambda_1|_{(3.10)} = \frac{-g+\zeta-b(1-\alpha)(1-g)}{2-3g-b(1-\alpha)(1-g)+\zeta}\not = 1 \ \mathrm{or}\ -1$ but $\lambda_2|_{(3.10)} = -1$. So, period-doubling bifurcation may take place if $\mathfrak{L}$ crosses through the curve

(e) if (3.10) holds, then roots of $J|_{\mathrm{TIFP}}$ satisfying $\lambda_1|_{(3.10)} = 1+\frac{2}{\zeta-1}(\zeta-g-b(1-\alpha)(1-\zeta))\not = 1 \ \mathrm{or}\ -1$ but $\lambda_2|_{(3.10)} = -1$. So, period-doubling bifurcation may take place if $\mathfrak{L}$ crosses through the curve

$(f)$ if (3.14) holds, then roots of $J|_{\mathrm{TIFP}}$ satisfying $\lambda_1|_{(3.14)} = -1$ but $\lambda_2|_{(3.14)} = \frac{g-\zeta+b(1-\alpha)(1-\zeta)}{2-3\zeta+b(1-\alpha)(1-\zeta)+g}$. So, period-doubling bifurcation may take place if $\mathfrak{L}$ crosses through the curve

(g) if (1.3) holds, then a pair of complex roots of $J|_{\mathrm{PIFP}}$ satisfying $\left|\lambda_{1, 2}\right|_{(1.3)} = 1$. So, Neimark-Sacker bifurcation may take place if $\mathfrak{L}$ crosses through the curve

(h) if (3.21) is a true, then the value of $\lambda_{1, 2}$ at (3.21) has been calculated and is depicted in (A.1). This implies that no period-doubling bifurcation takes place if $\mathfrak{L}$ crosses through the curve

(i) if (3.22) holds, then the value of $\lambda_{1, 2}$ at (3.22) has been calculated and is depicted in (A.2). This implies that no period-doubling bifurcation takes place if $\mathfrak{L}$ crosses through the curve

Now we will study detailed bifurcation analysis at fixed points based on the bifurcation sets. It is mentioned that no period-doubling bifurcation occurs at LFFP if (4.1) and (4.2) hold, and therefore, LFFP is degenerate. Now, we examine bifurcation analysis at DFFP, TIFP, and PIFP of HCV model (1.6) if (4.3)–(4.9) hold.

Lemma 4.1. The discrete HCV model (1.6) does not undergo period-doubling at respective fixed points if (4.3)–(4.6) hold.

Proof. The HCV model (1.6), takes the following form after being restricted to $v = 0$,

From (4.10), we define

Now, if $h = \hat{h} = \frac{2}{g-1}$ and $u = \hat{u} = 1-g$, then from (4.11), we get

and

So, from (4.14), no period-doubling bifurcation occurs at DFFP if (4.3) holds. Similarly, it is also easy to obtain that no period-doubling exists at the respective fixed points if (4.4)–(4.6) hold. □

Remark 4.2. In order to study the Neimark-Sacker bifurcation at PIFP of the HCV model (1.6), first it is imperative to highlight the seminal work of Rahman and Blackmore ^[20]. This paper presented significant advancements in understanding the complex behavior exhibited by discrete dynamical systems. The study meticulously examined the onset of chaos through Neimark-Sacker bifurcations, providing crucial insights into the mechanisms that drive such phenomena. Rahman and Blackmore ^[20] contributions to this research are particularly noteworthy for their depth and clarity in elucidating the mathematical underpinnings of chaotic behavior in these systems. In the broader context of modern dynamical systems analysis, Rahman and Blackmore's work has consistently demonstrated a pioneering spirit. Their investigations into the stability and bifurcation of nonlinear systems have laid a foundational framework that subsequent researchers have built upon. By systematically exploring the conditions under which chaotic dynamics emerge, Rahman and Blackmore have advanced our understanding of the intricate balance between order and disorder in mathematical models of physical phenomena. The meticulous approach adopted in their study underscores the critical role that rigorous mathematical analysis plays in unveiling the complexities of dynamical systems. Moreover, the implications of Rahman and Blackmore's research extend beyond theoretical exploration, influencing practical applications in various fields where understanding chaotic behavior is essential. From engineering to biological systems, the insights derived from their work on Neimark-Sacker bifurcations have provided valuable tools for predicting and managing chaotic dynamics. This integration of theoretical rigor with practical relevance is a hallmark of Rahman and Blackmore's legacy in the field, reflecting the profound impact on both the academic community and applied sciences. In conclusion, the work of Rahman and Blackmore on the study of Neimark-Sacker bifurcations and chaotic behavior stands as a testament to their pioneering contributions to dynamical systems analysis. By bridging the gap between abstract mathematical theory and real-world applications, Rahman and Blackmore left an indelible mark on the field, inspiring future generations of researchers to explore the rich and intricate world of nonlinear dynamics. This legacy continues to influence contemporary studies and will undoubtedly remain a cornerstone of dynamical systems research for years to come.

Theorem 4.3. If (4.7) holds, then the HCV model (1.6) undergoes a Neimark-Sacker bifurcation with $h$ being a bifurcation parameter.

Proof. By considering $h$ as a bifurcation parameter, the discrete HCV model (1.6) becomes

where $\epsilon < < 1$, and the characteristic roots of $J|_{\mathrm{PIFP}}$ at PIFP of $\epsilon$-dependence HCV model (4.15) are

where

From (4.16) and (4.17), one has

Now, (4.15) becomes

by

where $\hat{u} = \frac{\zeta-g-(1-\alpha)(1-\zeta)b}{(1-\alpha)^2 b^2}$ and $\hat{v} = \frac{g-\zeta+(1-\alpha)(1-g)b}{(1-\alpha)^2 b^2}$. Hereafter, following Kuznetsov's argument for a Neimark-Sacker bifurcation involves a detailed analysis of the model's dynamics as it undergoes this type of bifurcation, which is characterized by the emergence of a closed invariant curve from a fixed point as parameters are varied. The first step in this analysis is to bring the system into a normal form that simplifies the bifurcation process. This normal form representation allows for the computation of the first Lyapunov coefficients, which are essential in determining the stability and nature of the bifurcating solutions. The first Lyapunov coefficient, in particular, indicates whether the bifurcation will lead to stable or unstable quasiperiodic orbits. By computing these coefficients, one can predict the onset of chaotic dynamics and the subsequent behavior of the system. Kuznetsov's methodology provides a robust framework for analyzing the local behavior near the bifurcation point, ensuring a thorough understanding of the system's transition to chaos. This approach is crucial for identifying the critical parameters and conditions that govern the complex behavior in nonlinear dynamical systems. For more detail, we also refer the reader to the existing theory ^[17,18,20]. Now, in order to study the normal form of the HCV model (4.19) if $\epsilon = 0$, the Taylor's series expansion at $\left(0, 0\right)$ yields

where

Now, if

then (4.21) becomes

where

and

with

From (4.26), we get

Now, by the existing theory ^[21,22,23]. the map (4.24) undergoes Neimark-Sacker bifurcation if $\psi\not = 0$, that is,

where

After calculation, one gets

Using (4.31) along with (4.25) in (4.29) if one obtains $\psi\not = 0$, then at PIFP, the HCV model (1.6) undergoes N-S bifurcation. Additionally, it also requires that $\lambda_{1, 2}^m\neq 1, \ m = 1, \cdots, 4$ if $\epsilon = 0$ that corresponds to $\varrho_1(0)\neq -2, 0, 1, 2$. But if $h = \frac{(1-\alpha)^3b^3(\zeta-g)}{(2-b^2(1-\alpha)^2)(g+b(1-g)(1-\alpha)-\zeta)(g+b(1-\zeta)(1-\alpha)-\zeta)}$ holds, then from (4.17), we get $\varrho_2(0) = 1$, and so $\varrho_1(0)\not = -2, 2$. Therefore, it is only required that $\varrho_1(0)\not = 0, 1$, i.e.,

On the other hand, also by existing theory ^{[24,25,26,27]}, the Neimark-Sacker bifurcation at PIFP is also supercritical (resp. subcritical) if $\psi < 0$ (resp. $\psi > 0$). □

Remark 4.4. It is important to note that no period-doubling bifurcation occurs if (4.8) and (4.9) hold. In conclusion, at PIFP, no other bifurcations occur except the Neimark-sacker bifurcation.

5.   Chaos control

In this section, the OGY method, which was proposed by Ott et al. ^[28], is utilized for the HCV model (1.6). For this, one writes the HCV model (1.6) as follows

where the control parameter is $\zeta$ to acquire the desired chaos control by small perturbations. Due to this, we restrict it where $\zeta \in \left(\zeta_0- \sigma, \zeta_{0}+ \sigma \right)$ with $\sigma > 0$, and $\zeta_0$ indicates small value corresponding to chaotic region. The trajectory is moved in the direction of the target orbit using the stabilizing feedback control strategy. Assuming that PIFP is an unstable point in the chaotic region created by the occurrence of N-S bifurcation, then the following map can be used to approximate the model (5.1)

where

and

Moreover, the HCV model (5.1) is controlled provided that the following matrix

is of rank 2. So, to apply the OGY method for the HCV model (1.6), we choose $\zeta$ as a control parameter. Now, if $\zeta-\zeta_0 = -T\left(\begin{array}{cc} u_t-\hat{u}\\ v_t-\hat{v} \end{array} \right)$, where $T = (\hat{k_1} \ \ \hat{k_2})$, then model (5.2) can be written as

Thus, the corresponding controlled model of (1.6) is given by

By the theory of stability, PIFP is a sink iff roots of $J_{\left(\hat{u}, \hat{v}, \zeta_0\right)}-PT$ satisfying $|\lambda_{1, 2}| < 1$ where its Jacobian is

$J_{\left(\hat{u}, \hat{v}, \zeta_0\right)}-PT$ has the following characteristic equation

If one denotes roots of (5.9) are $\lambda_{1, 2}$, then

and

Now, to get lines of marginal stability for (5.7), one requires $\lambda_1 = \pm 1$ and $\lambda_1\lambda_2 = 1$, where these restrictions give characteristics roots satisfying $|\lambda_{1, 2}| < 1$. If $\lambda_1\lambda_2 = 1$, then from (5.11), we have

If $\lambda_1 = 1$, then from (5.10) and (5.11), we have

Finally, if $\lambda_1 = -1$, then from (5.10) and (5.11), we have

Therefore, the stable eigenvalues are located within the triangular region $(\hat{k_1}\hat{k_2})$-plane bounded by $L_i(i = 1, 2, 3)$ for parametric values $b$, $\alpha$, $g$, $\zeta$, and $h$.

6.   Numerical simulations

The following two cases are to be considered for the completeness of this section.

Case Ⅰ. If $b = 12.01, \alpha = 0.665, g = 0.1232, \zeta = 2.002, h = [0.55, 1]$ with the initial condition $(u_0, v_0) = (0.012352, 0.232)$, then at $h = 0.8850211364486024$, the HCV model (1.6) undergoes the Neimark-Sacker bifurcation, and Figure 3 shows the maximum Lyapunov exponent with bifurcation diagrams. Further, if $b = 12.01$, $\alpha = 0.665$, $g = 0.1232$, $\zeta = 2.002$, $h = 0.8850211364486024$, then model (1.6) has $\mathrm{PIFP} = (0.36511216748407743, 0.1018618715629854)$ and, moreover, from (3.15), one gets

with

The roots are $\lambda_{1, 2} = 0.8365437672560481\pm 0.5479001053705491 \iota$ with $|\lambda_{1, 2}| = 1$, which implies from the respective bifurcation curve (4.7) that $\left(b, \alpha, g, \zeta, h\right)$ = $(12.01, 0.665, 0.1232, 2.002, 0.8850211364486024)$$\in \mathcal{N}|_{\mathrm{PIFP} = (0.36511216748407743, 0.1018618715629854)}$. Furthermore, by fixing $b = 12.01$, $\alpha = 0.665$, $g = 0.1232$, $\zeta = 2.002$ and varying $h = 0.8850211364486024$, then one can conclude that the corresponding PIFP is a stable focus. To verify this, if $h = 0.795 < 0.8850211364486024$, then Figure 4(a) indicates that PIFP is a stable focus. Similarly, if $h = 0.805$, $0.835$, $0.855$, $0.865$, $0.875 < 0.8850211364486024$, then Figures 4(b)–(f) also show that the corresponding PIFP is stable focus. Conversely, if $h > 0.8850211364486024$, then respective PIFP is an unstable focus, and as a result, supercritical N-S bifurcation takes place for the HCV model (1.6). For instance, if $b = 12.01, \alpha = 0.665, g = 0.1232$, and $\zeta = 2.002$, then from (4.18), one gets $\frac{d|\lambda_{1, 2}|}{d \epsilon}\biggr|_{\epsilon = 0} = 0.6070024854770324 > 0$, and if $h = 0.895 > 0.8850211364486024$, then from (4.22) and (4.25), one gets

From (6.3), (6.4), and (4.27), one gets

From (6.5) and (4.31), one gets

Now, utilizing (6.6) with $\lambda, \bar{\lambda} = 0.8350788018880821 \pm 0.39028624560846026\iota$ into (4.29), one obtains $\psi = -0.6734309764490927 < 0$, indicating the existence of a closed invariant curve. Thus, at the PIFP, the HCV model undergoes a supercritical N-S bifurcation (see Figure 5(a). Similarly, for $h = 0.905$, $0.915$, $0.925$, $0.935$, $0.945 > 0.8850211364486024$, stable invariant curves also emerge, as depicted in Figures 5(b)–(f). It is worth highlighting that similar bifurcation behavior, as depicted in Figures 4 and 5, has been observed in another system ^[29].

Case Ⅱ. If one reconsiders the data from Case Ⅰ with $\zeta_0 = 0.0001$, then (5.7) becomes

where

with lines of marginal stability being

and

Therefore, the lines (6.9)–(6.11) determine a triangular region that gives $|\lambda_{1, 2}| < 1$ (see Figure 6).

7.   Conclusions

This work investigates the dynamical characteristics of a discrete two-dimensional Hepatitis C virus (HCV) model (1.6), specifically focusing on the local behavior at the fixed points of the discrete model. By examining these fixed points, we gain insights into the stability and potential bifurcations that may arise within the system. The study particularly emphasizes the absence of period-doubling bifurcations at various fixed points under specific conditions. For example, no period-doubling bifurcation occurs at liver-free fixed point (LFFP), disease-free fixed point (DFFP), and total infection fixed point (TIFP). Furthermore, the study proves that at the PIFP, the HCV model undergoes a N-S bifurcation. Additionally, no period-doubling bifurcation exists at the PIFP. To manage the chaos induced by the Neimark-Sacker bifurcation, the OGY control strategy is applied, showcasing the practical implications of controlling chaotic behavior in the HCV model (1.6). Numerical simulations confirm the theoretical results, reinforcing the validity and applicability of the discrete HCV model (1.6) in understanding the dynamics of Hepatitis C virus. The discrete HCV model is crucial for comprehending the complexities of HCV dynamics. Local stability analysis at various fixed points, such as the liver-free fixed point, disease-free fixed point, total infection fixed point, and partial infection fixed point, provides a comprehensive understanding of the virus's behavior under different conditions. The study's exploration of Neimark-Sacker bifurcation at the partial infection fixed point offers significant insights into the potential for chaotic dynamics in the virus's spread. Understanding the existence of chaos in the model underscores the importance of this research in devising effective control strategies for Hepatitis C. It is noted here that the existence of a Neimark-Sacker bifurcation at a partial infection fixed point in a discrete HCV model (1.6) indicates a transition from stable periodic behavior to quasi-periodic or chaotic dynamics as a key parameter crosses a critical threshold. Biologically, this implies that the HCV infection dynamics can shift from predictable oscillations of viral load and immune response to more complex and potentially erratic patterns. Such a bifurcation may correspond to scenarios where the immune system's response becomes less synchronized with the viral replication cycle, leading to irregular flare-ups and remissions of the infection. Understanding this transition is crucial for effective treatment strategies, as it suggests that minor changes in treatment or the patient's immune status could lead to significant changes in infection dynamics. Chaos control techniques, such as parameter adjustments or therapeutic interventions, are essential to stabilize the system, preventing chaotic outbreaks and maintaining the infection at manageable levels. This understanding provides insights into optimizing treatment regimens to avoid triggering complex and less controllable infection states. Future work can build upon this study by exploring more complex models and control strategies. Previous research, such as Blackmore et al.'s ^[29,30] work on discrete dynamical modeling and analysis of the R-S flip-flop circuit and their study on sigma map dynamics and bifurcations, provide a foundation for further investigations. Additionally, Rahman and Blackmore's ^[31] examination of bifurcations in walking droplet dynamics and Joshi et al.'s ^[32] exploration of generalized attracting horseshoes and chaotic strange attractors, suggest avenues for future research. Studies like Murthy et al. ^[33] on generalized attracting horseshoes in the Rössler attractor also offer valuable insights that can be integrated into future models of HCV dynamics. These directions can enhance our understanding of the complex behavior of the HCV model and improve strategies for managing and controlling the spread of Hepatitis C.

Author contributions

Abdul Qadeer Khan: Conceptualization, formal Analysis, investigation, methodology, resources, software, supervision, validation, visualization, writing-original draft, writing-review and editing; Ayesha Yaqoob: Formal Analysis, investigation, writing-original draft, writing-review and editing; Ateq Alsaadi: Conceptualization, funding acquisition, resources, writing-original draft, writing-review and editing. All authors have read and approved the final version of the manuscript for publication.

Acknowledgements

The authors extend their appreciation to Taif University, Saudi Arabia, for supporting this work through project number (TU-DSPP-2024-259).

Conflict of interest

The authors declare that they have no conflict of interest.

Appendix

A.   Values of $\lambda_{1,2}$ at (3.21) and (3.22)