Effect of DAA therapy in hepatitis C treatment — an impulsive control approach

1.   Introduction

Hepatitis C is an infectious disease caused by the Hepatitis C virus. According to the World Health Organisation (WHO), an estimated $71$ million people globally has been suffering from chronic Hepatitis C syndromes, resulting in cirrhosis and liver cancer. The fatality rate is approximately $39, 900$ every year ^[1]. It is highly blood contagious and at very low risk of sexual and vertical transmission ^[2]. Unhygienic clinical conditions and improper sterilization are the main reasons behind the Hepatitis C infection ^[3].

Hepatitis C syndromes are multiple and demographically manipulated. The virus generally spreads and affects between 2 weeks to $6$ months in the human body. Fever, fatigue, nausea, vomiting, abdominal pain, dark urine, grey face, joint pain, and jaundice are the symptoms of Hepatitis C affected patient. But the worst part of the disease is that the virus sometime remains undiagnosed for a long time and prolonged Hepatitis C infection leads to liver damage (fibrosis and cirrhosis) ^[4].

Chronic HCV infected patients have a risk of fibrosis, cirrhosis, and hepatocellular carcinoma. It is observed that 20–30$\%$ of patients with chronic HCV infection will develop cirrhosis. The ultimate end stage of HCV infection leads to Cirrhosis and Hepatocellular Carcinoma (HCC), which causes death or the need for transplantation. Direct-acting antivirals (DAAs) therapy plays an important role to prevent cirrhosis. The patients already have been existing cirrhosis in need of DAA therapy. DAAs were first approved by the Food and Drug Administration (FDA) in 2011 ^[5].

DAAs act to target specific steps in the HCV viral life cycle. DAAs try to shorten the length of therapy, minimize the side effects, target the virus and improve the virological responses rate. DAAs target one or more of these proteins and enzymes. This results in the delay in viral life cycle as well as it diminishes viral load. Many DAAs are taken in combination with one another or with other medications to improve efficacy and SVR rates ^[6,7].

In recent times mathematical models describing the pathogenic interaction between the human immune system and different kinds of the virus have been of enormous international importance. Appropriate mathematical models can be helpful in answering biologically important questions concerned with the dynamics of the immune response to persistence virus. The effectiveness of drug therapy has been modeled by several authors. Various theoretical studies have been carried out on the mathematical model of HCV transmission dynamics. Nowak and Banghum ^[8] used a mathematical model to explore the effect of individual variation in immune responses on virus load and diversity. They found better indications of CTL responses in the equilibrium virus load, rather than the abundance of virus-specific CTLs. Bonhoeffer et al. ^[9] analysed the virus populations' role of the immune system and resistance of the drug therapy for the HIV or Hepatitis B virus. Neumann et al. ^[3] used a mathematical model to analyse the efficacy of treatment with IFN-a therapy. Avendan et al. ^[10] formulated a mathematical model to describe HCV considering four population susceptible or healthy liver cell, infected liver cell, virus, and CTL responses whereas Zhao et al. ^[11] assumed the incidence rate of the virus model according to Beddington-DeAngelis functional responses. Numerous mathematical models describing the temporal dynamics of HCV have been proposed ^{[9,12,13,14,15,16]}. In all these articles, mathematical modeling plays a pivotal role in understanding and quantifying the biological mechanisms that govern HCV dynamics with or without therapies ^[17,18].

Various research groups have started or are on the way of starting research activities ^[19,11] in the allied fields considering theoretical control model of a different kind of infection in order to gain insights about optimal treatment strategies ^[20]. Different ways of optimal control of treatment are being explored currently by researchers ^{[21,22,23,24]}. Ahmed et al. ^[25] presented a fractional order generalization of Perelson et al. basic hepatitis C virus (HCV) model including an immune response term.

Mathematical models using impulsive differential equations ^[26] have got a lot of attention in the treatment policies of many diseases. For example, Lou et al. ^[27] have used impulsive differential equations to develop a rigorous approach to analyze the threshold behaviours of nonlinear virus dynamics models with impulsive drug effects and to examine the feasibility of virus clearance. Lou and Smith ^[28] have proposed a mathematical model to describe the interaction of HIV virus with CD4$^{+}$T cells in order to describe the fusion process. But in HCV treatment, in our knowledge, there is no such mathematical model with impulsive control therapeutic approach.

