Dynamics and control strategy for a delayed viral infection model

In this paper, we derive a delayed epidemic model to describe the characterization of cytotoxic T lymphocyte (CTL)-mediated immune response against virus infection. The stability of equilibria and the existence of Hopf bifurcation are analysed. Theoretical results reveal that if the basic reproductive number is greater than 1, the positive equilibrium may lose its stability and the bifurcated periodic solution occurs when time delay is taken as the bifurcation parameter. Furthermore, we investigate an optimal control problem according to the delayed model based on the available therapy for hepatitis B infection. With the aim of minimizing the infected cells and viral load with consideration for the treatment costs, the optimal solution is discussed analytically. For the case when periodic solution occurs, numerical simulations are performed to suggest the optimal therapeutic strategy and compare the model-predicted consequences.


Introduction
Over the past years, mathematical modelling has grown as an essential approach to exploring biological mechanism of within-host viral infection. One such model was proposed by Nowak et al. in 1996 to analyse the effect of antiviral treatment on reducing viral loads [19], known as the basic model of viral dynamics. Since then, a variety of dynamic models has been established to describe the dynamical process of viral infection for various infectious diseases such as HIV, HBV and influenza, among others. These models allow investigating the interaction between disease progression and immune response, predicting the outcome of different therapeutic alternatives, testing specific hypotheses based on clinical data and evaluating the effect of potential factors on dynamical behaviours of virus load. In general, the existing models provide useful frameworks for the management of viral infection and valuable reference for further research.
Different efforts have been undertaken to model individual aspects of the immune response and its interplay with the virus. It is observed that time delay needed to generate immune response during an infection period, such as cytotoxic T lymphocytes (CTLs) responsible for killing infected cells is one of the key factors in viral dynamics [1,7,12,22,23,25,26,28,31]. When a within-host viral model is characterized by the immune lag described by discrete or distributed delay, complex dynamics including stability switches, periodic oscillations and even chaotic behaviours may occur. For a model incorporating delayed CTL response to HTLV-I infection, the existence of multiple stable periodic oscillations, arising from an overlap of multiple stable Hopf branches, was investigated analytically in [12]. Global stability and Hopf bifurcation in a delayed viral model were analysed in [25], revealing that considering only two intracellular delays did not affect the dynamical behaviour of the model, but incorporating an immune delay greatly affects the dynamics, i.e. an immune delay may destabilize the immunity-activated equilibrium and lead to rich dynamics.
Currently, optimal control has been used to design different interventions and administrations for epidemic models described by systems with no delays [14,18,27,30] and systems with delays [2,3,17]. It is a powerful way to target epidemic outbreaks by defining strategies to control the system and obtaining the best possible outcome. In particular, it can serve for seeking the optimal response for control schedule while minimizing the cost involved in the management. There are many dynamical models of infectious diseases that can be treated with optimal control theory, like tuberculosis [24], hepatitis B [6,10], hepatitis C [29], cancer [4] and so on. Applying Pontryagin's maximum principle [20], optimal control models of delayed differential equations have been considered to investigate the optimal therapy in [3,4,6] by constructing the corresponding adjoint equations.
In this paper, we will study the effect of time delay in the proliferation of CTLs and discuss the corresponding optimality system based on the available therapy for hepatitis B infection. Motivated by [1,26], we build a delayed model capable to present the characterization of CTL-mediated immune response against virus infection. Our aim is to conduct qualitative analysis of the model that may provide useful insights about the dynamical properties, including stabilities of the equilibria and the existence of Hopf bifurcation, and to advise what is the best therapy that ensures the least amount of infected cells and viral load as well as the costs of control strategy over a period of time. In the literature mentioned above, most of them only discussed one of the two aspects.
The organization of the paper is as follows. In the next section, we introduce the delayed model and give the preliminaries. The dynamical results are studied in Section 3 by establishing stabilities and discussing conditions for the existence of Hopf bifurcation. In Section 4, the optimal control problem is formulated and the numerical simulations of optimality system are presented. A brief conclusion can be found in Section 5.

