key: cord-0054465-2wzjesm4
authors: Pan, Sonjoy; Chakrabarty, Siddhartha P.
title: Hopf Bifurcation and Stability Switches Induced by Humoral Immune Delay in Hepatitis C
date: 2021-01-05
journal: Indian J Pure Appl Math
DOI: 10.1007/s13226-020-0489-2
sha: 9d8d50481f5b6c3ebc165a58c12e3d0192f72e4f
doc_id: 54465
cord_uid: 2wzjesm4

The role of humoral immune delay on the dynamics of HCV infection incorporating both the modes of infection transmission, namely, viral and cellular transmissions with a non-cytolytic cure of infected hepatocytes is studied. The local and global asymptotic stability of the boundary equilibria, namely, infection-free and immune-free equilibrium are analyzed theoretically as well as numerically under the conditions on the basic reproduction number and the humoral immune reproduction number. The existence of Hopf bifurcation and consequent occurrence of bifurcating periodic orbits around the humoral immune activated equilibrium are illustrated. The findings show that Hopf bifurcation and stability switches occur under certain conditions as the bifurcation parameter crosses the critical values. Furthermore, the dynamical effect of the development rate of B cells is investigated numerically. The results obtained show that the system becomes unstable from stable and regains stability from instability depending on the development rate of B cells for a fixed delay value. Further, the results suggest that a high antigenic stimulation in humoral immunity is beneficial for uninfected hepatocytes with a significant reduction in virions density.

An estimated 71 million individuals are affected worldwide by chronic hepatitis C virus (HCV) infection (a blood-borne hepatological condition) resulting in about 399 thousand cases of fatalities happening in 2016, due to liver cirrhosis and hepatocellular carcinoma [1] . Among HCV infected individuals 60-80% of the cases become chronic, out of which about 15-30% face the risk of developing liver cirrhosis in the long run [1] . The transmission of HCV primarily happens through unscreened blood transfusions, injecting drug use and inadequately sterilized medical equipment [2] . HCV, which is a positive single-stranded viral RNA genome from the family of Flaviviridae [3] , replicates very rapidly resulting in difficulty in developing a vaccine for HCV [4] . The pioneering in-vivo model for HCV dynamics [5] based on similar models for human immunodeficiency virus (HIV) and hepatitis B virus (HBV) dynamics, incorporated the interferon-α (IFN-α) antiviral therapy for HCV infection with three model populations, namely, uninfected and infected hepatocytes, and virions. The analysis in [5] demonstrated that the antiviral therapy is more effective in the reduction of HCV RNA load with the efficacy increasing with an increase in the dosage of IFN-α, accompanied by a minor effect in blocking of production of infected hepatocytes. Dixit et al. [6] included ribavirin along with IFN in the antiviral therapy protocol and described the effect of this combination therapy in HCV infection. Due to the action of ribavirin, a fraction of the virions become non-infectious and hence are not involved in the production of infected hepatocytes. The analysis [6] demonstrated that the ribavirin in combination with IFN significantly improves the process of HCV RNA decline, but as a monotherapy, it produces only short-term early response in terms of decrease of viral load. The extended models [7, 8] of the models in [5, 6] included the argument of homeostatic mechanism by which liver can regenerate itself, thereby which the hepatocytes can proliferate up to a certain maximum level. A comparison of the extended models [7, 8] with the models in [5, 6] showed that the analysis of these models [5, 6] exhibits a biphasic decline in viral load (a rapid viral decline followed by a constant level of viral load), while the models in [7, 8] can better predict the kinetics of HCV RNA load in chronic stage and also explain the biphasic as well as triphasic viral decay. A fast viral decline in the third phase was observed in cases where the majority of hepatocytes were already infected before the initiation of antiviral therapy [8] . The effectiveness of pegylated IFN and ribavirin as antiviral therapy for HCV infection was analyzed by estimating the clearance rate of HCV RNA and the infected cells as well [9] .

