Optimal Control of Drug Therapy in a Hepatitis B Model

Combination antiviral drug therapy improves the survival rates of patients chronically infected with hepatitis B virus by controlling viral replication and enhancing immune responses. Some of these drugs have side effects that make them unsuitable for long-term administration. To address the trade-off between the positive and negative effects of the combination therapy, we investigated an optimal control problem for a delay differential equation model of immune responses to hepatitis virus B infection. Our optimal control problem investigates the interplay between virological and immunomodulatory effects of therapy, the control of viremia and the administration of the minimal dosage over a short period of time. Our numerical results show that the high drug levels that induce immune modulation rather than suppression of virological factors are essential for the clearance of hepatitis B virus.


Introduction
Hepatitis B virus (HBV) is the leading viral cause of liver disease, affecting 250-350 million people worldwide.Chronic HBV leads to the development of liver cirrhosis and liver cancer [1].Despite the availability of effective HBV vaccination [2], the prevalence of chronic HBV has only marginally declined [3,4].The natural course of chronic HBV includes an immune-tolerant phase, hepatitis B e-antigen positive immuno-active and -inactive phases and a hepatitis B e-antigen negative immuno-active phase [5].Understanding the virological and immunological characteristics of each of these stages can provide a useful framework for the management of chronic HBV [6,7].
Currently, seven drugs have been approved for treating chronic HBV disease: standard and PEGylated interferon (IFN-α) and five nucleo(t)side analogs (NAs) [8].These medications suppress HBV replication and liver inflammation, but do not lead to cure.The interferon-α treatments modulate immune responses that may lower viral levels.It is given for a finite time (usually 12 months) due to its toxic side effects [8,9].The nucleo(t)side analogues are administrated for many years and sometimes for life and are responsible for viral suppression.However, life-long therapy is difficult due to costs, side effects, compliance and, most importantly, the development of antiviral drug resistance [10].Combination therapy has not shown an increased effect on treatment response (but has reduced the rate of drug resistance) [11].This has made it difficult to establish a universal guideline for treatment start, duration, which type or combinations of drugs to use [7] and how to define the success of therapy (virologically, serologically and/or immunologically) [7,12].
To provide insight into the optimal combination therapy and, in particular, into the immune-mediated effects of interferon-α, we design an optimal control study for a mathematical model of hepatitis B infection.Mathematical models have been used to address transition from acute to chronic HBV infections [13][14][15][16] and to study the effects of drug therapy [17][18][19].Optimal control theory was developed by Pontryagin et al. for obtaining necessary conditions to characterize optimal controls for systems of ordinary differential equations (ODEs) [20].Optimal control models have been used previously to design treatment strategies for disease models described by systems of ordinary differential equations with no delays [21][22][23][24][25] and systems of delay differential equations (DDEs) [26].Kharatishvili [27] developed the extension of the Pontryagin's maximum principle for systems of DDEs with constant time delays.For approximation and numerical methods for such problems, see [28,29].
In this paper, we modify an in-host DDE model of immune responses to hepatitis B infection introduced in [14] by adding the effects of combination drug therapy.In particular, we consider the drug effects in blocking viral production, reducing viral infection, enhancing the killing of infected cells by immune responses and removing immune cell exhaustion.We will use this as a starting model for designing an optimal control problem that advises what is the best combination therapy to ensure viral clearance, immune activation and the least amount of liver damage.Time-varying rates for therapies will be the controls in the system.Optimal control has been applied to the study of hepatitis B therapy using ODE [30,31] and DDE [32] models that did not consider an immune system component and using an ODE model that considers the enhancement of immune responses following administration of traditional Chinese medicine [33].
Here, we introduce the DDE system with constant time delay from [14], derive its corresponding control formulation and use it to determine the temporal, quantitative and qualitative effects of the drugs that lead to hepatitis B virus control, balancing the goals of reducing viral load and minimizing the negative side effects of therapy.Our results indicate that early drug therapy that mainly modulates and restores the immune responses against the virus is mandatory for the success of the therapy.

