A FRACTIONAL-ORDER DELAY DIFFERENTIAL MODEL WITH OPTIMAL CONTROL FOR CANCER TREATMENT BASED ON SYNERGY BETWEEN ANTI-ANGIOGENIC AND IMMUNE CELL THERAPIES

. In this paper, we present an optimal control problem of fractional-order delay-diﬀerential model for cancer treatment based on the synergy be- tween anti-angiogenic and immune cells therapies. The governed model consists of eighteen diﬀerential equations. A discrete time-delay is incorporated to rep- resent the time required for the immune system to interact with the cancer cells, and fractional-order derivative is considered to reﬂect the memory and hereditary properties in the process. Two control variables for immunotherapy and anti-angiogenic therapy are considered to reduce the load of cancer cells. Necessary conditions that guarantee the existence and the uniqueness of the solution for the control problem have been considered. We approximate numer- ically the solution of the optimal control problem by solving the state system forward and adjoint system backward in time. Some numerical simulations are provided to validate the theoretical results.

1. Introduction. Nowadays, mathematical models can be considered as a successfully powerful tool to test hypotheses, confirm experiments, and simulate the dynamics of complex systems. Moreover, it has become an essential tool that is used to simulate tumor-immune cell interactions and to predict the therapeutic effect of various cancer treatments ( [2], [13], [15], [19], [23], [30], [39], [45]). Mathematical models, based on ordinary differential equations, delay differential equations, partial differential equations, have proven to be useful tools in analysing and understanding the interactions with viral, bacterial infections and cancerous cells [31]- [33].
It is well-known that Cancer can be considered as one of the biggest killer of humans over all the world ( [17], [18], [25] ). In spite of recent advancements in cancer treatment, the battle against this deadly disease still rages on. Immunotherapy (which is sometimes referred to biological therapy) is quickly becoming one of the most important components of cancer treatments, especially in multi-pronged approaches [26]. The goal of immunotherapy is to reinforce the body's own natural ability to combat cancer by enhancing the effectiveness of the immune system to act against cancer cells, which involves the use of cytokines with Adoptive Cellular Immunotherapy (ACI), derived from the body or laboratory-produced versions of such substances, to improve or restore immune system function [20,22]. Traditional treatment modalities include chemotherapy, radiotherapy and surgery.
Clearly, interactions between cancer (tumour) cells and immune system are very complex and need sophisticated models to describe such interactions. Several mathematical models have been proposed to describe the dynamics of tumour and immune system over time. Study of the effects of time-delays in the growth of tumors by using the methods of mathematical models was initiated by [6] in which the author introduced two time-delays describing two processes: proliferation and regulatory apoptosis. The authors in ( [28], [12], [35]) focus on delayed-induced oscillations in cancer-immune system dynamics to explain some periodicity observed during the progress of the disease. All of that efforts creates the basis for the more complex considerations focusing on at least control of cancer progression. In [38], the authors suggested a system of eighteen ordinary differential equations to predict synergy when anti-angiogenic treatment is administered first to reduce immuno-suppression, and is immediately followed by immunotherapy to boost and enhance the tumor-targeted immune response.
The main goal of this paper is to provide an optimal control problem governed by a system of fractional-order differential equations with time-delays for cancer treatment based on synergy between anti-angiogenic and immune cell therapies. Two control variables for immunotherapy and anti-angiogenic therapy are incorporated to eliminate or reduce the load of cancer cells. We use iterative optimal control method (IOCM) and generalized Euler method (GEM) to numerically solve the optimal control problem. Comparative studies are implemented. Up to the best of our knowledge, optimal control problem for the fractional mathematical model with time-lags of cancer treatment based on synergy between anti-angiogenic and immune cell therapies with the modified parameters was never explored before. This paper is organized as follows: In Section 2, the proposed model of fractionalorder derivatives with time-delays and two control variables is presented. Formulation of the optimal control problem and the necessary optimality conditions are derived in Section 3. Numerical simulations are provided in Section 4. In section 5, some conclusions are presented.
1.1. Preliminaries. Firstly, we provide some useful definitions and properties of fractional-order. Definition 1.1. [29] The left Caputo fractional derivative (LCFD) is defined as follows: the right Caputo fractional derivative (RCFD) is defined as follows : The Caputo operator is given by the Grünwald-Letnikov approach as follows: where, N n be a natural number and the coordinate of each mesh point is: , i = 1, 2, ..., n + 1.
Assume that [37]: 2. Fractional cancer treatment model with time-delay. The interactions between the immune system and tumor cells are highly complex, and understanding the interactions between a growing tumor and immune system is fundamental for optimizing current treatments and proposing suitable strategy. In this section, we present a delayed fractional-order differential model that represent the interactions of immune system with cancer cells in presence of treatment based on synergy between anti-angiogenic and immune cell therapies. The mathematical model presented in this paper examines the T-cell response to tumor cells. This model is an extension of the model of integer-order which is given in [38]. This mathematical model was designed to integrate tumor angiogenesis and tumor-targeted cytotoxicity by immune cells to identify the therapeutic window of two distinct modes to treat cancer: (i) an anti-angiogenesis treatment based on the monoclonal antibody bevacizumab that targets tumor vasculature, and (ii) immu-notherapy involving the injection of unlicensed dendritic cells to boost the anti-tumor adaptive response. The angiogenic cytokine Vascular Endothelial Growth Factor (VEGF) contributes to the immunosuppressive tumor microenvironment, which is responsible for the short-lived therapeutic effect of cancer-targeted immunotherapy. The model was used to determine treatment protocols involving the combination of anti-VEGF and unlicensed dendritic cell injections that enhance tumor regression. The model simulations predicted that the most effective method to treat tumors involves administering a series of biweekly anti-VEGF injections to disrupt angiogenic processes and limit tumor growth. The literature reveals that most mathematical models of tumor-immune interactions are based either on ordinary differential equations (ODEs) or delay differential equations (DDEs) with integer-order. Due to the fact that both time-delays and fractional-order play a vital role in biological systems with memory which gives more degree of freedom [3], [11], [16]. Herein, we develop the model of [38] to include fractional-order 0 < α ≤ 1 and time-lags. Moreover, the parameters of modified model will also depend on the fractional-order to unify the units of the differential equations. The model consists of eighteen time dependent variables. However, if the model is too simple, the results will be trivial and will fail to produce the complex dynamics observed in experimental and clinical settings. A discrete timedelay is incorporated to represent the time required for immune system to interact with tumor. Two control variables u A and u M are introduced to measure the antiangiogenic therapy and the immunotherapy, respectively. The model presented here uses18 biological variable, defined as follows: The developed system takes the form: ), where, the parameters ω α 1 and ω α 2 are the weight factors. It should be mention that the left sides of all equations of the system have dimension time −α . To make the system more consistent the reality, we must therefore make sure that the right-hand sides of these equations have the same dimensions. Therefore, we need to modify the right-hand sides to make the dimensions match. The most straightforward way of doing this is to put power α of each parameter in the right sides so that when α → 1 the system reduces to classical one.
The described of parameters given in Table1, ([7], [34], [38] ). Equation (2), describes the evalution of tumor cells with time-delay T (t−d τ ). The first term of the right hand side of equation (2), accounts for delay of tumor growth. Equations (3) and (4) describe the mature dendrtic cells. The first term of equation (3) describes the delay of the maturation of dendritic cells by an encounter with tumor antigen with fractional memory. The second term of equation (3) and first term of equation (4) account for the licensing of dendritic cells upon encounter with helper cells. The first term of equation (5), activates the effector memory cells, in the presence ofmature licensed dendritic cells. The second term describes the expiration of the activated phase of these cells. The first term of equation (6), gives the proliferation of CD8 + effector T cells. This depends on the presence of activated cells and the proliferative cytokine IL − 2(C). Equations (7) and (8), differ from those for the effector T cells in two ways. First, the activation of memory helper cells can be prompted by both licensed and unlicensed dendritic cells. Second, helper T cells can be converted in to Tregs in the presence of T GF − b, and this is represented by the second term in equation (8). The main difference in equations (9 ) and (10), is that Treg prolife ration is not suppressed by T GF −β, incontrast to the other T -cell populations. The first term of equations (11 ) represents the production of IL−2 by the activated helper cells. The first terms of each equations (12) and (13), represent the production by Tregs, the second terms represent the delay production by tumor cells, and the final terms represent the removal from the system with characteristic times τ α s and τ α . Equation (14), describes the behavior of Angiopoietin-1. Ang-1 is produced by vasculature at a constant rate, α A1 and degrades at a rate δ α A1 relative to the concentration of angiopoietin-1 in the system. Equation (15), describes the behavior of Angiopoietin-2. Angiopoietin-2 is produced by vasculature at a maximum rate α α A2 in response to the presence of cancer cells. It degrades at a rate δ α A2 relative to the concentration of Ang-2 in the system. Equation (16), describes the behavior of vascular endothelial growth factor. Vascular endothelial growth factor is produced by cancer cells both at a constant rate, α α v , and at a hypoxic rate, α α v2 , when there is not enough vasculature to support the growth of tumor cells (low V EGF per endothelial cell). V EGF degrades at a rate δ α v relative to the concentration of V EGF in the system. V EGF is also removed at a rate τ α when it binds with the anti-V EGF antibody. Equation (17), describes the behavior of endothelial cells. Endothelial cells proliferate at a maximum rate α α y in response to the amount of V EGF per endothelial cell. Endothelial cells mature into vasculature at a maximum rate ω α when there is a high enough Ang−1 Ang−2 ratio to mature vasculature, and when there is enough V EGF per endothelial cell to cause angiogenic sprouting. Endothelial cells die at a maximum rate δ α y when there is not enough V EGF per endothelial cell to continue survival. Equation (18), describes the behavior of mature blood vessels. Equation (19), describes the behavior of an anti-V EGF treatment. For more detiales see ( [7], [38]).
3. Formulation of the fractional optimal control problem. Consider the variables in the system (2)- (19), belongs to R 18 , let is the set of admissible control functions: The cost functional can be defined as follows: where, A, B 1 , and C 1 represent the weight constants of tumor cell numbers for immunotherapy and anti-angiogenic therapy, respectively. The goal now is to minimize the following cost functional: subjected to the constraints Table 1. The parameters of system (2)- (19) and their descriptions.
Parameters of fraction power Descriptions γ α The tumor growth rate a α The antigenicity for tumor The parameters for dendritic cell expansion α α Parameters of IL-2, TGF-β and IL-10 concentrations . . 18 represent the right hand side of the state system (2)- (19) which satisfying suitable initial initial conditions. According to [40] the modified cost functional takes the form where, H a is the Hamiltonian functional of the form The necessary conditions as in [1] that can be derived from (22) and (23) as follows: where, and it is also required that: Where λ i , i = 1, 2, 3, ..., 18, are the Lagrange multipliers. , 4. Numerical simulations. The fractional optimality system (48)-(65) and (27)- (42) with the transversality conditions λ i (T f ) = 0, i = 1, ..., 18, will be studied here using the methods IOCM and GEM for more details on these method see [42]. We consider that, the initial state of variable are T  Table 2 where, 0 < α ≤ 1. The approximate solutions of the proposed system are given in We noted that, the result of B(t), R(t), Y (t), C(t) and V a (t) are not affected by changing the value of d τ . Fig. 3, shows the behavior of the approximate solutions of the optimal controls u * A ,u * M , and different values of d τ , we noted that, the result of u * M is not affected by changing the value of d τ at α = 0.96. Fig. 4, shows the relation between the state variables T (t), E(t) and Y (t) when d τ = 3 and α = 0.98 using IOCM with two controls treatment cases. Fig. 5, shows the numerical simulations of the control variables with different values of α and d τ = 1, using IOCM. Fig. 6, shows the relationship between the variables T (t−d τ ), I(t−d τ ), V (t−d τ ) and R(t−d τ ) and T (t), I(t), V (t) and R(t) values of α = 0.92 and d τ = 2. Table 3, shows comperative studies between the value of objective function using IOCM with and without controls cases when T f = 100, d τ = 2 and different values of α. We noted that the values of objective functional with control are less than the values of objective functional without the controlled case for all value of 0.3 < α ≤ 1 by using IOCM, i. e., for small values of the α, the tumor already doesn't take a completed shape so the value of objective functional is expected to be small as well. Table 4, shows the values of objective functional J(u * A , u * M ), T f = 50 with different value of d τ and control case. We noted that, the values of objective functional increases by increasing the value of the d τ . Table 5, shows the comperative studies between the two proposed methods: IOCM and GEM to find the optimal values of the objective functional in controlled case at different values of α. We can conclude that IOCM is better than GEM, because, GEM can be failed to derive the optimal values for 0 < α ≤ 0.8. Additionally, in Table 5, increase in the value of objective functional indicates that as the time for tumor immune system interaction increases, control of the tumor becomes harder. All results were obtained by using MATLAB (R2013a), on a computer machine with intel(R) core i3-3110M @ 2.40GHz and 4GB RAM.

5.
Conclusions. In this paper, we provided an optimal control problem for tumorimmune interactions in presence of treatments. The mathematical evolution model is governed by a system of fractional-order differential equations with time-delays. The cancer treatment is based on synergy between anti-angiogenic and immune cell therapies. A discrete time-delay is incorporated to the model to represent the time required for immune system to interact with tumor. Two optimal control variables u A and u M are introduced to measure the anti-angiogenic therapy and the immunotherapy, respectively. Necessary optimality conditions have been derived. Iterative optimal control method and generalized Euler scheme have been used effectively to study fractional optimal control problems with delay. Numerical simulations support the theoretical findings and show the positive impact of the fractional-order and time-delays in the model. A combination of fractional-order and time-delay in the governing model enriches the dynamics and complexity of the model. The numerical simulations show also that combination of immunotherapy and anti-angiogenic therapy can lead to eliminate or reduce the load of cancer cells.
Sensitive analysis to investigate the impact of small changes in the model parameters in the stability of the tumor-free steady state can be investigated in future study.