The role of cytotoxic T lymphocyte (CTL) and humoral immune response in HCV infection was investigated through the mathematical modeling in [10, 11] . The role of CTL is to reduce the HCV infection with the antibody playing a role of neutralizing HCV RNA. The analysis found the correlation between the CTL and antibody when both are activated during chronic HCV infection. It was observed that the infection would be asymptotic if the CTL activates strongly, but persistent leading to pathology in case of weak CTL response, even in presence of antibody. A complex mechanism of the immune system for HCV patients was considered in the model formulation of [12] which includes the role of dendritic cell (DC) with CTL response in such a way that CTL is produced through the cross-presentation of activated DC and decays by a direct presentation of infected cells. It was also observed that the activation of immune response is dependent on initial DC and CTL levels. The consideration of cure of infected cells through the non-cytolytic process in which a part of the infected cells get converted to uninfected cells, was introduced in [13, 14] .

Timpe et al. [15] observed that the mode of HCV transmission can not only be virus-to-cell but also cell-to-cell. A HCV model considering cell-to-cell transmission was studied for the optimal antiviral treatment policy in [16] . Further, from the epidemiological point of view, the consideration of intracellular delay (the time needed for the hepatocytes to be infected or for replication of the virions) and immune delay (the time needed for antigenic stimulation for the development of B cells in case of antibody response or for the development of T cells in case of CTL response) are more realistic when analyzing HCV viral dynamics. The intracellular delay was incorporated in a HCV model in [17] which included the activation of CTL and antibody response. The role of the humoral immune delay in a HBV model with two stages of infected cells, namely, latently infected and actively infected was illustrated in [18] . A general viral dynamics model with the effect of the humoral immune delay was analyzed in [19] , where the existence of Hopf bifurcation was observed. The effect of the various intracellular delays on the dynamics of other viral infections (like HBV, HIV infection), incorporating both virus-to-cell and cell-to-cell transmissions, was studied in [20] [21] [22] . A general viral infection model considering the viral as well as cellular transmission of infection with cell-mediated immune response was proposed and analyzed with the effect of the intracellular delay as well as cellmediated immune delay [23] . The results showed that the delays could lead to stability switches and occurrence of bifurcating periodic solution, depending upon the intrinsic rate of logistic growth and infection transmission rate as well.

The limited number of models in the literature which illustrated the effect of humoral immune response in viral dynamics did not include either cell-to-cell transmission or cure of infected cells.

On the other hand, the models which considered cell-to-cell transmission did not take into account the humoral immune delay as well as cure of infected cells. However, in this work, we include all these factors in the proposed model.

The dynamics of HCV infection with both viral and cellular transmissions and cure rate in the presence of humoral immune response was analyzed in [24] . This model assumed that subsequent to the entry of the virions, the humoral immune response is stimulated to instantaneously generate B cells.

However, from a realistic epidemiological point of view, there is a time delay between the initial vi-ral entry into the liver and subsequent process of antigenic stimulation to generate B cells [19, 25] .

However, it is considered that the antibody neutralizing the virions is immediate [19] . Accordingly, we incorporate the humoral immune delay to propose the delay model as follows:

Here T (t), I(t), V (t) and Z(t) represent the densities of the uninfected hepatocytes, actively infected hepatocytes, virions and humoral immunity (B cell or antibody) at time t, respectively. It is assumed that the uninfected hepatocytes are being sourced at a constant rate λ (within the liver)

and are cleared at a natural death rate of d 1 . It is also considered that the virions and the infected hepatocytes infect the healthy hepatocytes at a rate of β 1 and β 2 respectively. The infected hepatocytes have a natural death rate of d 2 . Due to the non-cytolytic process, the infected hepatocytes are converted to uninfected ones at a rate α. The infected hepatocytes abet the production of free virions at a rate k, which in turn decay at a rate of 

The uniqueness and positivity of the solution to the system (2.1) with the initial condition (2.2) can be proved using the theory of functional differential equations [26] , and by Theorem A.4 in [27] and Lemma 2 in [28] , respectively. Further, the boundedness of the solution can be proved on the lines of [24] .