Model of Hepatitis B Virus (HBV) Drug Therapy
To understand the various modes of action of antiviral therapy, we modify the model of acute infection published in [14] by considering the combined effects of NAs and interferon-α drugs.As in [14], we consider five state variables, corresponding to uninfected hepatocytes (T), productively-infected hepatocytes (I), free virus (V), immune effector cells (E) and a population of refractory hepatocytes (R).
All hepatocyte populations, uninfected, infected and refractory to reinfection, are maintained by homeostasis described by a logistic equation, with carrying capacity K and maximal per capita growth rate r.Virus infects target cells at rate β.Infected cells are killed by the immune responses at rate µ.As in [14], we assume that infected cells can be lost due to the non-cytolytic response at rate ρ, dependent on the effector cell population E [34], and move into a refractory class R. Refractory cells will still be assayed as infected, since surface antigens persist for some time [35].We assume that they have lost their viral replicates and do not produce virus, which also makes them poor targets for cytotoxic T lymphocyte (CTL) responses.Therefore, the refractory population will be killed at a smaller rate, ν < µ.The refractory state is not permanent, and the R population may eventually become susceptible to reinfection.We model this by allowing R cells to move into the uninfected population at rate q.Free virus is produced at rate π and cleared at rate c.
In the absence of infection, we assume the immune effector cells E are at equilibrium value s/d, where s corresponds to a source of effector cells specific for HBV and 1/d is their average life-span.
Upon encountering antigen (on the surface of infected liver cells), these cells expand at rate α, with a constant time delay τ accounting for the lag between antigen encounter and effector cell expansion.
We model drug therapy as interference with virus production and infection and as immune modulation, as follows.The contributions of interferon-α are complex, with both direct antiviral and immunomodulatory effects.Here, we consider three effects of IFN-α: enhancement of effector cell killing rates, reduction of viral production rate and the delay of effector cell programmed death [36,37].We model this by changing rates µ, π and d to µ ), respectively.In each case, 0 ≤ ≤ 1 represents the efficacy of IFN-α, and a i ≥ 0 are scalar parameters representing the strength of the corresponding effect of interferon.a 1 can be any nonnegative number, while a 2 and a 3 are between 0 and 1.
The nucleos(t)ide analogues, on the other hand, interfere with both the ability of virus to infect a cell and with the generation of HBV DNA by an infected cell [38].We model this by assuming that the infectivity rate in the presence of NAs becomes and the viral production rate becomes Here, 0 ≤ η ≤ 1 represents NAs' efficacy; b i ≥ 0 are scalar parameters representing the relative strength of the corresponding effect of NAs; and f and 1 − f represent the relative contribution of interferon-α and NAs, respectively, to reducing viral production, for 0 ≤ f ≤ 1.
For t > t d , where t d represents the time of therapy onset, the dynamics of these populations is governed by the following differential equations As this is a system of delay differential equations with fixed delay τ, the initial conditions for the system are defined on the interval [t d − τ, t d ] by functions: In many cases, these will be assumed to be constant and equal to the value of the state variables at t = t d .

Model Analysis
We will consider that therapy leads to viral removal and only investigate the long-term behavior of the equilibrium where V = 0. Model (1) has such a disease free steady state, S 0 = (K, 0, 0, s/d 1 , 0).In the absence of delay (τ = 0), the Jacobian matrix associated with (1) is: We evaluate J(S 0 ) and obtain the Jacobian matrix: and the characteristic equation: where By the Routh-Hurwitz criteria, the characteristic equation has roots with negative real parts if and only if We can show that these conditions are satisfied when A 4 > 0.
Therefore, when τ = 0, S 0 is locally asymptotically stable when: and is unstable otherwise.This condition gives us a minimal drug efficacy for viral clearance in the absence of delay in immune activation.When τ > 0, the characteristic equation at state S 0 is: where: and: We see that the transcendental Equation ( 5) reduces to polynomial: as in the τ = 0 case.Therefore, the infection-free equilibrium S 0 is locally asymptotically stable for all τ when (4) holds.