Model formulation
Mathematical models for within-host infections typically include the variables of healthy target cells T(t), productively infected cells I(t), free virus particles V(t) and CTL cells Z(t). In recent years, models of viral infection incorporating virus-to-cell transmission and cell-to-cell transmission have been developed to describe the mechanisms underlying certain diseases [1,26]. Xu and Zhou [26] introduced an immune delay needed to activate the immune response in the model of HIV-1 infection, and used bilinear incidence to conduct the bifurcation analysis induced by the delay. Three time delays accounting for the period of chemical reaction in virus-to-cell infection, an intracellular incubation period in cellto-cell infection and the immune delay were considered in [1], and the bilinear incidence was also employed.
It is known that the incidence function in modelling viral infection plays an important role in determining qualitative behaviour of the proposed model and in giving a reasonable description of the dynamics. In this paper, we use general incidence function instead of bilinear form to describe the two routes of viral infection. It is assumed that the population of healthy target cells has a logistic growth function and there is only immune delay to investigate the possible effect of the latent period in production of the immune response. Our model is presented as follows: with initial conditions . In the logistic growth function, r is the intrinsic growth rate and the population of healthy target cells T(t) is limited by a carrying capacity of T max . The constant α (α ≥ 1) is the limitation coefficient of infected cells imposed on the growth of target cells. Parameters d 1 , d 2 and d 3 are the death rates of infected cells, free virus and CTL cells respectively. p denotes the killing rate of infected cells by CTL cells. Constant τ represents the time delay needed to activate CTLs. Free viral particles are released by infected cells at a rate kI(t), and the term cI(t − τ ) is used to present the delayed production rate of the effector cells.
In model (1), general incidence functions f 1 (V(t))T(t), f 2 (I(t))T(t) are used to describe virus-to-cell transmission and cell-to-cell transmission respectively. The sufficiently smooth function f i : R + → R + (i = 1, 2) is assumed to satisfy the following condition: which possesses biological meanings, for example, both of virus-to-cell and cell-to-cell transmission become faster as the populations of virus particles and infected cells increase.
From (A1), it is easy to check that ≤ 0 (similar for f 2 ), implying that the infections are possible to slow down due to some inhibition effect. In particular, the functions

Theorem 3.1: The solution Y(t, ϕ) of model (1) with initial condition (2) is nonnegative and ultimately bounded for all t ≥ 0.
Proof: Assume that there exists t 1 > 0 such that min{T(t 1 ), From the first two equations of (1), we have the following For sufficiently small > 0 and t ∈ (t 1 − 1 , t 1 ], it follows that Otherwise, if I(t 1 ) = 0, from the second equation of (1), we have From the equations of V(t) and Z(t), we have In terms of boundedness of the solutions, we have known that T + I ≤ T max . For V(t), from the third equation of (1), we have Therefore, the solutions of the system (1) are ultimately bounded for all t ≥ 0.
For the model, the basic reproductive number is defined by Note that R 0 is independent of time delay τ . In fact, R 0 can be biologically divided into two parts. R 1 0 measures the average value of secondary infections caused by an existing free virus, and R 2 0 gives that average number caused by an infected cell, which indicate the effect of virus-to-cell and cell-to-cell transmission on the secondary infected generation respectively.
By the equations of V(t) and Z(t) in model (1), we have V * = k d 2 I * and Z * = c d 3 I * . Substituting them into the equation of I(t) leads to . From the first equation in model (1), we obtain Denote and then there exists a positive I * such that G(I * ) = 1. Therefore, there is one infected steady state E * = (T * , I * , V * , Z * ) when R 0 > 1.

Dynamical results
In the following, we will conduct the dynamic properties of model (1) by establishing stabilities of the three equilibria E 0 , E 1 , E * and discussing conditions for the existence of Hopf bifurcation.  (1),

Stability analysis
Proof: Firstly, linearizing system (1) at E 0 , the characteristic equation is given by which has a positive eigenvalue λ = r, leading to the fact that E 0 is always unstable. The characteristic equation at E 1 is Except for the obvious eigenvalues λ 1 = −r, λ 2 = −d 3 , the remaining roots are determined by the following equation: Noticing 3) has only roots with negative real parts. Then the infection-free equilibrium E 1 is locally stable. If R 0 > 1, Equation (3) has at least one eigenvalue with positive real part, which implies that E 1 is unstable. Further, we claim that if R 0 < 1, lim t→∞ (T(t), I(t), V(t), Z(t)) = (T max , 0, 0, 0). In fact, by the fluctuation lemma in [9], we know that there exists a sequence {t n } with t n → ∞ such that By the second and third equations in model (1), we have which implies that I ∞ = 0 is valid when R 0 < 1. And it is easy to get V ∞ = 0 and Z ∞ = 0 when considering the equations for V(t) and Z(t) in model (1). As a result, by the first equation in (1), we obtain lim t→∞ T(t) = T max . This completes the proof.
To discuss the stability of E * , denote The characteristic equation at E * is which is equivalent to here Noticing that at the endemic equilibrium E * , it is valid that from which we can see that a 3 , a 2 and b i (i = 0, 1, 2) are all positive. Meanwhile, it is obvious that showing that a 0 and a 1 are both positive. Particularly, if τ = 0, the characteristic Equation (4) is reduced to By Routh-Hurwitz criterion, it is known that all solutions of (5) have negative real parts if and only if then we have the following result.
However, assuming (6) being satisfied, in the case of τ > 0, some roots of (4) may cross the imaginary axis to the right side in the complex plane, then it is necessary to analyse the possible complicated dynamic behaviours that will appear at E * when τ increases from 0.