The model system (2.1) has three equilibria, namely [24] ,

3. The humoral immune activated equilibrium, E * = (T * , I * , V * , Z * ), where

Further, the system (2.1) has two reproduction numbers, namely [24] , the basic reproduction number,

, and the humoral immune reproduction number,

Theorem 1 -The infection-free equilibrium E 0 is locally as well as globally asymptotically stable for any τ > 0 if R 0 < 1 and unstable if R 0 > 1.

PROOF : The proof for local stability easily follows from the Routh-Hurwitz criteria.

Further, the global stability can be proved using the Lyapunov-LaSalle invariance principle [26] by choosing the Lyapunov functional,

PROOF : The characteristic equation obtained by linearizing (2.1) at E 1 is given by

All the three roots of the equation

which has been proved in [24] . We now investigate whether any complex root with positive real part exists for the following equation:

3) and separating the real and imaginary parts, we obtain the following:

which is a contradiction. Therefore every root of (2.3) must have negative real part. Hence, E 1 is locally asymptotically stable for any Further, the global stability can be proved using the the Lyapunov-LaSalle invariance principle [26] by choosing the Lyapunov functional,

The characteristic equation obtained by linearizing (2.1) at E * is given by

where

When τ > 0, (2.6) becomes transcendental and therefore some of the roots may cross the imaginary axis to the right. We now investigate the existence of purely imaginary roots of (2.6). Substi-

) and separating the real and imaginary parts, we obtain the following,

Squaring and adding the equations (2.8) and (2.9), we obtain

Let γ = ω 2 . Then (2.10) becomes

Therefore (2.6) has a pair of purely imaginary roots ±iω if and only if (2.11) has a positive real root ω 2 . If (2.6) has no positive real root, then (2.11) has no purely imaginary root, in which case the existence of Hopf bifurcation is ruled out and hence E * is locally asymptotically stable for any τ > 0. Now, we suppose that (2.11) has m (1 ≤ m ≤ 4) positive roots, say, γ n , n = 1, 2, . . . , m.

Then (2.10) has m positive roots, say, ω n = √ γ n , n = 1, 2, . . . , m.

Solving (2.8) and (2.9) for cos (ωτ ), we obtain

When ω = ω n (n = 1, 2, . . . , m), we obtain the following from (2.12),

where n = 1, 2, . . . , m and j = 0, 1, 2, . . . . Hence (2.6) has a pair of purely imaginary roots ±iω n

n . Further, note that τ (j) n is a monotonically increasing sequence for every n = 1, 2, . . . , m and lim j→∞ τ (j) n = ∞. Therefore there exists a n 0 ∈ {1, 2, . . . , m} such that

Since E * is locally asymptotically stable for τ = 0 if R H > 1 [24] . Therefore, by Butler's Lemma [31] , E * remains locally asymptotically stable for τ < τ 0 if R H > 1. Let x(τ ) = ξ(τ ) + iω(τ ) be a root of (2.6) near τ = τ 0 with ξ(τ 0 ) = 0, ω(τ 0 ) = ω 0 . Therefore to prove the transversality condition for the existence of Hopf bifurcation [32] at τ = τ 0 , we establish the following Lemma. PROOF : Differentiating (2.6) with respect to τ , we obtain

Using (2.8) and (2.9), we obtain

.

.

.

Thus sign dRe(x) dτ

Thus the result about the existence of Hopf bifurcation is stated in the following theorem: 2

The interior equilibrium E * is locally asymptotically stable for any τ > 0, if (2.11) has no positive real root.

(ii) The interior equilibrium E * is locally asymptotically stable for τ ∈ (0, τ 0 ), if (2.11) has at least one positive real root.

(iii) The system (2.1) undergoes a Hopf bifurcation from the interior equilibrium E * as τ crosses the critical value τ 0 , if γ 0 is a simple root of (2.11), where

.

PROOF : (i) This case has already been proved in the preceding discussion.

