Mathematical Analysis and Treatment for a Delayed Hepatitis B Viral Infection Model with the Adaptive Immune Response and DNA-Containing Capsids

We model the transmission of the hepatitis B virus (HBV) by six differential equations that represent the reactions between HBV with DNA-containing capsids, the hepatocytes, the antibodies and the cytotoxic T-lymphocyte (CTL) cells. The intracellular delay and treatment are integrated into the model. The existence of the optimal control pair is supported and the characterization of this pair is given by the Pontryagin’s minimum principle. Note that one of them describes the effectiveness of medical treatment in restraining viral production, while the second stands for the success of drug treatment in blocking new infections. Using the finite difference approximation, the optimality system is derived and solved numerically. Finally, the numerical simulations are illustrated in order to determine the role of optimal treatment in preventing viral replication.


Introduction
Hepatitis B virus (HBV) infects hepatocytes and causes approximately one million deaths annually [1]. With more than 257 million infected persons; HBV is considered a global public health problem [2]. This dangerous epidemic can be easily transmitted through contact with infected body fluids [3]. After transmission of the infection, HBV can cause acute or chronic illness [4]. Many mathematical models have been developed in order to study and model the dynamics of this serious viral infection [5][6][7][8]. All these models include the interaction between HBV and both the healthy and infected liver cells. The models which include the adaptive immune response in fighting the free viruses and in reducing the infected cells have been studied [9][10][11][12]. This adaptive immunity is represented by cytotoxic T-lymphocytes (CTL) and antibody immune responses. The mathematical analysis of HBV viral infection with HBV DNA-containing capsids was determined [13][14][15][16][17]. It has also been noted that the infected liver cells release the HBV DNA-containing capsids under the form of mature viruses after being enveloped by both cellular membrane lipids and viral envelope proteins [18,19]. More recently, the optimal control of HBV infection including HBV DNA-containing capsids and CTL immune response was studied [17]. In this paper, we are interested in the same problem, but we will introduce antibodies into this model. This work is motivated by the role of antibodies in reducing the viral infection severity [20][21][22]; thus it will be important to consider such an interesting element in studying the HBV viral dynamics. The model under consideration will be stated under the form of the following nonlinear system of six differential equations: The uninfected cells X are produced with an average s, die with a rate µ, and become infected by the virus with a rate k. Infected cells Y die with an average δ and are killed by the CTL immune system with an average p. The constant λ is the death average of infected but still not virus-producing cells. The intracellular delay, τ, stands for the time needed for infected cells to produce new viruses after viral entry. The term e −λτ is the probability of surviving between t − τ and t. The capsids D are produced with a rate a, they are transmitted to blood with a rate β and die with a rate δ. The free viruses V grow with a rate β, decay at a rate u and are neutralized by antibodies with a rate q. Antibodies W expand in response to free virus with a rate g and decay at a rate h. CTLs Z develop in response to viral antigen derived from infected cells with a rate c and decay in the absence of antigenic stimulation with a rate b. Finally, η 1 and η 2 denote the efficiency of pegylated interferon (PEG IFN) and lamivudine (LMV) drugs, respectively. It is noteworthy to mention that the main function of the PEG IFN drug is to block new infections of the healthy hepatocytes in the liver, while the prime role of the second drug, LMV, is to stop viral production [16,17].
The organization of this paper is as follows. The next section is concerned to the analysis of the model. Section 3 is devoted to an optimization analysis of our suggested viral infection model. In Section 4, we construct an appropriate numerical algorithm and show some numerical simulations. The last section concludes the work.

Non-Negativity and Boundedness of Solutions
The model (1) represents a system of six delayed differential equations. For such kind of problems, initial functions have to be stated and the functional framework needs to be specified. Let X = C([−τ, 0]; R 6 ) be the Banach space of continuous mapping from [−τ, 0] to R 6 supplied by the sup-norm ϕ = sup −τ≤t≤0 ϕ(t). The initial functions of the problem verify the following: Also, for biological reasons, these six initial functions X(θ), Y(θ), D(θ), V(θ), W(θ) and Z(θ) have to be non-negative: We have the following result about the boundedness and the positivity of any solutions of the system (1): Theorem 1. For any initial functions X(θ), Y(θ), D(θ), V(θ), W(θ) and Z(θ)) verifying (2) and (3), the system (1) has an unique solution; moreover, this solution is non-negative and bounded for all t ≥ 0.
Proof. By the classical theory of the functional differential equations (see for instance [23], and the references therein), we know that there is a unique local solution (H(t), I(t), D(t), V(t), W(t), Z(t)) to system (1) in [0, t m ), where t m is a finite number. By using the system (1), we have By recursive argument we get that For the boundedness result of all the solutions, we will consider the following functional: Therefore, when we use (1), we have because of the fact η 2 ∈ [0, 1], we have 1 − η 2 ≤ 1, from which, it follows Assuming that ρ = min(µ, δ 2 , u, h, b), we obtain then, this shows that H(t) is bounded, and so are the other functions X(t), Y(t), D(t), V(t), W(t) and Z(t). Therefore, every local solution can be prolonged up to any time t m > 0, which means that the solution exists globally.