Mathematical models play a major role to provide biologically relevant explanations of HCV kinetics under therapy of DAAs. This treatment can reduce the drug resistance and toxicity, as well as having a pharmacokinetic profile that would allow once a day dosing of an all oral combination ^[12]. Impulsive mathematical models may offer a convenient method for rationally designing DAAs therapy based on the properties of single agents.

In this article, the objectives of the research work are, to develop a mathematical model associating the infection by HCV virus and the related treatment technique that deals with impulsive theory, analytical and numerical studies of the fundamental mathematical model from the control viewpoint and to develop the methods for determining pulse therapy based treatment on model predictive controls.

The article is organised as follows. In section 2, we have presented the mathematical model for HCV dynamics. Then we have modified the model incorporating impulsive DAA therapy. Equilibria and their stability of the system without DAA have been analysed in section 3. Dynamics of the impulsive system is studied in section 4. Section 5 contains the numerical results of the main outcomes. In section 6, a final discussion concludes the paper.

2.   Mathematical model formulation

The long term dynamics of HCV infection during the antiviral therapy needs more specific models which consist the complexity of HCV biology and its interplay with the host's immune system. Here we assume the following model which is extended the model as in ^[25]:

$H_s$ represents the uninfected liver cells population, $H_i$ represents the infected liver cells population, and $V$ is the HCV population. The sum of target and infected cells is assumed to remain roughly constant and equal to the total hepatocyte number in normal liver is denoted as $H_{T}$. Using saturated infection rate with maximum transmission rate $\beta$, we have the following model:

Here, the term $\mu_i H_i(1-\frac{H_i}{H_T})$ is considered as the immune response by considering high and low tolerance of the immune system. If $H_i\rightarrow 0$ or $H_i\rightarrow H_T$, then the term $\mu_i H_i(1-\frac{H_i}{H_T})$ vanishes in both cases which suggest that immune response is sufficient.

$\Lambda$ represents the constant production of liver cells, $\mu_s$ is the natural death rate of healthy liver cells and $k$ represents the half saturation constant for the infection. The clearance rate of infected liver cells is $\mu_i$ and HCV proliferates at a rate $p{, }$ and $\mu_v$ is the clearance rate of virus.

We want to study the effect of DAA therapy through impulsive mode. Therefore, we modify the model (2.1) and propose the following model using an impulsive differential equation as follows,

Here $D(t)$ denotes the concentration profile of drug in the human body and $g$ is the rate at which the drug is cleared. We assume that $D(t)$ follows the exponential decay curve ^[28]. $\mu_d$ is assumed as the removal rate of virus by drug therapy.

$D(t^-_k)$ denotes the drug dose concentration immediately before the impulse, $D(t^+_k)$ denotes the concentration after the impulse and $\omega$ is the dose that is taken at each impulse time $t_k$, $k \in \mathbb{N}$.

3.   Dynamics of the system without therapy

In this case we study the system (2.1). There exist two equilibria, namely

$(i)$ the disease-free equilibrium $\overline{E}(\overline{H}_s, \overline{H}_i, \overline{V})$, where $\overline{H}_s = \frac{\Lambda}{\mu_s}, \; \; \overline{H}_i = 0, \; \; \overline{V} = 0$, and

$(ii)$ the endemic equilibrium $E^*(H^*_s, H^*_i, V^*)$, where

and $V^*$ is the positive root of the quadratic equation,

where,

By applying the Descartes rule of sign on (3.1), we can draw the following proposition.

Proposition 1. (i) If $A_3 < 0$, there exists unique endemic equilibrium $E^*$.

(ii) If $A_3 = 0, \; A_2 < 0$, there exists unique endemic equilibrium $E^*$ with $V^* = -\frac{A_2}{A_1}$.

(iii) If $A_3 > 0, \; A_2\geq 0$ then there exists no positive endemic equilibrium.

(iv) Suppose that $A_3 > 0, \; A_2 < 0$.

3.1.   Stability analysis

For the stability analysis of the system (2.1), we need the Jacobian matrix at any equilibrium point $E(H_s, H_i, V)$,

The characteristic equation at the disease-free steady state is

Thus using Routh-Hurwitz condition, the disease-free equilibrium $\overline{E}$ is stable if

Now, we define the basic reproduction number $R_0$ as (a short description on the derivation of $R_0$ is given in Appendix A),

Remark 1. (a) From the condition (i) of Proposition 1, we can conclude that a sufficient condition for the existence of a unique endemic point is $A_3 < 0$, which implies $R_0 < 1$, i.e., when the disease-free equilibrium is stable. Numerically, we have checked that this endemic equilibrium is always unstable (See Figure 1).

