Bone remodeling and biological effects of mechanical stimulus

Chao Hu; Qing-Hua Qin; Chao Hu; Qing-Hua Qin

doi:10.3934/bioeng.2020002

AIMS Bioengineering

2020, Volume 7, Issue 1: 12-28. doi: 10.3934/bioeng.2020002

Previous Article Next Article

Review

Bone remodeling and biological effects of mechanical stimulus

Chao Hu ,
Qing-Hua Qin ^,

Research School of Electrical, Energy and Materials Engineering, The Australian National University, Canberra, ACT 2601, Australia

Received: 05 November 2019 Accepted: 17 January 2020 Published: 19 January 2020

This review describes the physiology of normal bone tissue as an introduction to the subsequent discussion on bone remodeling and biomechanical stimulus. As a complex architecture with heterogeneous and anisotropic hierarchy, the skeletal bone has been anatomically analysed with different levelling principles, extending from nano- to the whole bone scale. With the interpretation of basic bone histomorphology, the main compositions in bone are summarized, including various organic proteins in the bone matrix and inorganic minerals as the reinforcement. The cell populations that actively participate in the bone remodeling—osteoclasts, osteoblasts and osteocytes—have also been discussed since they are the main operators in bone resorption and formation. A variety of factors affect the bone remodeling, such as hormones, cytokines, mechanical stimulus and electromagnetic stimulus. As a particularly potent stimulus for bone cells, mechanical forces play a crucial role in enhancing bone strength and preventing bone loss with aging. By combing all these aspects together, the information lays the groundwork for systematically understanding the link between bone physiology and orchestrated process of mechanically mediated bone homoestasis.

Keywords:

Citation: Chao Hu, Qing-Hua Qin. Bone remodeling and biological effects of mechanical stimulus[J]. AIMS Bioengineering, 2020, 7(1): 12-28. doi: 10.3934/bioeng.2020002

Related Papers:

[1]	Abdul Qadeer Khan, Ayesha Yaqoob, Ateq Alsaadi . Discrete Hepatitis C virus model with local dynamics, chaos and bifurcations. AIMS Mathematics, 2024, 9(10): 28643-28670. doi: 10.3934/math.20241390
[2]	Yao Shi, Zhenyu Wang . Bifurcation analysis and chaos control of a discrete fractional-order Leslie-Gower model with fear factor. AIMS Mathematics, 2024, 9(11): 30298-30319. doi: 10.3934/math.20241462
[3]	Xiongxiong Du, Xiaoling Han, Ceyu Lei . Dynamics of a nonlinear discrete predator-prey system with fear effect. AIMS Mathematics, 2023, 8(10): 23953-23973. doi: 10.3934/math.20231221
[4]	A. Q. Khan, Ibraheem M. Alsulami, S. K. A. Hamdani . Controlling the chaos and bifurcations of a discrete prey-predator model. AIMS Mathematics, 2024, 9(1): 1783-1818. doi: 10.3934/math.2024087
[5]	Abdul Qadeer Khan, Zarqa Saleem, Tarek Fawzi Ibrahim, Khalid Osman, Fatima Mushyih Alshehri, Mohamed Abd El-Moneam . Bifurcation and chaos in a discrete activator-inhibitor system. AIMS Mathematics, 2023, 8(2): 4551-4574. doi: 10.3934/math.2023225
[6]	Yulong Li, Long Zhou, Fengjie Geng . Dynamics on semi-discrete Mackey-Glass model. AIMS Mathematics, 2025, 10(2): 2771-2807. doi: 10.3934/math.2025130
[7]	Ahmad Suleman, Rizwan Ahmed, Fehaid Salem Alshammari, Nehad Ali Shah . Dynamic complexity of a slow-fast predator-prey model with herd behavior. AIMS Mathematics, 2023, 8(10): 24446-24472. doi: 10.3934/math.20231247
[8]	A. Q. Khan, Ibraheem M. Alsulami . Discrete Leslie's model with bifurcations and control. AIMS Mathematics, 2023, 8(10): 22483-22506. doi: 10.3934/math.20231146
[9]	Xiaoming Su, Jiahui Wang, Adiya Bao . Stability analysis and chaos control in a discrete predator-prey system with Allee effect, fear effect, and refuge. AIMS Mathematics, 2024, 9(5): 13462-13491. doi: 10.3934/math.2024656
[10]	Figen Kangalgil, Seval Isșık . Effect of immigration in a predator-prey system: Stability, bifurcation and chaos. AIMS Mathematics, 2022, 7(8): 14354-14375. doi: 10.3934/math.2022791

Abstract

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

$\begin{equation} \begin{aligned} \dot{T}& = \lambda- \beta_1 T V-d_1 T-\beta_2T I+\alpha I, \\ \dot{I}& = \beta_1 T V-\alpha I-d_2I+\beta_2TI,\\ \dot{V}& = k I-p V Z-d_3V,\\ \dot{Z}& = cZ(t-\tau) V(t-\tau)-d_4 Z, \end{aligned} \end{equation}$

(1.1)

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

$\begin{equation} \begin{aligned} \dot{T}& = s-\beta T V-d_1 T, \\ \dot{I}& = \beta T V-d_2 I,\\ \dot{V}& = k I-u V, \end{aligned} \end{equation}$

(1.2)

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

$\begin{eqnarray} \begin{split} \dot{T}& = s+rT\left(1-\frac{T+I}{T_{\max}}\right)-d T-(1-\eta)\beta V T, \\ \dot{I}& = (1-\eta)\beta V T+rT\left(1-\frac{T+I}{T_{\max}}\right)-\delta I,\\ \dot{V}& = (1-\epsilon)p I-c V, \end{split} \end{eqnarray}$

(1.3)

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

$\begin{equation} \begin{aligned} \dot{u}& = u(1-u-v)-g u-(1-\alpha)b u v, \\ \dot{v}& = v(1-u-v)-\zeta u+(1-\alpha)b u v. \end{aligned} \end{equation}$

(1.4)

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

$\begin{equation} \begin{aligned} \frac{u_{t+1}-u_t}{h}& = u_t-u_tu_{t+1}-u_tv_t-g u_t-(1-\alpha)b u_t v_t, \\ \frac{v_{t+1}-v_t}{h}& = v_t-u_tv_t-v_tv_{t+1}-\zeta u_t+(1-\alpha)b u_t v_t, \end{aligned} \end{equation}$

(1.5)

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

$\begin{equation} \begin{aligned} u_{t+1}& = \frac{(1+h(1-g))u_t-h(1+(1-\alpha)b )u_t v_t}{1+h u_t }, \\ v_{t+1}& = \frac{(1+h(1-\zeta))v_t-h(1-(1-\alpha)b )u_t v_t}{1+h v_t}. \end{aligned} \end{equation}$

(1.6)

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

$\begin{equation} \begin{split} \hat{u}& = \frac{(1+h(1-g))\hat{u}-h(1+(1-\alpha)b )\hat{u} \hat{v}}{1+h \hat{u}}, \\ \hat{v}& = \frac{(1+h(1-\zeta))\hat{v}-h(1-(1-\alpha)b )\hat{u} \hat{v}}{1+h \hat{v}}. \end{split} \end{equation}$

(2.1)

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

$\begin{equation} \begin{split} -\hat{u}+(-1-(1-\alpha)b)\hat{v}+1-g = 0, \\ -\hat{v}+(-1+(1-\alpha)b)\hat{u}+1-\zeta = 0. \end{split} \end{equation}$

(2.2)

From the first equation of (2.2), one gets

$\begin{equation} \hat{v} = \frac{-\hat{u}+1-g}{1+(1-\alpha) b}. \end{equation}$

(2.3)

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

$\begin{equation} \hat{u} = \frac{\zeta-g-(1-\alpha)(1-\zeta)b}{(1-\alpha)^2 b^2}. \end{equation}$

(2.4)

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

$\begin{equation} \hat{v} = \frac{g-\zeta+(1-\alpha)(1-g)b}{(1-\alpha)^2 b^2}. \end{equation}$

(2.5)

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

$\begin{equation} J|_{(\hat{u},\hat{v})}: = \left( \begin{array}{cc} \frac{1-h(-1+g+(1+b(1-\alpha))\hat{v})}{(1+h\hat{u})^2} & \frac{-h(1+b(1-\alpha))\hat{u}}{1+h \hat{u}}\\ \frac{-h(1-b(1-\alpha))\hat{v}}{1+h \hat{v}} & \frac{1-h(-1+\zeta+(1-b(1-\alpha))\hat{u})}{(1+h\hat{v})^2} \end{array} \right), \end{equation}$

(3.1)

where

$\begin{equation} \begin{aligned} g_1 = \frac{(1+h(1-g))u_t-h(1+(1-\alpha)b )u_t v_t}{1+h u_t }, \\ g_2 = \frac{(1+h(1-\zeta))v_t-h(1-(1-\alpha)b )u_t v_t}{1+h v_t }. \end{aligned} \end{equation}$

(3.2)

Now, at LFFP, (3.1) becomes

$\begin{equation} J|_{\mathrm{LFFP}}: = \left( \begin{array}{cc} 1+h-g h & 0\\ 0 & 1+h-h\zeta \end{array} \right), \end{equation}$

(3.3)

with

$\begin{equation} \lambda_1 = 1+h-g h, \ \lambda_2 = 1+h-h\zeta. \end{equation}$

(3.4)

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

$\begin{equation} h = \frac{2}{g-1}, \end{equation}$

(3.5)

$\begin{equation} h = \frac{2}{\zeta-1}. \end{equation}$

(3.6)

At DFFP, (3.1) becomes

$\begin{equation} J|_{\mathrm{DFFP}}: = \left( \begin{array}{cc} \frac{1}{1+h(1-g)} & \frac{h(-1-(1-\alpha)b)(1-g)}{1+h(1-g)}\\ 0 & 1+h( g-\zeta +(1-\alpha)(1-g)b) \end{array} \right), \end{equation}$

(3.7)

with

$\begin{equation} \lambda_1 = \frac{1}{1+h-hg}, \ \lambda_2 = 1+h( g-\zeta +(1-\alpha)(1-g)b). \end{equation}$

(3.8)

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

$\begin{equation} h = \frac{2}{g-1}, \end{equation}$

(3.9)

$\begin{equation} h = \frac{2}{-g+\zeta-(1-\alpha)(1-g)b}. \end{equation}$

(3.10)

Furthermore, at TIFP, (3.1) becomes

$\begin{equation} J|_{\mathrm{TIFP}}: = \left( \begin{array}{cc} 1+h(\zeta-g -(1-\alpha)(1-\zeta)b) & 0\\ \frac{h(-1+(1-\alpha)b)(1-\zeta)}{1+h(1-\zeta)}& \frac{1}{1+h(1-\zeta)} \end{array} \right), \end{equation}$

(3.11)

with

$\begin{equation} \lambda_1 = 1+h(\zeta-g)-h(1-\alpha)(1-\zeta)b, \ \lambda_2 = \frac{1}{1+h(1-\zeta)}, \end{equation}$

(3.12)

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

$\begin{equation} h = \frac{2}{\zeta-1}, \end{equation}$

(3.13)

$\begin{equation} h = \frac{2}{-\zeta+g+(1-\alpha)(1-\zeta)b}. \end{equation}$

(3.14)

Finally, at PIFP, (3.1) becomes