Bifurcation analysis
When τ increases, suppose (4) has a purely imaginary root λ = iω(ω > 0), which is a critical case under small perturbation, then we obtain Separating the real and imaginary parts gives Squaring and adding the above two equalities lead to with If a 0 < b 0 , that is N 0 < 0, it is easy to see that F(ω) = 0 has at least one positive root. For simplicity, let z = ω 2 , then it is necessary to discuss the existence of positive root of equation F(z) = 0. To make use of the results in [13], we take the transformation for z as y = z + N 3 4 , then F (z) = 0 is equivalent to y 3 + n 1 y + n 0 = 0, here n 1 = N 2 2 −  Without loss of generality, assume that there exists four positive roots z * i (i = 1, 2, 3, 4) for equation F(z) = 0, then we have ω i = z * i . By (7), trigonometric functions cos ωτ and sin ωτ can be expressed as from which time delay τ can be represented by ω i as For some i = 1, 2, 3, 4, let τ * = τ i j , then correspondingly ω * = ω i . From the above lemma, we can see that Hopf bifurcation can possibly occur at the endemic equilibrium E * when τ = τ * . We have the following theorem.
which completes the proof.
The critical delay τ = τ * does cause stability switching at E * and the resulting Hopf bifurcation occurs, which is presented numerically in Figure 1. The parameter values in the simulation are listed in Table 1. Two types of functional incidence are taken in the simulations, i.e. bilinear incidence with f 1 = β 1 V, f 2 = β 2 I (the first row) and saturation incidence with f 1 = β 1 V 1+aV , f 2 = β 2 I 1+bI (the second row). With these values, Figure 1(a) plots the solution of model (1) for τ = 0.5, which shows that the endemic equilibrium is asymptotically stable. Conversely, when τ = 1.2, a periodic solution is bifurcated from the endemic equilibrium, as shown in Figure 1(b). For the mentioned case of saturation incidence, take a = b = 0.3, then the dynamic behaviours of model (1) are displayed for τ = 3 in Figure 1(c) and τ = 5 in Figure 1(d) respectively.

Control strategy
In this section, we will dedicate to dealing with the optimal control problem by Pontryagin Maximun Principle with delay in state. The interpretation of model (1) is typical of a disease such as hepatitis B, which is a viral infection that targets the liver. Currently, the therapy for hepatitis B infection focuses on the impairment of viral processes and on the use of immunomodulators, i.e. inhibiting viral production by the antiretroviral nucleside analogues such as tenofovir, entecavir or lamivudine and blocking new infection of the healthy hepatocyte by interferon.