(b) Form condition (iv) of Proposition 1, two different endemic equilibrium points exist if $A_3 > 0, \; A_2 < 0$ and $D = A^2_2-4 A_1 A_3 > 0$. Here, $A_3 > 0$ implies $R_0 > 1$. But one endemic equilibrium point (with $V^* = \frac{-A_2+\sqrt{D}}{2A_1}$) is unstable and another one (with $V^* = \frac{-A_2-\sqrt{D}}{2A_1}$) is stable (See Figure 1).

Now, at the endemic equilibrium $E^*$, the characteristic equation is

where,

According to Routh-Hurwitz criteria, the endemic steady state is stable if

We have the following theorem.

Theorem 1. Disease-free equilibrium is stable if $R_0 < 1$ and unstable otherwise. Forward transcritical bifurcation occurs at $R_0 = 1$. The endemic steady state $E^*$ is stable if the condition $(3.7)$ holds.

Remark 2. In light of the fact that $R_0$ is monotonically decreasing with increasing $\mu_v$, this suggests that eradication of disease, as represented by a stable disease-free steady state $E_0$ is possible if $R_0 < 1$. The available means to achieve this is by increasing the clearance rate of virus, $\mu_v$. This can be done using DAA therapy ^[5].

4.   Dynamics of the system with impulsive therapy

Now, we shall analyse the dynamics of the drug (or combined drug) and its effects on the system population. The aim is to find a better treatment strategy which can suppress the viral load and also inhibit viral entry into the host cell. To study the effect of therapy in regular intervals, we study the system of impulsive differential equations given in (2.2).

Remark 3. In this article, our main aim is to justify the effect of antibody therapy in impulsive modes. We have not carried out the stability analysis for the endemic state in presence of impulsive therapy. Actually there will not exist any equilibria, rather equilibria like periodic orbits ^[28].

4.1.   Dynamics of the drug

Dynamics of the drug is governed by the following impulsive differential equation,

To study the dynamics of the perfect drug adherence, we assume the dosing interval is $\tau$ defined as

The solution of the impulsive differential equation (4.1) is

Calculating the least value of the concentration of $D(t)$ for the perfect drug adherence with fixed interval length $\tau > 0$, we get

which is the required concentration of dosage to control the virus. If $D(0) = 0, \; D(t^+_1) = \omega$, then

Hence,

Thus the start and end point of periodic trajectories are $D_u = \frac{\omega }{1-e^{-g\tau}}$ and $D_l = \frac{\omega e^{-g\tau}}{1-e^{-g\tau}}$.

For perfect therapy, the antibody response after the $n^{th}$ dosage is

To control the virus and avoid resistance, the minimum value $\bar{D}$ of the periodic orbit must satisfy

Remark 4. If we can restrict the dosing interval of $\tau$ satisfying the condition $0\leq \tau \leq \tau_{max}$, then the disease can be controlled. For $\tau > \tau_{max}$, the disease progression continues.

The lower and upper limits for $D(t)$ are $D_l$ and $D_u$. The critical dosage $D_c$ must satisfy, $D_l > D_c$ which holds true iff $\tau < \frac{1}{g} \ln (\frac{\omega +D_c}{D_c})$. Further $D_u < D_c$ implies that $\tau > \frac{1}{g} \ln (\frac{D_c}{D_c-\omega })$. By helping this, we can conclude the result by the following theorem.

Theorem 2. A treatment regimen $(\omega, \tau)$ is successful if $\tau < \tau_s = \frac{1}{g}\ln (\frac{\omega +D_c}{D_c})$. A treatment is unsuccessful if $\tau > \tau_u = \frac{1}{g}\ln (\frac{D_c}{D_c-\omega })$. For a fixed dose $\omega$, $\tau_s$ is the longest dosing interval that guarantees a successful treatment regimen.

Remark 5. For a fixed $\tau$, we can find out the safe dose $\omega_s$ and unsafe dose $\omega_u$. Solving $\frac{\omega e^{-g\tau}}{1-e^{-g\tau}} = D_c$ for $\omega$ leads to

and solving $\frac{\omega }{1-e^{-g\tau}} = D_c$ for $\omega$ leads to

If $\tau_s < \tau < \tau_u$ or $\omega_u < \omega < \omega_s$, a decision/conclusion can not be reached because $R_0$ will fluctuate around 1.

4.2.   Stability of disease-free periodic orbit