$\begin{equation} J|_{\mathrm{PIFP}}: = \left( \begin{array}{cc} \frac{(1-\alpha)^2b^2}{b^2(1-\alpha)^2-h(g-\zeta)-hb(1-\alpha)(1-\zeta)}& \frac{h(b(1-\alpha)+1)(g+b(1-\alpha)(1-\zeta)-\zeta)}{(1-\alpha)^2b^2-h(g-\zeta)-hb(1-\alpha)(1-\zeta)}\\ \frac{h(b(1-\alpha)-1)(g+b(1-\alpha)(1-g)-\zeta)}{(1-\alpha)^2b^2+h(g-\zeta)+hb(1-\alpha)(1-g)}& \frac{(1-\alpha)^2b^2}{(1-\alpha)^2b^2+h(g-\zeta)+hb(1-\alpha)(1-g)} \end{array} \right), \end{equation}$

(3.15)

with

$\begin{equation} \lambda^2-\varrho_1\lambda+\varrho_2 = 0, \end{equation}$

(3.16)

where

$\begin{eqnarray} \begin{aligned} \varrho_1& = \frac{b^3(1-\alpha)^3(2b(1-\alpha)-h(g-\zeta))}{((1-\alpha)^2b^2-h(g-\zeta)-hb(1-\alpha)(1-\zeta))((1-\alpha)^2b^2+h(g-\zeta)+hb(1-\alpha)(1-g))},\\ \varrho_2& = \frac{b^4(1-\alpha)^4+h^2(1-b^2(1-\alpha)^2)(g+b(1-g)(1-\alpha)-\zeta)(g+b(1-\zeta) (1-\alpha)-\zeta)}{((1-\alpha)^2b^2-h(g-\zeta)-hb(1-\alpha)(1-\zeta))((1-\alpha)^2b^2+h(g-\zeta)+hb(1-\alpha)(1-g))}. \end{aligned} \end{eqnarray}$

(3.17)

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

$\begin{eqnarray} \lambda_{1,2}& = \frac{\varrho_1\pm\sqrt{\Delta}}{2}, \end{eqnarray}$

(3.18)

where

$\begin{eqnarray} \begin{aligned} \Delta = &\varrho^2_1-4\varrho_2\\ = &\left(\frac{b^3(1-\alpha)^3(2b(1-\alpha)-h(g-\zeta))}{\left\{{((1-\alpha)^2b^2-h(g-\zeta)- h b(1-\zeta)(1-\alpha))\times}{((1-\alpha)^2b^2+h(g-\zeta)+h b(1-\alpha)(1-g))}\right\}} \right)^2-\\&4\left(\frac{\left\{{b^4(1-\alpha)^4+h^2(1-b^2(1-\alpha)^2)(g+b(1-g) (1-\alpha)-\zeta)\times}{(g+b(1-\alpha)(1-\zeta)-\zeta)}\right\}}{\left\{{ ((1-\alpha)^2b^2-h(g-\zeta)-h b(1-\alpha)(1-\zeta))\times}{(b^2(1-\alpha)^2+h(g-\zeta)+ h b(1-\alpha)(1-g))}\right\}}\right). \end{aligned} \end{eqnarray}$

(3.19)

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

$\begin{equation} 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)}. \end{equation}$

(3.20)

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

$\begin{eqnarray*} \begin{aligned} h > \max\left\{\frac{4b\left(1-\alpha\right)}{g-\zeta-\sqrt{\left\{{5\left(g-\zeta\right)^2+4b^2\left(1-\alpha\right)^2\left(1-g\right) \left(1-\zeta\right)-}{4b\left(1-\alpha\right)\left(g-\zeta\right)\left(-2+g+\zeta\right)}\right\}}},\right.\\ \left. {}\frac{4b\left(1-\alpha\right)}{g-\zeta+\sqrt{\left\{{5\left(g-\zeta\right)^2+4b^2\left(1-\alpha\right)^2\left(1-g\right) \left(1-\zeta\right)-}{4b\left(1-\alpha\right)\left(g-\zeta\right)\left(-2+g+\zeta\right)}\right\}}}\right\}, \end{aligned} \end{eqnarray*}$

a stable node if

$\begin{eqnarray*} \begin{aligned} 0 < h < \min\left\{\frac{4b\left(1-\alpha\right)}{g-\zeta-\sqrt{\left\{{5\left(g-\zeta\right)^2+ 4b^2\left(1-\alpha\right)^2\left(1-g\right)\left(1-\zeta\right)-}{4b\left(1-\alpha\right)^2\left(g-\zeta\right)\left(-2+g+\zeta\right)}\right\}}},\right.\\ \left. {} \frac{4b\left(1-\alpha\right)}{g-\zeta+\sqrt{\left\{{5\left(g-\zeta\right)^2+4b^2\left(1-\alpha\right)^2\left(1-g\right) \left(1-\zeta\right)-}{4b\left(1-\alpha\right)\left(g-\zeta\right)\left(-2+g+\zeta\right)}\right\}}}\right\}, \end{aligned} \end{eqnarray*}$

and non-hyperbolic if

$\begin{equation} h = \frac{4b(1-\alpha)}{g-\zeta-\sqrt{\left\{{5(g-\zeta)^2+4b^2(1-\alpha)^2(1-g)(1-\zeta)-}{4b(1-\alpha)(g-\zeta)(-2+g+\zeta)}\right\}}}, \end{equation}$

(3.21)

$\begin{equation} h = \frac{4b(1-\alpha)}{g-\zeta+\sqrt{\left\{{5(g-\zeta)^2+4b^2(1-\alpha)^2(1-g)(1-\zeta)-}{4b(1-\alpha)(g-\zeta)(-2+g+\zeta)}\right\}}}. \end{equation}$

(3.22)

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.

Figure 1. Dynamics of PIFP with

$\alpha\in(0, 0.5)$ and

$h\in(0, 0.9)$ .

DownLoad: Full-Size Img PowerPoint

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.

Figure 2. Dynamics of PIFP with

$\alpha\in(0, 1.8)$ and

$h\in(0, 2.9)$ .

DownLoad: Full-Size Img PowerPoint

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

$\begin{equation} \mathcal{F}|_{\mathrm{LFFP}} = \left\{\mathfrak{L}: \ h = \frac{2}{g-1} \right\}; \end{equation}$

(4.1)

