www.biomathforum.org/biomath/index.php/biomath ORIGINAL ARTICLE Modelling and analysis of a within-host model of hepatitis B and D co-infections Plaire Tchinda Mouofo∗, Jean Jules Tewa†, Samuel Bowong‡ ∗ Department of Mathematics, University of Yaounde I, P.O. Box 812 Yaounde, Cameroon tchindaplaire@yahoo.fr †Department of Mathematics and Physics, National Advanced School of Engineering (Polytechnic), University of Yaounde I, P.O. Box 8390 Yaounde, Cameroon tewajules@gmail.com ‡Department of Mathematics and Computer Science, Faculty of Science, University of Douala, PO Box 24157 Douala, Cameroon sbowong@gmail.com Received: 17 July 2017, accepted: 21 July 2018, published: 7 August 2018 Abstract—The Hepatitis delta virus (HDV) is a defect RNA virus that requires the presence of the hepatitis B virus (HBV) for cellular infection. A co- infection may result in a more severe acute disease and a higher risk of developing acute liver failure compared with those infected with HBV alone. At the present time, there has been very little to the modeling of HDV. The derivation and analysis of such a mathematical model poses difficulty as it requires the inclusion of (HBV). In this paper, a within-host model for the co-interaction of HDV and HBV is presented and rigorously analyzed. We calculate the basic reproduction number (R0), the disease-free equilibrium, boundary equilibrium, which we define as the existence of one disease along with the complete eradication of the other disease, and the co-infection equilibrium. We determine stability criteria for the disease-free and boundary equilibrium. We also use the optimal control theory to assess the disease control. Numerical simulations have been presented to illustrate analytical results. Keywords-Hepatitis D, Hepatitis B, Immune sys- tem, Basic reproduction number, Optimal control. AMS subject classifications: 34A34, 34D23, 34D40, 92D30 I. INTRODUCTION Hepatitis D is a liver disease caused by the hepatitis D virus (HDV), a defective virus that needs the hepatitis B virus to exist. The hepatitis D virus requires the outer coating of the hepatitis B virus called the surface antigen in order to re- produce itself in a human host. The virus currently Copyright: c©2018 Mouofo et al. This article is distributed under the terms of the Creative Commons Attribution License (CC BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Citation: Plaire Tchinda Mouofo, Jean Jules Tewa, Samuel Bowong, Modelling and analysis of a within-host model of hepatitis B and D co-infections, Biomath 7 (2018), 1807219, http://dx.doi.org/10.11145/j.biomath.2018.07.219 Page 1 of 17 http://www.biomathforum.org/biomath/index.php/biomath https://creativecommons.org/licenses/by/4.0/ https://creativecommons.org/licenses/by/4.0/ http://dx.doi.org/10.11145/j.biomath.2018.07.219 P. T. Mouofo, J. J. Tewa, S. Bowong, Modelling and analysis of a within-host model of hepatitis ... infects 15 million worldwide, nearly all adults, and it is most common among injecting drug user populations and in countries bordering the Mediterranean Sea. HDV is transmitted through blood and body fluids, similar to the hepatitis B virus. There are two types of HDV infection, co- infection and super-infection. Co-infection oc- curs when a patient is simultaneously infected with HDV and HBV. The majority of these pa- tients completely recover but there is a higher rate of fulminant hepatic failure and death than with HBV infection alone. Super-infection occurs when someone with an existing chronic HBV infection becomes infected with HDV. These pa- tients usually experience a sudden worsening of liver disease. Patients with hepatitis B who be- come chronically infected with HDV experience a very high rate of cirrhosis and end stage liver disease, which makes this super-infection a very dangerous disease. There is no specific treatment for HDV infec- tion. The most common therapeutic approach is based on the administration of interferon-α. How- ever, the clinical response is variable, and in most cases reversible upon interruption of treatment [1], [2]. The concomitant use of antiviral drugs like ribavirin or lamivudine, showed no significant benefits in the treatment of hepatitis delta patients [3], [4]. Although these drugs may have some inhibitory effect on HBV replication, they do not suppress HDV replication probably due to the fact that HBsAgs expression, at least in part, seems not to be affected. The use of mathematical models to study dy- namics of virus infections may represent a pow- erful approach to simulate the course of infection and predict the potential response to different therapies. They have been previously developed for a number of pathologies including HBV and HCV [5], [6], [7], [8], [9], [11], [12], [34]. the humoral immune response is universal and nec- essary to eliminate or control the disease after viral infection [23]. Therefore, several mathemat- ical models have been proposed to describe the virus dynamics with humoral immunity [24], [25], [26], [27], [28], [30], [31], [32]. Mostafa et al. introduce an improved HBV model with stan- dard incidence function, cytotoxic T lymphocytes (CTL) immune response, and take into account the effect of the export of precursor CTL cells from the thymus and the role of cytolytic and non- cytolytic mechanisms [33]. Hattaf et al. [29] study the global stability of a generalized model of a viral dynamic that includes the adaptive immune response, represented by Cytotoxic Lymphocyte T-cell (CTL-cell ). So, their work does not take in consideration the role of the innate immune response. Noura Yousfi et al. [38] investigate a new mathematical model that describes the interactions between Hepatitis B virus (HBV), liver cells (hepatocytes), and the adaptive immune response. More recently, the numerical simulation of the spread of HDV and HBV in a population was reported [13]. However, a mathematical model to study HDV and HBV dynamics in infected individuals is still lacking. Moreover, to our knowledge, none of the existing models take into account the reaction of immune CTL cells. The main interest in studying HBV and HDV model is to understand the long and short term behavior of the dynamics of both diseases and to predict whether the diseases will die out or will persist. In this paper, the dynamical behavior of a co- interaction of HDV and HBV virus model with CTL immune responses is studied. The main objective is in carrying out a detailed qualitative analysis of the resulting model. The existence and the uniqueness results of the solution are discussed. We compute the basic reproduction number. We also investigate the existence of equilibria and study their stability. The model is used to determine the optimal methodology for administering anti-viral medication therapies to fight HBV and HDV infection. In particular, we investigates the fundamental role of chemotherapy treatment in controlling the virus reproduction. A characterization of the optimal control via adjoint variables is also established. We obtain an opti- Biomath 7 (2018), 1807219, http://dx.doi.org/10.11145/j.biomath.2018.07.219 Page 2 of 17 http://dx.doi.org/10.11145/j.biomath.2018.07.219 P. T. Mouofo, J. J. Tewa, S. Bowong, Modelling and analysis of a within-host model of hepatitis ... mality system that we seek to solve numerically. It is our view that this study represents the very first modelling work that provides an in-depth analysis of the qualitative dynamics of HBV-HDV co-infection and its control. The paper is organized as follows. In the next section, we propose a new mathematical model for HBV infection alone that takes into account the CTL immune system. To be more realistic, we have assumed that the infection rate is given by the standard incidence function [16]. The model is rigorously analyzed. In Section 3, we formulate and analyze a realistic mathematical model for HBV-HDV co-infection, which incorporates the key epidemiological and biological features of each of the two diseases. Optimal control ap- proach is applied in Section 4, in order to find the best way to fight the co-infection between the two diseases. We end this paper with a brief discussion and conclusion. II. THE HBV MODEL A. The model description We consider the HBV model given by the following differential equations:  ẋ = λ−dx x− β(1 −η)xv x + y ẏ = β(1 −η) xv x + y −ay y −δ y I v̇ = k(1 −ε) y −dv v, İ = ρy I + pI −q I2 (1) where x is the number of uninfected liver cells, y is the number of infected liver cells, v is the number of free virus, and I is the number of CTL cells. All the parameters λ, dx, η, β, ay, δ, ε, k, dv, ρ, p, and q are positive. dx, ay and dv are the death rates of uninfected liver cells, the infected cells and free virus, respectively. The constant parameter λ represents the production of the liver cells. β is the contact rate between uninfected cells and free virus. Free virus is produced from infected cells at rate ky. Infected cells are removed at rate δI by CTL immune responses. The virus-specific CTL cells proliferate at rate ρy by contacts with infected cells. The parameter η is the efficacy of inhibiting new virus infections as a consequence of virus clearance and ε the efficacy of inhibiting viral production from infected cells. Of course the behaviour when there is no treatment is obtained by setting η = ε = 0. The parameter p denotes the proliferation rate of immune cells and q the density-dependent rate of immune cells suppression. more precisely, we suppose that the immune cells expand at a net rate p, which encapsulate the positive feedback upon the immune system. The parameter q comes from the fact that we assume a regulatory negative feedback force such as the effect of cell density, inhibitory cytokines or natural apoptosis, which oparates to suppress immune population growth. we suppose in this case that immune population is suppressed at a net rate q which is proportional to the square of its density (qI2). So, the term pI − qI2 can be written as pI ( 1 − Ip q ) which express a logistic law for evolution of immune population in the absence of infected cells. B. Analysis of the model Herein, we present some basic results, such as the positive invariance of model system (1), the boundedness of solutions, the existence of equilibria and and its stability analysis. 1) Positivity and boundedness of solutions: The following result guarantees that model system (1) is biologically well behaved and its dynamics is concentrated on a bounded region of R4+. More precisely, the following result holds. Theorem 1. Let R4+ = {(x,y,v,I) ∈ R4 : x ≥ 0, y ≥ 0,v ≥ 0,I ≥ 0}. Then, R4+ is positively invariant under the flow induced by model system (1). Moreover, the region ∆ = { (x,y,v,I) ∈ R4 : x + y ≤ λ dx , v ≤ kλ(1 −ε) dxdv , p q ≤ I ≤ p + ρλ dx q } is positively invariant and absorbing with respect to model system (1). Biomath 7 (2018), 1807219, http://dx.doi.org/10.11145/j.biomath.2018.07.219 Page 3 of 17 http://dx.doi.org/10.11145/j.biomath.2018.07.219 P. T. Mouofo, J. J. Tewa, S. Bowong, Modelling and analysis of a within-host model of hepatitis ... Proof: No solution of model system (1) with initial conditions (x(0),y(0),v(0),I(0)) ∈ R4+ is negative. In fact, for (x(t),y(t),v(t),I(t)) ∈ R4+, we have ẋ |x=0= λ > 0, ẏ |y=0= (1 − η)βv ≥ 0, v̇ |v=0= (1 − ε)k ≥ 0, İ |I=0= 0 ≥ 0, this immediately implies that all solu- tions of model system (1) with initial condition (x(0),y(0),v(0),I(0)) ∈ R4+ stay in the first quadrant. For the invariance property of ∆, it suffices to show that the vector field, on the boundary, does not point to the exterior. Adding the first and second equations of model system(1) yields on the boundary of ∆: d(x + y) dt ∣∣∣∣ x+y= λ dx = λ−dxx−ayy −δyI |x+y= λ dx ≤ (λ−dx(x + y))|x+y= λ dx = 0. Similarly, we get dv dt ∣∣∣∣ v= kλ(1−ε) dxdv ≤ kλ(1 −ε) dx −dvv ∣∣∣∣ v= kλ(1−ε) dxdv = 0, dI dt ∣∣∣∣ I= p q ≥ (p−qI) I|I= p q = 0, i.e I(t) ≥ p q ∀t ∈ [0, +∞), and dI dt ∣∣∣∣ I= p+ ρλ µ q ≤ ( ρλ dx + p−qI ) I ∣∣∣∣ I= p+ ρλ µ q = 0. Therefore, solutions starting in ∆ will remain there for t ≥ 0. Now, we prove the attractiveness of the trajec- tories of model system (1). To do so, from model system (1), one has d(x + y) dt ≤ λ−dx(x + y). Therefore, lim sup t→∞ (x + y)(t) ≤ λ dx . Similarly, since dv dt ≤ kλ(1 −ε) dx −dvv, one has lim sup t→∞ v(t) ≤ kλ(1 −ε) dxdv . Concerning the variable I, we have dI dt ≤ ( ρλ dx + p ) I −qI2, which implies that I(t) ≤ 1 q p+ ρλ dx + ( 1 I(0) − q p+ ρλ dx ) e − ( p+ ρλ dx ) t . So, I is bounded and hence, ∆ is attracting, that is, all solutions of model system (1) eventually enters ∆. This concludes the proof. 2) Basic reproduction number and equilibria: The basic reproduction number is defined as the average number of secondary infections produced by one infected cell during the period of infection when all cells are uninfected. This threshold is obtained at the virus free equilibrium. The virus free equilibrium is obtained by setting v = 0 in all equations of model system (1) at the equilibrium. Thus, when v = 0, we have P0 = (x∗, 0, 0, 0) and P1 = (x ∗, 0, 0,I∗) where x∗ = λ dx and I∗ = p q . Since I(t) ≥ p q , the virus free equilibrium is P1. We use the method of van den Driessche[17] to compute the basic reproduction number. Us- ing the notations of van den Driessche and Watmough[17], for model system (1), we have F = ( 0 β(1 −η) k(1 −ε) 0 ) and V = ( ay + δI ∗ 0 0 dv ) . Thus, the basic reproduction number is given by: R0 = k(1 −ε)β(1 −η) dv(ay + δI∗) . (2) Biomath 7 (2018), 1807219, http://dx.doi.org/10.11145/j.biomath.2018.07.219 Page 4 of 17 http://dx.doi.org/10.11145/j.biomath.2018.07.219 P. T. Mouofo, J. J. Tewa, S. Bowong, Modelling and analysis of a within-host model of hepatitis ... From theorem 2 of Van Den Driessche[17], we have the following result. Lemma 1. The virus-free equilibrium P1 of the model system (1) is locally asymptotically stable (LAS) if R0 < 1, and unstable if R0 > 1. We now study the existence of equilibria of model system (1). Setting the right-hand sides of model system (1) equals to zero gives λ−dx x− β(1 −η)xv x + y = 0, (3) β(1 −η) xv x + y −ay y −δ y I = 0, (4) k(1 −ε) y −dv v = 0, (5) ρy I + pI −q I2 = 0. (6) Model system (1) has always equilibrium P1 = (x∗, 0, 0,I∗) which is the virus free equilibrium and represents the state when the viruses are absent. To find the endemic equilibrium, we look for P2 = (x̄, ȳ, v̄, Ī) that represents the state in which both the viruses and CTL cells are present. Let γ = dv k(1 −ε) , a = ρdv qk(1 −ε) , b = p q and assume that R0 > 1. From Eqs.(5) and (6), one obtains ȳ = γv̄ and Ī = av̄ + b. (7) Substituting the above expression of ȳ in Eq.(4) yields v̄ ( β(1 −η) x̄ x̄ + γv̄ −ayγ −δγ(av̄ + b) ) = 0 (8) Since v̄ 6= 0, Eq.(8) leads to x̄ = γ2v̄(ay + δb + δav̄) β(1 −η) −γ(ay + δb + δav̄) . (9) Since x̄ > 0, one can deduce that v̄ < β(1 −η) −γ(ay + δb) γδa . Note that β(1 −η) −γ(ay + δb) = dv(ay + δb) k(1 −ε) [R0 − 1]. Thus, x̄ > 0 implies that v ∈]0, ṽ1[, where ṽ1 = dv(ay + δb) γδak(1 −ε) [R0 − 1]. (10) Moreover, since x̄ < λ dx , the above relation implies that γ2δadxv̄ 2 + [ dxγ 2(ay + δb) + λγδa ] v̄ −λ [ β(1 −η) −γ(ay + δb) ] < 0. (11) Since R0 > 1, the discriminant of (11) is delta = [ dxγ 2(ay + δb) + λγδa ]2 +4λγ2δadx [ β(1 −η) −γ(ay + δb) ] > 0. Thus, the condition x̄ < λ dx implies that v̄ ∈]0, ṽ2[ where ṽ2 = − [ dxγ 2(ay + δb) + λγδa ] + √ delta 2γ2δadx (12) We have the following result. Lemma 2. Let ṽ1 and ṽ2 given in Eqs. (10) and (12), respectively. Then, ṽ2 < ṽ1 whenever R0 > 1. Proof: Since R0 > 1, applying some techni- cal manipulation one can show that γdx [ β(1−η)−γ(ay + δb) ] + dxγ 2[ay + δb] > 0 (13) is equivalent to ṽ1 > ṽ2. Remark 1. 1) From Lemma 2, conditions v̄ ∈ ]0, ṽ1[ and v̄ ∈]0, ṽ2[ can be limited to v̄ ∈ ]0, ṽ2[. Thus, when R0 > 1, model system (1) may have an endemic equilibrium P2 = (x̄, ȳ, v̄, Ī) with v̄ ∈]0, ṽ2[. 2) We point out that v̄ = ṽ2 implies that x̄ = λ dx . Substituting Eq. (7) and Eqs. (9) into (3), we obtain the following cubic equation in v̄: a3v̄ 3 + a2v̄ 2 + a1v̄ + a0 = 0, (14) Biomath 7 (2018), 1807219, http://dx.doi.org/10.11145/j.biomath.2018.07.219 Page 5 of 17 http://dx.doi.org/10.11145/j.biomath.2018.07.219 P. T. Mouofo, J. J. Tewa, S. Bowong, Modelling and analysis of a within-host model of hepatitis ... where a3 = [ d2δρ qk2(1 −ε)2 ]2 , a2 = d2v ( ay + δ p q ) k2(1 −ε)2 − dxδρd 3 v qk3(1 −ε)3 − δρd3v ( ay + δ p q ) qk3(1 −ε)3 [ R0 − 1 ] , a1 = − λd2vδρ qk2(1 −ε)2 − dxd 2 v ( ay + δ p q ) k2(1 −ε)2 − γdv ( ay + δ p q )2 k(1 −ε) [ R0 − 1 ] , a0 = λdv ( ay + δ p q ) k(1 −ε) [ R0 − 1 ] . The coefficient a3 is always positive. Also, if R0 > 1, then a0 > 0 and a1 < 0. Using the Descartes’ rules of sign, if R0 > 1, the above equation given in (14) has exactly nought or two positive solutions. Let us consider the following polynomial P : [0, ṽ2]→R v̄ 7→P(v̄) =a3v̄3 +a2v̄2 +a1v̄+a0. (15) We have • P(0) = a0 > 0 since R0 > 1 • P(ṽ2) = − β(1 −η)k(1 −ε)λṽ2 λk(1 −ε) + dxdvṽ2 < 0 Since P(0) > 0 and P(ṽ2) < 0, Eq.(14) has only one or three solutions on ]0,v2[. By the Descartes’ rules of sign, Eq.(14) cannot have three solutions on ]0, ṽ2[. We have proved the following result. Theorem 2. If R0 > 1, model system (1) has exactly one endemic equilibrium P2 = (x̄, ȳ, v̄, Ī) where x̄, ȳ and Ī are defined in Eqs. (7) and (9) with v̄ ∈]0, ṽ2[ satisfying Eq. (14). Remark 2. When R0 < 1, we have x̄ < 0. So, there is no endemic equilibrium. 3) Stability of equilibria: Here, we study the stability of equilibria. We have the following result about the global stability of the virus free equilibrium. Theorem 3. When R0 < 1, then the virus free equilibrium P1 = (x∗, 0, 0,I∗) is globally asymp- totically stable in ∆. Proof: We consider the following LaSalle- Lyapunov candidate function: L(t) = b1y + b2v, (16) where b1 and b2 are positive constants to be determined later. Using the fact that x ≤ x∗ and I ≥ I∗, its time derivative along the trajectories of (1) satisfies L̇ = ( b1β(1 −η)x x + y − b2dv ) v +[b2k(1 −ε) − b1(ay + δI)]y, ≤ ( b1β(1 −η)x∗ x∗ + y − b2dv ) v +[b2k(1 −ε) − b1(ay + δI∗)]y. The constants b1 and b2 are chosen such that the coefficient of y are equal to zero, that is, b1 = k(1 −ε) and b2 = ay + δI∗. Then, we obtain L̇ ≤ (ay + δI∗)dv(R0 − 1)v. (17) Thus, L̇(t) ≤ 0 when R0 < 1. By LaSalle’s invariance principle, the largest invariant set in {(x,y,v,I) ∈ R4+, L̇(t) = 0} ⊂ ∆ is reduced to the virus free equilibrium. This proves the global asymptotic stability of P1 = (x∗, 0, 0,I∗) on ∆ [18] (Theorem 3.7.11, page 346). This achieves the proof. Now, we study the stability of the endemic equilibrium. The Jacobian matrix of model system (1) at the endemic equilibrium P2 = (x̄, ȳ, v̄, Ī) is J3 =   −dx−U W −V 0 U −W−ay−δĪ V −δȳ 0 k(1 −ε) −dv 0 0 ρĪ 0 −qĪ   , Biomath 7 (2018), 1807219, http://dx.doi.org/10.11145/j.biomath.2018.07.219 Page 6 of 17 http://dx.doi.org/10.11145/j.biomath.2018.07.219 P. T. Mouofo, J. J. Tewa, S. Bowong, Modelling and analysis of a within-host model of hepatitis ... where U = (1 −η)βv̄ȳ (x̄ + ȳ)2 , V = (1 −η)βx̄ x̄ + ȳ and W = (1 −η)βv̄x̄ (x̄ + ȳ)2 . The characteristic equation of J3 is ξ4 + c1ξ 3 + c2ξ 2 + c3ξ + c4 = 0, (18) where c1 = dx + U + W + ay + δ Ī + dv + qĪ, c2 = dxW + dxay + dxδ Ī + dxdv + dxqĪ + Uay + Uδ Ī + Udv + UqĪ + Wdv + qĪW + qĪay + qĪ 2δ + qĪdv + δ ȳρ Ī, c3 = Uaydv + Uδ Īdv + dxWdv + dxqĪW + dxqĪay + dxqĪ 2δ + dxqĪdv + dxδ ȳρ Ī + UqĪay + UqĪ 2δ + UqĪdv + Uδ ȳρĪ + WdvqĪ + δ ȳρ Īdv c4 = dxWdvqĪ + dxδ ȳρ Īdv + Uδ ȳρĪdv + UaydvqĪ + Uδ Ī 2dvq. It is clear that c1 > 0, c3 > 0 and c4 > 0. After some technical computation, we obtain that c1c2c3 − c23 − c 2 1c4 is positive. Hence, we have proved the following result. Theorem 4. The endemic equilibrium P2 of model system (1) is locally asymptotically stable if R0 > 1. III. THE HBV-HDV CO-INFECTION MODEL A. Model construction In this section, we incorporate the HDV in- fection into the previous model in order to ob- tain a mathematical models that can describe the dynamics of HDV and HBV co-infection. We add three state variables, namely the HDV viral load, hepatocyte infected with HDV and those co- infected with both HBV and HDV. So, we use the following state variables: • x(t) the number of uninfected cells at time t, • y(t) the number of HBV infected cells at time t, • z(t) the number of HDV infected cells at time t, • w(t) the number of infected cells with both HBV and HDV at time t, • v1(t) the HBV viral load at time t, • v2(t) the HDV viral load at time t, • I(t) the number of CTL cells at time t. We make the following hypothesis: 1) Uninfected liver cells (x) can only be in- fected by HBV virions alone and become y, or by HDV virions alone and become z. 2) Infected cells by HBV virions can be super- infected by HDV and become w. 3) Infected cells by HDV virions can be super- infected by HBV and become w. Putting the above formulations and assumptions together gives the following system of differential equations:  ẋ=λ−dxx− β1(1 −η)xv1 P − β2(1−η)xv2 P , ẏ = β1(1−η)xv1 P − β2(1−η)yv2 P −ayy−δ1yI, ż = β2(1−η)xv2 P − β1(1−η)zv1 P −azz −δ2zI, ẇ= β2(1−η)yv2 P + β1(1−η)zv1 P −aww−δ3wI, v̇1 =k1(1−ε)y+k3(1−ε)w−d1v1, v̇2 =k2(1−ε)w−d2v2, İ =ρ1yI+ρ2zI+ρ3wI+pI−qI2, (19) where P = x + y + z + w. (20) Susceptible host (healthy hepatocytes) cells are produced at a rate λ, died at a rate dx. β1 and β2 are respectively the contact rates between unin- fected cells with free HBV virions, and uninfected cells with free HDV virions. ay, az and aw are respectively death rate of HBV only infected cells, HDV only infected cells and coinfected cells by HDV and HBV. In these equations, all the parameters are positive and we assume that the death rate of uninfected cells is not greater than the death rate of infected cells, that is, min{dx,ay,az,aw} = dx. HBV virions v1 are Biomath 7 (2018), 1807219, http://dx.doi.org/10.11145/j.biomath.2018.07.219 Page 7 of 17 http://dx.doi.org/10.11145/j.biomath.2018.07.219 P. T. Mouofo, J. J. Tewa, S. Bowong, Modelling and analysis of a within-host model of hepatitis ... TABLE I NUMERICAL VALUES FOR THE PARAMETERS OF MODEL SYSTEM (19) Symbols Definition value Source λ Production rate of hepathocyte 252666.6667 day−1 [10] dx Normal death rate of healthy liver cells 0.0039 days−1 [15] ay Infected cell death rate (HBV) 0.0693 − 0.00693 day−1 [9] az Infected cell death rate (HDV) 0.009 day−1 Assumed aw Infected cell death rate (HBV-HDV) 0.0091 day−1 Assumed d1 death rate of free HBV 0.693 day−1 [14] d2 death rate of free HDV 0.6 day−1 Assumed β1 Infection rate of HBV 3.6 × 10−5 − 1.8 × 10−3 [14] β2 Infection rate of HDV 0.002 day−1 Assumed k1 HBV production per infected liver cells 200 − 1000 days−1 [8] k2 HDV production per infected liver cells 300 days−1 Assumed k3 HBV production per co-infected HBV& HDV cells 300 days−1 Assumed δ1 rate of CTL elimination 0.02 day−1 Assumed in infected HBV cells δ2 rate of CTL elimination in infected HDV cells 0.02 day−1 Assumed δ3 rate of CTL elimination in infected HBV-HDV cells 0.02 day−1 Assumed p the proliferation rate of immune cells 0.5 day−1 Assumed q density-dependent rate of immune 0.03 day−1 Assumed cells suppression ρi HBV-specific CTL stimulation rate 0.02 day−1 Assumed produced by HBV only infected cells y and by coinfected cells w. In fact, in order for HDV to successfully complete its replication cycle, a hep- atocyte must be coinfected with HBV and HDV. In these coinfected cells, the replication of HBV is suppressed by HDV, although not completely abolished [19], [20]. So, v1 can be produced by w. HBV only infected cells y produce HBV virions particles v1 at a rate k1(1 − ε)y and are killed by the CTL immune response at a rate δ1yI. For the same reason, HDV are killed by the CTL immune response at a rate δ1yI. Finally, CTL cells increase at a rates ρ1yI, ρ2zI and ρ3wI as a result of stimulation by the viral antigen of the infected cells. p and q have the same meaning as in the previous model. The parameter values used for numerical sim- ulation are given in Table 1. B. Mathematical analysis of the model 1) Positivity and boundedness of trajectories: We have the following result. Theorem 5. The nonegative orthant R7+ is posi- tively invariant for model system (19). Moreover, system (19) is pointwise dissipative and the ab- Biomath 7 (2018), 1807219, http://dx.doi.org/10.11145/j.biomath.2018.07.219 Page 8 of 17 http://dx.doi.org/10.11145/j.biomath.2018.07.219 P. T. Mouofo, J. J. Tewa, S. Bowong, Modelling and analysis of a within-host model of hepatitis ... sorbing set is given by Σ = { (x,y,z,w,v1,v2,I) ∈ R7+ : T ≤ λ dx , v1 ≤ (k1 +k3)(1 −ε)λ dxd1 , v2 ≤ λk2(1−ε) dxd2 , p q ≤ I ≤ ρ0 q } , where T = x+y+z+w and ρ0 = (ρ1+ρ2+ρ3)λ dx +p. The proof is similar to the proof of Theorem 1. 2) Basic reproduction number and equilibria: The disease-free equilibrium of model system (19) is given by E0 = (x∗, 0, 0, 0, 0, 0,I∗). Using the notations in van den Driessche and Watmough[17] for model system (19), the matri- ces F and V for the new infection terms and the remaining transfer terms are, respectively, given by F =   0 0 0 (1−η)β1 0 0 0 0 0 (1 −η)β2 0 0 0 0 0 k1(1−ε) 0 k3(1−ε) 0 0 0 0 k2(1−ε) 0 0   and V =   ay +δ1I ∗ 0 0 0 0 0 az +δ2I ∗ 0 0 0 0 0 aw +δ3I ∗ 0 0 0 0 0 d1 0 0 0 0 0 d2   . It then follows that the basic reproduction number is given by R0 = k1(1 −ε)β1(1 −η) d1(ay + δ1I∗) . (21) Thus, using Theorem 2 of van den Driessche and Watmough[17], we have the following result. Lemma 3. : The virus free equilibrium E0 of model system (19) is locally asymptotically stable (LAS) if R0 < 1, and unstable if R0 > 1. Remark 3. Note that the basic reproduction num- ber of the co-infection model of HBV and HDV is the same than the basic the basic reproduction of the HBV model alone. This suggests that to control the co-infection of HBV and HDV within the body of a host, one only needs to control the HBV infection. We now process with the existence of steady states. The steady states of model system (19) satisfy the following equations:  λ−dxx̄− β1(1−η)x̄v̄1 P̄ − β2(1−η)x̄v̄2 P̄ = 0, β1(1−η)x̄v̄1 P̄ − β2(1−η)ȳv̄2 P̄ −ayȳ−δ1ȳĪ = 0, β2(1−η)x̄v̄2 P̄ − β1(1−η)z̄v̄1 P −azz̄−δ2z̄Ī = 0, β2(1−η)ȳv̄2 P̄ + β1(1−η)z̄v̄1 P̄ −aww̄−δ3w̄Ī = 0, k1(1 −ε)ȳ + k3(1 −ε)w̄ −d1v̄1 = 0, k2(1 −ε)w̄ −d2v̄2, ρ1ȳĪ + ρ2z̄Ī + ρ3w̄Ī + pĪ −qĪ2 = 0, (22) where P̄ = x̄ + ȳ + z̄ + w̄. From the last equation of (22), we have (ρ1ȳ + ρ2z̄ + ρ3w̄ + p−qĪ)Ī = 0, (23) which has two possible solutions: Ī = 0 and ρ1ȳ+ ρ2z̄ + ρ3w̄ + p−qĪ = 0. Since Ī ≥ I∗ and Ī 6= 0, one has that Ī = 1 q (ρ1ȳ+ρ2z̄+ρ3w̄+p). From the sixth equation of (22), one has w̄ = d2 k2(1 −ε) v̄2. Adding the third and fourth equations of (23) and using the expression of w̄ yields z̄ = [ β2(1−η)(x̄+ȳ) (az +δ2Ī)P̄ − (aw + δ3Ī)d2 k2(1−ε)(az +δ2Ī) ] v̄2. (24) Substituting Eq.(24) into the fourth equation of (22) gives{ β2(1−η)ȳ P̄ + β1(1−η)v̄1 P̄ [ β2(1−η)(x̄+ȳ) (az +δ2Ī)P̄ − (aw +δ3Ī)d2 k2(1−ε)(az +δ2Ī) ]} v̄2 = 0. (25) Eq. (25) has two possible solutions. If v̄2 = 0, then from Eq. (22) one has that v̄1 = 0 or z̄ = 0. Biomath 7 (2018), 1807219, http://dx.doi.org/10.11145/j.biomath.2018.07.219 Page 9 of 17 http://dx.doi.org/10.11145/j.biomath.2018.07.219 P. T. Mouofo, J. J. Tewa, S. Bowong, Modelling and analysis of a within-host model of hepatitis ... Note that v̄1 = 0 leads to the virus free equi- librium E0 = (x∗, 0, 0, 0, 0, 0,I∗). If z̄ = 0, we obtain the HBV-only persistence equilibrium Ē = (x̄, ȳ, 0, 0, v̄1, 0, Ī) which has biological signifi- cance if R0 > 1. If v̄2 6= 0, the model admits a co- infection equilibrium Ê = (x̂, ŷ, ẑ, ŵ, v̂1, v̂2, Î) where the existence is verified numerically. 3) Stability of equilibria: We have the follow- ing result. Theorem 6. If R0 < 1, then the infection free equilibrium E0 of model system (19) is globally asymptotically stable in Σ whenever k3(1−ε)(ay +δ1I∗)+ k1(1−ε)k2(1−ε)β2(1−η) d2 −k1(1−ε)(aw +δ3I∗) ≥ 0. (26) Proof: Consider the following LaSalle func- tion candidate: L(t) = ay + az + aw + bv1 + cv2, (27) where a, b and c are positive constants to be determined later. Using the fact that x ≤ x∗ and I ≥ I∗, the derivative of L along the solution of (19) satisfies L̇ = aβ1(1−η)xv1 P + aβ2(1−η)xv2 P −a(ay+δ1I)y −a(az + δ2I)z −a(aw + δ3I)w+bk1(1−ε)y + bk3(1−ε)w−bd1v1 + ck2(1−ε)w−cd2v2 ≤ [ aβ1(1 −η) − bd1 ] v1 + [ bk1(1 −ε) −a(ay + δ1I∗) ] y + [ aβ2(1 −η) − cd2 ] v2 + [ bk3(1−ε)+ck2(1−ε)−a(aw +δ3I∗) ] w. We choose a = k1(1 − ε), b = ay + δ1I∗ and c = k1(1 −ε)β2(1 −η) d2 so that the coefficients of y, v2 and w are equal to zero. In this case, if the condition (26) holds, we obtain L̇ ≤ (ay + δ1I∗) [ R0 − 1 ] v1. Thus, if R0 ≤ 1 then L̇ ≤ 0 ∀x, y, z, w, v1, v2, I ≥ 0 and L̇ = 0 if only if (x, y, z, w, v1, v2) = (x∗, 0, 0, 0, 0, 0,I∗). Then the globally asymptotically attractivity of E0 follows from Lyapunov LaSalle Invariance Principle [18]. This completes the proof. Now, we study the stability of the HBV-only persistence equilibrium Ē = (x̄, ȳ, 0, 0, v̄1, 0, Ī). The Jacobian matrix of model system (19) at Ē is JĒ=   −dx−S2 S1 S1 S2 −S1−ay−δ1Ī −S1 0 0 −S1−az−δ2Ī 0 0 S1 0 k1(1−ε) 0 0 0 0 0 ρ1Ī ρ2Ī S1 −S3 −S4 0 −S1 S3 −S5 −δ1ȳ 0 0 S4 0 −aw−δ3Ī 0 S5 0 k3(1−ε) −d1 0 0 k2(1−ε) 0 −d2 0 ρ3Ī 0 0 −qĪ   , where S1 = β1(1 −η)x̄v̄1 (x̄ + ȳ)2 , S2 = β1(1 −η)ȳv̄1 (x̄ + ȳ)2 , S3 = β1(1 −η)x̄ x̄ + ȳ , S4 = β2(1 −η)x̄ x̄ + ȳ and S5 = β2(1 −η)ȳ x̄ + ȳ . The local stability of Ē is governed by the eigenvalues of the JĒ. Hence, conditions for local stability of Ē have been derived by applying the Routh-Hurwitz criterion to the char- acteristic equation of JE. The expresions are com- plicated and are not presented here, but available from the authors on request. Importantly, the set of parameters satisfying these conditions is not empty. Figure 1 presents the time evolution of model system (19) when β1 = 0.02 β2 = 0.07, ε = η = 0, k1 = 600 and ρ1 = 0.00002 (so that R0 > 1). All other parameter values are as in Table 1. In this case, the co-infection equilibrium is Ê(6.5× 104, 1.8 × 104, 490.93, 310.3, 1.572 × 107, 1.55 × 105, 670.9). Initial conditions have been chosen to be x(0) = 2 × 107, y(0) = 104, z(0) = 4 × 104, w(0) = 3 × 103, v1(0) = 2 × 103, v2(0) = 1.5 × Biomath 7 (2018), 1807219, http://dx.doi.org/10.11145/j.biomath.2018.07.219 Page 10 of 17 http://dx.doi.org/10.11145/j.biomath.2018.07.219 P. T. Mouofo, J. J. Tewa, S. Bowong, Modelling and analysis of a within-host model of hepatitis ... 0 10 20 30 40 50 60 70 80 90 100 0 0.5 1 1.5 2 2.5 x 10 7 Time(Days) U n in fe ct e d c e lls (a1) 0 10 20 30 40 50 60 70 80 90 100 0 2 4 6 8 10 12 14 x 10 4 Time(Days) In fe ct e d c e lls b y H B V (b1) 0 10 20 30 40 50 60 70 80 90 100 0 500 1000 1500 2000 2500 3000 3500 4000 Time(Days) In fe ct e d c e lls b y H D V (c1) 0 10 20 30 40 50 60 70 80 90 100 0 500 1000 1500 2000 2500 3000 Time(Days) C o − in fe ct e d c e lls (d1) 0 10 20 30 40 50 60 70 80 90 100 0 2 4 6 8 10 12 x 10 7 Time(Days) F re e H B V v ir u s (e1) 0 10 20 30 40 50 60 70 80 90 100 0 0.5 1 1.5 2 2.5 3 x 10 5 Time(Days) F re e H D V v ir u s (f1) 0 10 20 30 40 50 60 70 80 90 100 0 500 1000 1500 2000 2500 Time(Days) C T L c e lls (g1) Fig. 1. Time evolution of model system (19) when β1 = 0.02 β2 = 0.07, ε = η = 0, k1 = 600 and ρ1 = 0.00002 (so that R0 > 1). All other parameter values are as in Table 1. 103 and I(0) = 4 × 102. From this figure, it is evident that the trajectories of model system (19) converge to the co-infection equilibrium. IV. OPTIMAL CONTROL OF TREATMENT IN THE HBV-HDV CO-INFECTION MODEL This section deals with the problem of optimal control of the co-infection of HBV and HDV. More precisely, we are concerned with the prob- lem of adopting the best strategy of treatment to fight the HBV-HDV co-infection. We seek to search a maximum count of healthy cells with a minimum dose of the administered drugs. Hence, if we denote η(t) the first control variable which is efficacy of inhibiting new virus infections as a consequence of virus clearance and ε(t) the sec- ond control variable which represents the efficacy of inhibiting viral production from infected cells, model system (19) can be written, to accommo- date control actions or chemotherapy treatment, as follows:  ẋ=λ−dxx− β1(1−η(t))xv1 P − β2(1−η(t))xv2 P , ẏ = β1(1−η(t))xv1 P − β2(1−η(t))yv2 P −ayy −δ1yI, ż = β2(1−η(t))xv2 P − β1(1−η(t))zv1 P −azz −δ2zI, ẇ= β2(1−η(t))yv2 P + β1(1−η(t))zv1 P −aww −δ3wI, v̇1 =k1(1−ε(t))y+k3(1−ε(t))w−d1v1, v̇2 =k2(1−ε(t))w−d2v2, İ = ρ1yI + ρ2zI + ρ3wI + pI −qI2, (28) Biomath 7 (2018), 1807219, http://dx.doi.org/10.11145/j.biomath.2018.07.219 Page 11 of 17 http://dx.doi.org/10.11145/j.biomath.2018.07.219 P. T. Mouofo, J. J. Tewa, S. Bowong, Modelling and analysis of a within-host model of hepatitis ... with x0, y0, z0, w0, v01 , v 0 and I0 the given initial values of x, y, w, z, v1, v2, I at t0 = 0 respectively. The control functions, η(t) and ε(t), are bounded, Lebesgue integrable functions. Our objective functional to be maximized is J(η,ε) = ∫ tf 0 [ x(t) − ( A 2 η2 + B 2 ε2 )] dt. (29) In other words, we want to maximize the benefit based on the healthy liver cells count and min- imizing the cost based on the percentage effect chemotherapy given (i.e. η and ε). The parameters A and B are positive and represent the weights on the benefit and cost. The goal is to seek an optimal control pair (η∗,ε∗) such that J(η∗,ε∗) = max{J(η,ε) : (η,ε) ∈U}, (30) where U is the control set defined by: U = {u = (η,ε), η,εmeasurable, 0 ≤ η(t) ≤ 1, 0 ≤ ε(t) ≤ 1,∀t ∈ [0, tf ]}. The basic framework of this problem is to prove the existence and the uniqueness of the optimal control and to characterize it. A. Analysis of Optimal Controls The necessary conditions that an optimal pair must satisfy come from Pontryagin’s Maximum Principle [21]. This principle converts (28) - (30) into a problem of minimizing pointwise a Hamil- tonian, H, with respect to η and ε: H = x(t) − ( A 2 η2 + B 2 ε2 ) +λ1 [ λ−dxx− β1(1−η(t))xv1 P − β2(1−η(t))xv2 P ] +λ2 [ β1(1−η(t))xv1 P − β2(1−η(t))yv2 P −ayy −δ1yI ] +λ3 [ β2(1−η(t))xv2 P − β1(1−η(t))zv1 P −azz−δ2zI ] + λ4 [ β2(1−η(t))yv2 P + β1(1−η(t))zv1 P −aww−δ3wI ] +λ5 [k1(1 −ε(t))y + k3(1 −ε(t))w −d1v1] +λ6 [k2(1 −ε(t))w −d2v2] +λ7 [ ρ1yI + ρ2zI + ρ3wI + pI −qI2 ] . (31) By applying Pontryagins Maximum Principle [21] and the existence result for the optimal control pairs from [22], we have the following result. Theorem 7. There exists an optimal control pair η∗, ε∗ and corresponding solution x∗, y∗, z∗, w∗, v∗1 , v ∗ 2 and I ∗, such that J(η∗,ε∗) = max U J(η,ε). Furthermore, there exist adjoint functions λ1, λ2, λ3, λ4, λ5, λ6, and λ7 such that equation (32) holds with the transversality conditions λi(tf ) = 0, for i = 1, ..., 7 (33) and P = x∗ + y∗ + z∗ + w∗. The following characterization holds: η∗ = min{max{0, R1(t)}, 1}, ε∗ = min{max{0, R2(t)}, 1} (34) where R1(t) = 1 A [ λ1 β1x ∗v∗1 +β2x ∗v∗2 P +λ2 β2y ∗v∗2 −β1x ∗v∗1 P +λ3 β1z ∗v∗1−β2x ∗v∗2 P −λ4 β1z ∗v∗1 + β2y ∗v∗2 P ] R2(t) = −λ5(k1y∗ + k3w∗) −λ6k2w∗ B Proof: Corollary 4.1 in Fleming and Rishel [22] gives the existence of an optimal control pair due to the concavity of the integrand of J with respect to (η,ε), a priori boundedness of the state solutions, and the Lipschitz property of the state system with respect to the state variables. Applying Pontryagin’s Maximum Principle, we obtain Biomath 7 (2018), 1807219, http://dx.doi.org/10.11145/j.biomath.2018.07.219 Page 12 of 17 http://dx.doi.org/10.11145/j.biomath.2018.07.219 P. T. Mouofo, J. J. Tewa, S. Bowong, Modelling and analysis of a within-host model of hepatitis ... λ̇1 = −1 + λ1 ( dx + β1(1 −η)v∗1 (y ∗ + z∗ + w∗) P 2 + β2(1 −η)v∗2 (y ∗ + z∗ + w∗) P 2 ) , −λ2 ( β1(1 −η)v∗1 (y ∗ + z∗ + w∗) P 2 + β2(1 −η)y∗v∗2 P 2 ) −λ3 ( β2(1 −η)v∗2 (y ∗ + z∗ + w∗) P 2 + β1(1 −η)z∗v∗1 P 2 ) −λ4 ( β2(1 −η)y∗v∗2 P 2 + β1(1 −η)z∗v∗1 P 2 ) λ̇2 = −λ1 ( β1(1 −η)x∗v∗1 P 2 + β2(1 −η)x∗v∗2 P 2 ) + λ3 ( β2(1 −η)x∗v∗2 P 2 − β1(1 −η)z∗v∗1 P 2 ) +λ2 ( β1(1 −η)x∗v∗1 P 2 + β2(1 −η)v∗2 (x ∗ + z∗ + w∗) P 2 + ay + δ1I ∗ ) +λ4 ( β1(1 −η)z∗v∗1 P 2 − β2(1 −η)v∗2 (x ∗ + z∗ + w∗) P 2 ) −k1(1 −ε)λ5 −ρ1I∗λ7, λ̇3 = −λ1 ( β1(1 −η)x∗v∗1 P 2 + β2(1 −η)x∗v∗2 P 2 ) + λ2 ( β1(1 −η)x∗v∗1 P 2 − β2(1 −η)y∗v∗2 P 2 ) +λ3 ( β2(1 −η)x∗v∗2 P 2 + β1(1 −η)v∗1 (x ∗ + y∗ + w∗) P 2 + az + δ2I ∗ ) +λ4 ( β2(1 −η)y∗v∗2 P 2 − β1(1 −η)v∗1 (x ∗ + y∗ + w∗) P 2 ) −ρ2I∗λ7, λ̇4 = −λ1 ( β1(1 −η)x∗v∗1 P 2 + β2(1 −η)x∗v∗2 P 2 ) + λ2 ( β1(1 −η)x∗v∗1 P 2 − β2(1 −η)y∗v∗2 P 2 ) +λ3 ( β2(1 −η)x∗v∗2 P 2 − β1(1 −η)z∗v∗1 P 2 ) −k3(1 −ε)λ5 −k2(1 −ε)λ6 −ρ3I∗λ7 +λ4 ( β2(1 −η)y∗v∗2 P 2 + β1(1 −η)z∗v∗1 P 2 + aw + δ3I ∗ ) , (32) λ̇5 = β1(1 −η)x∗ P (λ1 −λ2) + β1(1 −η)z∗ P (λ3 −λ4) + d1λ5, λ̇6 = β2(1 −η)x∗ P (λ1 −λ3) + β2(1 −η)y∗ P (λ2 −λ4) + d2λ6, λ̇7 = δ1y ∗λ2 + δ2z ∗λ3 + δ3w ∗λ4 − (ρ1y∗ + ρ2z∗ + ρ3w∗ + p− 2qI∗)λ7, Biomath 7 (2018), 1807219, http://dx.doi.org/10.11145/j.biomath.2018.07.219 Page 13 of 17 http://dx.doi.org/10.11145/j.biomath.2018.07.219 P. T. Mouofo, J. J. Tewa, S. Bowong, Modelling and analysis of a within-host model of hepatitis ... λ̇1 = − ∂H ∂x , λ̇2 = − ∂H ∂y , ..., λ̇7 = − ∂H ∂I , evaluated at the optimal control pair and corre- sponding states, which results in the stated adjoint system (32) and (33) [35]. By considering the optimality conditions, ∂H ∂η = 0 and ∂H ∂ε = 0 and solving for η∗, ε∗, subject to the constraints, the characterizations (34) can be derived. To il- lustrate the characterization of ε∗, we have ∂H ∂ε = −Bε−λ5(k1y + k3w) −λ6k2w = 0 at ε∗ on the set {t | 0 < ε∗(t) < 1}. On this set, we have ε∗(t) = −λ5(k1y + k3w) −λ6k2w B . Taking into account the bounds on ε, we obtain the character- ization of ε in (34). Next, we discuss the numerical solutions of the optimality system and the corresponding optimal control pairs, the parameter choices, and the in- terpretations from various cases. B. Numerical results In this section, we study numerically an optimal treatment strategy of our HBV-HDV co-infection model. The optimal treatment strategy is obtained by solving the optimality system, consisting of 14 ODEs from the state and adjoint equations. An iterative method is used for solving the optimality system. We start to solve the state equations with a guess for the controls over the simulated time using a forward fourth order Runge-Kutta scheme. Because of the transversality conditions (33), the adjoint equations are solved by a back- ward fourth order Runge-Kutta scheme using the current iteration solution of the state equations. Then, the controls are updated by using a convex combination of the previous controls and the value from the characterizations (34). This process is repeated and iteration is stopped if the values of unknowns at the previous iteration are very close to the ones at the present iteration. In order to evaluate our control strategy, we consider six co-infected patients by HBV and HDV: A1, A2, B1, B2, C1, and C2 divided in three groups such that The first group consists of A1 and A2 and the viral HDV load of patients A1 and A2 are the same and are 1.16×103 copies/ml. The second group are patients B1 and B2 with the same HDV viral loads of 1.16 × 105 copies/ml; the end group are patients C1, and C2 with high levels of HDV viremia which are 1.16 × 107 copies/ml. We suppose that patients A1, B1 and C1 are not under treatment. Otherwise, patients A2, B2 and C2 are under our treatment control strategy. We want to evaluate in 100 days the evolution of viral HDV load in every patient. We suppose that the viral load of HBV in those groups are 108 copies/ml. Note that some studies report improvements in patients, with interferon-α (IFN) efficacy as high as 90% (η = 0.9) [36]. Although lamivudine (LMV) does not have a direct effect on HDV viral production k2, its effect on the viral production of HBV will also have an effect on the HDV viral dynamics. Studies such as Lewin et al [37] show efficacy levels of therapies based on LMV varying between 90% to 99% (ε > 0.9) for the HBV infection. With this in mind, we consider combination of interferon-α and lamivudine in our simulations. In Fig.(2), parameter values are the same as in Fig.(1). Figure (2) (a1) − (c1) presents the time evolution of viral load of patients of group 1. In this figure, we observe that when there is not strategy of treatment, the disease persists in the host. But, with our strategy of treatment, the disease disappear in the host. To have this result, the control ε(t) is at it maximal value almost all the period of our control strategy. This means that, to control HDV, we need medication with hight efficacy. We have the same observation on the other figures for group 2 and group 3. Fig.(3) illustrate the importance of immune system in our model. We observe that without CTL cells, the viral load of patients A1, B1 and C1 increase. Otherwise, the viral load of patients A2, B2 and C2 converge to zero if only if ε(t) and η(t) are to their maximal values during all the period of control strategy. So, it is difficult Biomath 7 (2018), 1807219, http://dx.doi.org/10.11145/j.biomath.2018.07.219 Page 14 of 17 http://dx.doi.org/10.11145/j.biomath.2018.07.219 P. T. Mouofo, J. J. Tewa, S. Bowong, Modelling and analysis of a within-host model of hepatitis ... 0 10 20 30 40 50 60 70 80 90 100 0 2 4 6 8 10 12 x 10 4 Time (Days) F re e H D V v ir u s HDV viral load of patient A 2 with control treatment HDV viral load for patient A 1 without treatment (a1) 0 10 20 30 40 50 60 70 80 90 100 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Time (Days) T re a tm e n t co n tr o l ε (t ) fo r p a tie n t A 2 (b1) 0 10 20 30 40 50 60 70 80 90 100 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 Time (Days) T re a tm e n t co n tr o l η (t ) o f p a tie n t A 2 (c1) 0 10 20 30 40 50 60 70 80 90 100 0 2 4 6 8 10 12 x 10 4 Time (Days) F re e H D V v ir u s HDV viral load of patient B 2 with treatment strategy HDV viral load of patient B 1 without treatment (a2) 0 10 20 30 40 50 60 70 80 90 100 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Time (Days) T re a tm e n t co n tr o l ε( t) f o r p a tie n t B 2 (b2) 0 10 20 30 40 50 60 70 80 90 100 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 Time (Days) T re a tm e n t co n tr o l η (t ) fo r p a tie n t B 2 (c2) 0 10 20 30 40 50 60 70 80 90 100 0 2 4 6 8 10 12 x 10 4 Time (Days) F re e H D V v ir u s HDV viral load of patient C 2 with control strategy HDV viral load of patient C 1 without treatment (a3) 0 10 20 30 40 50 60 70 80 90 100 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Time (Days) T re a tm e n t co n tr o l ε (t ) fo r p a tie n t C 2 (b3) 0 10 20 30 40 50 60 70 80 90 100 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 Time (Days) T re a tm e n t co n tr o l η (t ) fo r p a tie n t C 2 (c3) Fig. 2. Time evolution of viral load of patients A1, B1, C1 without treatment, and patients A2, B2 and c2 with our treatment strategy. to control the disease in this case. This show the importance of immune system in our model. V. DISCUSSION AND CONCLUSION Hepatitis D virus infection has a worldwide distribution. Studying the transmission, epidemi- ology and dynamic virus of HDV infection is an important topic. It is a unique virus for which many open questions remain. For example, HDV- specific treatment protocols still do not exist. The investigation of the inter-related dynamics of chronic HBV and HDV infections are important to understanding how treatment may affect this complex system. To this end, the development of biologically realistic mathematical models is an important tool. The main objective of this paper was to shed light on the interaction between HBV and HDV. A realistic deterministic ODE based compart- mental model for the transmission of HBV and HDV co-dynamics within the body of a host has been proposed and analyzed. The HBV-only model was qualitatively examined, first of all. The mathematical analysis results show that the basic reproductive number R0 of the co- interaction HBV-HDV model is the same than the basic reproduction number of the HBV model alone. This suggests that the eradication of HDV is conditioned by the eradication of HBV. The epidemiological implication of this is that for Biomath 7 (2018), 1807219, http://dx.doi.org/10.11145/j.biomath.2018.07.219 Page 15 of 17 http://dx.doi.org/10.11145/j.biomath.2018.07.219 P. T. Mouofo, J. J. Tewa, S. Bowong, Modelling and analysis of a within-host model of hepatitis ... 0 20 40 60 80 100 0 1 2 3 4 5 6 7 8 x 10 9 Time (Days) F re e H D V v ir u s HDV viral load of patient C 2 without immune system reaction HDV viral load of patient C 1 without immune system reaction (d1) 0 20 40 60 80 100 0 1 2 3 4 5 6 7 8 x 10 9 Time (Days) F re e H D V v ir u s HDV viral load of patient B 2 without immune system reaction HDV viral load of patient B 1 without immune system reaction (d2) 0 20 40 60 80 100 0 1 2 3 4 5 6 7 8 x 10 9 Time (Days) F re e H D V v ir u s HDV viral load of patient A 2 without immune system reaction HDV viral load of patient A 1 without immune system reaction (d3) Fig. 3. Time evolution of viral load different groups without CTL cells . the effective eradication and control of HDV, R0 should be less than one. Moreover, achieving this may be too costly, because it means that for constant controls, one needs to keep treating for infinite time. Therefore, we considered time dependent controls as a way out, to ensure the eradication of the disease in a finite time. In this light, we addressed the optimal control by deriving and analyzing the conditions for optimal eradication of the disease. From our numerical results, we can conclude that immune responses play a significant role in eradication of disease. Moreover, to eradicate the disease, it is important to manufacture a drug treatment with hight effi- cacy which can block the viral production. REFERENCES [1] Hoofnagle J, Mullen K, Peters M (1987);“ Treatment of chronic delta hepatitis with recombinant alpha inter- feron”. In: Rizzetto M, Gerin JL, Purcell RH, editors. The hepatitis delta virus and its infection. New York: Alan R. Liss. [2] Lau DT, Kleiner DE, Park Y, Bisceglie AMD, Hoofnagle JH (1999); “Resolution of chronic delta hepatitis after 12 years of interferon alfa therapy”. Gastroenterology 117(5): 12291233. [3] Lau DT, Doo E, Park Y, Kleiner DE, Schmid P, et al. (1999); “Lamivudine for chronic delta hepatitis. Hepa- tology”. 30(2): 546549. [4] Yurdaydin C, Bozkaya H, Gurel S, Tillman HL, Aslan N, et al. (2002); “Famciclovir treatment of chronic delta hepatitis”. Journal Hepatology 37(2): 266271 [5] R. Qesmi, J. Wu, J. M. Heffernan, (2010); “Influence of backward bifurcation in a model of hepatitis B and C viruses”. Math. Biosci. 224, 118125. [6] K. Wang, W. Wang,(2007); “Propagation of HBV with spatial dependence”. Math. Biosci. 210, 7895 [7] F. Brauer, C. Castillo-Chvez,(2001); “Mathematical Models in Population Biology and Epidemiology, Texts in Applied Mathematics”. Springer, New York, NY, . [8] Murray J. M., Purcell R. H., Wieland S. F. (2006); “The half-life of hepatitis B virions”. Hepatology 44, 1117- 1121. [9] M. A. Nowak, R. M. May,, (2000); “Virus Dynamics”. Oxford University Press, Oxford. [10] L. Min, Y. Su, Y. Kuang, (2008); “Mathematical Anal- ysis of a Basic Virus Infection Model With Application to HBV Infection”. Rocky Mountain Journal of Mathe- matics, 38, 113. [11] K. Wang, W. Wang, S. Song,(2008); “Dynamics of an HBV model with diffusion and delay”, J. Theor. Bio. 253 3644. [12] S. Eikenberry, S. Hews, J. Nagy, Y. Kuang, (2009); “The dynamics of a delay model of hepatitis B virus infection with logistic hepatocyte growth”, Math. Biosci. Eng. 6 283299, . [13] Xiridou M, Borkent-Raven B, Hulshof J, Wallinga J (2009);“ How hepatitis D virus can hinder the con- trol of hepatitis B virus”. PLoS ONE 4(4): e5247. doi:10.1371/journal.pone.0005247. [14] Whalley S. A., Murray J. M., Brown D., Webster G. J. M., Emery V. C., Dusheiko G. M., Perelson A. S. (2001); “Kinetics of Acute Hepatitis B Virus Infection in Humans”, J. Exp. Med. 193, 847-853. [15] Bralet M., Branchereau S., Brechot C. and Ferry N. (1994);“ Cell Lineage Study in the Liver Using Retroviral Mediated Gene Transfer”; Am. J. Pathol. 144, 896-905. [16] S. A. Gourley, Y. Kuang, and J. D. Nagy, “Dynamics of a delay differential equation model of hepatitis B virus infection”, Journal of Biological Dynamics, vol. 2, no. 2, pp. 140153, 2008. [17] P. van den Driessche, J. Watmough,(2002); “Reproduc- tion numbers and sub-threshold endemic equilibria for compartmental models of disease transmission”, Math. Biosci. 180 29-48. Biomath 7 (2018), 1807219, http://dx.doi.org/10.11145/j.biomath.2018.07.219 Page 16 of 17 http://dx.doi.org/10.11145/j.biomath.2018.07.219 P. T. Mouofo, J. J. Tewa, S. Bowong, Modelling and analysis of a within-host model of hepatitis ... [18] Bhatia NP, Szeg GP.1970; “Stability theory of dynami- cal systems”. Springer-Verlag. [19] M. Schaper, F. Rodriguez-Frias, R. Jardi, D. Tabernero, M. Homs, G. Ruiz, J. Quer, R. Esteban, M. Buti,(2010); Quantitative longitudinal evaluations of hepatitis delta virus rna and hepatitis B virus DNA shows a dynamic, complex replicative profile in chronic hepatitis B and D, Journal of Hepatology 52 658. [20] R. Jardi, F. Rodriguez, M. Buti, X. Costa, M. Cotrina, R. Galimany, R. Esteban, J. Guardia, (2001); Role of hepatitis B, C, and D viruses in dual and triple infection: Influence of viral genotypes and hepatitis B precore and basal core promoter mutations on viral replicative interference, Hepatology 34 404. [21] L. S. Pontryagin, V. G. Boltyanskii, R. V. Gamkrelidze, and E. F. Mishchenko, (1962) “The mathematical theory of optimal processes”, Wiley, New York. [22] W. H. Fleming and R. W. Rishel, (1975); “Deterministic and Stochastic Optimal Control”, Springer Verlag, New York. [23] J. A. Deans, S. Cohen, (1983); ”Immunology of malaria. Ann. Rev. Microbiol., 37 , 25-50. [24] A. M. Elaiw,(2015) “Global stability analysis of hu- moral immunity virus dynamics model including latently infected cells”, J. Biol. Dyn., 9 , 215-228. [25] A. M. Elaiw,(2015); “Global threshold dynamics in humoral immunity viral infection models including an eclipse stage of infected cells”, J. Korean Soc. Ind. Appl. Math., 19 , 137-170. [26] A. M. Elaiw, N. H. AlShamrani,(2015); “Global sta- bility of humoral immunity virus dynamics models with nonlinear infection rate and removal”, Nonlinear Anal. Real World Appl., 26 , 161-190. [27] A. Murase, T. Sasaki, T. Kajiwara,(2005); “Stability analysis of pathogen-immune interaction dynamics”, J. Math. Biol., 51 , 247-267. [28] M. A. Obaid, A. M. Elaiw,(2014); “Stability of virus infection models with antibodies and chronically infected infection models with antibodies and chronically infected cells”, Abstr. Appl. Anal., 12 pages. [29] K. HATTAF, N. YOUSFI AND A. TRIDANE; (2012)“ Global stability analysis of a generalized virus dynamics model with the immune response”, Canadian Applied Mathematics Quarterly 20 (4) 499518. [30] T. Wang, Z. Hu, F. Liao,(2014) “Stability and Hopf bifurcation for a virus infection model with delayed humoral immunity response”, J. Math. Anal. Appl., 411 63-74. [31] T. Wang, Z. Hu, F. Liao, W. Ma, “Global stability analysis for delayed virus infection model with general incidence rate and humoral immunity”, Math. Comput. Simulation, 89 (2013), 13-22. [32] S. Wang, D. Zou,(2012); “Global stability of in-host viral models with humoral immunity and intracellular delays”, J. Appl. Math. Model, 36 , 1313-1322. [33] Mostafa Khabouze, Khalid Hattaf, and Noura Yousfi1;(2014) “Stability Analysis of an Improved HBV Model with CTL Immune Response, ”,International Scholarly Research Notices. [34] Hattaf, K.; Yousfi, N. (2015); “A generalized HBV model with diffusion and two delays”. Comput. Math. Appl., 69, 3140. [35] M. I. Kamien and N. L. Schwarz,(1991); “Dynamic optimization: the calculus of variations and optimal control”, North Holland, Amsterdam. [36] Chien RN, (2008); “Current therapy for hepatitis C or D or immunodeficiency virus concurrent infection with chronic hepatitis B. Hepatol Int 2: 296303. [37] Lewin SR, Ribeiro RM, Walters T, Lau GK (2001) “Analysis of hepatitis B viral load decline under potent therapy: complex decay profiles observed”. Hepatology 34: 10121020. [38] Noura Yousfi Khalid Hattaf Abdessamad Tridane, (2011); Modeling the adaptive immune response in HBV infection, J. Mathematical Biology 63 933957. Biomath 7 (2018), 1807219, http://dx.doi.org/10.11145/j.biomath.2018.07.219 Page 17 of 17 http://dx.doi.org/10.11145/j.biomath.2018.07.219 Introduction The HBV model The model description Analysis of the model Positivity and boundedness of solutions Basic reproduction number and equilibria Stability of equilibria The HBV-HDV co-infection model Model construction Mathematical analysis of the model Positivity and boundedness of trajectories Basic reproduction number and equilibria Stability of equilibria Optimal control of treatment in the HBV-HDV co-infection model Analysis of Optimal Controls Numerical results Discussion and conclusion References