Numerical Results
We assume that in the absence of drug therapy, the virus will persist, and the patient will experience chronic infection.We therefore set all parameters values and initial conditions to the values for the chronically-infected Patient 7 in [14] (see Table 1).The dynamics of the five variables in the absence of drug therapy, = η = 0, throughout acute infection and transition to chronic disease is presented in Figure 1.Effector cell clearance rate 0.5 day −1 q Waning of refractory cell immunity 2 × 10 −5 day −1 We next investigate the change in the model dynamics in the presence of treatment.Initially, we set time-constant drug efficacy, η = 0.5 and = 0.9 as in [17]; assume scalars a i and b i to be either 1 or 0, representing an effect or lack of effect in therapy; and set f = 0.5.We define HBV DNA clearance as the presence of less than one HBV DNA in the host serum.Since we assume that HBV DNA can distribute throughout the 3 liters of serum in an average 70-kg person, the viral extinction becomes V ≤ V ext = 3 × 10 −4 copies per mL.We notice that the speed of viral extinction is dependent on both the timing of therapy initiation and on the type of effects considered.
Indeed, when the therapy is started at t d = T peak = 96 days post-infection (corresponding to the peak HBV DNA), our model predicts HBV DNA clearance 145 days after the start of therapy when 2a, dashed line) and 2226 days (6.1 years) after the start of therapy when a 1 = a 2 = b 1 = b 2 = 1, a 3 = 0 and f = 0.5 (see Figure 2a, dotted line).This suggests that the immune modulation effect of interferon is important in clearing the HBV DNA in a short amount of time.
Similar immune-mediated effects of interferon are observed when drug therapy is started at the chronic state of Model (1), t d = T ss = 20 years post-infection.Indeed, viremia clearance occurs faster when 3a, dotted lines).However, viral clearance is delayed by 12 days compared to the case in which therapy is initiated at peak HBV DNA.The rapidity of clearance during peak therapy is due to the transient effects of immune cells, which are activated and expanding at the peak HBV DNA (see Figure 2b, solid lines) and tolerant at the chronic HBV infection (see Figure 3b, solid lines).Temporal evolution for the variables in Model (1) without drug therapy, i.e., = η = 0, and the parameters in Table 1.The circles represent patient data.effector cells E per mL given by Model (1) and the parameters in Table 1 for: η = = 0 (solid lines); η = 0.5, = 0.9, a 1 = a 2 = a 3 = b = b 2 = 1, f = 0.5 (dashed lines); and η = 0.5, = 0.9, Here, time t = 0 corresponds to both the start of therapy and the peak viral load, i.e., the initial conditions are T(0) = 1.32 × 10 5 cells per mL, I(0) = 1.23 × 10 5 cells per mL, V(0) = 3 × HBV DNA per mL, E(0) = 20.23 cells per mL and R(0) = 1.23 × 10 6 cells per mL.  1) and parameters in Table 1 for η = = 0 (solid lines); η = 0.5, = 0.9, Here, t = 0 corresponds to both the start of therapy and the steady state of the viral load, i.e., the initial conditions are T(0) = 3 × 10 5 cells per mL, I(0) = 3.8 × 10 4 cells per mL, V(0) = 9.33 × 10 6 HBV DNA per mL, E(0) = 20.3 cells per mL and R(0) = 1.33 × 10 7 cells per mL.
These results are supported by relative sensitivity curves.Briefly, sensitivity functions are numerical solutions of the following system: where , the function g represents the right-hand side of Model (1) and: The partial derivatives ∂x i /∂ξ j , for i = 1, . . ., 5 and j = 1, . . ., 20, are time functions denoting the rate of change in a state variable with respect to variations in model parameters (see [39] and the references therein).The functions ξ j /x i × ∂x i /∂ξ j are the relative sensitivity curves (similar to relative error), which allow for the comparison between the sensitivity of two variables x 1 i and x 2 i with respect to the same parameter ξ j .Here, we are interested in the sensitivity of V(t, ξ) and E(t, ξ) with respect to model parameters ξ = {a 2 , a 3 , b 2 }.As the numerical solutions displayed in Figure 4 show, the a 3 effect is the most important drug effect in both reducing virus load and increasing immune response.Our model suggests that the timing and the type of interferon effects are important in the success of the treatment.In the next sections, we will consider a time-dependent effect of both interferon and nucleo(s)tide analogues, and we will formulate an optimal control problem for our model to determine the connection between successful therapy and the optimal temporal efficacy of the effects considered.