(ii) By the definition of τ 0 , (2.11) has no positive real roots for τ ∈ (0, τ 0 ). Hence all the roots of (2.6) have negative real parts. Thus E * is locally asymptotically stable for τ ∈ (0, τ 0 ).

(iii) Suppose γ 0 is a simple root of (2.11). Then we have G (ω 2 0 ) = 0. If G (ω 2 0 ) < 0, then (2.6) has at least one root with positive real part when τ is slightly less than τ 0 , which contradicts conclusion (ii) of Theorem 3. Therefore, we have G (ω 2 0 ) > 0. Hence there exists a Hopf bifurcation for the system (2.1) when τ crosses the critical value τ 0 .

Next, we determine the conditions in terms of the model parameters, for which Hopf bifurcation occurs around the interior equilibrium E * . Accordingly, we define,

Therefore, using Lemma 2.1 and Lemma 2.2 of [33] and Lemma 3.1 of [34] , we can restate Theorem 3 as follows:

Theorem 4 -Suppose R 1 > 1 with τ 0 and ω 0 being already defined in (2.14) . Then (i) The humoral immune activated equilibrium E * is locally asymptotically stable for any τ > 0, if P 4 ≥ 0 and one of the following conditions is satisfied:

(a) ∆ > 0 and y 1 < 0.

(b) ∆ = 0 and y 2 < 0.

(c) ∆ < 0 and y 3 < 0.

(ii) The humoral immune activated equilibrium E * is locally asymptotically stable for τ ∈ (0, τ 0 ), if one of the following conditions is satisfied:

(c) P 4 ≥ 0, ∆ < 0 and there exists at least one y ∈ {y 1 , y 2 , y 3 } such that y > 0 and G(y) ≤ 0.

the humoral immune activated equilibrium E * when τ crosses the critical value τ 0 provided G (ω 2 0 ) > 0.

In the previous section, the existence of bifurcating periodic orbits was investigated. The occurrence of the periodic orbits in a small neighborhood of τ 0 happens either for τ < τ 0 or τ > τ 0 . In order to analyze the stability of the periodic orbits, we estimate (following the approach in [35] ) the maximum length of the delay to preserve the stability of the bifurcating limit cycle.

Taking Laplace transform of (3.1), we obtain Combining all the equations of (3.2) and using (2.7), we obtain,

A necessary and sufficient condition for E * to be locally asymptotically stable is that all the poles ofP (s) must have negative real parts [35] . Therefore, by using the Nyquist criterion [36] , we obtain the sufficient conditions for local asymptotic stability of E * as follows

where H(s) = s 4 + a 1 s 3 + a 2 s 2 + a 3 s + a 4 − cV * e −sτ (s 3 + b 1 s 2 + b 2 s + b 3 ), and η 0 is the smallest positive root of ReH(iη 0 ) = 0, which satisfies ImH(iη 0 ) > 0 as well. Now, (3.3) and (3.4) give

Therefore the sufficient condition for stability of E * is that (3.5) and (3.6) hold simultaneously.

In order to estimate the length of the delay to preserve the stability of E * , we have to find an upper bound of η 0 (independent of τ ). From (3.6), we have,

Using the bounds | sin(η 0 τ )| ≤ 1 and | cos(η 0 τ )| ≤ 1, we obtain from (3.7),

Let η + be the smallest positive root of (3.8) when the equality holds. Then η 0 ≤ η + . Inequality (3.5) can be written as

Adding (3.6) and (3.9), we obtain,

We have sin(η 0 τ ) ≤ η + τ and 1 − cos(η 0 τ ) = 2 sin 2 ( η 0 τ 2 ) ≤ η + 2 τ 2 2 . We now suppose that 