There does not exist particular equilibrium point of an impulsive system but equilibrium-like periodic orbits can be evaluated. Using the following analysis, we can show that there are two periodic orbits namely the disease-free periodic orbit and the endemic periodic orbit for the system (2.2) with impulses. Using the results from ^[31,32,33], we have the following result.

Lemma 1. The system (4.1) has unique positive periodic globally asymptotically stable solution $D_1(t)$ with period $\tau = t_{k+1}-t_k$

The general solution of the system (4.1), $D(t)$ can be written as

On the above basis, we study the stability of periodic orbits. We only focus on the disease-free periodic orbit deriving the following theorem.

Theorem 3. The disease free periodic solution $(\tilde{H_s}, 0, 0, \tilde{D})$ of the system (2.2) is locally asymptotically stable if

where,

Proof. Variational matrix at $(\tilde{H_s}, 0, 0, \tilde{D})$ is given by,

The monodromy matrix $\mathbb{D}$ of the variational matrix $J_v(t)$ is

where $I$ is the identity matrix. We thus have: $\mathbb{D}(\tau) = \rm{diag} (\lambda_1, \lambda_2, \lambda_3, \lambda_4)$. Here, $\lambda_i, i = 1, 2, 3, 4$, are the Floquet multipliers given by

Here $A = \mu_i+\mu_s+\mu_d \tilde{D}$ and $B = \mu_i(\mu_v+\mu_d\tilde{D})-p\beta \tilde{H_s}$. Clearly $\lambda_{1}$ and $\lambda_4 < 1$. It is easy to check that $A^2-4B > 0$. Further if $B > 0$, then $\lambda_{2, 3} < 1$. Thus, according to Floquet theory, the periodic solution $(\tilde{H_s}(t), 0, 0, \tilde{D}(t))$ of the system (2.2) is locally asymptotically stable if the condition (4.9) holds. Note that $\tilde{R_0} < 1$ if $B > 0$.

5.   Numerical simulations

In this section, we observe the dynamical behaviours of system (2.1) and impulsive effect of the therapy using the model (2.2) through numerical simulations taking the parameters from Table 1. We have assumed $H_s(0) = 200\; mm^{-3}, \; H_i(0) = 20\; mm^{-3}, \; V(0) = 20\; mm^{-3}$ as initial biological conditions of the model populations.

In Figure 1, we have seen the transcritical bifurcation at $R_0 = 1$. When $R_0 < 1$, the disease-free equilibrium $E_0$ is stable and a unique endemic equilibrium $E^*$ with $V^* = \frac{-A_2+\sqrt{D}}{2A_1}$ exists but it is unstable (Figure 1a, red dotted line). There exists two endemic equilibria exist for $R_0 > 1$. Hence Proposition 1, condition (iv) is verified. Using the conditions in Theorem 1, we have seen that the endemic equilibrium $E^*$ with $V^* = \frac{-A_2+\sqrt{D}}{2A_1}$ is unstable but the endemic equilibrium $E^*$ with $V^* = \frac{-A_2-\sqrt{D}}{2A_1}$ is stable. From this figure, we have also observed the effect of the infection rate $\beta$. For lower infection rate, $R_0 < 1$, the disease-free state is stable (Figure 1a, green dotted line). In Figure 1b, region of stability of the equilibria are presented in $\beta-\Lambda$ parameter plane. It can be observed that when the product of $\beta$ and $\Lambda$ crosses a threshold value (which corresponds to $R_0 > 1$), the disease-free equilibrium $E_0$ becomes unstable and the endemic equilibrium $E^*$ (with $V^* = \frac{-A_2-\sqrt{D}}{2A_1}$) exists and stable.

Numerical solution of system (2.1) without any control measures is plotted in Figure 2. In this case, we have chosen the set of parameters so that $R_0$ is greater than unity, i.e., the system is endemic. It can be observed that the system trajectories converge to the endemic equilibrium $E^*$ with $V^* = \frac{-A_2-\sqrt{D}}{2A_1}$. According to Theorem 1, this $E^*$ is stable. But there exists another endemic equilibrium $E^*$ (see Figure 1) which is unstable.

In Figure 3, we have plotted the numerical solution of system (2.2) taking two different doses for fixed dosing interval $\tau$. The system becomes free of disease/infection as infected liver cell $H_i$ and virus population are not appearing for $\omega = 10$ (see the red line). But for lower dosing rate ($\omega = 5$), the system is endemic in nature as both infected liver cell and virus are present in the system (indicated by blue line). In Figure 4, we have plotted the numerical solution of system (2.1) taking two different dosing interval $\tau$ for a fixed dosing rate $\omega$. System becomes free of disease/infection as infected liver cell $H_i$ dies out for $\tau = 7$. But for wider dosing interval ($\tau = 12$), the system is endemic in nature as indicated by infected liver cell $H_i$.