Optimal Control Problem
To better understand the interactions between drug therapy and the host immune reaction, we allow the effects of the two drug types (NAs and IFN-α) to vary with time.We replace the NAs' efficacy parameter η with u 1 (t) to represent the time-varying effective drug dosage needed for optimal therapy.Similarly, u 2 (t) will replace the parameter to represent the time-varying dosage of IFN-α.As before, a 1 , a 2 and a 3 indicate the relative strength of the three effects of IFN-α under consideration; b 1 and b 2 represent the relative strengths of the two effects of the NA; and f and 1 − f are the relative contribution of interferon and NAs, respectively, in reducing viral production.
For computational convenience and consistency of notation, we will relabel the state variables as in Table 2.With these adjustments, our model becomes as follows.1) and ( 8).
Our goal is to find the regime of drug therapy that minimizes the objective functional: subject to the system of delay differential Equation ( 8) with initial conditions (2).Here, T represents the duration of therapy.We will study two different scenarios for drug treatment by considering two different sets of initial conditions.For the case of chronic HBV infection, the initial values of the state variables x i will be constant and equal to their chronic infection equilibrium values for the entire interval [t d − τ, t d ].We will also study the case of drug therapy initiation at the peak of viral load, during acute infection.In this case, the initial functions for the state variables x i will be set by their trajectories during acute infection in the absence of treatment.This is described in more detail in Section 3.3.
Given upper bounds M i on u i (for i = 1, 2), we seek to find an optimal pair in the control set: The integral portion of the objective functional encapsulates the goal of minimizing total virus concentration (x 3 ), infected cell concentration (x 2 ) and the amount of drug used (u 1 and u 2 terms) over the entire treatment period.The final term c 5 x 3 (T) represents the goal of minimizing the final viral concentrations at time T, which ideally would be below clearance levels.The parameters c i > 0 and ε j > 0 give the relative weight of each of these factors.

Analysis of the Optimal Control Problem
Pontryagin's key idea was to use adjoint functions to attach the state dynamics to the objective functional.This converts the problem into minimizing the Hamiltonian and generates the adjoint differential equations and final-time transversality conditions.In this paper, we use a generalization of Pontryagin's maximum principle to systems of delay differential equations (developed by [27]).
Letting f denote the integrand of the objective functional (9) and g 1 , . . ., g 5 be the right-hand sides of System (8), the Hamiltonian for our optimal control problem is: where the adjoint functions λ i correspond to the states x i , for i ∈ {1, . . ., 5}.
From the Hamiltonian, we derive the adjoint equations.For i = 1, 3 and 5, we have: on [0, T].For i = 2 and 4, on the interval [0, T − τ] we have where z i (t) = x i (t − τ) represent the delayed state variables, while on the interval [T − τ, T], we have: just as in [11].
The equations containing delayed state variables produce adjoint equations with a "forward" delay, due to the opposite time orientation of the adjoint differential equation.The adjoint equations are subject to the final condition that λ 3 (T) = c 5 , λ i (T) = 0 and i ∈ {1, 2, 4, 5}.On the interior of the control set, minimizing the Hamiltonian gives: , and then, the optimal control pair becomes: Using the bounds on the controls, we obtain the following characterization of the optimal control:

