Optimal Control of a Cell-to-Cell Fractional-Order Model with Periodic Immune Response for HCV

: In this paper, a Caputo fractional-order HCV Periodic immune response model with saturation incidence, cell-to-cell and drug control was proposed. We derive two different basic reproductive numbers and their relation with infection-free equilibrium and the immune-exhausted equilibrium. Moreover, there exists some symmetry in the relationship between the two equilibria and the basic reproduction numbers. We obtain the global stability of the infection-free equilibrium by using Lyapunov function and the local stability of the immune-exhausted equilibrium. The optimal control problem is also considered and two control strategies are given; one is for ITX5061 monotherapy, the other is for ITX5061 and DAAs combination therapy. Matlab numerical simulation shows that combination therapy has lower objective function value; therefore, it is worth trying to use combination therapy to treat HCV infection.


Introduction
Hepatitis C virus (HCV), which is an enveloped flavivirus [1], causes both acute and chronic infection. It is reported that about 290,000 people have died of HCV infection so far; most of them from liver cirrhosis and hepatocellular carcinoma [2]. Hence, understanding the mechanism of HCV virus infection is very important for the treatment of HCV.
Mathematical models are widely used to study the dynamic behavior of HCV. At the same time, the development of mathematical theory also makes it possible for us to study the treatment of HCV. In 1998, Neumann et al. [3] proposed a simple integer-order HCV model with antiviral treatment. Furthermore, in 2005, Dahari et al. [4] extended the model in paper [3] to predict long-term response rates to interferon monotherapy and combination therapy. However, the development of fractional order generalizes the integer order. The reason is that, compared with integer order, it is non-local, so it is more suitable for dynamical systems with memory influences on their state variables [5], such as the immune response. So, some fractional-order models related to viruses were proposed by [6][7][8][9][10][11]. However, these models didn't consider immune system.
With the development of molecular techniques, there have been fundamental insights into the molecular mechanisms of the immune system for HCV infection [12]. In 1996, Nowak et al. [13] presented a model to explore the relationship between antiviral immune responses, virus load, and virus diversity. In 2003, Wodarz [14] added antibody responses to investigate the role of cytotoxic T lymphocytes (CTLs) and antibody responses in HCV dynamics and pathology. However, the human immune system has high complexity. The immune would be influenced by taking some drugs regularly. Papers [15][16][17] considered periodic lytic immune response, whose main purposes were to study the effect of oscillation of the immune system on the viral dynamic behaviors. The model in paper [17] is presented as follows: where x, y, z denote uninfected cells, infected cells, CTLs respectively. Uninfected cells are assumed to be produced at the constant rate λ and become infected by infected cells at the rate of βxy. Infected cells are thus producted at the rate of βxy and become lysed by CTLs at a rate p(t)yz. p(t) stands for the strength of the lytic component. T(days) presents the period of oscillation of human immune system. β 0 is the basic strength of the lytic component, the amplitude parameter β 1 is the degree of oscillation of immune system, and ϕ is the acrophase. CTLs are produced at a rate cyz. The mortality rates of uninfected cells, infected cells and CTLs are d 1 x, ay and bz respectively. Except for p(t), which varies with time, all parameters of system (1) are constants. However, [17] did not consider viruses. It is known that cellular and humoral responses are of limited effectiveness in achieving viral responses in most patients with acute infection HCV [18]. It is generally believed that, due to factors such as delays in the life cycle of the virus and the evasion of immune responses, the cell-to-cell transmission of viral particles is far greater than a cell-free spread. In 2014, Fei Xiao et al. [19] showed that cell-cell transmission is the route of viral spread. Papers [20][21][22] proposed some cell-to-cell models, and in 2018, a mathematical model of HCV infection with cell-to-cell and cure rate was proposed by Sonjoy et al. [23]. This model is given by where x denotes uninfected hepatocytes, y represents infected hepatocytes, v and z denote free virus and CTLs respectively. The parameters of system (2) are all constants. αy denotes the cure rate of infected cells. Uninfected hepatocytes are produced at the constant rate λ, become infected hepatocytes at the rate of β 2 xy, and are infected by free virus at the rate of β 1 xv. Infected hepatocytes are produced at the rate of β 1 xv + β 2 xy. Free viruses are produced by infected hepatocytes at a rate ky and become lysed by CTLs at a rate pvz. CTLs are producted at the rate of cvz. The mortality rate of x, y, v and z are d 1 x, ay, γv and bz, respectively. However, this model is not fractional-order and p is a constant. Meanwhile, this system did not consider the control of drugs. After that, some authors extended the system (2). In 2020, a cell-to-cell model with humoral immunity is considered by Mausumi et al. [24]. After this, paper [25] further considered the time delay of humoral immunity. In addition, a stochastic HCV infection model with therapy was discussed by paper [26]. However, paper [24][25][26] did not consider fractional order.
It should be noted that the model (2) described the process of free virus infecting normal cells using bilinear incidence βxv, and so did paper [17,[20][21][22][23]. Paper [27] uses Holling type-II incidence instead of bilinear incidence, that is βxv 1+mx . Yang et al. [28] changed it to saturation incidence function βxv x+v , and the reproductive number is R 0 = βk µ(a+δ) , which seems more reasonable because it is independent of the number of total cells of liver. Since infected liver cells can also infect normal liver cells, it might also be better to replace βxy with βxy x+y .
HCV can be controlled in many ways. In 2016, Ian A et al. [29] proposed that ITX5061 blocks HCV entry and infection in vitro. The recent introduction of new direct-acting antivirals (DAAs), targeting HCV replication, allowed us to achieve sustained virological response (SVR) rates in over 90% of treated patients [30]. In 2021 [31], efficacy and safety of pangenotypic DAAs against hepatitis C virus infection in Taiwan are proven.
Based on the above discussion, in order to better explain the phenomenon of HCV clinical treatment, we consider an HCV cell-to-cell fractional-order system with the saturation incidence rate. The model is as follows: where D α denotes the Caputo fractional derivative. The variable x, y, v and z and the parameter d 1 , k, a, β 1 , β 2 , γ and r have the same meaning as that in model (2). and η denote the efficacy of the drugs, namely, ITX5061 and DAAs. m 0 and m 1 have the same meaning as β 0 and β 1 , respectively, in model (1). Infected hepatocytes are lysed by CTLs at a rate p(t)yz. CTLs are producted by Infected hepatocytes at the rate of qyz. This paper consists of the following seven parts: in part 2, some related definitions and lemmas were presented. The existence and uniqueness of positive solution of system (3) were proven in part 3. The stability of the equilibrium point of the system (3) was analyzed in part 4. The fifth part presents the theoretical analysis of the optimal control of the model. Numerical simulation and discussion are carried out in part 6, and in part 7, we give the conclusion.

Preliminaries
In this section, we give some basic definitions and lemmas of the Caputo fractionalorder. We denote by AC n [a, b] the space of complex-valued functions f (x) which have continuous derivatives up to order Definition 1 (see [32]). Let α > 0 and n = [α] + 1. If y(x) ∈ AC n [a, b], then the Caputo fractional derivative (D α y)(x) exists almost everywhere on [a, b], which is defined by Lemma 1 (see [33]). Assume that the vector function f : R + × R n → R n satisfies the following conditions: is Lebesgue measurable with respect to t on R + ; (ii) Function f (t, X(t)) is continuous with respect to t on R n ; (iii) (∂ f (t, X(t))/∂X) is continuous with respect to X(t) on R n ; (iv) f (t, X) ≤ ω + λ X , ∀t ∈ R + and X ∈ R n . Here, ω and λ are two positive constants; then, initial value problems then, there exists a unique solution.

The Existence and Uniqueness of Positive Solutions
Firstly, we analyze the system (3) without control; that is, = η = 0. (3) is always non-negative.

Theorem 1. The solution of system
Proof of Theorem 1. For system (3), we can obtain Based on Lemma 2, we can prove that the solution of the system (3) will remain nonnegative for all t ≥ 0.
here µ = min d 1 , a 2 , γ, r , so we obtain it is easy to see that D is a positively invariant region for system (3).
Theorem 2. For any initial condition in D, system (3) has a unique solution X = (x, y, v, z) .
Proof of Theorem 2. We denote Obviously f (t, X(t)) satisfies conditions (i), (ii) and (iii) of Lemma 1, and we only prove that system (3) satisfies the last condition (iv) of Lemma 1.
In the case of HCV infection, let the whole hepatocytes population T(t) = x(t) + y(t), then add the first two formulas of the system (3), and we get from the biological significance, the death of infected hepatocytes is greater than the natural death of uninfected hepatocytes, so a ≥ d, then then, we have Therefore, we obtain where Vector field function f (t, X(t)) satisfies Lemma 1, so the solution of the system (3) exists and is unique.

Proof of Theorem 3. Define a global Lyapunov function
Then, we have The solution is So, we obtain by the positive invariance of E, we can learn E = {E 0 }, and by the Lvapunov-Lasalle theorem [36], all solutions in the set D are attracted to E 0 . Noting that E 0 is locally asymptotically stable, therefore E 0 is global asymptotically stable.

Optimal Control Problem
In this section, we consider the optimal control problem of system (3). Here, and η denote the effects of the drugs ITX5061 and DAAs, respectively. Suppose the curative effect is positive proportional to the dosage. We will give the optimal control strategy at the minimum expense within a supposed time. The objective function of the optimal control was given as follows [38], where t f is the endpoint of treatment.: Here a 2 , a 3 , b 2 , b 3 , c 1 , c 2 represent the corresponding weight of each variable. Let X(t) = (x(t), y(t), v(t), z(t)) , u = ( , η) , X(0) = (x 0 , y 0 , v 0 , z 0 ) , then we can obtain Then, the Hamilton function of J(u) was given as follows: According to Pontryagin's minimum principle, the necessary conditions for the existence of the optimal solution of system (3) are shown as follows: (iv) So, we obtain the optimal control as:

Simulation and Discussion
In this part, we will present some numerical simulations. First of all, our model is not based on actual data, but only qualitative discussion, so the parameters can be taken within the range of biological significance. Secondly, the parameters in the literature that we refer to are also given on the basis of biological significance. Moreover, the parameter values were presented in Table 1. We give the simulation with different orders (α = 0.95, 0.9, 0.8), and the initial conditions are x 0 = 66; y 0 = 1000; v 0 = 2000; z 0 = 30. Let the weight a 2 = 5000; a 3 = 0.01; b 2 = 0.1; b 3 = 0.1; c 1 = 0.01; c 2 = 0.1. T Period of oscillation of human immune system 1 [17] ϕ acrophase π/12 [17] • Strategy A: monotherapy of control (i.e., = 0 and η = 0); In this section, we simulated the optimal procedure without considering DAAs drug, that is η = 0. We assume that the maximum and minimum treatment effects were 0.98 and 0.5, respectively. The simulation results are shown in Figures 1 and 2. Figure 1 gives the optimal control strategy with different values of α; we can see that the maximun dosage only should be used at the former stage, and then after a period of time, taking the dosage should be gradually reduced to a minimum. Obviously, the time points of minimum dosage are different with different α.
As for the therapy effect, we can see from Figure 2 that with the same initial conditions, the change rates and terminal values of infected hepatocytes, viruses and CTLs vary with α, which reflect individual differences. Therefore, it is worth trying to use fractional order model to describe the treatment of HCV infection.  Table 2 shows the objective functions values of optimal therapy, maximum effect and minimum effect. Obviously, compared with the objective function values of maximum effect and minimum effect, we find that the objective functions values of optimal therapy are lowest for different α values.

•
Strategy B: combination therapy of control (i.e., = 0 and η = 0); Next, we will present the optimal strategy for combination therapy of ITX5061 and DAAs, which means = 0, η = 0. The numerical simulation results are shown in Figures 3 and 4.
We can see from Figure 3 that the maximum dose of ITX5061 and DAAs only needs to be used at the former stage. However, the time points of lowering the drug dose are different with different medicine when α is the same. In addition, the time when DAAs reach the lowest point first is related to the value of α, and so is ITX5061.
We also compared the objective function values of constant control = η = 0.98 and = η = 0.5 with optimal control. The results are shown in Table 3. We also find from Table 3 that the value of optimal control is lower than maximum effect and minimum effect when α is the same.   We compare the terminal value of infected cells and viruses between strategy A and B in Tables 4 and 5. The results show that the terminal values of infected cells and viruses in strategy B are lower than that in strategy A. At the same time, by comparing Tables 2 and 3, we can see that the objective function values of strategy A are higher than for strategy B. Therefore, the combination of ITX5061 and DAAs is more suitable for treating HCV. From the simulation of strategy A and strategy B, we can also see that even if the treatment effect can be achieved in the end, the therapy strategy will vary with the order α, which also reflects the clinical personalized differences.

Conclusions
In this paper, we discussed a Caputo fractional-order HCV periodic immune response model with saturation incidence, cell-to-cell, and drug control. We proved the existence and uniqueness of positive solutions. We derived two different equilibrium E 0 and E 1 , and two different basic reproductive numbers R 0 and R 1 . And then, we dissused their relations between equilibrium and basic reproductive numbers. The result shows that R 0 > 1 is equivalent to the existence of E 1 . When R 0 < 1, the equilibrium point E 1 is globally asymptotically stable, and when R 0 > 1 and R 1 < 1, the immune-exhausted equilibrium is locally asymptotically stable. We condsidered the optimal control problem, and then gave two optimal control strategies. We compared the objective function of two strategies by simulating. The result shows that the combination of two drugs is better than any single drug in treating those infected with HCV. Even if patients take the same dose of drugs, different people will have different curative effects. The order of exactly fractional order can reflect individual differences in clinical treatment. Therefore, it is worth being used to study clinical phenomena.
On the other hand, we give an idealized result based on the assumption that the dosage is proportional to the effect, but unfortunately, we lack clinical experimental data.