In Figures 5 and 6, we plotted the numerical solution of system (2.2) for varying dose and dosing interval respectively. Figure 5 is plotted by changing dose at fixed interval namely $\omega = 5, \; 6, \; 7, \; 8\; \mu$ mg/L etc.; that means during successive dosing attempts, the dosing rate should be increased by $1\; \mu$ mg/L from the previous dosing. Figure 6 is plotted by changing dosing interval ($\tau$) at increasing order such as 1, 3, 5, 7 days etc. which means in every time dosing interval time is increased by $2$ days from the previous dosing. In this way we can find the suitable dose and dosing intervals. This enables us to reduce the disease with minimum side effects.

6.   Discussion and conclusions

In this article, we have formulated a basic mathematical model to study the infection for HCV on human population. In this model, we have considered uninfected liver cells, infected liver cells and virus population. Then we formulate the impulsive differential equations as the therapy is given periodically. By using the impulsive differential equation, we mainly worked on the suitable interval of therapy period and size of the dosage so that disease can be controlled.

The model without impulse is an extension of the model in ^[10]. The authors of ^[10] have proposed a model for HCV taking four populations. Here we also assume four populations and the model is more novel as our model can describe more phenomenon of HCV dynamics. For example, we have assumed the infection term as $\frac{\beta H_s V}{1+kV}$ in place of ${\beta H_s V}$ and also incorporated impulse control in the model. Our model contains impulsive differential equations.

We have observed that the model system has two equilibria, one is disease free and another is endemic equilibrium. The disease free equilibrium is asymptotically stable if the basic reproduction number is below unity, but when the basic reproduction number is grater than unity, then the disease free equilibrium becomes unstable.

From the above analysis it is clear that infection is minimum at considerably higher dose $\omega$, as well as lower dosing interval $\tau$. But, since infection rate depends on densities of virus and infected liver cells, the therapy with a constant dose and at constant time interval is not appropriate to maintain disease-free situation for the whole duration. Variations in these two parameters should be studied for maximum incidence to maximize eradication of the disease. During maximum infection, the dosing rate should be higher and time interval should be lower but the interval time should be increased at ascending order during disease eradication process. From this study it can be seen that at higher time interval ($\tau = 7$ days) and also at comparatively lower therapy rate (e.g., $\omega = 12\; \mu$ mg/L) the system reaches a disease-free stable periodic state.

In this study, the analytical and numerical findings reveal the theoretical eradication of the disease. However, in reality it is quite difficult to conclude the complete eradication of the disease. Because HCV has other reservoirs like spleen, intestine, pancreas, heart, kidneys, brain, lymph nodes, dendritic cells, B and T lymphocytes ^[34,35], from which the virus can come back once again before it is controlled below the level of detection. However, our model results in large-scale declines in viral load, which can be practically estimated by using impulsive differential equations. From the mathematical point of view "theoretical eradication" should be understood to mean a substantial drop of viral load below its threshold value but not complete eradication ^[36].

From this research, we can conclude that DAA therapy with regular adherence to HCV can be effective at controlling the virus. Also for a perfect adherence of drug dose interval and drug dosage, cellular infection can be controlled and immune system performs accurately. However, as of adherence delays, it results in extreme variations in the system. The critical therapy period suggests that careful follow-up must be taken. Hence optimal level of therapy period as well as size of pulse therapy affect the disease progression and disease replication. Hence, the patients should be advised on the significance of adherence to this DAA therapy against HCV.

Acknowledgements

YT's research is supported by Japan Society for the Promotion of Science "Grant-in-Aid 20K03755".

Conflict of interest

All authors declare no conflicts of interest in this paper.

Appendix A

We have followed the method used by Hefferman et al. ^[38] for the derivation of $R_0$ of the system (2.1). We consider the next generation matrix M which comprised of two parts namely F and V, where

Here, ${\bf F}_i$ are the new infections and ${\bf V}_i$ are for the transfer of infections from one compartment to another, and $\bar{E}$ is the disease-free equilibrium. The basic reproduction number is the dominant eigenvalue of the matrix ${\bf M} = {\bf F}{\bf V}^{-1}$ and is denoted as $R_0$.