In this section, we present several numerical illustrations to analyze the effect of time delay in the generation of B cells, in addition to investigating the effect of the development rate of B cells, on the dynamical behavior of the model system (2.1). In order to perform the numerical simulation, we chose the parameter values given in Table 1 In order to illustrate the case R 0 < 1, we choose the parameter values listed in Table 1 , except λ = 1 [10, 11] , β 2 = 0.001 [24] , d 3 = 6 [7, 40] , which correspond to R 0 = 0.5775. This scenario is presented in Fig. 1 , which shows that the uninfected hepatocytes (Fig. 1a) with three different initial levels, increase gradually and then finally stabilize at the level T = 100, whereas the infected hepatocytes (Fig. 1b) as well as virions (Fig. 1c) gradually decrease and eventually converge to zero.

The behavior of B cells (Fig. 1d) depends on the initial condition at the beginning, but after a period of time, it also converges to zero. This simulation also indicates that the variation of immune delay does not in any way affect the trajectories of uninfected and infected hepatocytes as well as virion population (solid, dashed and dotted lines merged on same path in Figs. 1a, 1b and 1c) , but the slow convergence of B cells (solid lines of all colors converging to zero earlier than other patterned lines of corresponding colors in Fig. 1d ) occurs due to increase in value of the immune delay. This suggests that the infection-free equilibrium E 0 (100, 0, 0, 0) is globally asymptotically stable for any time delay in the generation of B cells, which supports the theoretical result in Theorem 1.

In order to study the case R H < 1 < R 0 , we choose λ = 1 [10, 11] with the other parameter values being the ones as in Table 1 , which correspond to R 0 = 3.8613 < 1 + and G (ω 2 0 ) = 1.7072 > 0, which satisfies the existence condition of Hopf bifurcation in Theorem 3. We plot the densities (Figs. 3a-8a ) of the four model populations against time and the phase portraits (Figs. 3b-8b) of the trajectories of uninfected and infected hepatocytes as well as virion population for various time delays (τ ). It is clearly noticed in Fig. 3 that E  *  (76.2402, 9 .2375, 3, 79.2967) is locally asymptotically stable when τ = 5 < τ 0 . From the numerical simulation, it can be observed that E * is locally asymptotically stable when τ ∈ [0, τ 0 ), which is obtained theoretically in Theorem 3. Moreover, by Theorem 3, when τ is increased past the critical value τ 0 , a Hopf bifurcation occurs at τ = τ 0 . The result for τ = 6 > τ 0 is shown in Fig. 4 , which implies that E * becomes unstable and consequently a bifurcating periodic solution exists for τ ∈ (τ 0 , τ

2 ). Again, when τ crosses another critical value τ (0) 2 , E * regains local asymptotical stability which is presented for a particular value τ = 14 in Fig. 5 . This shows that E * is asymptotically stable again for τ ∈ (τ

2 ), E * losses stability and a periodic oscillation exists around the equilibrium E * which is exhibited for τ = 25 in Fig. 6 . If τ is increased and exceeds τ (1) 2 , then the interior equilibrium E * becomes stable again, which is illustrated with τ = 35 ∈ (τ

1 ) in Fig. 7 . Further, we simulate the system for a large value of τ , namely, τ = 43 ∈ (τ

2 ), which is demonstrated in Fig. 8 , showing that the stability switch occurs again with a periodic oscillation of the populations.

In the same way, when the bifurcation parameter τ is increased and passes the critical bifurcation values τ 0 , τ

1 , τ (2) 2 and so on; stability of the interior equilibrium E * of the system (2.1) changes from stable to unstable, Hopf bifurcation occurs at these critical values and then regains asymptotic stability from unstable behavior successively. From this numerical discussion, we can conclude that Hopf bifurcation and stability switches occur at the critical values of the bifurcation parameter, which are τ (j) 1 and τ (j) 2 , j = 0, 1, 2, . . . . Therefore, combining the numerical results obtained, one can finally conclude that the interior equilibrium E * (76.2402, 9.2375, 3, 79.2967) is locally asymp-