The optimization problem
To evaluate the efficacy of the two treatments, model (1) is rewritten as with u 1 (t) and u 2 (t) denoting the functions of time appearing in the time-dependent optimality system, which refer to the therapeutic effects of reducing the risk of transmission and restraining the replication of the virus. Supposing u i (t) ∈ [0, 1](i = 1, 2) be Lebesgue integrable on a finite time interval [t 0 , t f ], our optimal control problem is to find u 1 (t), u 2 (t) and the associated state variables such that the following given objective functional can be minimized The above integrand consists of the number of infected hepatocytes, free virus load and the cost of two treatments, here B i (i = 1, 2) denoting the weight constant that reflects the value and the importance of control cost. In addition to be more tractable, the quadratic form of costs is used to model the decreasing efficacy of treatment.
In order to investigate the optimal control of system (9), we first prove that there exists solution of the optimality system such that being the control set, By the results in [5,15], we have the following theorem.
Proof: It is necessary to verify that the following conditions are satisfied.
Based on [15], the boundedness of the state system with strategy of treatment ensures the existence of a solution, which implies that condition (i) is valid. For condition (ii), it is obvious by definition that the control set is bounded and convex. Using the boundedness of the solution and the property of f i (i = 1, 2), we know that condition (iii) is satisfied for the form of u 1 and u 2 in the state system. In order to show the concavity of the integrand L on , its Hessian matrix is given as then for any (u 1 , u 2 ) ∈ , the determinant of H L is det(H L ) = B 1 B 2 ≥ 0 which leads to the condition (iv). Finally, for condition (v), when c 1 is taken to be dependent on lower bound of I and V, c 2 = min 1 2 B 1 , 1 2 B 2 , we have Accordingly, the existence of optimal control for the minimization problem has been established.
In the following, by optimal control theory and method in [11,20], the therapeutic strategy is analysed to achieve the desired goal of minimizing objective function in (10). In fact, the necessary condition of this optimal problem can be generated by dealing with the Hamilton function, which is defined as with g = (g 1 , g 2 , g 3 , g 4 ) denoting the right-hand side of system (9) and λ = (λ 1 , λ 2 , λ 3 , λ 4 ) being the associated adjoint variable satisfying The problem of exploring u * = (u * 1 , u * 2 ) that minimizes the objective function (10) subjected to (9) is converted to minimizing the Hamiltonian with respect to the control, then we have the necessary optimality in the following theorem.

Theorem 5.2:
Let u * = (u * 1 , u * 2 ) be an optimal control and (T,Ī,V,Z) be the corresponding solution of (9), then there exists adjoint variableλ = (λ 1 ,λ 2 ,λ 3 ,λ 4 ) satisfying (12) such that Moreover, the optimal control is given as Proof: Computing for the adjoint system (12) leads to By using Pontryagin Maximum Principle, the optimal control u * = (u * 1 , u * 2 ) can be deduced from the optimality condition Figure 2. Numerical solution of model (9) with bilinear incidence, the same parameter values as in Figure 1(b) and 0 ≤ u 1 (t), u 2 (t) ≤ 0.6. from which u = (u 1 , u 2 ) can be solved in terms of state variables and adjoint variables, denoted asū = (ū 1 ,ū 2 ), here (9) and (14), we can obtain the corresponding state variables and adjoint variables, denoted as (T,Ī,V,Z) andλ = (λ 1 ,λ 2 ,λ 3 ,λ 4 ) respectively, then we replace them intoū again, which is still denoted asū for convenience. Meanwhile, we need to consider the constrains on the control and the sign of ∂H/∂u, i.e. component u * 1 of the optimal control u * can be expressed as follows: For u * 2 , it can be conducted in the similar way, which leads to the expression of the optimal control u * = (u * 1 , u * 2 ) given in (13). The proof is completed. Figure 3. Implementation of the optimal control strategy with time.