Implementation of the Optimal Control Problem
We wish to find optimal controls numerically by applying a forward-backward iterative method [25].Initially, a constant control is assumed, and the state equations are solved in the forward-time direction from our standard set of initial conditions.Given this solution to the state equations, the adjoint equations are then solved in the backwards-time direction, beginning with the final time T and final condition λ 3 (T) = c 5 and λ i (T) = 0 for i ∈ {1, 2, 4, 5}.The controls to be used for the next forward run are then updated using the characterization of the optimal control given in (15), and the forward-backward solution process is repeated with the updated control functions.This process is iterated until the controls and the solutions to all of the differential equations converge to within acceptable numerical tolerances.See [28,29,40] for the background on this procedure and the approximation of delay equations.
To use the forward-backward sweep, we rewrite the adjoint equations for our system as: Notice that in the time interval [T − τ, T], the advance terms (those with argument t + τ, found in the second and fourth equations) drop out, so we have five ordinary differential equations.On the interval [0, T − τ], we once again have advance equations, but the solutions to the ODEs on [T − τ, T] provide the required initial data to solve these equations.Thus, the adjoint equations are advance differential equations on [0, T − τ] and ordinary differential equations on [T − τ, T], subject to the final condition λ 3 (T) = c 5 and λ i (T) = 0 for i ∈ {1, 2, 4, 5}.
For our numerical simulations, we used the built-in MATLAB numerical delay differential equation solver, dde23.This tool is not capable of solving advance equations directly, so we made a change of variables to convert the advance differential equations to a system of delay differential equations.Specifically, we define a new reversed-time variable σ = T − t and new adjoint variables In terms of these new variables, adjoint equations for L i (σ) are ordinary differential equations on [0, τ] and delay differential equations on [τ, T], subject to the initial condition L 3 (0) = c 5 and L i (0) = 0 for i ∈ {1, 2, 4, 5}.The new system was solved numerically using dde23 as part of the implementation of the forward-backward method to find the optimal control.