In the preceding discussion, we have shown that the time delay in generation of B cells has a significant influence on the dynamical behavior of the system. Moreover, as B cells are directly involved in the neutralization of the virions, so the viral infection can be affected depending on the development rate of B cells as well. We therefore investigate (in Fig. 9 ), the effect of the development rate (c) of B cells on the dynamics of uninfected hepatocytes (Fig. 9a) , infected hepatocytes (Fig.   9b ) and viral load (Fig. 9c) as well as antibody response (Fig. 9d) . For this, the numerical simulation is performed for five different values of c, namely, c = 0.01, 0.03, 0.07, 0.1, 0.5 with the other parameters values being retained as in Table 1 . It is observed from Fig. 9 that when B cells are generated very slowly (in case of c = 0.01), the viral load as well as infected hepatocytes increase to a very high level, which results in a decrease of the uninfected hepatocytes to a very low level resulting in the antibody not responding anymore. Moreover, a slight increase in the development rate of B cells, raises the uninfected hepatocytes highly with a significant decrease in virions density. The simulation showed that the development rate of B cells drives the system from stable to unstable and then from unstable to stable again. The findings suggested that a small increment in the development rate of B cells significantly increases the uninfected hepatocytes being neutralizing the virions. Therefore, a high antigenic stimulation in the generation of B cells is beneficial for uninfected hepatocytes. However complete cure from infection not possible by only magnifying the development rate of B cells.

The first author is grateful to Indian Institute of Technology Guwahati for the financial support provided to pursue his Ph.D.

World Health Organization, Global hepatitis report 2017

Ethnic and cultural determinants influence risk assessment for hepatitis C acquisition

Chronic hepatitis C unfection

Vaccine development for hepatitis C, Seminars in liver disease

Hepatitis C viral dynamics in vivo and the antiviral efficacy of interferon-α therapy

Modelling how ribavirin improves interferon response rates in hepatitis C virus infection

Modeling hepatitis C virus dynamics: Liver regeneration and critical drug efficacy

Triphasic decline of hepatitis C virus RNA during antiviral therapy

New kinetic models for the hepatitis C Virus

Hepatitis C virus dynamics and pathology: the role of CTL and antibody responses

Mathematical models of immune effector responses to viral infections: Virus control versus the development of pathology

Dynamical analysis on a chronic hepatitis C virus infection model with immune response

Mathematical modeling of primary hepatitis C infection: Noncytolytic clearance and early blockage of virion production

Analysis of hepatitis C virus infection models with hepatocyte homeostasis

Hepatitis C virus cell-cell transmission in hepatoma cells in the presence of neutralizing antibodies

Dynamical analysis of a class of hepatitis C virus infection models with application of optimal control

Global dynamics for a delayed hepatitis C virus infection model

Stability properties and Hopf bifurcation for a Hepatitis B infection model with exposed state and humoral immunity-response delay

Stability and Hopf bifurcation for a virus infection model with delayed humoral immunity response

Global dynamics of a delayed within-host viral infection model with both virus-to-cell and cell-to-cell transmissions

Stability and Hopf bifurcation in a HIV-1 infection model with delays and logistic growth

Threshold dynamics of HIV-1 virus model with cell-to-cell transmission, cell-mediated immune responses and distributed delay

Stability analysis in delayed within-host viral dynamics with both viral and cellular infections

Threshold dynamics of HCV model with cell-to-cell transmission and a non-cytolytic cure in the presence of humoral immunity

Complex dynamic behavior in a viral model with delayed immune response

Introduction to functional differential equations

Mathematics in population biology

Permanence and positive periodic solution for single-species nonautonomous delay diffusive model

Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission

Perspectives on the basic reproductive ratio

The trade-off between mutual interference and time lags in predatorprey systems

The Hopf bifurcation and its applications

On the zeros of a fourth degree exponential polynomial with applications to a neural network model with delays

Stability analysis of a virus infection model with humoral immunity response and two time delays

Three-species food-chain models with mutual interference and time delays

Regeneration theory

Modeling immune response and its effect on infectious disease outbreak dynamics

Modeling the intracellular pathogen-immune interaction with cure rate

Stability analysis of pathogen-immune interaction dynamics

Modeling the dynamics of hepatitis C virus with combined antiviral drug therapy: Interferon and ribavirin