Simulations of optimality system
In order to perform the numerical simulations, a discretized scheme on the basis of forward and backward finite-difference approximation method is used to solve the optimization system numerically. The bilinear incidence is taken and the parameter values are given the same as in Figure 1(b). We consider the simulation period for 200 days. To compare the numerical results, the simulated procedure is implemented without treatment for the first period of 100 days and with treatment for the second period of 100 days. For the weight coefficients B 1 and B 2 , we take them as B 1 = 1000 and B 2 = 1000, which means that the therapeutic strategy of inhibiting viral production and blocking new infection costs the same. In the following, we mainly consider the implementing of optimal control when there exists periodic solution for model (1). Figure 2 presents the comparison of simulation results in two different cases, i.e. with and without optimal control effect on the bifurcated periodic solution. From the comparison, it can be seen that the application of optimal treatment with 60% efficacy causes significant quantitative change in the model-simulated consequences of state variables, which specifically results in an increased concentration of target cells and decreased concentrations of infected cells and free virus as well as CTL cells correspondingly. As shown in Figure 2, the solution oscillates in both cases, but the magnitude decreases obviously in the case of implementing control. Meanwhile, with the aim of minimizing the objective functional subject, optimal controls u 1 (t) and u 2 (t), corresponding to blocking new infections and inhibiting viral production, must be manipulated properly and economically as time evolves, which is plotted in Figure 3. The periodic control curves demonstrate the therapy administration during the given period of 100 days, corresponding to the optimal control in Figure 3 from the start of treatment.
To reveal the effectiveness of optimal control strategy, 90% efficacy of treatment is also considered, with the same incidence function and parameter values as in Figure 2, which is shown in Figure 4. From the simulation, it can be noted that the application of treatment results in qualitative change in the system dynamics, comparing with that when 60% efficacy is adopted in Figure 2. As desired, the concentrations of infected cells and free virus are reduced to zero, and the system approaches to a stable disease-free steady state, implying the successful controlling of infection. This figure suggests that the disease clearance can be achieved if the therapeutic efficacy is improved to some high value and (u 1 (t), u 2 (t)) is administrated according to the optimal control solution. However, it is still currently challenging since the available drugs do not directly target cccDNA in chronic hepatitis B and the long-term treatment can lead to the problem of drug resistance.

Conclusion
In this paper, we investigate an in-host DDE model with immune response by incorporating the lag between infection and activation of CTLs. First, dynamical properties of our system are studied analytically. In particular, stability analysis of the steady states has shown that the trivial equilibrium with no biological implication is always unstable, and the endemic equilibrium becomes feasible and probably stable once the infection-free equilibrium loses its stability. When constant time delay is further taken as the bifurcation parameter, the existence of pure imaginary roots of the characteristic equation is verified and the endemic steady state can lose its stability via Hopf bifurcation, resulting in periodic solution, which is presented intuitively in Figure 1.
Furthermore, considering the drug effects in blocking viral production and reducing viral infection for hepatitis B infection, we design an optimal control problem based on the DDE model. Allowing the efficacy of drug combination to vary with time, two control variables needed for optimal therapy are introduced in the system, with one acting as reducing the appearance of new infections in the incidence term, the other one inhibiting viral production in the releasing term for free virus particles. By the optimal control theory and method, we derive the control formulation and use it to search for the best control strategy, i.e. achieving the least amount of liver damage and the lowest cost for the drug combination that lead to hepatitis B virus control. When the bifurcated periodic solution occurs, the comparisons of model-predicted consequences of the state variables are presented numerically with different efficacy (60% and 90%), from which the quantitative effects of the control show that the sustained infection can be reduced to a lower level or even to viral clearance, if the control strategy is implemented correspondingly and the rapid progress in the therapy can be achieved.
However, it is known that long-term antiviral therapy results in side effects, compliance and, most importantly, rapid emergence of drug resistance. This has made it difficult to design an optimal control for treatment start, duration, which type of drugs to use and how to assess the effectiveness of treatment measures. Therefore, the limitation in our research can be considered and more accurate models can be developed to characterize the dynamic process and address the optimal control strategy in the presence of these events.