(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

$\begin{equation} \mathcal{F}|_{\mathrm{LFSS}} = \left\{\mathfrak{L}: \ h = \frac{2}{\zeta-1} \right\}; \end{equation}$

(4.2)

(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

$\begin{equation} \mathcal{F}|_{\mathrm{DFFP}} = \left\{\mathfrak{L}: \ h = \frac{2}{g-1} \right\}; \end{equation}$

(4.3)

(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

$\begin{equation} \mathcal{F}|_{\mathrm{DFFP}} = \left\{\mathfrak{L}: \ h = \frac{2}{-g+\zeta-(1-\alpha)(1-g)b} \right\}; \end{equation}$

(4.4)

(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

$\begin{equation} \mathcal{F}|_{\mathrm{TIFP}} = \left\{\mathfrak{L}: \ h = \frac{2}{\zeta-1} \right\}; \end{equation}$

(4.5)

$(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

$\begin{equation} \mathcal{F}|_{\mathrm{TIFP}} = \left\{\mathfrak{L}: \ h = \frac{2}{-\zeta+g+(1-\zeta)(1-\alpha)b} \right\}; \end{equation}$

(4.6)

(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

$\begin{equation} \begin{split} \mathcal{N}|_{\mathrm{PIFP}} = \left\{\mathfrak{L}: \Delta < 0\ \mathrm{and}\ h = \frac{(1-\alpha)^3b^3(\zeta-g)}{\left\{{(2-b^2(1-\alpha)^2)(g+b(1-g)(1-\alpha)-\zeta)}{(g+b(1-\zeta)(1-\alpha)-\zeta)}\right\}}\right\}; \end{split} \end{equation}$

(4.7)

(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

$\begin{equation} \mathcal{F}|_{\mathrm{PIFP}} = \left\{\mathfrak{L}: \Delta > 0\ \mathrm{and}\ h = \frac{4b(1-\alpha)}{g-\zeta- \sqrt{\left\{{5(g-\zeta)^2+4b^2(1-\alpha)^2(1-g)(1-\zeta)-}{4b(1-\alpha)(g-\zeta)(-2+g+\zeta)}\right\}}}\right\}; \end{equation}$

(4.8)

(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

$\begin{equation} \mathcal{F}|_{\mathrm{PIFP}} = \left\{\mathfrak{L}: \Delta > 0\ \mathrm{and}\ h = \frac{4b(1-\alpha)}{g-\zeta+\sqrt{\left\{{5(g-\zeta)^2+4b^2(1-\alpha)^2(1-g)(1-\zeta)-}{4b(1-\alpha)(g-\zeta)(-2+g+\zeta)}\right\}}}\right\}. \end{equation}$

(4.9)

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$ ,

$\begin{equation} u_{t+1} = \frac{(1+h(1-g))u_t}{1+hu_t}. \end{equation}$

(4.10)

From (4.10), we define

$\begin{equation} g_1: = \frac{(1+h(1-g))u}{1+hu}. \end{equation}$

(4.11)

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

$\begin{equation} \frac{\partial g_1}{\partial u}\biggr|_{u = \hat{u} = 1-g, \ h = \hat{h} = \frac{2}{g-1}} = -1, \end{equation}$

(4.12)

$\begin{equation} \frac{\partial^2 g_1}{\partial u^2}\biggr|_{u = \hat{u} = 1-g,\ h = \hat{h} = \frac{2}{g-1}} = \frac{-4}{g-1}, \end{equation}$

(4.13)

and

$\begin{equation} \frac{\partial g_1}{\partial h}\biggr|_{u = \hat{u} = 1-g,\ h = \hat{h} = \frac{2}{g-1}} = 0. \end{equation}$

(4.14)

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

$\begin{equation} \begin{aligned} u_{t+1} = \frac{(1+(h^*+\epsilon)(1-g))u_t-(h^*+\epsilon)(1+(1-\alpha)b )u_t v_t}{1+(h^*+\epsilon) u_t}, \\ v_{t+1} = \frac{(1+(h^*+\epsilon)(1-\zeta))v_t-(h^*+\epsilon)(1-(1-\alpha)b )u_t v_t}{1+(h^*+\epsilon) v_t}, \end{aligned} \end{equation}$

(4.15)

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

$\begin{eqnarray} \lambda_{1,2} = \frac{\varrho_1(\epsilon)\pm\iota\sqrt{4\varrho_2(\epsilon)-\varrho_1^2(\epsilon)}}{2}, \end{eqnarray}$

(4.16)

where

$\begin{eqnarray} \begin{split} \varrho_1(\epsilon)& = \frac{b^3(1-\alpha)^3(2b(1-\alpha)-(h^*+\epsilon)(g-\zeta))}{\left\{{((1-\alpha)^2b^2-(h^*+\epsilon)(g-\zeta) -(h^*+\epsilon)b(1-\alpha)(1-\zeta))}{((1-\alpha)^2b^2+(h^*+\epsilon)(g-\zeta)+(h^*+\epsilon)b(1-\alpha)(1-g))}\right\}},\\ \varrho_2(\epsilon)& = \frac{\left\{{b^4(1-\alpha)^4+(h^*+\epsilon)^2(1-b^2(1-\alpha)^2)}{(g+b(1-g)(1-\alpha)-\zeta)(g+b(1-\zeta) (1-\alpha)-\zeta)}\right\}}{\left\{{((1-\alpha)^2b^2-(h^*+\epsilon)(g-\zeta)-(h^*+\epsilon)b(1-\alpha)(1-\zeta))}{((1-\alpha)^2b^2+ (h^*+\epsilon)(g-\zeta)+(h^*+\epsilon)b(1-\alpha)(1-g))}\right\}}. \end{split} \end{eqnarray}$

(4.17)

From (4.16) and (4.17), one has

$\begin{eqnarray} \begin{aligned} \frac{d|\lambda_{1,2}|}{d \epsilon}\biggr|_{\epsilon = 0} = &\frac{\left\{{-b^4 h (1 - \alpha)^4 (g + b (1 - \alpha) (1 - \zeta) - \zeta) (-2 g -b (1 -\alpha)\times}{(2 - g - b (1 - \alpha) (g + b (1- g) (1 - \alpha) - \zeta) - \zeta) + 2 \zeta)}\right\}} {\left\{ \begin{array}{l} {{\left( {{b^4}{{(1 - \alpha )}^4} - {h^2}{{(g - \zeta )}^2} - 2b{h^2}(1 - \alpha )(g - \zeta )(1 - \zeta ) - {b^2}{h^2}{{(1 - \alpha )}^2}{{(1 - \zeta )}^2}} \right)}^{3/2} \times } \\ \sqrt {\left\{ {{b^4}{{(1 - \alpha )}^4} - {h^2}( - 1 + {{(1 - \alpha )}^2}{b^2})(g + b(1 - \alpha )(1 - g) - \zeta )(g + (1 - \alpha )b(1 - \zeta ) - \zeta )} \right\}} \end{array} \right\}} \\&\not = 0. \end{aligned} \end{eqnarray}$

(4.18)

Now, (4.15) becomes

$\begin{eqnarray} \begin{aligned} I_{t+1} = &\frac{\left(1+(h^*+\epsilon)(1-g)\right)(I_{t}+\hat{u})-(h^*+\epsilon)(1+(1-\alpha)b) (I_t+\hat{u})(S_t+\hat{v})}{1+(h^*+\epsilon)(I_t+\hat{u})}-\hat{u},\\ S_{t+1} = &\frac{\left(1+(h^*+\epsilon)(1-\zeta)\right)(S_{t}+\hat{v})-(h^*+\epsilon) (1-(1-\alpha)b)(I_t+\hat{u})(S_t+\hat{v})}{1+(h^*+\epsilon)(S_t+\hat{v})}-\hat{v}, \end{aligned} \end{eqnarray}$

(4.19)

$\begin{equation} \begin{aligned} I_t = &u_t-\hat{u},\\ S_t = &v_t-\hat{v}, \end{aligned} \end{equation}$

(4.20)

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

$\begin{eqnarray} \begin{aligned} I_{t+1} = &l_{11}I_t+l_{12}S_t+l_{13} I_t^2+l_{14} I_t S_t, \\ S_{t+1} = &l_{21}I_t+l_{22}S_t+l_{23}I_t S_t+l_{24}S_t^2, \end{aligned} \end{eqnarray}$

(4.21)

where

$\begin{eqnarray} \begin{aligned} l_{11} = &\frac{1+h(1-g)-h(1+(1-\alpha)b)\hat{v}}{(1+h\hat{u})^2}, \\ l_{12} = &\frac{-h(1+(1-\alpha)b)\hat{u}}{1+h\hat{u}}, \\ l_{13} = &\frac{-2h(1+h(1-g)-h(1+(1-\alpha)b)\hat{v})}{(1+h\hat{u})^3}, \\ l_{14} = &\frac{h(-1-(1-\alpha)b)}{(1+h\hat{u})^2} ,\\ l_{21} = &\frac{h(-1+(1-\alpha)b)\hat{v}}{1+h\hat{v}}, \\ l_{22} = &\frac{1+h(1-\zeta)-h(1-(1-\alpha)b)\hat{u}}{(1+h\hat{v})^2}, \\ l_{23} = &\frac{h(-1+(1-\alpha)b)}{(1+h\hat{v})^2},\\ l_{24} = &\frac{-2h(1+h(1-\zeta)-h(1-(1-\alpha)b)\hat{u})}{(1+h\hat{v})^3}. \end{aligned} \end{eqnarray}$

(4.22)

Now, if

$\begin{equation} \left( \begin{array}{c} I_t \\ S_t \\ \end{array} \right): = \left( \begin{array}{cc} l_{12}&0 \\ \eta-l_{11}&-\rho \\ \end{array} \right)\left( \begin{array}{c} {u_t} \\ {v_t} \\ \end{array} \right), \end{equation}$

(4.23)

then (4.21) becomes

$\begin{eqnarray} \begin{split} \left( \begin{array}{c} u_{t+1} \\ v_{t+1} \\ \end{array} \right): = \left( \begin{array}{cc} \eta&-\rho \\ \rho &\eta \\ \end{array} \right)\left( \begin{array}{c} u_{t} \\ v_{t} \\ \end{array} \right)+\left( \begin{array}{c} \Upsilon_1(u_{t},v_{t}) \\ \Upsilon_2(u_{t},v_{t}) \\ \end{array} \right), \end{split} \end{eqnarray}$

(4.24)

where

$\begin{eqnarray} \begin{aligned} \eta = &\frac{b^3(1-\alpha)^3(2b(1-\alpha)-h(g-\zeta))}{2((1-\alpha)^2b^2-h(g-\zeta)- h b(1-\alpha)(1-\zeta))((1-\alpha)^2b^2+h(g-\zeta)+hb(1-\alpha)(1-g))},\\ \rho = &\frac{{\sqrt{\left\{ \begin{array}{l} -b^6 (1 - \alpha)^6 (-2(1 - \alpha) b + h (g - \zeta))^2 + 4 ((1 - \alpha)^2b^2 + h (g +\\ \;\;\;\; b(1 - \alpha) (1 - g) - \zeta)) (b^4 (1 - \alpha)^4 - h^2 (-1 + b^2 (1 -\alpha)^2)\times\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;(g + b (1 - g) (1 - \alpha) - \zeta) (g + b (1 - \alpha) (1 - \zeta) -\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\zeta)) (b^2 (1 - \alpha)^2 + \zeta)) \end{array} \right\}}}} {2((1-\alpha)^2b^2-h(g-\zeta)-hb(1-\zeta)(1-\alpha))((1-\alpha)^2b^2+h(g-\zeta)+hb(1-\alpha)(1-g))}, \end{aligned} \end{eqnarray}$

(4.25)

and

$\begin{eqnarray} \begin{aligned} \Upsilon_1& = \sigma_{11} u_t^2+\sigma_{12}u_t v_t,\\ \Upsilon_2& = \sigma_{21}u_t^2+\sigma_{22}u_t v_t+\sigma_{23}v_t^2, \end{aligned} \end{eqnarray}$

(4.26)

with

$\begin{eqnarray} \begin{aligned} \sigma_{11}& = l_{13}l_{12}+l_{14}\eta-l_{14}l_{11},\\ \sigma_{12}& = -\rho l_{14},\\ \sigma_{21}& = \frac{(\eta-l_{11})l_{12}}{\rho}\left(l_{13}-l_{23} \right)+\frac{(\eta-l_{11})^2}{\rho}\left(l_{14}-l_{24} \right),\\ \sigma_{22}& = l_{23}l_{12}+(\eta-l_{11})(2l_{24}-l_{14}), \\ \sigma_{23}& = -l_{24}\rho. \end{aligned} \end{eqnarray}$

(4.27)

From (4.26), we get

$\begin{eqnarray} \begin{aligned} \frac{\partial^2\Upsilon_1}{\partial u_t^2}\biggr|_{(0,0)}& = 2 \sigma_{11},\ \frac{\partial^2\Upsilon_1}{\partial{u_t}{\partial v_t}}\biggr|_{(0,0)} = \sigma_{12},\\ \frac{\partial^2\Upsilon_1}{\partial{v_t^2}}\biggr|_{(0,0)}& = \frac{\partial^3\Upsilon_1}{\partial{u_t^3}}\biggr|_{(0,0)} = \frac{\partial^3\Upsilon_1}{\partial{u_t^2}\partial{v_t}}\biggr|_{(0,0)} = \frac{\partial^3\Upsilon_1}{\partial{u_t}\partial{v_t^2}}\biggr|_{(0,0)} = \frac{\partial^3\Upsilon}{\partial{v_t^3}}\biggr|_{(0,0)} = 0,\\ \frac{\partial^2\Upsilon_2}{\partial{u_t^2}}\biggr|_{(0,0)}& = 2 \sigma_{21},\ \frac{\partial^2\Upsilon_2}{\partial{u_t}\partial{v_t}}\biggr|_{(0,0)} = \sigma_{22}, \frac{\partial^2\Upsilon_2}{\partial{v_t^2}}\biggr|_{(0,0)} = 2\sigma_{23},\\ \frac{\partial^3\Upsilon_2}{\partial{u_t^3}}\biggr|_{(0,0)}& = \frac{\partial^3\Upsilon_2}{\partial{u_t^2}\partial{v_t}}\biggr|_{(0,0)} = \frac{\partial^3\Upsilon_2}{\partial{u_t}\partial{v_t^2}}\biggr|_{(0,0)} = \frac{\partial^3\Upsilon_2}{\partial{v_t^3}}\biggr|_{(0,0)} = 0. \end{aligned} \end{eqnarray}$

(4.28)

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

$\begin{eqnarray} \psi = -\Re\left(\frac{(1-2\bar{\lambda})\bar{\lambda}^2}{1-\lambda}\rho_{11}\rho_{20}\right)-\frac{1}{2}\|\rho_{11}\|^2-\|\rho_{02}\|^2+ \Re\left(\bar{\lambda}\rho_{21}\right), \end{eqnarray}$

(4.29)

where

$\begin{eqnarray} \begin{aligned} \rho_{02}& = \frac{1}{8}\left(\frac{\partial^2\Upsilon_1}{\partial{u_t^2}}-\frac{\partial^2\Upsilon_1}{\partial{v_t^2}}+ 2\frac{\partial^2\Upsilon_2}{\partial{u_t}\partial{v_t}}+\iota\left(\frac{\partial^2\Upsilon_2}{\partial{u_t^2}}-\frac{\partial^2\Upsilon_2}{\partial{v_t^2}}+ 2\frac{\partial^2\Upsilon_1}{\partial{u_t}\partial{v_t}}\right)\right)\biggr|_{\left(0,0\right)}, \\ \rho_{11}& = \frac{1}{4}\left(\frac{\partial^2\Upsilon_1}{\partial{u_t^2}}+\frac{\partial^2\Upsilon_1}{\partial{v_t^2}}+ \iota\left(\frac{\partial^2\Upsilon_2}{\partial{u_t^2}}+\frac{\partial^2\Upsilon_2}{\partial{v_t^2}}\right)\right)\biggr|_{\left(0,0\right)},\\ \rho_{20}& = \frac{1}{8}\left(\frac{\partial^2\Upsilon_1}{\partial{u_t^2}}-\frac{\partial^2\Upsilon_1}{\partial{v_t^2}}+ 2\frac{\partial^2\Upsilon_2}{\partial{u_t}\partial{v_t}}+\iota\left(\frac{\partial^2\Upsilon_2}{\partial{u_t^2}}-\frac{\partial^2\Upsilon_2}{\partial{v_t^2}}- 2\frac{\partial^2\Upsilon_1}{\partial{u_t}\partial{v_t}}\right)\right)\biggr|_{\left(0,0\right)},\\ \rho_{21}& = \frac{1}{16}\left(\frac{\partial^3\Upsilon_1}{\partial{u_t^3}}+\frac{\partial^3\Upsilon_1}{\partial{v_t^3}}+ \frac{\partial^3\Upsilon_2}{\partial{u_t^2}\partial{v_t}}+\frac{\partial^3\Upsilon_2}{\partial{v_t^3}}+ \right.\\ &\qquad \left. {} \iota\left(\frac{\partial^3\Upsilon_2}{\partial{u_t^3}}+\frac{\partial^3\Upsilon_2}{\partial{u_t} \partial{v_t^2}}-\frac{\partial^3\Upsilon_1}{\partial{u_t^2}\partial{v_t}}-\frac{\partial^3\Upsilon_1}{\partial{v_t^3}}\right) \right)\biggr|_{\left(0,0\right)}. \end{aligned} \end{eqnarray}$

(4.30)

After calculation, one gets

$\begin{eqnarray} \begin{aligned} \rho_{02}& = \frac{1}{4}\left(\sigma_{11}+\sigma_{22}+\iota (\sigma_{21}-\sigma_{23}+\sigma_{12})\right),\\ \rho_{11}& = \frac{1}{2}\left(\sigma_{11}+\iota (\sigma_{21}+\sigma_{23}) \right),\\ \rho_{20}& = \frac{1}{4}\left(\sigma_{11}+\sigma_{22}+\iota(\sigma_{21}-\sigma_{12}-\sigma_{23})\right),\\ \rho_{21}& = 0. \end{aligned} \end{eqnarray}$

(4.31)

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.,

$\begin{eqnarray} g \not = \frac{{\left\{ \begin{array}{l} {8 b (-1 + \alpha) + 8 b^3 (1 - \alpha)^3 -2 b^5 (1 - \alpha)^5 - b^6 (1 - \alpha)^6 (-1 + \zeta)+ 8 \zeta + b^4 (1 - \alpha)^4}\\ {(-4 + 5 \zeta) -2 b^2 (-1 + \alpha)^2 (-2 + 5 \zeta)\pm \sqrt{b^6 (-2 + b^2 (1 - \alpha)^2)^3 (1 - \alpha)^6 (1 - \zeta)^2}} \end{array} \right\}}} {8 - 8 b (1 - \alpha) - 6 b^2 (1 - \alpha)^2 + 8 b^3 (1 -\alpha)^3 + b^4 (1 -\alpha)^4 - 2 b^5 (1 - \alpha)^5}. \end{eqnarray}$

(4.32)

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

$\begin{eqnarray} \begin{aligned} u_{t+1}& = \frac{(1+h(1-g))u_t-h(1+(1-\alpha)b )u_t v_t}{1+h u_t} = g_1(u_t, v_t, \zeta), \\ v_{t+1}& = \frac{(1+h(1-\zeta))v_t-h(1-(1-\alpha)b )u_t v_t}{1+h v_t} = g_2(u_t, v_t, \zeta), \end{aligned} \end{eqnarray}$

(5.1)

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)

$\begin{equation} \left( \begin{array}{c} u_{t+1}-\hat{u} \\ v_{t+1}-\hat{v} \\ \end{array} \right)\approx J_{\left(\hat{u},\hat{v},\zeta_0\right)}\left( \begin{array}{c} {u_t-\hat{u}} \\ {{v_t-\hat{v}}} \\ \end{array} \right)+P(\zeta-\zeta_0), \end{equation}$

(5.2)

where

$\begin{equation} J_{\left(\hat{u},\hat{v},\zeta_0\right)} = \left( \begin{array}{cc} \frac{(1-\alpha)^2b^2}{(1-\alpha)^2b^2+h(\zeta_0-g)+hb(1-\alpha)(\zeta_0-1)}& \frac{h(1+b(1-\alpha))(g+b(1-\alpha)(1-\zeta_0)-\zeta_0)}{(1-\alpha)^2b^2+h(\zeta_0-g)+hb(1-\alpha)(\zeta_0-1)}\\ \frac{h(-1+b(1-\alpha))(g+b(1-\alpha)(1-g)-\zeta_0)}{(1-\alpha)^2b^2+h(g-\zeta_0)+ hb(1-\alpha)(1-g)}&\frac{(1-\alpha)^2b^2}{(1-\alpha)^2b^2+h(g-\zeta_0)+hb(1-\alpha)(1-g)} \end{array} \right), \end{equation}$

(5.3)

and

$\begin{equation} P = \left( \begin{array}{cc} \frac{\partial g_1(\hat{u},\hat{v},\zeta_0)}{\partial \zeta }\\ \frac{\partial g_2(\hat{u},\hat{v},\zeta_0)}{\partial \zeta} \end{array} \right) = \left( \begin{array}{cc} 0\\ \frac{h(-g+\zeta_0-(1-g)(1-\alpha)b)}{(1-\alpha)^2b^2+h(g-\zeta_0+(1-g)(1-\alpha)b)} \end{array} \right). \end{equation}$

(5.4)

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

$\begin{equation} C = \left( \begin{array}{cc} P:J_{\left(\hat{u},\hat{v},\zeta_0\right)}P\\ \end{array} \right) = \ \begin{aligned} & \left(\begin{matrix} 0\\ \\ \frac{h(-g+\zeta_0-(1-\alpha)(1-g)b)}{b^2(1-\alpha)^2+h(g-\zeta_0+(1-g)(1-\alpha)b)}\\ \\ \end{matrix}\right.\\ \\ &\qquad\qquad \left.\begin{matrix} \frac{h^2(g-\zeta_0+b(1-\zeta_0)(1-\alpha))(g-\zeta_0+(1-\alpha)(1-g)b)(-1-b(1-\alpha))}{((1-\alpha)^2b^2+ h(\zeta_0-g)+hb(1-\zeta_0)(1-\alpha))((1-\alpha)^2b^2+h(g-\zeta_0)+hb(1-\alpha)(1-g))}\\ \\ \frac{hb^2(1-\alpha)^2(-g+\zeta_0-(1-\alpha)(1-g)b)}{((1-\alpha)^2b^2+h(g-\zeta_0+(1-\alpha)(1-g)b))^2} \end{matrix}\right), \end{aligned} \end{equation}$

(5.5)

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

$\begin{equation} \left( \begin{array}{c} u_{t+1}-\hat{u} \\ v_{t+1}-\hat{v} \\ \end{array} \right)\approx \left(J_{\left(\hat{u},\hat{v},\zeta_0\right)}-PT\right)\left( \begin{array}{c} {u_t-\hat{u}} \\ {{v_t-\hat{v}}} \\ \end{array} \right). \end{equation}$

(5.6)

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

$\begin{eqnarray} \begin{aligned} u_{t+1}& = \frac{(1+h(1-g))u_t-h(1+(1-\alpha)b )u_t v_t}{1+h u_t}, \\ v_{t+1}& = \frac{(1+h(1-(\zeta_0-\hat{k_{1}}(u_{t}-\hat{u})-\hat{k_{2}}(v_{t}-\hat{v}))))v_t-h(1-(1-\alpha)b )u_t v_t}{1+h v_t}. \end{aligned} \end{eqnarray}$

(5.7)

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

$\begin{equation} J_{\left(\hat{u},\hat{v},\zeta_0\right)}-PT = \ \begin{aligned} & \left(\begin{matrix} \frac{(1-\alpha)^2b^2}{(1-\alpha)^2b^2+h(\zeta_0-g)+hb(1-\alpha)(\zeta_0-1)} \\ \\ \frac{-h(g+b(1-\alpha)(1-g)-\zeta_0)(-\hat{k_1}+1-b(1-\alpha))}{(1-\alpha)^2b^2+h(g-\zeta_0)+hb(1-\alpha)(1-g)} \\ \\ \end{matrix}\right.\\ \\ &\qquad\qquad \left.\begin{matrix} \frac{h(1+b(1-\alpha))(g+b(1-\alpha)(1-\zeta_0)-\zeta_0)}{(1-\alpha)^2b^2+h(\zeta_0-g)+hb(1-\alpha)(\zeta_0-1)}\\ \\ \frac{(1-\alpha)^2b^2+h(g-\zeta_0+b(1-\alpha)(1-g))\hat{k_2}}{(1-\alpha)^2b^2+h(g-\zeta_0)+hb(1-\alpha)(1-g)} \end{matrix}\right). \end{aligned} \end{equation}$

(5.8)

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

$\begin{eqnarray} \begin{aligned} \lambda^2-\frac{{\left\{ \begin{array}{l} {(1-\alpha)^3b^3(2b(1-\alpha)-h(g-\zeta))+h(g+(1-\alpha)b(1-g)-\zeta_0) ((1-\alpha)^2b^2}\;\;\;\\ {\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;+h(\zeta_0-g)+h b(1-\alpha)(\zeta_0-1))\hat{k_2}} \end{array} \right\}}} {{\left\{ \begin{array}{l} ((1-\alpha)^2b^2+h(\zeta_0-g)+h b(\zeta_0-1)(1-\alpha))\\ ((1-\alpha)^2b^2+h(g-\zeta_0)+h b(1-g)(1-\alpha)) \end{array} \right\}}} \lambda +\\ \frac{\left\{ \begin{array}{l} (1-\alpha)^4b^4-h^2((1-\alpha)b+1)(g+b(1-\zeta_0)(1-\alpha)-\zeta_0) (g+b(1-\alpha)\times\\ (1-g)-\zeta_0)(-1+b(1-\alpha)+\hat{k_1})+h(1-\alpha)^2b^2(g-\zeta_0+ (1-g)(1-\alpha)b)\hat{k_2} \end{array} \right\}} {((1-\alpha)^2b^2+h(\zeta_0-g)+h b(1-\alpha)(\zeta_0-1)) ((1-\alpha)^2b^2+h(g-\zeta_0)+h b(1-\alpha)(1-g))} = 0. \end{aligned} \end{eqnarray}$

(5.9)

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

$\begin{equation} \lambda_1+\lambda_2 = \frac{\left\{{(1-\alpha)^3b^3(2b(1-\alpha)-h(g-\zeta))+h(g+b(1-\alpha)(1-g)-\zeta_0)\times}{((1-\alpha)^2b^2+h(\zeta_0-g)+hb(1-\alpha)(\zeta_0-1))\hat{k_2}}\right\}}{\left\{{((1-\alpha)^2b^2+h(\zeta_0-g)+hb(1-\alpha)(\zeta_0-1))\times}{((1-\alpha)^2b^2+h(g-\zeta_0)+hb(1-\alpha)(1-g))}\right\}}, \end{equation}$

(5.10)

and

$\begin{eqnarray} \begin{split} \lambda_1\lambda_2 = & \frac {\left\{ \begin{array}{l} (1-\alpha)^4b^4-h^2(b(1-\alpha)+1)(g+b(1-\zeta_0)(1-\alpha)-\zeta_0)(g+b(1-\alpha)\times\\ \;\;\;(1-g)-\zeta_0)(-1+b(1-\alpha)+\hat{k_1})+h(1-\alpha)^2b^2(g-\zeta_0+(1-g)(1-\alpha)b)\hat{k_2} \end{array} \right\}} {\left\{ \begin{array}{l} ((1-\alpha)^2b^2+h(\zeta_0-g)+hb(1-\alpha)(\zeta_0-1))\times\\ \;\;\;\;\;\;((1-\alpha)^2b^2+h(g-\zeta_0)+hb(1-\alpha)(1-g)) \end{array} \right\}}. \end{split} \end{eqnarray}$

(5.11)

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

$\begin{eqnarray} \begin{split} L_1: \ &h^2((1-\alpha)b+1)(g+b(1-\zeta_0)(1-\alpha)-\zeta_0)(g+b(1-g)(1-\alpha)- \zeta_0)\hat{k_1}-\\&h(1-\alpha)^2b^2(g-\zeta_0+(1-g)(1-\alpha)b)\hat{k_2}+ ((1-\alpha)^2b^2+h(\zeta_0-g)+hb(\zeta_0-1)(1-\alpha))\times\\&((1-\alpha)^2 b^2+h(g-\zeta_0)+hb(1-g)(1-\alpha))+b^4(1-\alpha)^4-h^2(1-(1-\alpha)^2b^2) \times \\&(g+b(1-\zeta_0)(1-\alpha)-\zeta_0)(g+b(1-g)(1-\alpha)-\zeta_0) = 0. \end{split} \end{eqnarray}$

(5.12)

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

$\begin{equation} \begin{split} L_2: \ &h^2((1-\alpha)b+1)(g+(1-\alpha)b(1-\zeta_0)-\zeta_0)(g+(1-\alpha)b(1-g)- \zeta_0)\hat{k_1}-\\&h^2(g-\zeta_0+(1-g)(1-\alpha)b)(g+b(1-\alpha)(1-\zeta_0) -\zeta_0)\hat{k_2}+\\&h^2(1-\alpha)^2b^2(g-\zeta_0+(1-g)(1-\alpha)b)(g+b(1-\alpha)(1-\zeta_0)-\zeta_0) = 0. \end{split} \end{equation}$

(5.13)

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

$\begin{eqnarray} \begin{split} L_3:\ &h^2((1-\alpha)b+1)(g+(1-\alpha)b(1-\zeta_0)-\zeta_0)(g+(1-\alpha)b(1-g)- \zeta_0)\hat{k_1}-\\&h(g-\zeta_0+(1-g)(1-\alpha)b)\Big(2(1-\alpha)^2b^2-h(g+b(1-\alpha)(1-\zeta_0)- \zeta_0)\Bigr)\hat{k_2}+\\&b^2(1-\alpha)^2(-4b^2(1-\alpha)^2+hb^2(1-\alpha)^2(1-g)(1-\zeta_0)+ \\&h^2(g-\zeta_0)^2-bh(1-\alpha)(g-\zeta_0)(-2+h(-2+g+\zeta_0))) = 0. \end{split} \end{eqnarray}$

(5.14)

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 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

$\begin{equation} J|_{\mathrm{PIFP} = (0.36511216748407743,0.1018618715629854)} = \left( \begin{array}{cc} 0.7557825001858909& -1.2267899776912048\\\\ 0.2500158244698668& 0.9173050343262055 \end{array} \right), \end{equation}$

(6.1)

with

$\begin{equation} \lambda^2-1.6730875345120961\lambda+1 = 0. \end{equation}$

(6.2)

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 indicates that PIFP is a stable focus. Similarly, if $h = 0.805$ , $0.835$ , $0.855$ , $0.865$ , $0.875 < 0.8850211364486024$ , then – 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

$\begin{eqnarray} \begin{aligned} l_{11}& = 0.7537070762796483, l_{12} = -1.2372155583706284,\\ l_{13}& = -1.0168530987328883, l_{14} = -2.5540045067599264,\\ l_{21}& = 0.252599297693408, l_{22} = 0.916450527496516,\\ l_{23}& = 2.272634068707745, l_{24} = -1.50338800913407, \end{aligned} \end{eqnarray}$

(6.3)

$\begin{eqnarray} \eta = 0.8350788018880821, \zeta = &0.39028624560846026. \end{eqnarray}$

(6.4)

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

$\begin{eqnarray} \begin{aligned} \sigma_{11}& = 1.0502427204029419, \sigma_{12} = 0.9967928302104191, \sigma_{21} = 0.8306993329208056,\\ \sigma_{22}& = -2.848581027486128, \sigma_{23} = 0.5867516617777137. \end{aligned} \end{eqnarray}$

(6.5)

From (6.5) and (4.31), one gets

$\begin{eqnarray} \begin{aligned} \rho_{02}& = -0.4495845767707965 + 0.6035609562272346 \iota,\\ \rho_{11}& = 0.5251213602014709 + 0.7087254973492596\iota,\\ \rho_{20}& = -0.4495845767707965 - 0.1882112897668318\iota, \rho_{21} = 0. \end{aligned} \end{eqnarray}$

(6.6)

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 . 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].

Figure 3. N-S bifurcation diagrams at PIFP of the HCV model (1.6) for 3a

$u_t$ ; 3b

$v_t$ ; 3c

$u_t$ and

$v_t$ ; 3d

$b$ and

$u_t$ ; 3e

$g$ and

$v_t$ ; 3f MLEs.

DownLoad: Full-Size Img PowerPoint

Figure 4. Phase portrait of the HCV model (1.6) at PIFP if

$h$ 4a 0.7064; 4b 0.7164; 4c 0.7264; 4d 0.7464; 4d 0.7564; 4f 0.7664.

DownLoad: Full-Size Img PowerPoint

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

$\begin{eqnarray} \begin{aligned} u_{t+1} = &\frac{2.652824992 u_t-9.46934126775 u_t v_t}{1+0.885065 u_t},\\ v_{t+1} = &\frac {{\left\{ \begin{array}{l} 5.699211267749999 u_t v_t+\Big(1 + 0.885065\Bigr(0.9999+\\ \;\;\;\;\hat{k_1}(0.25612894566782324+u_t) +\hat{k_2} (-0.22553255211518677+ v_t)\Big)\Bigr)v_t \end{array} \right\}}}{1+0.885065 v_t}, \end{aligned} \end{eqnarray}$

(6.7)

Figure 5. Invariant closed curves of the HCV model (1.6) at PIFP if

$h$ 5a 0.7864; 5b 0.7964; 5c 0.8064; 5d0.8164; 5e 0.8364; 5f 0.8464.

DownLoad: Full-Size Img PowerPoint

where

$\begin{equation} J_{\mathrm{PIFP}}-PT = \ \begin{aligned} & \left(\begin{matrix} 1.933561702180291 \\ \\ 0.9019145397657196 + 0.29831628483824885\hat{k_1}\\ \\ \end{matrix}\right.\\ \\ &\qquad\qquad \left.\begin{matrix} 4.689607176647366\\ \\ 1. + 0.4251435203530145\hat{k_2} \end{matrix}\right), \end{aligned} \end{equation}$

(6.8)

with lines of marginal stability being

$\begin{eqnarray} \begin{split} L_1:1272.0324493907451+270.18749966926106\hat{k_1}-111.40041843372863\hat{k_2} = 0, \end{split} \end{eqnarray}$

(6.9)

$\begin{equation} \begin{split} L_2: 870.6576949546483+270.18749966926106\hat{k_1}-53.786317829588036\hat{k_2} = 0, \end{split} \end{equation}$

(6.10)

and

$\begin{eqnarray} \begin{split} L_3:-530.511530819259+270.18749966926106\hat{k_1}-169.01451903786923\hat{k_2} = 0. \end{split} \end{eqnarray}$

(6.11)

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

Figure 6. Region of stability where

$|\lambda_{1, 2}| < 1$ .

DownLoad: Full-Size Img PowerPoint

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)

$\begin{aligned} \lambda_{1,2}|_{(3.21)}=& \frac{b^3(1-\alpha)^3(2b(1-\alpha)-h(g-\zeta))}{2((1-\alpha)^2b^2-h(g-\zeta)-hb(1-\alpha)(1-\zeta))((1-\alpha)^2b^2+h(g-\zeta)+hb(1-\alpha)(1-g))}\pm \\ & 2\Big(\frac{\sqrt{{64(1 - b^2 (1 - \alpha)^2) (g - \zeta)^2 (g + b (1 - g) (1 - \alpha) - \zeta)^2}\\ {{\;\;\;(g + b (1 - \alpha) (1 - \zeta) - \zeta)^2 - 16 b^2 (-1 + b^2 (1 - \alpha)^2)}\\ {{\;\;\;\;\;\; (1 - \alpha)^2 (g - \zeta) (g + b (1 - g) (1 - \alpha) - \zeta) (g +b (1 - \alpha)}\\ {{\;\;\;\;\;\;\;\;\;(1 - \zeta) - \zeta) (-g +\zeta+\sqrt{{\left\{ \begin{array}{l} 5(g-\zeta)^2+4b^2(1-\alpha)^2(1-g)(1-\zeta)-\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;4b(1-\alpha)(g-\zeta)(-2+g+\zeta) \end{array} \right\}}})}\\ {{\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;+b^4 (1 - \alpha)^4 (9 (g - \zeta)^2 +8 b^2 (1 - g) (1 - \alpha)^2 (1 - \zeta) +}\\ {{\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;8 b (-1 + \alpha) (g - \zeta) (-2 + g + \zeta)) (-g + \zeta +}\\ {\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\sqrt{{\left\{ \begin{array}{l} 5(g-\zeta)^2+4b^2(1-\alpha)^2(1-g)(1-\zeta)-\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;4b(1-\alpha)(g-\zeta)(-2+g+\zeta) \end{array} \right\}}})^2 +}\\ {\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{2 b^4 (1 - \alpha)^4 (g - \zeta) (-g + \zeta +}\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\sqrt{\left\{\begin{array}{l} {5(g-\zeta)^2+4b^2(1-\alpha)^2(1-g)(1-\zeta)-}\\ {\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;4b(1-\alpha)(g-\zeta)(-2+g+\zeta)}\end{array}\right\}})^3}}}}}}}}}} {{\left\{ \begin{array}{l} {(-g + \zeta + \sqrt{{\left\{ \begin{array}{l} 5(g-\zeta)^2+4b^2(1-g)(1-\alpha)^2(1-\zeta)-\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;4b(g-\zeta)(1-\alpha)(-2+g+\zeta) \end{array} \right\}}})^4}\\ \;\;\;\;\;\;(16 (g - \zeta)^2 +16 b^2 (1 - g) (1 -\alpha)^2 (1 - \zeta) +\\ \;\;\;\;\;\;\;\;\;16 b (-1 + \alpha) (g - \zeta) (-2 + g + \zeta) - 4 b^2 (1 - \alpha)^2 (g - \zeta)\\ \;\;\;\;\;\;\;\;\;\;\;\;{(g - \zeta - \sqrt{{\left\{ \begin{array}{l} 5(g-\zeta)^2+4b^2(1-\alpha)^2(1-g)(1-\zeta)-\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;4b(1-\alpha)^2(g-\zeta)(-2+g+\zeta) \end{array} \right\}}}) -}\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{b^2 (1 - \alpha)^2 (-g + \zeta +\sqrt{{\left\{ \begin{array}{l} 5(g-\zeta)^2+4b^2(1-\alpha)^2(1-g)(1-\zeta)-\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;4b(1-\alpha)(g-\zeta)(-2+g+\zeta) \end{array} \right\}}}} \end{array} \right\}}}\Bigr)\\ &\not=1\ \mathrm{or} \ -1 \end{aligned}$

(A.1)

$\begin{aligned} \lambda_{1,2}|_{(3.22)}=& \frac{(1-\alpha)^3b^3(2b(1-\alpha)-h(g-\zeta))}{2((1-\alpha)^2b^2-h(g-\zeta)-hb(1-\alpha)(1-\zeta))((1-\alpha)^2b^2+h(g-\zeta)+hb(1-\alpha)(1-g))}\pm \\& 2\Big(\frac{\sqrt{{64(1 - b^2 (1 - \alpha)^2) (g - \zeta)^2 (g +b (1 - g) (1 - \alpha) - \zeta)^2}\\ {{\;\;\;(g +b (1 - \alpha) (1 - \zeta) - \zeta)^2 +16 b^2 (-1 + b^2 (1 - \alpha)^2)}\\ {{\;\;\;\;\;\;(1 - \alpha)^2 (g - \zeta) (g +b (1 - g) (1 - \alpha) - \zeta) (g +b (1 - \alpha)}\\ {{\;\;\;\;\;\;\;\;(1 - \zeta) -\zeta)(-g+\zeta+\sqrt{\left\{\begin{array}{l}{5(g-\zeta)^2+4b^2(1-\alpha)^2(1-g)(1-\zeta)-}\\ {\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;4b(1-\alpha)(g-\zeta)(-2+g+\zeta)}\end{array}\right\}})}\\ {{\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;+b^4 (1 - \alpha)^4 (9 (g - \zeta)^2 +8 b^2 (1 - g) (1 - \alpha)^2 (1 - \zeta) +}\\ {{\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;8 b (-1 + \alpha) (g - \zeta) (-2 + g + \zeta)) (-g + \zeta+}\\ {\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\sqrt{\left\{\begin{array}{l}{5(g-\zeta)^2+4b^2(1-g)(1-\alpha)^2(1-\zeta)-}\\ {\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;4b(g-\zeta)(1-\alpha)(-2+g+\zeta)}\end{array}\right\}})^2 -}\\ {{\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;2 b^4(1-\alpha)^4(g-\zeta)(-g+\zeta+}\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\sqrt{\left\{\begin{array}{l}{5(g-\zeta)^2+4b^2(1-\alpha)^2(1-g)(1-\zeta)-}\\ {\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;4b(1-\alpha)(g-\zeta)(-2+g+\zeta)}\end{array}\right\}})^3}}}}}}}}}} {\left\{\begin{array}{l}{(-g + \zeta + \sqrt{\left\{\begin{array}{l}{5(g-\zeta)^2+4b^2(1-\alpha)^2(1-g)(1-\zeta)-}\\ {\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;4b(1-\alpha)(g-\zeta)(-2+g+\zeta)}\end{array}\right\}})^4}\\ {{\;\;\;(16 (g - \zeta)^2 +16 b^2 (1 - g) (1 -\alpha)^2 (1 - \zeta) + }\\ {\;\;\;\;\;\;{16 b (-1 + \alpha) (g - \zeta) (-2 + g + \zeta) - 4 b^2 (1 - \alpha)^2 (g - \zeta) }\\ {{\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;(g - \zeta-\sqrt{\left\{\begin{array}{l}{5(g-\zeta)^2+4b^2(1-\alpha)^2(1-g)(1-\zeta)-}{4b(1-\alpha)(g-\zeta)(-2+g+\zeta)}\end{array}\right\}}) -}\\ {\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;b^2 (1 - \alpha)^2 (-g + \zeta +\sqrt{\left\{\begin{array}{l}{5(g-\zeta)^2+4b^2(1-g)(1-\alpha)^2(1-\zeta)-}\\ {\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;4b(g-\zeta)(1-\alpha)(-2+g+\zeta)}\end{array}\right\}}}}}}\end{array}\right\}}\Bigr)\\ &\not=1\ \mathrm{or} \ -1 \end{aligned}$

(A.2)

Acknowledgments

The authors would like to acknowledge the financial support from Australian Research Council (Grant No. DP160102491).

Conflict of interest

The authors declare no conflict of interests.

References

[1]	Hu C, Ashok D, Nisbet DR, et al. (2019) Bioinspired surface modification of orthopedic implants for bone tissue engineering. Biomaterials 119366. doi: 10.1016/j.biomaterials.2019.119366
[2]	Karsenty G, Olson EN (2016) Bone and muscle endocrine functions: unexpected paradigms of inter-organ communication. Cell 164: 1248-1256. doi: 10.1016/j.cell.2016.02.043
[3]	Rossi M, Battafarano G, Pepe J, et al. (2019) The endocrine function of osteocalcin regulated by bone resorption: A lesson from reduced and increased bone mass diseases. Int J Mol Sci 20: 4502. doi: 10.3390/ijms20184502
[4]	Loebel C, Burdick JA (2018) Engineering stem and stromal cell therapies for musculoskeletal tissue repair. Cell Stem Cell 22: 325-339. doi: 10.1016/j.stem.2018.01.014
[5]	Dimitriou R, Jones E, McGonagle D, et al. (2011) Bone regeneration: current concepts and future directions. BMC Med 9: 66. doi: 10.1186/1741-7015-9-66
[6]	Nordin M, Frankel VH (2001) Basic Biomechanics of the Musculoskeletal System, 3 Eds USA: Lippincott Williams & Wilkins.
[7]	Kobayashi S, Takahashi HE, Ito A, et al. (2003) Trabecular minimodeling in human iliac bone. Bone 32: 163-169. doi: 10.1016/S8756-3282(02)00947-X
[8]	Bartl R, Bartl C (2019) Control and regulation of bone remodelling. The Osteoporosis Manual Cham: Springer, 31-39. doi: 10.1007/978-3-030-00731-7_4
[9]	Kenkre JS, Bassett JHD (2018) The bone remodelling cycle. Ann Clin Biochem 55: 308-327. doi: 10.1177/0004563218759371
[10]	Prendergast PJ, Huiskes R (1995) The biomechanics of Wolff's law: recent advances. Irish J Med Sci 164: 152-154. doi: 10.1007/BF02973285
[11]	Wegst UGK, Bai H, Saiz E, et al. (2015) Bioinspired structural materials. Nat Mater 14: 23-36. doi: 10.1038/nmat4089
[12]	Reznikov N, Shahar R, Weiner S (2014) Bone hierarchical structure in three dimensions. Acta Biomater 10: 3815-3826. doi: 10.1016/j.actbio.2014.05.024
[13]	Weiner S, Wagner HD (1998) The material bone: structure-mechanical function relations. Annu Rev Mater Sci 28: 271-298. doi: 10.1146/annurev.matsci.28.1.271
[14]	Recker RR, Kimmel DB, Dempster D, et al. (2011) Issues in modern bone histomorphometry. Bone 49: 955-964. doi: 10.1016/j.bone.2011.07.017
[15]	Eriksen EF, Vesterby A, Kassem M, et al. (1993) Bone remodeling and bone structure. Physiology and Pharmacology of Bone Heidelberg: Springer, 67-109. doi: 10.1007/978-3-642-77991-6_2
[16]	Augat P, Schorlemmer S (2006) The role of cortical bone and its microstructure in bone strength. Age Ageing 35: ii27-ii31. doi: 10.1093/ageing/afl081
[17]	Kozielski M, Buchwald T, Szybowicz M, et al. (2011) Determination of composition and structure of spongy bone tissue in human head of femur by Raman spectral mapping. J Mater Sci: Mater Med 22: 1653-1661. doi: 10.1007/s10856-011-4353-0
[18]	Cross LM, Thakur A, Jalili NA, et al. (2016) Nanoengineered biomaterials for repair and regeneration of orthopedic tissue interfaces. Acta Biomater 42: 2-17. doi: 10.1016/j.actbio.2016.06.023
[19]	Zebaze R, Seeman E (2015) Cortical bone: a challenging geography. J Bone Miner Res 30: 24-29. doi: 10.1002/jbmr.2419
[20]	Liu Y, Luo D, Wang T (2016) Hierarchical structures of bone and bioinspired bone tissue engineering. Small 12: 4611-4632. doi: 10.1002/smll.201600626
[21]	Brodsky B, Persikov AV (2005) Molecular structure of the collagen triple helix. Adv Protein Chem 70: 301-339. doi: 10.1016/S0065-3233(05)70009-7
[22]	Cui FZ, Li Y, Ge J (2007) Self-assembly of mineralized collagen composites. Mater Sci Eng R Rep 57: 1-27. doi: 10.1016/j.mser.2007.04.001
[23]	Wang Y, Azaïs T, Robin M, et al. (2012) The predominant role of collagen in the nucleation, growth, structure and orientation of bone apatite. Nat Mater 11: 724-733. doi: 10.1038/nmat3362
[24]	Bentmann A, Kawelke N, Moss D, et al. (2010) Circulating fibronectin affects bone matrix, whereas osteoblast fibronectin modulates osteoblast function. J Bone Miner Res 25: 706-715.
[25]	Szweras M, Liu D, Partridge EA, et al. (2002) α2-HS glycoprotein/fetuin, a transforming growth factor-β/bone morphogenetic protein antagonist, regulates postnatal bone growth and remodeling. J Biol Chem 277: 19991-19997. doi: 10.1074/jbc.M112234200
[26]	Boskey AL, Robey PG (2013) The regulatory role of matrix proteins in mineralization of bone. Osteoporosis, 4 Eds Academic Press, 235-255. doi: 10.1016/B978-0-12-415853-5.00011-X
[27]	Boskey AL (2013) Bone composition: relationship to bone fragility and antiosteoporotic drug effects. Bonekey Rep 2: 447. doi: 10.1038/bonekey.2013.181
[28]	Stock SR (2015) The mineral–collagen interface in bone. Calcified Tissue Int 97: 262-280. doi: 10.1007/s00223-015-9984-6
[29]	Nikel O, Laurencin D, McCallum SA, et al. (2013) NMR investigation of the role of osteocalcin and osteopontin at the organic–inorganic interface in bone. Langmuir 29: 13873-13882. doi: 10.1021/la403203w
[30]	He G, Dahl T, Veis A, et al. (2003) Nucleation of apatite crystals in vitro by self-assembled dentin matrix protein 1. Nat Mater 2: 552-558. doi: 10.1038/nmat945
[31]	Clarke B (2008) Normal bone anatomy and physiology. Clin J Am Soc Nephro 3: S131-S139. doi: 10.2215/CJN.04151206
[32]	Olszta MJ, Cheng X, Jee SS, et al. (2007) Bone structure and formation: A new perspective. Mater Sci Eng R Rep 58: 77-116. doi: 10.1016/j.mser.2007.05.001
[33]	Nair AK, Gautieri A, Chang SW, et al. (2013) Molecular mechanics of mineralized collagen fibrils in bone. Nature Commun 4: 1724. doi: 10.1038/ncomms2720
[34]	Landis WJ (1995) The strength of a calcified tissue depends in part on the molecular structure and organization of its constituent mineral crystals in their organic matrix. Bone 16: 533-544. doi: 10.1016/8756-3282(95)00076-P
[35]	Hunter GK, Hauschka PV, POOLE RA, et al. (1996) Nucleation and inhibition of hydroxyapatite formation by mineralized tissue proteins. Biochem J 317: 59-64. doi: 10.1042/bj3170059
[36]	Oikeh I, Sakkas P, Blake D P, et al. (2019) Interactions between dietary calcium and phosphorus level, and vitamin D source on bone mineralization, performance, and intestinal morphology of coccidia-infected broilers. Poult Sci 11: 5679-5690. doi: 10.3382/ps/pez350
[37]	Boyce BF, Rosenberg E, de Papp AE, et al. (2012) The osteoclast, bone remodelling and treatment of metabolic bone disease. Eur J Clin Invest 42: 1332-1341. doi: 10.1111/j.1365-2362.2012.02717.x
[38]	Teitelbaum SL (2000) Bone resorption by osteoclasts. Science 289: 1504-1508. doi: 10.1126/science.289.5484.1504
[39]	Yoshida H, Hayashi SI, Kunisada T, et al. (1990) The murine mutation osteopetrosis is in the coding region of the macrophage colony stimulating factor gene. Nature 345: 442-444. doi: 10.1038/345442a0
[40]	Roodman GD (2006) Regulation of osteoclast differentiation. Ann NY Acad Sci 1068: 100-109. doi: 10.1196/annals.1346.013
[41]	Martin TJ (2013) Historically significant events in the discovery of RANK/RANKL/OPG. World J Orthop 4: 186-197. doi: 10.5312/wjo.v4.i4.186
[42]	Coetzee M, Haag M, Kruger MC (2007) Effects of arachidonic acid, docosahexaenoic acid, prostaglandin E2 and parathyroid hormone on osteoprotegerin and RANKL secretion by MC3T3-E1 osteoblast-like cells. J Nutr Biochem 18: 54-63. doi: 10.1016/j.jnutbio.2006.03.002
[43]	Steeve KT, Marc P, Sandrine T, et al. (2004) IL-6, RANKL, TNF-alpha/IL-1: interrelations in bone resorption pathophysiology. Cytokine Growth F R 15: 49-60. doi: 10.1016/j.cytogfr.2003.10.005
[44]	Mellis DJ, Itzstein C, Helfrich M, et al. (2011) The skeleton: a multi-functional complex organ. The role of key signalling pathways in osteoclast differentiation and in bone resorption. J Endocrinol 211: 131-143. doi: 10.1530/JOE-11-0212
[45]	Silva I, Branco J (2011) Rank/Rankl/opg: literature review. Acta Reumatol Port 36: 209-218.
[46]	Martin TJ, Sims NA (2015) RANKL/OPG; Critical role in bone physiology. Rev Endocr Metab Dis 16: 131-139. doi: 10.1007/s11154-014-9308-6
[47]	Wang Y, Qin QH (2012) A theoretical study of bone remodelling under PEMF at cellular level. Comput Method Biomec 15: 885-897. doi: 10.1080/10255842.2011.565752
[48]	Weitzmann MN, Pacifici R (2007) T cells: unexpected players in the bone loss induced by estrogen deficiency and in basal bone homeostasis. Ann NY Acad Sci 1116: 360-375. doi: 10.1196/annals.1402.068
[49]	Duong LT, Lakkakorpi P, Nakamura I, et al. (2000) Integrins and signaling in osteoclast function. Matrix Biol 19: 97-105. doi: 10.1016/S0945-053X(00)00051-2
[50]	Stenbeck G (2002) Formation and function of the ruffled border in osteoclasts. Semin Cell Dev Biol 13: 285-292. doi: 10.1016/S1084952102000587
[51]	Jurdic P, Saltel F, Chabadel A, et al. (2006) Podosome and sealing zone: specificity of the osteoclast model. Eur J Cell Biol 85: 195-202. doi: 10.1016/j.ejcb.2005.09.008
[52]	Väänänen HK, Laitala-Leinonen T (2008) Osteoclast lineage and function. Arch Biochem Biophys 473: 132-138. doi: 10.1016/j.abb.2008.03.037
[53]	Vaananen HK, Zhao H, Mulari M, et al. (2000) The cell biology of osteoclast function. J cell Sci 113: 377-381.
[54]	Sabolová V, Brinek A, Sládek V (2018) The effect of hydrochloric acid on microstructure of porcine (Sus scrofa domesticus) cortical bone tissue. Forensic Sci Int 291: 260-271. doi: 10.1016/j.forsciint.2018.08.030
[55]	Delaissé JM, Engsig MT, Everts V, et al. (2000) Proteinases in bone resorption: obvious and less obvious roles. Clin Chim Acta 291: 223-234. doi: 10.1016/S0009-8981(99)00230-2
[56]	Logar DB, Komadina R, Preželj J, et al. (2007) Expression of bone resorption genes in osteoarthritis and in osteoporosis. J Bone Miner Metab 25: 219-225. doi: 10.1007/s00774-007-0753-0
[57]	Lorget F, Kamel S, Mentaverri R, et al. (2000) High extracellular calcium concentrations directly stimulate osteoclast apoptosis. Biochem Bioph Res Co 268: 899-903. doi: 10.1006/bbrc.2000.2229
[58]	Nesbitt SA, Horton MA (1997) Trafficking of matrix collagens through bone-resorbing osteoclasts. Science 276: 266-269. doi: 10.1126/science.276.5310.266
[59]	Xing L, Boyce BF (2005) Regulation of apoptosis in osteoclasts and osteoblastic cells. Biochem Bioph Res Co 328: 709-720. doi: 10.1016/j.bbrc.2004.11.072
[60]	Hughes DE, Wright KR, Uy HL, et al. (1995) Bisphosphonates promote apoptosis in murine osteoclasts in vitro and in vivo. J Bone Miner Res 10: 1478-1487. doi: 10.1002/jbmr.5650101008
[61]	Choi Y, Arron JR, Townsend MJ (2009) Promising bone-related therapeutic targets for rheumatoid arthritis. Nat Rev Rheumatol 5: 543-548. doi: 10.1038/nrrheum.2009.175
[62]	Harvey NC, McCloskey E, Kanis JA, et al. (2017) Bisphosphonates in osteoporosis: NICE and easy? Lancet 390: 2243-2244. doi: 10.1016/S0140-6736(17)32850-7
[63]	Ducy P, Schinke T, Karsenty G (2000) The osteoblast: a sophisticated fibroblast under central surveillance. Science 289: 1501-1504. doi: 10.1126/science.289.5484.1501
[64]	Katagiri T, Takahashi N (2002) Regulatory mechanisms of osteoblast and osteoclast differentiation. Oral dis 8: 147-159. doi: 10.1034/j.1601-0825.2002.01829.x
[65]	Kretzschmar M, Liu F, Hata A, et al. (1997) The TGF-beta family mediator Smad1 is phosphorylated directly and activated functionally by the BMP receptor kinase. Gene Dev 11: 984-995. doi: 10.1101/gad.11.8.984
[66]	Bennett CN, Longo KA, Wright WS, et al. (2005) Regulation of osteoblastogenesis and bone mass by Wnt10b. P Natl A Sci 102: 3324-3329. doi: 10.1073/pnas.0408742102
[67]	Wang Y, Qin QH, Kalyanasundaram S (2009) A theoretical model for simulating effect of parathyroid hormone on bone metabolism at cellular level. Mol Cell Biomech 6: 101-112.
[68]	Elefteriou F, Ahn JD, Takeda S, et al. (2005) Leptin regulation of bone resorption by the sympathetic nervous system and CART. Nature 434: 514-520. doi: 10.1038/nature03398
[69]	Proff P, Römer P (2009) The molecular mechanism behind bone remodelling: a review. Clin Oral Invest 13: 355-362. doi: 10.1007/s00784-009-0268-2
[70]	Katsimbri P (2017) The biology of normal bone remodelling. Eur J Cancer Care 26: e12740. doi: 10.1111/ecc.12740
[71]	Fratzl P, Weinkamer R (2007) Nature's hierarchical materials. Prog Mater Sci 52: 1263-1334. doi: 10.1016/j.pmatsci.2007.06.001
[72]	Athanasiou KA, Zhu CF, Lanctot DR, et al. (2000) Fundamentals of biomechanics in tissue engineering of bone. Tissue Eng 6: 361-381. doi: 10.1089/107632700418083
[73]	Takahashi N, Udagawa N, Suda T (1999) A new member of tumor necrosis factor ligand family, ODF/OPGL/TRANCE/RANKL, regulates osteoclast differentiation and function. Biocheml Bioph Res Co 256: 449-455. doi: 10.1006/bbrc.1999.0252
[74]	Nakashima T, Hayashi M, Fukunaga T, et al. (2011) Evidence for osteocyte regulation of bone homeostasis through RANKL expression. Nat Med 17: 1231-1234. doi: 10.1038/nm.2452
[75]	Prideaux M, Findlay DM, Atkins GJ (2016) Osteocytes: the master cells in bone remodelling. Curr Opin Pharmacol 28: 24-30. doi: 10.1016/j.coph.2016.02.003
[76]	Dallas SL, Prideaux M, Bonewald LF (2013) The osteocyte: an endocrine cell… and more. Endocr Rev 34: 658-690. doi: 10.1210/er.2012-1026
[77]	Rochefort GY, Pallu S, Benhamou CL (2010) Osteocyte: the unrecognized side of bone tissue. Osteoporosis Int 21: 1457-1469. doi: 10.1007/s00198-010-1194-5
[78]	Rowe PSN (2012) Regulation of bone–renal mineral and energy metabolism: The PHEX, FGF23, DMP1, MEPE ASARM pathway. Crit Rev Eukaryot Gene Expr 22: 61-86. doi: 10.1615/CritRevEukarGeneExpr.v22.i1.50
[79]	Pajevic PD, Krause DS (2019) Osteocyte regulation of bone and blood. Bone 119: 13-18. doi: 10.1016/j.bone.2018.02.012
[80]	Frost HM (1987) The mechanostat: a proposed pathogenic mechanism of osteoporoses and the bone mass effects of mechanical and nonmechanical agents. Bone Miner 2: 73-85.
[81]	Tate MLK, Adamson JR, Tami AE, et al. (2004) The osteocyte. Int J Biochem Cell Biol 36: 1-8. doi: 10.1016/S1357-2725(03)00241-3
[82]	Bonewald LF, Johnson ML (2008) Osteocytes, mechanosensing and Wnt signaling. Bone 42: 606-615. doi: 10.1016/j.bone.2007.12.224
[83]	Manolagas SC, Parfitt AM (2010) What old means to bone. Trends Endocrinol Metab 21: 369-374. doi: 10.1016/j.tem.2010.01.010
[84]	Wang Y, Qin QH (2010) Parametric study of control mechanism of cortical bone remodeling under mechanical stimulus. Acta Mech Sinica 26: 37-44. doi: 10.1007/s10409-009-0313-z
[85]	Qu C, Qin QH, Kang Y (2006) A hypothetical mechanism of bone remodeling and modeling under electromagnetic loads. Biomaterials 27: 4050-4057. doi: 10.1016/j.biomaterials.2006.03.015
[86]	Parfitt AM (2002) Targeted and nontargeted bone remodeling: relationship to basic multicellular unit origination and progression. Bone 1: 5-7. doi: 10.1016/S8756-3282(01)00642-1
[87]	Hadjidakis DJ, Androulakis II (2006) Bone remodeling. Ann NYAcad Sci 1092: 385-396. doi: 10.1196/annals.1365.035
[88]	Vaananen HK, Zhao H, Mulari M, et al. (2000) The cell biology of osteoclast function. J cell Sci 113: 377-381.
[89]	Goldring SR (2015) The osteocyte: key player in regulating bone turnover. RMD Open 1: e000049. doi: 10.1136/rmdopen-2015-000049
[90]	Silver IA, Murrills RJ, Etherington DJ (1988) Microelectrode studies on the acid microenvironment beneath adherent macrophages and osteoclasts. Exp Cell Res 175: 266-276. doi: 10.1016/0014-4827(88)90191-7
[91]	Delaissé JM, Andersen TL, Engsig MT, et al. (2003) Matrix metalloproteinases (MMP) and cathepsin K contribute differently to osteoclastic activities. Microsc Res Techniq 61: 504-513. doi: 10.1002/jemt.10374
[92]	Delaisse JM (2014) The reversal phase of the bone-remodeling cycle: cellular prerequisites for coupling resorption and formation. Bonekey Rep 3: 561. doi: 10.1038/bonekey.2014.56
[93]	Bonewald LF, Mundy GR (1990) Role of transforming growth factor-beta in bone remodeling. Clin Orthop Relat R 250: 261-276.
[94]	Locklin RM, Oreffo ROC, Triffitt JT (1999) Effects of TGFβ and bFGF on the differentiation of human bone marrow stromal fibroblasts. Cell Biol Int 23: 185-194. doi: 10.1006/cbir.1998.0338
[95]	Lee B, Oh Y, Jo S, et al. (2019) A dual role of TGF-β in human osteoclast differentiation mediated by Smad1 versus Smad3 signaling. Immunol Lett 206: 33-40. doi: 10.1016/j.imlet.2018.12.003
[96]	Koseki T, Gao Y, Okahashi N, et al. (2002) Role of TGF-β family in osteoclastogenesis induced by RANKL. Cell Signal 14: 31-36. doi: 10.1016/S0898-6568(01)00221-2
[97]	Anderson HC (2003) Matrix vesicles and calcification. Curr Rheumatol Rep 5: 222-226. doi: 10.1007/s11926-003-0071-z
[98]	Bellido T, Plotkin LI, Bruzzaniti A (2019) Bone cells. Basic and Applied Bone Biology, 2 Eds Elsevier, 37-55. doi: 10.1016/B978-0-12-813259-3.00003-8
[99]	Weinstein RS, Jilka RL, Parfitt AM, et al. (1998) Inhibition of osteoblastogenesis and promotion of apoptosis of osteoblasts and osteocytes by glucocorticoids. Potential mechanisms of their deleterious effects on bone. J Clin Invest 102: 274-282. doi: 10.1172/JCI2799
[100]	Vezeridis PS, Semeins CM, Chen Q, et al. (2006) Osteocytes subjected to pulsating fluid flow regulate osteoblast proliferation and differentiation. Biochem Bioph Res Co 348: 1082-1088. doi: 10.1016/j.bbrc.2006.07.146
[101]	Lind M, Deleuran B, Thestrup-Pedersen K, et al. (1995) Chemotaxis of human osteoblasts: Effects of osteotropic growth factors. Apmis 103: 140-146. doi: 10.1111/j.1699-0463.1995.tb01089.x
[102]	Russo CR, Lauretani F, Seeman E, et al. (2006) Structural adaptations to bone loss in aging men and women. Bone 38: 112-118. doi: 10.1016/j.bone.2005.07.025
[103]	Ozcivici E, Luu YK, Adler B, et al. (2010) Mechanical signals as anabolic agents in bone. Nat Rev Rheumatol 6: 50-59. doi: 10.1038/nrrheum.2009.239
[104]	Rosa N, Simoes R, Magalhães FD, et al. (2015) From mechanical stimulus to bone formation: a review. Med Eng Phys 37: 719-728. doi: 10.1016/j.medengphy.2015.05.015
[105]	Noble BS, Peet N, Stevens HY, et al. (2003) Mechanical loading: biphasic osteocyte survival and targeting of osteoclasts for bone destruction in rat cortical bone. Am J Physiol-Cell Ph 284: C934-C943. doi: 10.1152/ajpcell.00234.2002
[106]	Robling AG, Castillo AB, Turner CH (2006) Biomechanical and molecular regulation of bone remodeling. Annu Rev Biomed Eng 8: 455-498. doi: 10.1146/annurev.bioeng.8.061505.095721
[107]	Qin QH, Mai YW (1999) A closed crack tip model for interface cracks inthermopiezoelectric materials. Int J Solids Struct 36: 2463-2479. doi: 10.1016/S0020-7683(98)00115-2
[108]	Yu SW, Qin QH (1996) Damage analysis of thermopiezoelectric properties: Part I—crack tip singularities. Theor Appl Fract Mec 25: 263-277. doi: 10.1016/S0167-8442(96)00026-2
[109]	Qin QH, Mai YW, Yu SW (1998) Effective moduli for thermopiezoelectric materials with microcracks. Int J Fracture 91: 359-371. doi: 10.1023/A:1007423508650
[110]	Jirousek J, Qin QH (1996) Application of hybrid-Trefftz element approach to transient heat conduction analysis. Comput Struct 58: 195-201. doi: 10.1016/0045-7949(95)00115-W
[111]	Qin QH (1995) Hybrid-Trefftz finite element method for Reissner plates on an elastic foundation. Comput Method Appl M 122: 379-392. doi: 10.1016/0045-7825(94)00730-B
[112]	Qin QH (1994) Hybrid Trefftz finite-element approach for plate bending on an elastic foundation. Appl Math Model 18: 334-339. doi: 10.1016/0307-904X(94)90357-3
[113]	Qin QH (2013) Mechanics of Cellular Bone Remodeling: Coupled Thermal, Electrical, and Mechanical Field Effects CRC Press. doi: 10.1201/b13728
[114]	Wang H, Qin QH (2010) FE approach with Green's function as internal trial function for simulating bioheat transfer in the human eye. Arch Mech 62: 493-510.
[115]	Qin QH (2003) Fracture analysis of cracked thermopiezoelectric materials by BEM. Electronic J Boundary Elem 1: 283-301.
[116]	Qin QH, Ye JQ (2004) Thermoelectroelastic solutions for internal bone remodeling under axial and transverse loads. Int J Solids Struct 41: 2447-2460. doi: 10.1016/j.ijsolstr.2003.12.026
[117]	Qin QH, Qu C, Ye J (2005) Thermoelectroelastic solutions for surface bone remodeling under axial and transverse loads. Biomaterials 26: 6798-6810. doi: 10.1016/j.biomaterials.2005.03.042
[118]	Ducher G, Jaffré C, Arlettaz A, et al. (2005) Effects of long-term tennis playing on the muscle-bone relationship in the dominant and nondominant forearms. Can J Appl Physiol 30: 3-17. doi: 10.1139/h05-101
[119]	Robling AG, Hinant FM, Burr DB, et al. (2002) Improved bone structure and strength after long-term mechanical loading is greatest if loading is separated into short bouts. J Bone Miner Res 17: 1545-1554. doi: 10.1359/jbmr.2002.17.8.1545
[120]	Rubin J, Rubin C, Jacobs CR (2006) Molecular pathways mediating mechanical signaling in bone. Gene 367: 1-16. doi: 10.1016/j.gene.2005.10.028
[121]	Tatsumi S, Ishii K, Amizuka N, et al. (2007) Targeted ablation of osteocytes induces osteoporosis with defective mechanotransduction. Cell Metab 5: 464-475. doi: 10.1016/j.cmet.2007.05.001
[122]	Robling AG, Turner CH (2009) Mechanical signaling for bone modeling and remodeling. Crit Rev Eukar Gene 19: 319-338. doi: 10.1615/CritRevEukarGeneExpr.v19.i4.50
[123]	Galli C, Passeri G, Macaluso GM (2010) Osteocytes and WNT: the mechanical control of bone formation. J Dent Res 89: 331-343. doi: 10.1177/0022034510363963
[124]	Robling AG, Duijvelaar KM, Geevers JV, et al. (2001) Modulation of appositional and longitudinal bone growth in the rat ulna by applied static and dynamic force. Bone 29: 105-113. doi: 10.1016/S8756-3282(01)00488-4
[125]	Burr DB, Milgrom C, Fyhrie D, et al. (1996) In vivo measurement of human tibial strains during vigorous activity. Bone 18: 405-410. doi: 10.1016/8756-3282(96)00028-2
[126]	Sun W, Chi S, Li Y, et al. (2019) The mechanosensitive Piezo1 channel is required for bone formation. Elife 8: e47454. doi: 10.7554/eLife.47454
[127]	Goda I, Ganghoffer JF, Czarnecki S, et al. (2019) Topology optimization of bone using cubic material design and evolutionary methods based on internal remodeling. Mech Res Commun 95: 52-60. doi: 10.1016/j.mechrescom.2018.12.003
[128]	Goda I, Ganghoffer JF (2018) Modeling of anisotropic remodeling of trabecular bone coupled to fracture. Arch Appl Mech 88: 2101-2121. doi: 10.1007/s00419-018-1438-y
[129]	Louna Z, Goda I, Ganghoffer JF, et al. (2017) Formulation of an effective growth response of trabecular bone based on micromechanical analyses at the trabecular level. Arch Appl Mech 87: 457-477. doi: 10.1007/s00419-016-1204-y
[130]	Goda I, Ganghoffer JF (2017) Construction of the effective plastic yield surfaces of vertebral trabecular bone under twisting and bending moments stresses using a 3D microstructural model. ZAMM Z Angew Math Mech 97: 254-272. doi: 10.1002/zamm.201600141
[131]	Qin QH, Wang YN (2012) A mathematical model of cortical bone remodeling at cellular level under mechanical stimulus. Acta Mech Sinica-Prc 28: 1678-1692. doi: 10.1007/s10409-012-0154-z

Reader Comments

Your name:*

Email:*
© 2020 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)

通讯作者: 陈斌, bchen63@163.com

1.
沈阳化工大学材料科学与工程学院沈阳 110142

AIMS Bioengineering

1.2

Metrics

Article views(11852) PDF downloads(1269) Cited by(18)

Preview PDF

Download XML

Export Citation

Article outline

Show full outline

Figures and Tables

Figures(2)

AIMS Bioengineering

Bone remodeling and biological effects of mechanical stimulus

Related Papers:

Abstract

1. Introduction

2. Existence of fixed points

3. Local stability

4. Bifurcation analysis

5. Chaos control

6. Numerical simulations

7. Conclusions

Author contributions

Acknowledgements

Conflict of interest

Appendix

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

References

Reader Comments

通讯作者: 陈斌, bchen63@163.com

Metrics

Figures and Tables

Other Articles By Authors

Catalog

AIMS Bioengineering

Bone remodeling and biological effects of mechanical stimulus

Related Papers:

Abstract

1. Introduction

2. Existence of fixed points

3. Local stability

4. Bifurcation analysis

5. Chaos control

6. Numerical simulations

7. Conclusions

Author contributions

Acknowledgements

Conflict of interest

Appendix

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

References

Reader Comments

通讯作者: 陈斌, bchen63@163.com

Metrics

Figures and Tables

Other Articles By Authors

Related pages

Tools

Export File

Citation

Format

Content

Catalog

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