Numerical Results with Optimal Controls
For optimal therapy during chronic HBV infection, we assume that at the time of treatment initiation, the patient has reached chronic steady state values x1 = 3 × 10 5 cells per mL, x2 = 3.8 × 10 4 cells per mL, x3 = 9.33 × 10 6 HBV DNA per mL, x4 = 20.3 cells per mL and x5 = 1.33 × 10 7 cells per mL.The terms in the objective functional (9) have different scales.In order to weigh each term equally, we choose parameters c i , such that c i x i (for i = 1, 2) and c i u i (for i = 3, 4) are one at the start of the optimal control.When the optimal control is started at the viral steady state, we normalize variables x 1 and x 2 by c 1 = 1 x3 and c 2 = 1 x2 = π c x3 .As a result, the factors c 1 x 1 and c 2 x 2 are one at the beginning of treatment and between zero and one from there on.Since we assume 0 ≤ u 1 (t), u 2 (t) ≤ 0.9 for all t, their corresponding weights are c 3 = c 4 = 1.Lastly, we choose the effects for the convex terms to be 1 = 2 = 0.1.The treatment will be given for T 1 = 400 days, which is sufficient for virus clearance and is still within the timeline of the HBV therapy guidelines [7].
The ideal clinical outcome of HBV drug therapy is to achieve viral suppression, such that x 3 ≤ V ext = 3 × 10 −4 HBV DNA per mL, corresponding to one HBV DNA in the body (as in Section 2.3).To account for this final condition, we set c 5 = 3 × 10 3 so that the final time condition is normalized to one.Under these assumptions, we run the optimal control problem over T 1 = 400 days with the aim of finding the best temporal drug usage at each time point that ensures viral clearance while minimizing the drug levels.
We will be considering several scenarios by varying the levels of control effects, described by constants a i and b i .Our simulations show that for high virus levels (above extinction), the optimal control will always be the maximal drug dosage for the scenario considered.For example, for a case with any a i > 0 and b 1 = b 2 = 0, the optimal control is u 1 (t) = 0 and u 2 (t) = 0.9 for all time points.Similarly, if a i > 0 and b j > 0, then the optimal control is u 1 (t) = u 2 (t) = 0.9 for all time points.Therefore, for the first 325 days, we assume the maximal drug dosage for a given scenario and only run the optimal control with variable u 1 (t) and u 2 (t) for the final 75 of the 400 days of treatment.
Initially, we include all effects of nucleos(t)ide analogues u 1 and interferon u 2 , i.e., a 1 = a 2 = a 3 = b 1 = b 2 = 1 and f = 0.5.For this choice of parameters, optimal control requires: (i) maximal NA drug dosage u 1 = 0.9 for 0 ≤ t ≤ 347 days and zero drug dosage for t ≥ 347 days; and (ii) maximal u 2 = 0.9 interferon dosage for 0 ≤ t ≤ 366.5 days, zero dosage for 366.5 ≤ t ≤ 399 days and maximal dosage after that (see Figure 5a, top row).Under these drug regimes, HBV DNA is cleared 331 days after the start of treatment (see Figure 5b, top row).(a) Optimal evolution for controls u 1 (dashed lines) and u 2 (solid lines) as given by ( 15); (b) virus V per mL (solid lines) and infected cells I per mL (dashed lines) over time; (c) effector cells E per mL over time when drugs are introduced at equilibrium virus concentration and: We next assume monotherapy with interferon-α.We set b 1 = b 2 = 0, a 1 = a 2 = a 3 = 1 and change the anti-viral production effect to f = 1.The optimal control problem predicts that the virus is cleared 338 days after the start of therapy, seven days later than under combination therapy (see Figure 5b, second row).Another downside of this monotherapy is the need for interferon dosage to be maximal (u 2 = 0.9) for 0 ≤ t ≤ 377 and for t ≥ 396.5 in order to compensate for the lack of NAs (see Figure 5a, second row).
In combination therapy, when u 2 neither increases CTLkilling abilities nor reduces virus production, i.e., a 1 = a 2 = 0, and all other effects are maximal a 3 = b 1 = b 2 = 1, f = 0 (so that NAs have maximal effect in blocking viral production), the HBV DNA will clear 341 days after the start of treatment (see Figure 5b, third row).The controls are: (i) u 1 = 0.9 for 0 ≤ t ≤ 358 and t ≥ 396, and zero otherwise; and (ii) u 2 = 0.9 for 0 ≤ t ≤ 369.5, and zero otherwise (see Figure 5a, third row).
Lastly, we consider monotherapy with NAs alone.We set a 1 = a 2 = a 3 = 0, b 1 = b 2 = 1 and f = 0.The optimal control model does not result in viral clearance even when u 1 = 0.9 throughout the duration of the treatment (see Figure 5a,b, bottom row).
For the first three cases, we find that an increase in the effector cells' lifespan is needed for viral suppression.Indeed, when a 3 = 1 and u 2 > 0, the increase of CTL lifespan to 1/d(1 − a 3 u 2 ) leads to elevated CTL concentrations (see Figure 5c, top three rows) and subsequent HBV DNA removal (see Figure 5b, top three rows).By contrast, when a 3 = 0, HBV DNA does not reach extinction during the T 1 = 400 days of therapy (see Figure 5b, bottom row) due to low HBV-tolerant CTL concentrations, E = s/d 1 = s/d, which do not expand in the presence of HBV (see Figure 5c, bottom row).
For optimal therapy during acute HBV infection, to further determine the relationship between CTL activation and the possible success of short-term anti-HBV therapy, we run the optimal control problem during acute HBV disease, where CTL activation has been reported [41].We start by running the DDE system mainly to the peak virus concentration occurring at time T peak , record the values of all variables for the times −τ + T peak ≤ t ≤ T peak and start the optimal control problem at t d = T peak .The treatment will be given for T 2 = 95 days and, as before, we include the following weights As in the chronic HBV case, we consider four scenarios.If we include all effects of nucleos(t)ide analogues u 1 and all effects of interferon u 2 , i.e., a 1 = a 2 = b 1 = b 2 = 1, f = 0.5 and a 3 = 0.9, we predict that HBV DNA is cleared 90 days after the start of treatment when interferon dosage is maximal at each time step and NAs are zero for 0 ≤ t ≤ 90 days and maximal from t ≥ 90 (see Figure 6a,b, top row).This result implies that NAs do not have a role in viral clearance.This is corroborated by the optimal control solution for interferon monotherapy.Indeed, when we set a 1 = a 2 = 1, a 3 = 0.9, b 1 = b 2 = 0 and f = 1, virus is cleared even faster, 87 days after the start of therapy (see Figure 6b, second row) when u 2 dosage is maximal u 2 = 0.9 at each time step (see Figure 6a, second row).This is due to the higher efficacy of interferon-α at blocking viral production assumed in this scenario.The results in both of these cases are due to fast expansion and transient maintenance of CTLs to high levels of 2500 cells per mL (see Figure 6c, top two rows).
To further determine which of the interferon effects are the most influential, we remove the first two interferon effects a 1 = a 2 = 0, keep a 3 = 0.9, b 1 = b 2 = 1 and set f = 0.Under these conditions, virus is cleared 91 days after the start of therapy, one and four days later than the previous two cases (see Figure 6b, third row) when: (i) NA dosage (u 1 ) is zero for 0 ≤ t ≤ 87 and maximal (u 1 = 0.9) afterwards; and (ii) interferon dosage is maximal (u 2 = 0.9) for 0 ≤ t ≤ 91 and zero afterwards (see Figure 6a, third row).This result implies, again, that interferon is needed for virus clearance.
Lastly, when we consider monotherapy with NAs, i.e., b 1 = b 2 = 1, a 1 = a 2 = a 3 = 0 and f = 0, the HBV DNA will not clear (see Figure 6b, bottom row) in spite of CTL levels being higher than the base of E = s/d 1 = s/d (see Figure 6c, bottom row).Interestingly, our numerical results show that the optimal NA dosage for the first 42 days is zero (see Figure 6a, bottom row).
Our study has used a delay of τ = 33.4days, since that represented the delay in the CTL activation in the only patient that developed chronic disease in [14].To check whether the size of the delay affects the results, we have investigated the optimal control problem for a shorter delay of τ = 15.2 days, which is the smallest delay among the patients in [14].We found that the length of the delay does not influence the results when the therapy is started at the viral steady state (not shown).When the therapy is started at the peak viral load, we find, not surprisingly, that a shorter time lag speeds viral decay and CTL expansion (see Figure 7b,c).Moreover, the optimal control problem for the shorter delay of τ = 15.2 days suggests no therapy (u 1 = u 2 = 0) for the first three days post-peak, maximal interferon therapy for 3 ≤ t ≤ 75 and t ≥ 90 and no interferon u 2 = 0 otherwise.The role of NAs is transient, with maximal dosages u 1 = 0.9 for 20 ≤ t ≤ 27, 38 ≤ t ≤ 51 and t > 90 days and u 1 = 0 elsewhere (see Figure 7a, bottom row).
Together, these results suggest that the immunomodulatory effects of interferon are needed for HBV DNA control and that immune modulation rather than suppression of virological factors is essential for inducing HBV clearance.