Steady States
By simple calculation the system (1) has the following disease free equilibrium Indeed, the system (1) has four steady states other than E f :

The Optimization Problem
To state the optimization problem, we first suppose that η 1 and η 2 vary in time. The problem (1) becomes then For this problem, we will have the following result of the boundedness and the positivity of any solutions: Theorem 2. For any initial conditions X(θ), Y(θ), D(θ), V(θ), W(θ) and Z(θ)) verifying (2) and (3), the system (6) has a unique solution; moreover, this solution is non-negative and bounded for all t ≥ 0.
Proof. Using the classical theory of the functional differential equations [23], it is clear to see that there is a unique local solution ( From the system (6), we have from all these previous equalities, we deduce that all solutions are non-negative in t ∈ [0, t m ). About the boundedness, we will prove that the solutions are bounded in each interval [nτ, (n + 1)τ] such that n ∈ N.
We will begin with n = 0. Let t ∈ [0, τ], from the first equation of system (6), we obtain so, this means that X is bounded. From the second equation of (6), we obtain therefore, (2) and (3), we have the fact that From the third equation of (6), we obtain this inequality implies that from the boundedness result of I, one can conclude that D is bounded. From the fourth equation of (6), we obtain then, from the boundedness result of D, we conclude that V is also bounded. From both the fourth and the fifth equations of (6), we obtain from the boundedness results of D and V, we deduce that W is bounded. From the second and the last equation of system (6), we obtain from the boundedness results of X, Y and V, it follows the result that Z is bounded. By following the same analysis as before, for each single interval [nτ, (n + 1)τ] with n ≥ 1, one can conclude that all the solutions are bounded for all t ≥ 0. Therefore, every local solution can be prolonged up to any time t m > 0, which means that the solution exists globally.
Let us consider the following objective functional: where t f is the time period of therapy and the two positive constants A 1 and A 2 are based on the benefit-cost of the therapy η 1 and η 2 , respectively. The two control functions, η 1 (t) and η 2 (t) are supposed to be bounded and also Lebesgue integrable. Our main purpose is to maximize the objective functional defined in the Equation (7) by maximizing the number of the uninfected cells, increasing the CTL immune responses and the antibodies, decreasing the viral load and also decreasing the cost of treatment. That means, we are seeking an optimal control pair (η * 1 , η * 2 ) such that where U is the control set given by

An Optimal Control Existence Result
The two optimal control pair existence result can be obtained via the results [24,25]. Indeed, we have the following result: Theorem 3. There exists an optimal control (η * 1 , η * 2 ) ∈ U such that Proof. To use the existence result [24], we should first check the following properties (C 1 ) The set of the corresponding state variables and controls is nonempty.
(C 2 ) The set U is closed and convex.
(C 3 ) The right hand side of the state system is bounded by a linear function in the state and control variables. (C 4 ) The integrand of the objective functional is concave on U. (C 5 ) There exists an α > 1 and two constants c 1 , c 2 > 0, such that the integrand I(X, W, Z, η 1 , η 2 ) of the objective functional satisfies where The boundedness of the state system equations with the two controls (6) ensures us the existence of a solution. We can therefore deduce that the set of controls and the corresponding state variables are non-empty, this gives us the condition (C 1 ). The control set is convex and closed by definition, which ensures the condition (C 2 ). Moreover, since the system of state is bi-linear in η 1 , η 2 , the right hand-side of (6) verifies condition (C 3 ), using the fact that the solutions are bounded. About the condition (C 4 ), we obtain the Hessian matrix for I as follows, its determinant is stated as follows, then I is concave on U.

Numerical Results
To illustrate the numerical simulations, we implement and solved numerically our optimality system by the finite difference approximation method [27][28][29]. We obtain the following algorithm (Algorithm 1):

Algorithm 1:
The forward-backward finite difference numerical scheme.

. end for
Using values of parameters from [12,17]; i.e., s = 2.6 × 10 7 , k = 1.67 × 10 −12 , µ = 0.01, δ = 0.053, a = 150, β = 0.87, u = 3.8, τ = 5, λ = 1.1 × 10 −2 , q = 10 −12 , g = 10 −4 , h = 0.1, p = 0.01, b = 0.2, c = 0.03, A 1 = 50, 000 and A 2 = 5000. The role of the two parameters A 1 and A 2 is to calibrate the terms size in the system equations. Figure 1 depicts the evolution of the uninfected cells as function of time for both cases with and without control therapy. It is shown that with control the number of the uninfected cells is higher than those observed for the case without control. This result support the fact that the control strategy is to maximize the number of the healthy cells. From Figure 2, one can observe that the plot representing the infected cells with control strategy converges towards 0.002, however without any control therapy it converges towards 7.69, which proves that administrating this treatment will help the patient by a significant reduction of the infected cells.  Figure 3 illustrates the evolution of the capsids during the period of observation. It was established that with a control strategy, the amount of capsids vanishes after the first weeks of the administrated therapy. Meanwhile, without any control strategy this number remains at very high positive level, 4.37 × 10 3 .
The role of the two administrated therapy controls is also remarked in Figure 4. It was shown that with therapy control, the number of virions dies out after the first weeks of therapy, while without any control strategy it remains equal to 11.53. This indicates clearly the impact of the administrated therapy in controlling the HBV viral replication.
The antibody immune response is clearly affected by the control. This is illustrated in Figure 5; indeed, with control, the curve of antibodies converges towards zero; however, without any control strategy it converges towards 33.09 which clearly indicates the importance of adding the antibody component to HBV viral dynamics.
The CTL cell dynamics are also affected by the optimal control. This is shown in Figure 6; indeed, with control, the curve of CTL cells converges towards zero; however, without any control it converges towards a high level of 3.41 × 10 3 ; this reveals the role of the CTL component in blocking HBV viral infection.   The two optimal controls η 1 and η 2 representing inhibiting new infections and blocking viral replication are represented in Figure 7. The two curves present the treatment administration schedule for the period of observation. Both of the controls start from zero and when the first immune boosting drug is administered at full scale, the second drug is at its lowest for the first days of observation. During the last days of observation both the therapies should be administrated at their full scales. In this case the new infection is totally blocked. It is worthy to notice that if we compare the numerical simulation of this paper and the recent work [17], we remark that the antibodies have a clear effect in maximizing the level of the healthy cells and reducing the amount of the free viruses. Also, we remark that the presence of antibodies improves the effectiveness of both treatments after the first 60 days by reaching their maximum value. However, without antibodies, only the second treatment presents efficacy (see [17]) which leads to the elimination of HBV viruses. The obtained numerical results show that with antibodies and the optimal control, we observe a significant reduce of the HBV infection. All our numerical simulations will performed during the acute HBV infection period. This period is also known as the early stage of the infection [12,30]; however it will be very useful to predict infection for a chronic type of the infection.  Figure 7. The optimal control η 1 (left) and the optimal control η 2 (right) versus time.

Conclusions
We have modeled HBV in regards to intracellular behavior, capsids and adaptive immunity. The considered adaptive immune system is represented by the cytotoxic T-lymphocyte cells and the antibodies. The model under consideration includes six nonlinear differential equations that describe the dynamics that occur between hepatitis B free viruses (HBV), HBV DNA-containing capsids, hepatocytes, the antibodies and the CTL cells. An intracellular infection time delay and the effect of two drugs are incorporated into the suggested model. We have established the existence and uniqueness of the optimal controls via Pontryagin's maximum principle. The problem was implemented and solved numerically using backward and forward finite numerical difference schemes. It was established that with the two administrated optimal therapies, the amount of the healthy hepatocytes increases considerably while the number of infected hepatocytes decreases remarkably. Moreover, it was also shown that, with the control strategy, the viral load decreases significantly comparing with the model without control case, and this may boost the patient's life quality. Finally, we would like to mention that the used optimal controllers are given by open-loop, it will be useful to test other feedback control methods as a new predictive control model.