Discussion
The management of chronic hepatitis B is based on guidelines regarding screening, diagnosis, duration of treatment, adherence and monitoring of immune and virological markers [12].Combined antiviral therapy improves the survival rates of patients chronically infected with hepatitis B virus [42].While the first goal of therapy is to reduce HBV replication, the ultimate goal is the removal of hepatitis B surface and e-antigens and of covalently-closed circular DNA (cccDNA).The removal of the hepatitis B e-antigen is associated with immunological factors, such as removal of the tolerant status of the hepatitis B-specific T cells [36].Such effects have been demonstrated by interferon therapy, which leads to improved immune function, sustained response rates with combination therapy and improved overall prognosis [37].These drugs have limitations, such as side effects, the use of injection and poor response in patients with compromised liver function.That makes them unsuitable for long-term administration [43].
To address the positive effects of interferon therapy and to account for its limitations due to adverse effects and limited time usage, we developed a control problem that accounts for: (i) three interferon functions: increased CTL life-span, reduced viremia and increased CTL killing of infected hepatocytes; (ii) two nucleos(t)ide analogue effects: blocking of viral infection and of viral production; (iii) hepatitis B DNA decay below the limits of detection; and (iv) minimal dosage administration over a short time period.
The presence of significant side effects and the need for lengthy treatment make the question of the optimal therapy strategy relevant for the case of chronic hepatitis B. Little work has been previously done on the optimal control of HBV treatment, and the existing models do not include immune system involvement [30].Ciupe et al. demonstrated that consideration of delayed cytotoxic and non-cytotoxic immune reactions and the presence of cells refractory to infection was necessary to properly understand the dynamics of HBV acute infection and progression to chronic disease [14].In this study, we derived an optimality system associated with this model, and we constructed the corresponding adjoint equations, which differed from the construction of adjoint equations for systems of ODE [44].In particular, we found that for a delay τ and a control period [0, T], the adjoint equations are ODEs on [T − τ, T] and advance differential equations on [0, T − τ], meaning that there are terms with time argument t + τ.After deriving the proper form of the adjoint equations, we have investigated the optimal control system numerically using a forward-backward sweep method, as in [25].
The control problem indicates the need for high immunomodulatory effects of IFN-α until HBV DNA is cleared.Most importantly, the immunomodulatory effects that increase the survival of effector cells are essential for timely reduction in viremia, which is needed to limit the IFN-α-induced side effects [43].We predict that starting interferon therapy at the peak viral load, rather than at viral equilibrium, shortens the time to HBV DNA removal.This is due to enhancement of an already activated T cell response.Most interestingly, monotherapy with interferon-α is sufficient for virus control, while the effects of nucleos(t)ide analogues emerge only at the end stages of combination therapy.This result suggests that sequential single therapy (interferon followed by nucleos(t)ide analogues) may be the optimal course of action for both viral suppression and the limitation of drug effects.
Our results imply that an increase in T cell response (in acute infections) and reversal of T cell inactivation (in chronic infections) are essential for fast control of viremia.However, experimental studies have shown that only a small percentage of patients on INF-α therapy experience loss of surface-antigen, e-antigen loss and reduction in virus replication [45].Moreover, interferon-α is responsible for only partial reversal of T cell inactivation in chronic hepatitis B infections [37] and is not induced naturally during acute hepatitis B infections [46].However, interferon-α therapy induces cccDNA degradation in cell culture [47] and enhances the innate immune response mediated by natural killer (NK) cells in e-antigen negative patients [45].Our optimal control study predicts that enhanced cccDNA degradation, which can be incorporated in our model as an effect on the term ρ, has limited effect on the timing of HBV DNA removal (not shown).The NK activation would increase infected cell removal in a complementary fashion to the CTL effect we consider now.Similar to our current results, increased survival of NK cells would be needed for fast HBV DNA removal.
Our study has considered a linear combination effect of interferon and nucleos(t)ide analogues on the blocking of viral production f a 2 u 2 + (1 − f )b 2 u 1 with f = 0.5.To determine the effect of f on the results, we investigated two control therapies, (i) strong interferon influence on blocking viral production f = 0.9 and (ii) strong NA influence on blocking viral production f = 0.1, and found that the size of f has little influence one the timing of viral clearance for both peak and equilibrium therapy (not shown).
One limitation in our research comes from ignoring the consequences that prolonged treatment exerts on the evolution of the viral population.Studies have shown that life-long therapy and lack of compliance leads to the development of antiviral drug resistance [10], even though combination therapy helps reduce drug resistance [11].Our model has not considered the emergence of HBV variants or mutation in the presence of combination therapy.Further work is needed to address the optimal therapy in the presence of these events.
In conclusion, we have designed an optimal control study that shows that a successful short-term anti-HBV therapy requires the modulation of strong innate and/or adaptive immune responses, rather than induction of anti-virological effects.Such therapy needs to be active and elevated for the entire period of viremia.

Figure 1 .
Figure 1.Temporal evolution for the variables in Model (1) without drug therapy, i.e., = η = 0, and the parameters in Table1.The circles represent patient data.
1 and c 5 = 3 × 10 3 , so that all factors are normalized to one.

Table 1 .
Variables, parameters and values used in simulations with mililiter (mL) and day (d).

Table 2 .
Correspondence between labeling of the state variables in Models (