Review of Optimization Methods for Cancer Chemotherapy Treatment Planning

Cancer is one of the major diseases that limited the human life; it is treated with surgery, radiation, chemotherapy, hormones, and immunotherapy. According to estimates from the International Agency for Research on Cancer (IARC), there were 12.7 million new cancer cases in 2008 worldwide. The corresponding estimates for total cancer deaths in 2008 were 7.6 million (about 21,000 cancer deaths a day). By 2030, the global burden is expected to grow to 21.4 million new cancer cases and 13.2 million ancer deaths simply due to the growth and aging of the population, as well as reductions in childhood mortality and deaths from infectious diseases in developing countries [1].


Introduction
Cancer is one of the major diseases that limited the human life; it is treated with surgery, radiation, chemotherapy, hormones, and immunotherapy. According to estimates from the International Agency for Research on Cancer (IARC), there were 12.7 million new cancer cases in 2008 worldwide. The corresponding estimates for total cancer deaths in 2008 were 7.6 million (about 21,000 cancer deaths a day). By 2030, the global burden is expected to grow to 21.4 million new cancer cases and 13.2 million ancer deaths simply due to the growth and aging of the population, as well as reductions in childhood mortality and deaths from infectious diseases in developing countries [1].
The development work in cancer prevention, detection, treatment, and management is recently very advanced. Mathematical modeling of cancer is one of the important methodology that have contributed in this domain. With their basis in clinical studies, multiple discussions, revisions and evaluations, mathematical models help researchers understand the effects that various factors, such as tumor growth and drug infusion rates, have on optimal treatment plans. The models are used to spotlight the importance of looking at cancer treatment as a formal optimization problem. The number of cancerous cells, toxicity, and drug resistance are the key factors in planning a chemotherapy treatment. Most of the studies examining chemotherapy treatment optimizations use various ways to model the interactions among these key factors. Certain authors consider the minimal number of cancer cells at the end of therapy to be a sufficient indicator of treatment quality (e.g., Costa and Boldrini [2]), but the negative cumulative effect of the administered drug(s) should also be explicitly included in the performance index.
The mathematical modeling of the various phases of solid-tumor growth has itself been developing and expanding over the years. However, over the last 20 years or so, models of cancer invasion have begun to appear in the research literature [3][4][5][6][7]. In [8], a hybrid, multiscale cellular-automaton model was presented with the goal of studying the effects of various combinations of cell-cycle specific chemotherapy drugs in the presence of internal and external heterogeneity.
Actually, surgery and radiation therapy are the most common direct therapies for curing solid tumors [9]. However, when cancer reaches the metastasis stage, when cancerous cells from the primary tumor are transmitted to other parts of the body, a systemic treatment such as chemotherapy must be applied to the diffused cancerous cells. Surgery is performed when a tumor is found in only one area, and it is likely that all of the tumor can be removed; occasionally radiation therapy is also used during or after an operation. When the tumor is entirely excised, we do not need to study the treatment optimization problem. However, when chemotherapy is used as treatment, there is an evolution of the number of cancer cells with time, as shown in Figure  1. Because chemotherapy works like a two-sided sword, destroying the normal cells while annihilating the cancerous cells, we cannot inject the drug frequently. Therefore, the administration of chemotherapy involves a tradeoff between cancerous cell reduction and tissue toxicity, which is a function of the size of the destroyed normal cell population and the drug dosage limits. For this reason, we are interested in optimizing cancer chemotherapy.
Randomized clinical trials are the standard method for the evaluation of chemotherapy treatment plans. For example, in reference [10], Andre et al. conduct a trial to compare the effectiveness of twodrug chemotherapy with various dosages, frequencies and treatment durations. Apart from the dosages and time duration, the choice of drug combination is another essential factor considered in clinical trials [11].
The choice of treatment depends on the type of cancer, its stage and its grade. The oncologist will also consider the overall health of the patient and their medical history, personal preferences and other relevant factors. Typically, chemotherapy is administered into the veins through an injection-this mode of injecting the drug into the vein is known as intravenous chemotherapy. Chemotherapy can also be administered orally with the help of tablets and capsules; this method is known as oral chemotherapy. Chemotherapy drugs can be injected directly into the muscles, known as intramuscular injection, or beneath the skin's surface, called subcutaneous injection. In intravenous chemotherapy, the drug reaches the blood stream directly [12].
Although clinical trials have been used to determine the most reliable and efficient chemotherapy treatment plans, they are limited by high costs, long trial times, and the difficulty of having to test multiple options. In addition to the treatment cost and the effectiveness of the chemotherapy plan, its feasibility should also be evaluated. All of these steps multiply the cost; for this reason, we are studying chemotherapy treatment planning as an optimization problem using a mathematical model. As several researchers note, further refinement of chemotherapy will require attention to rigorously derived models because clinical empiricism can be an inefficient method of understanding and developing a treatment strategy [13,14].
To this end, Mathematical modeling provides a low-cost method to evaluate different treatment strategies more efficiently, quantifying the relationships among several important factors, such as the population of cancerous cells, toxicity, and drug resistance. Mathematical models also aid researchers in understanding the effects of other variables, such as the tumor growth and drug infusion rates, on the performance of the optimal treatment plan. Therefore, there is a growing interest among researchers on the problem of chemotherapy treatment optimization. A close collaboration with an oncologist will improve the model, making the research more practical and meaningful.
In general, many researchers develop mathematical models to simulate the pharmacokinetic and pharmacodynamic processes. Pharmacokinetic processes describe the distribution and metabolism of the drug, and pharmacodynamic processes characterize the effects of a drug on cancerous cells and normal cells [15]. The researchers then apply computational methods to solve these models and compute the optimal chemotherapy plan, which generally specifies the combination, frequency, and dose of drug administration.
In real cases, the oncologist uses for a given cancer a standard chemotherapy treatment defined from the results of many clinical trials [16,17]. For example, for a patient with breast cancer, we have many standard protocol treatments. The oncologist uses the predetermined treatment dosage.
Most cancer modeling studies consider treatment in the form of a generic efficacy term that represents the effectiveness of the drug. Various treatment types have been used, such as chemotherapy as a means to reduce cancer cells [18][19][20][21][22][23] in an attempt to simulate clinical practice as closely as possible. In [18], the authors demonstrate that, in the chemotherapy treatment, the effect of the drugs was shown to reduce not only cancer cells but also normal and immune cells, as is the case in reality.
The concept of applying optimal control to various disease states began by the mid-1970s, and since then, it has become the subject of various publications [24][25][26][27]. In [28], engineering optimal control theory is applied to investigate the drug regimen for reducing an experimental tumor cell population. In [29], several optimal-control problems resulting from the simplest models of cancer chemotherapy leading to singular control solutions was proven.
Swan [30] presents a study that has used engineering optimal control theory for a chemotherapy problem. It involves a human tumor and minimizes the total amount of used drug for a specified value of tumor cell population. The study by Swan is critical for the basic understanding and comprehension of the early mathematical modeling approaches on chemotherapy treatment planning problem. The first published review of optimal control problems in the general area of cancer research appeared in [31]. In later papers, for example [32], Swan provides evidence for the use of continuous delivery of drugs. An excellent general reference for this whole topic is [26]. But in [33], Zietz and Nicolini attempted a compromise between toxicity and cell kill by using an objective function that is a combination of tumor cell final and normal population. As a different methodology, like in [34] and [35], the application of drugs is matched up with the progression of the cells through the cell cycle. In [34] the optimal period for drug application corresponds approximately to the normal cell cycle time.
In the literature, one can find several modeling approaches related to cancer chemotherapy problems and to optimal drug regimens. Several of these are mentioned in [36], and we can find many different approaches and models for the job of reducing the tumor burden while keeping the side effects of the drugs sufficiently in check.
References in recent years, such as Jinghua Shi et al. [37] and Nanda et al. [38], represent various optimal control models of cancer chemotherapy. Most of these references concern the treatment of solid tumors or cancer treatment in general.
In [39], the authors review many topics about i) the modelization of the cell cycle as an object control under the action of anticancer chemotherapy, which includes modeling of drug resistance; ii) the modelization of tumor angiogenesis and antiangiogenic therapy; and iii) the modelization of therapies targeting specific cellular regulatory networks, involved in carcinogenesis and cancer growth and progression. All of these topics exist in the literature and are discussed to help improve the situation and realize the potential for cancer treatment. This paper is rich with references on mathematical models and does not take into account the optimization part.
In another work [23], the authors confirm that the two major obstacles against successful chemotherapy of cancer are the cell-cyclephase dependence of treatment and the emergence of resistance of they consider models based on amplification of the resistance gene up to a very large number of copies. This approach is to study basic mathematical properties of the models to help in understanding the control problem.
The simplest mathematical models that describe the optimal control of cancer chemotherapy treat the entire cell cycle as one compartment [51]. In reality, each drug affects the cell evolution in a particular phase. Therefore, it makes sense to combine these drugs to produce the greatest cumulative effect on the cancer population. However, as 20 years have already passed, updating the level of knowledge in chemotherapy-treatment planning research is necessary. Clare et al. [13] introduce several models in their review of the application of mathematical models on breast cancer and discuss the mathematical modeling of adjuvant chemotherapy as one of the subsections. Parker and Doyle [52] provide a comprehensive review of the articles using mathematical modeling for drug delivery with only a subsection on the optimal cancer chemotherapy. This review describes the primary factors considered in the chemotherapy treatment models without surveying the literature in mathematical modeling systematically.
In the previous discourse, we discussed and appraised some mathematical models of cancer treatment. We noted the benefits of using the theory of optimal control to determine chemotherapy scheduling, and we underscored the need for cancer treatment to be viewed as a formal optimization problem.

The optimization model
Optimal control theory is widely used to model chemotherapy treatment planning problems. In this paper, we first explain the problem formulation. Then, we discuss the model describing the tumor growth and give an explication of several conditions that affect the optimal control results, such as the objective function and the treatment's function, duration and dosage.
To properly pose the optimal-control problem, we must define the goal we wish to minimize (i.e., the objective functional, the treatment dosage and duration). Hence, an optimization problem consists of a set of constraints that must be fulfilled and an objective function to be minimized (or maximized) with respect to a set of independent variables.
Let us describe the continuous optimal problem as follows: Tf Subject to: Where j w  represents the different type of cells (i.e., cancer, normal or immune cells) described in the model studied, u is the dosage of the drug and f f is the final value function. This paper will be organized according to the definition of the cancer cells to cytotoxic agents. These two problems are understood by applying optimal control theory to mathematical models of cell dynamics. The authors have defined a mathematical model that can be used to pose and solve an optimal chemotherapy problem under evolving resistance and estimated the parameters of the constructed models.
In [40], Clairambault presents a review concerning several problems encountered by biologists and physicians who address natural cell proliferation and disruptions of its physiological control in cancer disease. He concludes that the modeling of cell proliferation and drug disposition by continuous time evolution equations can help by i) studying the theoretical effects of the control functions; ii) proving the feasibility or unfeasibility of proposed therapeutic strategies, under accurate biological theoretical assumptions on the biological framework used; and iii) proposing optimized multidrug multi-targeted therapies according to criteria and constraints defined by oncologists. In this paper, the author describes in detail the biological phenomena related to cancer and treatment without going into the mathematical details of in which case they are mentioned and whether they are mentioned briefly or in a general and simple way.
Finally, in [41], Billy et al. review i) some models that describe the evolution of cancer cells under the influence of anticancer drugs, which have been used or may be used to tackle the general problem of therapeutic optimization in oncology; and ii) theoretical therapeutic optimization methods that can be used in the context of various models of cell-population growth. Then, they present several techniques used for the identification of the parameters of population-dynamics models used in chemotherapy, and they focus on a novel method of optimization under unwanted toxicity constraints. This method is based on the optimization of eigenvalues in an age-structured model of cell-population dynamics.
The following is another set of references where drug therapies for cancers that involve tumors, are studied: Chakrabarty and Hanson [42], Duda [43], Fister and Panetta [19], Murray [43,36] and Swan [27] on cancer chemotherapy. For alternate types of objective functions, we can cite the work of Swierniak, Polanski, and Kimmel [44], Murray [43] and Ledzewicz and Schättler [20]. The books by Sethi and Thompson [45], Cohen [46], Clark [47] and Eisen [48] provide the background for the optimal control theory we use and include some simple biological applications. The book of Martin and Teo [49] applies optimal control to several detailed models of cancer tumors, and provides a survey of many research results of optimal control applied to cancer.
Mathematical work that has been performed in the optimal control setting includes two non-cell-cycle-specific (drugs that are effective in all of the phases of the cell cycle) models by Murray [43,36]. Swan [27] provides a good review of the role of optimal control in noncell-cycle-specific cancer chemotherapy: he not only elaborates on miscellaneous-growth kinetic models and cell-cycle models, but he also briefly describes important issues in chemotherapy treatment planning, such as drug resistance that requires the use of quantitative research. Swierniak, Polanski, and Kimmel [44] use optimal control theory on a cell-cycle-specific chemotherapeutic model. Swierniak,Polanski,and Kimmel [44], along with Swan [27], investigate the effects of the drugs on the normal tissue and use this to limit the drug strength, but only for non-cell-cycle-specific treatments. Therefore, we develop an optimal-control problem that will directly determine the effects of the cell-cycle-specific treatments on the normal tissue, and we attempt to relate the mathematical results to known clinical information. In another work [50], the authors stress the aspect of drug resistance; mathematical optimization problem. For this reason, we will review i) a set of deterministic models that describe the evolution of cancer cells; ii) a set of mathematical models that describe the pharmacokinetics and pharmacodynamics of injected chemotherapy drugs; iii) a set of objective functions, considered in the literature, and affect the results of the optimization method; and iv) various optimization methods found in the literature. All these points are reviewed to optimize the theoretical drug delivery and the number of cancer cells in chemotherapy treatment.
The cancer models presented include deterministic, stochastic, compartmental, spatial, and hybrid models. They provide the reader with a state of the art overview of the rapidly evolving field of cancer modeling and treatment optimization.
The remainder of this paper is organized into sections. In the section of review of Cancer Modeling, we review the different deterministic mathematical models used to describe tumor growth under the influence of chemotherapy. In the section of review of treatments, we review the treatment type, the administration method, the drug combination and the treatment protocols for the chemotherapy. In the section of Review of the Objective Function, we review the various objective functions that affect the results of the optimization method. In the section of Solution Methods, we summarize the role of optimal control and the optimal control models applied to the chemotherapy treatment planning and group them according to the solution methods used. Finally, in the section of Discussion and Conclusions, we describe the challenges of using mathematical modeling and discuss the gap between theoretical research and their clinical applications. We also provide several future research directions.

Review of Cancer Modeling
Because there are three distinct stages (avascular, vascular, and metastatic) to cancer development, researchers often concentrate their efforts on answering specific questions on each of these stages [53]. This review aims to describe the current state of the mathematical modeling of avascular tumor growth, i.e., tumors without blood vessels. This is because, when attempting to model any complex system, it is wise to attempt to understand each of the components as well as possible before they are all put together. Avascular tumor growth is much simpler to model mathematically, and yet it contains many of the phenomena that we will need to address in a general model of tumor growth.
Cancer modeling has a wide variety of forms. Indeed, it can involve almost any type of applied mathematics, as shown in Figure 2. We use probability models to understand how genetic mutations lead to cancer progression, metastasis and resistance to therapy. Ordinary differential equations can be used to study the growth of tumor cell populations, often leading to a conclusion of Gompertzian growth [54].
Partial differential equation models (PDE) using cell densities and nutrient concentrations as state variables can be used to analyze various spatio-temporal phenomena [23].
The agent-based model (ABM) that treats cells as discrete objects with predefined rules of interaction can offer an improvement over PDE methods in some situations, such as the study of angiogenesis [55,56], the development of new blood vessels to bring nutrients to a growing tumor [57]. In [58], an extensive theoretical investigation of the process of tumor-induced angiogenesis is presented, and the results from computational simulations of the mathematical model have highlighted a number of important new targets for therapeutic intervention.
In general, PDE models are used with hyperthermia cancer treatment [59], and the ABM model is used with immunotherapy treatment [60]. The ABM model often contains components of multiscale models, which have the ability to simulate tumor properties across multiple scales in space and time [61,62].
Another phenomenon is incorporated into the process of cancer development: metastasis [63,64]. In case of the formation of another cancer type via metastasis phenomena, we resolve this problem using two different protocol treatments for each cancer to minimize the number of cancer cells with the minimum treatment dosage.
In this paper, we are interested in reviewing the models used to describe tumor growth with chemotherapy treatment, which can be viewed as a formal optimization problem.
The mathematical models that describe tumor growth are divided into five categories:

•
Deterministic model with non-linear differential equations: In this model, we used a system of ordinary differential equations (ODE) to describe the evolution of cells reflected by a welldefined curve.
Where ( ) f t may take different non-linear forms.
Compartmental model: The cell cycle is modeled as compartments, each of which either describe different cell phases or combine phases of the cell cycle into clusters. The function ( ) f t is linear, and the optimization solutions are analytical.
Stochastic model: Observations of cancer cell cycles are often presented in an erratic fashion. For this reason, we are trying to introduce probabilities, such as those used in stochastic models.
The various stochastic models that are able to describe biological processes, include the following:
• Spatio-temporal model [70][71][72]: In this model, the evolution of cancer cells is studied in space as well as time. We can add Laplacian terms into the different equations of the deterministic models to study the spatial evolution of cancer cells.

Deterministic model
Mathematical models describing the cancer-evolution phenomena can be used as an important tool in therapy planning. Several models of cancer at various stages have been formulated over many years of research. In this paper, we focus on models in the form of ordinary differential equations (ODEs), controlling cancer growth on a cellpopulation level. Many references, such as [25,28,30], studied cancer growth using Gompertzian growth [26]. Other ODE models have also been used to describe the exchange between two cells populations, proliferating and quiescence ones [25]. Other ODE models have considered some form of tumor immune interaction [31,36]. Normal cells have also been modeled in many forms to model toxicity disadvantage [25,74].
We must divide the deterministic model into three types: • Non-linear differential equations • Compartmental model

• Hybrid model
Non-linear differential equations: This model must be further divided based on the number of differential equations, which are used to describe the evolution of several cell types (including cancer, normal and immune cells). Historically, these models are cited as below in the various articles but scientifically they are not related.
We start with the simplest model: Model with one differential equation: The most frequently used model in oncology for growth retardation in malignant tumors is the Gompertz equation [75], developed by the 19th century English actuary, Benjamin Gompertz [54,76]. The equation is characterized by two parameters, the initial specific rate of exponential growth and the rate of exponential fall in the initial growth rate. Unlike the logistic function, which is symmetrical about the inflexion, the Gompertz function is asymmetrical about the inflexion, which is always 0.37 of the asymptotic value [76].
Typically a tumor will grow very rapidly with the growth rate proportional to the tumor size. However, as tumors increase, the growth rate decreases as tumor size increases. In order to reflect this, differential equation models have been proposed. It is customary to assume that tumor volume and number of tumor cells are proportional and so the model which follows describes ( ) dC t dt , where ( ) C t represents the number of tumor cells at time t. Three models have been widely used in the study of tumor growth: exponential, logistic and Gompertz growth. The effect of chemotherapy (pharmacodynamics) is incorporated into tumor growth model by adding a kill term to the differential equation by assuming that there is a concentration ( ) u t of a cytotoxic drug at time t and the presence of the drug will cause a decline in tumor cell population jointly proportional to concentration ( ) u t and the population size ( ) C t at any given instant. We will also assume that there is a threshold drug concentration level, ε below which no tumor cells are killed. The Gompertz growth equation is the most commonly used model of tumor growth and takes with treatment the form: Where λ is a positive constant related to the growth function and θ is the largest possible size of the tumor. k is a constant of proportionality, and H is the Heaviside unit function [77].
Among the proposed models, those based on Gompertz growth are frequently proposed [78]. It models the growth of a population consisting of a group of individuals of one or more similar species in the absence of migration and interaction with other species, by an ( ) that is the solution of: Within the context of tumor growth, C denotes the number of tumor cells at time t , 0 C is the size of the tumor at the start of treatment, ( ) u t denotes the drug concentration at the site of action.
( ) C t is assumed to be continuous and differentiable [79]. Parameters α and β (measured in 1 t − ), denote growth and decay rates, respectively, they characterize the evolution of different tumor types. During chemotherapy it is assumed that the basic growth kinetics of the tumor will be perturbed by the action of the cytotoxic drug ( ) u t .
Swan and Vincent [12] used assumptions about Gompertzian growth of immunoglobulin G (IgG) multiple myeloma cells and developed a single differential equation of drug action on such cells. The equation describing the drug action is as follows: Where, the parameter θ represents the greatest size of the tumor; and 1 k and 2 k are the positive constants in the saturation type loss term in the equation.

Model with two differential equations: The cell growth of various types of cells, whether they are immune, cancer or normal cells have certain features in common. Therefore, a mathematical model of a particular cell growth may be used to model other types of cell growth, provided some basic changes can be made.
A few models of acute myeloblastic leukemia (AML) were proposed by Afenya [80]. Normal and leukemic cells were assumed side-by-side with the two cell populations obeying Gompertzian dynamics but with the leukemic cells exercising inhibition over the normal cells. The kinetic equations and steady-state properties of one of the models are analytically obtained as: where A N and θ are the asymptotic bounds (limits) on, or carrying capacities of, the normal and leukemic cell populations, respectively, and ( ) , f C N is a function that can take many form: (10) a is the intrinsic growth rate or growth speed of the normal cells, and b is their death rate. c is a measure of the degree of inhibition of the normal cells with respect to the leukemic cells, k is the fraction of leukemic cells that are killed due to the drug effect, and h is the fraction of the normal cells that are destroyed by the lethal binding effects of the drugs (it is assumed that k h > ). r is the maximum fraction of normal cells that should be injected per unit time to ensure normal cell regrowth. It can also be interpreted to be the maximal rate per unit time at which growth factors are infused. We assume that, as the normal cell population decreases because of the drug effect, the regenerative process through which normal cells are injected or growth factors are infused could be directly related to this normal population. This is modeled by the term rN . The negative effects of the drugs on the normal cells are modeled by the term ( ) hu t N [78]. m is the maximum growth coefficient of the normal cells and ( ) 1 (1 , which could take a number of forms, is the growth fraction of the normal cells, which is dependent on the leukemic population, and p can be considered as an inhibitive parameter. Model with three differential equations: In this model, we let ( ) n T t denote the native T cell population, and ( ) e T t the effector T cell population at time t [38,81]. We assume that the effector T cells are specific to one type of cancer (e.g., chronic myelogenous leukemia, CML), activated by the presence of the CML antigen. The state system for our model is given by: This model consists of a system of three nonlinear ordinary differential equations. Each equation represents the rate of change with respect to time of one of the cell populations. All parameter descriptions are in the cited reference [38].
Another model with three differential equations is described as follows by Villasana and Ochoa [82]: This model divides the population of the tumor cells into interphase cells (pre-mitotic, phase period comprising G 1 , S and G 2 ) and mitosis cells, which are represented by ( ) Model with four differential equations: Immunotherapies are quickly becoming an important component in the multi-pronged approaches being developed to treat certain forms of cancer [84]. The goal of immunotherapy is to strengthen the body's own natural ability to combat cancer by enhancing the effectiveness of the immune system. The importance of the immune system in fighting cancer has been verified in the laboratory as well as with clinical experiments [85][86][87][88]. Additionally, it is known that those with weakened immune systems, such as those suffering from AIDS, are more likely to contract certain rare forms of cancer. This phenomenon can be interpreted as providing further evidence that the role played by the immune response in battling cancer is critical [89,90]. Through the mathematical modeling of tumor growth, the presence of an immune component has been shown to be essential for producing clinically observed phenomena such as tumor dormancy, oscillations in tumor size, and spontaneous tumor regression.
The mathematical modeling of the entire immune system can be an enormously intricate task, as demonstrated in [91], so models that describe the immune system response to a tumor challenge must necessarily focus on those elements of the immune system that are known to be significant in controlling tumor growth. In the work of de Boer and Hogeweg [92], a mathematical model of the cellular immune response was used to investigate such an immune reaction to tumors. It was found that, initially, small doses of antigens do lead to tumor dormancy.
This model describes the kinetics of four cell populations (tumor cells and three types of immune cells), as well as the concentrations of two drugs in the bloodstream, using a series of coupled ordinary differential equations that are based on the model developed by de Pillis and Radunskaya [93]. The populations at time t are denoted by the following terms [84]: The equations governing the population kinetics must take into account a net growth term for each population, the fractional cell kill, the per cell recruitment, the cell inactivation and the external intervention with medication. We attempt to use the simplest expressions for each term that still accurately reflect the experimental data and population interactions.
Bringing together the specific forms for each cell growth and interaction term leads to the full system of equations: All parameter descriptions are in Table 1 of the cited reference [84]. The term I represents the immunotherapy drug concentration in the bloodstream, and in our case we set I equal zero because we are interested only in the chemotherapy drug concentration.
Finally, as stated above, Villasana and Ocho's model does not include quiescent tumor cells. In [94], we show how Yafia [95] develops a delayed differential equation model for the interactions of proliferating and quiescent tumor cells. The concept of cancer cells in the quiescent phase ( ) Q T is introduced, and the model by Villasana, cited above, is corrected: (14) Compartmental model: In the compartmental model, the cell cycle is modeled in the form of compartments that describe different cell phases or combine phases of the cell cycle into clusters. Each cell passes through a sequence of phases from cell birth to cell division. The starting point is the growth phase G 1 after which the cell enters the phase S, in which the DNA synthesis occurs. Then, the second growth phase G 2 occurs, in which the cell prepares for mitosis or phase M, in which cell division occurs.
Each of the two offspring cells can either reenter phase G 1 or may simply lie dormant for some time in a separate phase G 0 until reentering G 1 , thus starting the entire process all over again.
The simplest mathematical models that describe the optimal control of cancer chemotherapy treat the entire cell cycle as one compartment, but solutions to these single-compartment models are  For the compartmental model, the problem of finding an optimal cancer chemotherapy protocol is formulated as an optimal control problem over the finite time interval of the fixed therapy horizon. The state variable is given by the average number of cancer cells, and the control is the effect of the drug dosages on the respective subpopulation. The goal is to maximize the number of cancer cells that the agent kills and to appropriately minimize the number of cancer cells at the end of the therapy session, while keeping the toxicity to the normal tissues acceptable. The last aspect is modeled implicitly by including an integral of the control over the therapy interval in the objective so that minimizing controls will balance the quantity of drugs given with the conflicting objective to kill cancer cells. The analytical approaches to these models are based on applications of the Pontryagin Maximum Principle.
The general compartmental model can be described uniformly by the following state equation [22]: compartment, i j ≠ . This is satisfied for each of the models described below. More generally, if the condition that the first orthant of the control system is positively invariant were violated, this would be a strong indication that the modeling is inconsistent.
We start from mathematical models describing the dynamics of the cancer cell population under a single drug treatment, from the simplest, two-compartmental, to infinite-dimensional ones. Problems of multidrug treatment and multidrug resistance are addressed afterward. A model with an infinite number of compartments was considered in a series of papers by Kimmel, Swierniak, and coworkers (e.g., references [22,23,42,51,96,97]), in the context of drug resistance caused by the gene amplification dynamical process.
The infinite-dimensional compartmental model of drug-resistance evolution has been proposed to approximate the dynamic branching random walk related to the process of gene amplification [98] and has led to results not achievable in simplified finite-dimensional models. These results in turn have made it possible to formulate several recommendations regarding treatment protocols. Optimization of the model has become possible after its transformation to the integrodifferential form.
We briefly recall some compartment models that fit into this general class. For a more detailed description of the models, we refer the reader to Swierniak et al. [44].
In a two-compartment model, the phases G 0 , G 1 and S are clustered into the first compartment, G 2 and M are combined into the second compartment, and only a killing agent = are given by: The i a values are positive coefficients related to the mean transit times of cells through the i th − compartment.
The differential equations of this model are: (17) In this three-compartment model, a blocking agent 2 u ν = is additionally considered which is active in the synthesis phase S; thus, S is modeled as a separate compartment. Now, 3, 2 n m = = , and the matrices are given by: Thus, the differential equations of this model are: (1 ) (19) In both models, the control 1 u u = represents the dose of the killing agent administered with the value 0 = u corresponding to no treatment and 1 u = corresponding to a maximum dose. It is assumed that the dose stands in direct relation to the fraction of cells that are being killed in the G 2 /M phase. Therefore, only the fraction 1 u − of the outflow of cells from the last compartment undergoes cell division and reenters the first compartment. However, all cells leave compartment G 2 /M. In the second model, the blocking agent 2 u ν = is additionally applied to slow the transit times of cancer cells during the synthesis phase S. As a result, the flow of cancer cells from the second into the third compartment is reduced by a factor of 1 ν − of its original flow to Here the control ( ) 0 t ν = corresponds to no drug being applied, and a maximal reduction occurs with a full dose max ν .
In the models mentioned above, the compartments are divided based on the cell cycle. Next, we must divide the compartments depending on the number of cells sensible or resistant to the drug injected [74]. The most basic model consists of only two compartments representing the sensitive and drug resistant subpopulations. It has been a subject analyzed by many researchers but arguably, the most recent comprehensive studies can be found in Ledzewicz and Schättler [20,21].
The model based on the assumption of the mutation being the basis for drug resistance caused the cell to acquire one additional copy of a certain gene. Therefore, a sensitive mother cell produces two daughter cells, one of which remains sensitive, while the other changes into a resistant one with a probability of γ . Similarly, if a resistant cell undergoes cell division, one of the offspring remains sensitive. The other one may mutate back into a sensitive cell with a probability d (it is assumed that no additional gene copies can be acquired).
Let us denote by 0 C and 1 C the average number of sensitive and drug-resistant cells, respectively Then, the system is described by the following set of equations: This model represent the evolution of cancer cells using a singleagent treatment; in Panetta [99], the simplest model possible was analyzed using two treatment agents, and the mutation that makes cells resistant to the first drug is assumed to be irreversible. Therefore, the second sub population resistant to this first drug 1 u is assumed to be sensitive to the second drug 2 u , so the resulting mathematical model is as follows: where γ is a probability of a mutational event leading to drug resistance. For this model, no formal optimization has been performed.
If we consider killing and blocking agent actions (denoted by 1 u and 2 u , respectively), the issue of resistance is much more complex.
However, if only the resistance to the killing agent is considered, then, the compartments of cells in the phase G 1 , S, G 2 and M, denoted by 0,1, 2 i = , respectively, and by the 3 i ≥ compartments of resistant cells, the following description is obtained: (1 ) The Another form of the compartmental model with two compartments is denote as the model of Gyllenberg and Webb [79,100]; the creators used this model to show that it yields Gompertz or logistic growth under special parameter selection. The model consists of proliferating and quiescent cell compartments; it allows for transition between the compartments and cell death from either compartment. The transition rate functions are considered to be functions of the total number of cells. The following set of ordinary differential equations describes the model: This idea was briefly introduced in a recent paper [101]; we start by defining the net transition rate function: When ( ) 0 C Φ > , the net transition rate is from the proliferating compartment into the quiescent compartment.
Analytical approaches to these models are based on applications of the Pontryagin Maximum Principle [102].
A mathematical model of two ordinary differential equations for the interactions between the cancer cell growth and the activity of the immune system during the development of cancer was proposed by Stepanova [103]. This model is described as follows [104]: Where T represents the immunocompetent cell densities related to various types of T cells activated during the immune reaction, x C denotes the tumor volume, 1 δ corresponds to a threshold beyond which the immunological system becomes depressed by the growing tumor, and the coefficients I µ and δ are used to calibrate the interactions between the cancer and the immune cell, and in the product with T collectively describe the state-dependent influence of the cancer cells on the stimulation of the immune system, The coefficients x k and Y k are chosen to normalize the control set, i.e., we assume that 0 1 u < < . In Stepanova's original formulation, F is simply given by ( ) 1 E x F C = , in Gompertzian growth, F is given by , and in the logistic growth models, F is given Hybrid model: Discrete or continuous time equations are used for the deterministic model. It is represented by ordinary differential equations (ODEs) with delay or partial differential equations (PDE) [79] or the agent-based model (ABM) [105]. Among the suggested models, the simplest ones are those of one or several differential equations. Despite their simplicity, these models must predict the evolution of many biological phenomena [106]. Models based on the deterministic Gompertz law appear to be particularly compatible with the evidence of tumor growth.
However, it is quite common that discrepancies are found to exist between clinical data and theoretical predictions due to intense environmental fluctuations. For instance, Ferreira et al. [107] analyzed the effect of various chemotherapeutic strategies on the vascular tumor growth. They confirmed that an environment such as chemotherapy would affect tumor growth behavior and lead to morphological transitions under certain conditions. Therefore, a better model is necessary to reflect the external randomness that affects the tumorgrowth behavior.
Ferrante et al. [108] proposed a stochastic version of the Gompertz law in which random fluctuations of the model parameters are considered. They assume that the growth deceleration factor β (death rate) does not change, whereas the variability of the environmental conditions induces fluctuations in the intrinsic growth rate α . The intrinsic growth rate is assumed to vary in time according to the following equation: The parameter α is the constant mean value of ( ) t θ , σ is the diffusion coefficient, and ( ) t ε is a Gaussian white noise process.
By adding the expression for ( ) t θ into Eq. (7) without treatment, we have: Briefly, all deterministic equations could then be written by adding the term to the equation, which describes the evolution of cancer cells, transforming them from deterministic to hybrid models.
Finally, the hybrid model is a combination between deterministic and stochastic models. It appears in three forms: , or λ and k are constants, and a white noise term is added. Then, Eq. (6) becomes: For the model described by Eq. (8), α is also a random number that can be represented by the same three distribution options described above, if we add a white noise to this model, the equation becomes: For a model with three differential equations, described by Eq. (11), we obtain the following equation, when α and β are set to constants and white noise is added:

Review of Treatments
Standard chemotherapies are typically administered with a constant dose scheduling [22,[109][110][111][112]. In the administration of cancer treatments, it is conventional that strategic dosing is used to maximize anticancer-drug effects while minimizing host toxicity [113]. Accordingly, many therapy schedules employ intensive therapy initially, when the tumor is the largest, and then the dose is decreased as the tumor is reduced.

Treatment type and administration method
Drug delivery is the method or process of administering a pharmaceutical compound to achieve a therapeutic effect; drugdelivery technologies modify the drug-release profile, absorption, distribution and elimination for the benefit of improving the product efficacy and safety, as well as patient convenience and compliance. The most common routes of administration are the oral, topical (skin), transmucosal and inhalation routes. Many medications may not be delivered using these routes because they might be susceptible to enzymatic degradation or cannot be absorbed into the systemic circulation efficiently due to their molecular size. For this reason, the drugs must be delivered by injection. Current efforts in the area of drug delivery include the development of targeted delivery in which the drug is only active in the target area of the body and of sustainedrelease formulations in which the drug is released over a period of time in a controlled manner from a formulation. To achieve efficient targeted delivery, the designed system must avoid the host's defense mechanisms and circulate to its intended site of action [114]. One type of sustained-release formulation includes liposomes.
We are interested i) in the different methods of chemotherapydrug administration that can be used in cancer treatments, and ii) for each type of drug administration to describe the pharmacokinetic (PK) models that predict the hematic-drug concentration after their administration.

Continuous-infusion drugs:
An infusion is a method used to put fluids, including saline and drugs, directly into the bloodstream as a body-wide way to fight cancer.
Because infusional chemotherapy is administered directly into the blood, every cell in the body is exposed to the drugs. Cancer cells as well as certain healthy cells may be affected. Blood counts may change after each treatment depending on the drugs given, so a test called a complete blood count (CBC) will be administered to check the levels of white and red cells as well as other elements in the blood.
Mathematical models used to describe the diffusion throughout the body of drugs that are administered via continuous infusion are described below.
For comparison, we first consider a therapy in which the dose linearly increases with time [115]. This can be implemented by simply replacing the constant dosage treatment u with the following: Next, we also examine the late logarithmic intensification therapy proposed by Gonzalez et al. [116]. To model such a logarithmic therapy, we replace the constant dosage u with the below term: where e is the 'Neper constant' equal to 2.718. We choose Γ and β to be two adjustable negative parameters that control the rate because clinically, a therapy with a linearly increasing intensity could be fatal to the patient.
Another form in the administration of drugs is used, the sizedependent therapy [117]. For this, we replace the constant dosage u by the following equation: ψ is the initial tumor size, 1 A is the intrinsic growth rate of the tumor (related to the initial mitosis rate) and 2 A is the growth deceleration factor (related to the antiangiogenic processes) and σ is the diffusion coefficient. Determining the schedule of treatment depends mainly on the initial state of the cancerous cell population.
The function cited above describes the best way to administer a treatment to the body. Next, we need to describe the pharmacokinetic (PK) models that predict the hematic drug concentration after their administration. In the literature, the majority of therapeutic studies did not include the pharmacokinetic action of the drug injected but a general efficacy term is introduced, in the form min max u u u < < , to represent the percentage effectiveness of the drug.
Several different approaches have been used in the past to model the fate of a drug within the human body [118]; several of them are summarized in this paper. First, we address the mechanistic compartmental models, the one-compartment model and the twocompartment model. The one-compartment model is based on the schematization of the human body as a single reactor; in this case, the evolution of the drug concentration in the tumor tissues over time [119] is given by: Here, 0 u is a control function characterizing the rate at which the drug is introduced into the tumor, and γ is the dissipation coefficient.
The generalizations of this model (25) were formulated in [51,120] to demonstrate its flexibility.
Another form is used to describe the PK in a one or two compartmental model [121]; this form is explained as follows: In the case of the one-compartment model, the human body is represented as a single reactor, as shown in Figure 3. The twocompartment model, shown in Figure 4, is based on the representation of the body by two reactors, the first of which is a central compartment, involved in absorption, distribution, metabolism and elimination, which represents the plasma, and the second represents tissues, involved only in distribution.
In these two models, ( ) g t represents the mass rate with which the drug reaches the reactor and depends on the route of administration,  Another approach is used to distinguish the plasma and the active drug concentrations [122]; most mathematical models assume that the drug is instantaneously delivered to the cancer site. We avoided this undesirable simplification by considering the dynamic relationship between the kinetic behavior of the drug administered and its corresponding concentration profile at the cellular level. We propose a compartment model, as shown in Figure 5, to express the   V is the effect compartment of the active concentrations 2 ( ) u t [118,123]. From the plasma compartment, the drug reaches the effect compartment, and then it is eliminated from it: ( ) u t represents the drug-input function, which is the intermittent or continuous drug delivery by intravenous infusion, the amount i d , denoting infused function between the start-and the end-of-infusion times 1 i t − and i t , respectively.
Another approach to determine the drug concentration in the body is described below [124]: The distribution of the drug in a tumor is represented in terms of its concentrations in three compartments: intracellular, ( ) The plasma concentration for an infusion time T is, for t T < [126]: And for t>T: Then, the plasma concentration ( ) v u t for bolus injection is given by: (41) Here, u is the total dose injected, β γ and C are given in [124].
The equation used for a liposomal drug in the plasma is a biexponential fit made by Gabizon et al. [128] for their clinical data: (42) For the liposomal drug in the extracellular space: Lv Le re du u P S u u dt (43) Therefore, the plasma concentration of the free drug v u is governed by: The parameters appearing in these equations are defined in Table  1 of reference [124].
The above model for liposomes has been adapted for thermo sensitive liposomes [128,129], which are designed to release their contents rapidly upon heating. First, re τ the time constant for drug release from the liposomes in the tumor extracellular space, is replaced by the following time function: if t t t t (45) Here, h t is the time at which hyperthermia begins (after the initial drug injection), and d t is the duration of the hyperthermia. A new constant, h re τ is the time of drug release from the heated liposomes [130]. The permeability of the vessel to liposomes is assumed to increase during heating, so the parameter L P is replaced by a function of time: E is an enhancement factor for L P at 45°C. Gaber et al. [130] observed a 76-fold increase in the liposome extravasation upon heating to 45°C. It is difficult to deduce from this the increase in the value of L P itself. The range of E=1 to 100 is considered here.
bolus drug injection: The injection of a drug (or drugs) in a high quantity is called a bolus [12]. It is a relatively large dose of medication administered into a vein in a short period, usually within 1 to 30 minutes. The intravenous (IV) bolus is commonly used when rapid administration of a medication is needed, such as in an emergency; when drugs that cannot be diluted, such as many cancer chemotherapeutic drugs, are administered; and when the therapeutic purpose is to achieve a peak drug level in the bloodstream of the patient. The IV bolus is not used when the medication must be diluted in a large-volume parenteral fluid before entering the bloodstream or when the rapid administration of a medication, such as potassium chloride, may be life threatening. Bolus injections allow medication to become useful to a patient faster, which can be the difference between life and death in some situations.
The mathematical model used to describe the diffusion throughout the body of drugs that are administered via bolus drug injection is similar to the model (Eq. (36)). However, in this case, ( ) g t is equal to zero because of the instantaneous (limited to the initial instant) drug absorption.
To describe the evolution in the time of the drug concentration with the one-and two-compartment models, the initial conditions that are required to solve the balanced equations can be found in the recent Orally Administered Drugs: Oral chemotherapy doses are set up so that the patients will have constant levels of the drugs in their bodies to kill the cancer cells. Not taking chemotherapy drugs as they should be taken can affect how well the treatment works, and it can even allow the cancer to grow, so sometimes changes may be needed. Even after starting to feel better, they may still have cancer cells in the body that must be kept under control with chemotherapy [12].
Oral chemotherapy drugs may be taken every day, every week, or once or twice a month; sometimes, the frequency is as much as several times a day. It may necessary to use chemotherapy for several months or longer. Chemotherapy is often given in cycles. This means that the chemotherapy will be used for a period of time, and then there will be a break. This allows the body of the patient to grow new, healthy cells.
The mathematical model used to describe the diffusion throughout the body of drugs that are administered via the oral method is similar to the model (Eq. (36)). However, in this case, ( ) g t is equal to t is the drug mass available at the time t at the site of absorption, the profile of which is defined by the pharmaceutical system-dissolution characteristics and by the drug solubility.
For the oral method, the drug mass balance at the site of absorption is [121,131]: The function ( ) f t corresponds to immediate availability of the entire dose u , and we define ( ) f t as follows: For case a, the ingestion is of a system with a high solubility, or a high-permeability drug; in this case, the drug's mass that lies at the site of absorption changes over time only because of absorption proceeding. For case b, the system ingested is a low solubility or low permeability drug; in this case, the drug's mass that lies at the site of absorption additionally changes between the instant t and the following instant because of the change in the fraction of the drug that is actually available for absorption. A k is the absorption constant, and the function ( ) r t expresses the in vitro release. The relationship between in vitro and in vivo release , where fu is the maximum fraction of the pharmacological dose ( u ) that could actually be absorbed, taking into account the liver metabolism.
To describe the evolution in the time of drug concentration for the one-and two-compartment models, the initial conditions that are required to solve the balanced equations are mentioned in the recent article [122] and are observed below: In oral assumption In these three modes of administration, the drug can be taken alone with dilution or in a liposome or a thermo liposomal capsule. Liposomes are artificially prepared vesicles made of a lipid bilayer, which can be filled with drugs and used to deliver drugs for cancer and other diseases. They are composite structures made of phospholipids and may contain small amounts of other molecules [132]. Thermo liposomal capsules have the same definition as liposomal capsules with the difference that the thermo liposomal capsule can release its drug content within minutes of heating [124].

Drug combinations and treatment protocols
Chemotherapy was first introduced in the 1940s [133]. For the next 20 years, it was considered an investigational treatment. In the last 30 years, chemotherapy information has evolved, and many more effective drugs have been developed. During this time, doctors have documented responses and conducted clinical trials comparing standard treatments to new treatments. This process of gathering chemotherapy information has helped to establish specific protocols regarding the types, doses and dosing schedule of drugs that are based on the type, stage, and other specifics of a person's cancer. There is no one correct choice in chemotherapy treatment. Each treatment protocol has advantages and disadvantages, and there may be more than one good option. In addition, treatment choices can change over time. A good chemotherapy treatment choice at one time may not be the right choice later.
The development of drug resistance is one reason that drugs are often given in combination. It is thought that this may reduce the incidence of a resistance developing to any one drug. Often, if a cancer becomes resistant to one drug or group of drugs, it is more likely that the cancer will also be resistant to other drugs. Thus, it is very important to select the best possible treatment protocol at the outset. Table 1 presents the different combinations of treatment used for some cancers. For some of the combinations, we cited treatment protocols: • CAF treatment consists of the cyclic administration of [134] drug "C" orally for 14 days, while "A" and "F" are given together, intravenously (i.v.) into the hand or arm, on days 1 and 8 of that 2-week period. This schedule will repeat four to six times, once every 4 weeks. The entire process takes a total of approximately 4 to 5 months, barring any complications that slow it down. Another option is that C, A, and F are all given simultaneously, via a drip into the arm or hand. This treatment is repeated every three weeks, four to six times, barring any complications. The standard dosages of CAF are as follows: 600 mg/ m 2 cyclophosphamide (19.9 mg/kg), 60 mg/m 2 adriamycin (1.9 mg/kg) and 600 mg/m 2 fluorouracil. Doxorubicin, bleomycin, vinblastine, and dacarbazine regimen (ABVD) [135] : Combined modality consisting of a doxorubicin 25mg/m 2 IV plus a bleomycin 10-IU/m 2 IV plus a vinblastine 6-mg/m 2 IV plus a dacarbazine 375-mg/m 2 IV on days 1 and 15; every 28ad for 2-4 cycles; followed by involved field radiation therapy (IFRT) at a dose of approximately 20 Gy.

Review of the Objective Function
Optimal control, which has been widely used in the cancer literature, is the process of determining the control and state trajectories for a dynamic system over a period of time to minimize (or maximize) the final value of a single variable J . This can be written mathematically as, Eq. (1): The objective function may be composed of any number of functions, combined to produce a single number, usually through a weighted sum. It is difficult to convey the knowledge and experience of a clinician to an optimizer through one or more simple functions. The use of the drug is subjected to a constraint represented by the limit of toxicity. Therefore, numerical solutions are found, considering various objective functions. It is assumed that the drug is administered by continuous infusion and that the number of tumor cells decreases during treatment. The toxicity of the treatment causes a limitation in the dosing [142]; then, u must be within the range max In this section, we give examples of objective functions considered in the literature on the treatment of cancers. Starting with a quadratic function: Where 2 ( ) C t denotes the number of cancer cells sensitive to treatment and 2 ( ) u t represents the nonlinear cost of the treatment. The quadratic forms do not make much sense, as they do not represent anything biological. The idea behind them came from control theory and closed-loop control, in which the quadratic error is minimized and the quadratic control terms arise from the energy interpretation.
The linear function (Eq. (52)), is used to minimize the number of tumor cells at the end of the treatment period while limiting the side effects of drugs [36].
The logarithmic function (Eq. (53)) is used to prevent, for example, the development of new subpopulations that are resistant to drugs. On the other hand, the logarithmic function (Eq. (52)) was mainly used because of the logarithmic transformation of state variables, which is especially useful in the case of Gompertz-type growth equations.
These functions have been used for the optimal control of the noncell-cycle-specific model (drugs that are effective in all of the phases of the cell cycle) introduced by Murray [36]. Additionally, Swan [27] provides a good review of the role of optimal control in non-cell-cyclespecific cancer chemotherapy. However, in another study we focus on optimal control problems as applied to cell-cycle-specific chemotherapy. First, Eisen [143] developed a system of linear differential equations describing the growth dynamics of the proliferating (drug-sensitive phase) and quiescent (drug-resistant phase) cells. In this work, the control over a given interval reduces the cancer to a fixed level while minimizing the total drug use. Another work by Swierniak, Polanski, and Kimmel [44] uses optimal control theory on a cell-cycle-specific chemotherapeutic model. In the stochastic model, different objectives can be used that combine to maximize the probability of cure while minimizing toxicity. Here, we achieve this goal by maximizing the probability of uncomplicated control [144]: where ( ) T CUMP t is the cumulative probability of a toxic event with: tox U denotes any explicit limits on the toxicity, and { ( ) 0} N P R t = is the probability that the number of cancer cells equals zero.
A patient can fail treatment for two reasons. Either the treatment may not be able to eliminate the tumor or the treatment itself may create toxicity for a patient, to the point that the regimen cannot be completed. Both of these results lead to treatment failure. This function expresses the probability that neither of these circumstances occur. It determines the probability for a given regimen that a patient can complete treatment (i.e., that it will not be prematurely stopped due to toxicity) and have the tumor eliminated.
Several studies using a compartmental model are able to solve their optimal control model and find an analytical expression for the optimal chemotherapy treatment plan.
In most optimal control studies, the sufficient and necessary conditions for optimality, optimality conditions, can be obtained using frameworks such as the Pontryagin Maximum Principle [102]. On the other hand, several other studies solve their optimal control model by using sequential quadratic programming [145]. The objective function to be minimized is given by [120]: Where m is the number of administered drugs and r and 1 r are weighing factors. This general model even allows for the inclusion of special cases, such as multidrug therapy. However, this model is valid only in cases when two of the following conditions are met: • Each drug affects cells of a different type (true in the basic model of a killing and a blocking agent treatment).

•
Either the molecular source of resistance to each drug is identical (as in multidrug resistance [146]) or the infinite subsystem representing gene amplification is required for only one type of a drug (the basis for resistance to other drugs requires only a single mutation and there is only one level of resistance for each of them).
The objective function of the compartmental model in the case of a single drug is described as: where i C denotes the average number of cancer cells in the i th − compartment, T denotes the time at the end of therapy, and r and i r are weighting factors. Similarly, Eq. (57) was used in the context of phase-specific chemotherapy in a number of papers, such as [20,22,44].
For a deterministic model with the three differential equations described above, (Eq. (11)), we have the following objective function: minimized. It is expected that the effects of the drugs are nonlinear, and we choose the quadratic cost terms 2 1 ( ) u t and 2 2 ( ) u t to reflect these effects, as in [38]. The coefficients B 1 and B 2 are weight constants for the controls and include a measure of toxicity of the drugs to the body. The salvage term 3 ( ) f B C t is included to counter the effect of using a fixed treatment time. If this term were not present, the controls could taper off earlier and allow a rise in the cancer-cell count at the end of the treatment period. The salvage term is included to create a penalty for low values of n T because this affects the patient's ability to fight off other diseases.
For the Stepanova model [103], we have the following objective function: Where a and b are positives coefficients determined by a stable eigenvector, and ε is a positive constant.
For the delay differential equation model [83], we have the following objective function: And: Another objective function is described in a recent paper [134], as follows: The first objective of treating a patient with stage IIB breast cancer, who is in adjuvant chemotherapy for four months, is to find the optimal value of ( ) u t in a way that the number of cancer cells, ( ) C t , resulting from this treatment, track ( ) r C t . Thus, the optimal value of ( ) u t can minimize the cost function, described as follows: Where, ( ) r C t is a reference trajectory, and is the period of treatment in this case. Because biological systems usually respond in a sigmoidal fashion to inputs, ( ) r C t is derived by the nonlinear rescaling [147,148].
The second objective function is used to preserving the normal cell population ( ) N t in the best way by minimizing the following cost function: where N ∞ is an asymptotic number of normal cells, and ( ) N i − for i = 21, 42, 63 is the normal cell population at the times just before the second, third and fourth cycles of the chemotherapy. (84) N is the normal cell population 21 days after the fourth cycle of chemotherapy.
Different beneficial cases can be found between these two extreme situations by minimizing the cost function, described as follows: (64) where, the weighting coefficients 1 W and 2 W are selected by the treating physician according to the particular conditions of each patient to indicate the cases for which the cancer-cell reduction has more priority than preserving the normal cells 1 2 ( ) W W > or vice versa.

Solution Methods
As we have already cited, optimal-control theory is usually used to model chemotherapy-treatment planning. Optimal-control theory uses a system of differential equations that are solved to determine the optimal choice of the chemotherapy-treatment plans. These problems are difficult to solve optimally. For these reasons, the optimal control theory is applied with simplifying assumptions that reduces its clinical validity.
In the first part, this paragraph presents a simple review of the ways in which optimal-control theory interacts with cancer chemotherapy. In the second part, we review the optimal-control models based on the solution methods used: the analytical solution, approximation with analytical intuition, and heuristics [37].

Role of optimal-control theory
There are three main areas of study: the first studies the diverse growth-kinetics models, the second studies cell-cycle models and the third is a classification of 'other models'.
Diverse growth-kinetics models: In this section, the authors in [27] use the Gompertz equation to describe the evolution of cancer cells under the influence of chemotherapy treatment. The optimal-control problem is formulated using the Hamiltonian function to determine the positive continuous time optimal controller u that minimizes the toxicity of the treatment subject to the perturbed growth kinetics of the tumor.
Another model used to describe the evolution of cancer cells is the Verhulst-Pearl equation [149]. For this model, several numerical results to the solution are presented. The optimal control u increases, but as time increases, its rate of increase slows down rapidly. For the Cox-Woodbury-Meyers equation [150], there is a nonlinear algebraic feedback relationship that connects the control and the state, and the optimal controller is monotonically decreasing and drives the tumor population to a more desirable level. It is hoped to see closed-loop control.
Finally, Abulesz and Lyberatos use a model to describe the evolution of cancer and the normal cells [151], with the introduction of a pharmacokinetic equation. They give an optimal solution for the drug administered, but it is not clear what the biological significance of the steady state means, and it is not clear how this applies to the dosage rate.
Cell-cycle models: The first application of optimal-control theory in a cancer-chemotherapy problem is to a cell-cycle model analyzed by Bahrami and Kim [28], using a discrete-time state vector approach. The optimal-control problem is to determine the controls u(k) such that the size of the tumor is minimized at the end of the treatment interval while minimizing the toxicity.
Creasey et al. [152] present a general overview of the optimal control and cancer chemotherapy, but no particular problem is solved. Long and Fegley [153] consider a discrete dosage program, which is cast into the framework of dynamic programming. They wish to determine the doses and treatment times to maximize the surviving fraction of normal cells at the end of the recovery period. They do not apply their procedure to any illustrative problem.
Biran an dMcKinnis [154] considered a single closed loop of the cell cycle in three phases. The time rate of change of each number of cells in one phase is related to the flux of cells into and out of each phase. Numerical results are presented for the case of use of the drug Melphalan. The optimal control results suggest that cultured cells treated with Melphalan accumulate and arrest their progress in the mitotic phase.
In another paper, Dibrov et al. [155] examine a dynamic model of a population of continuously dividing cells under multiple periodic treatments with a phase-specific agent. They introduce a model of the cell cycle that has both deterministic and what they refer to as "probabilistic" discrete compartments.
Other models: In this section, several models are presented that, are not directly dependent on the growth kinetics of the tumor but could be combined with them if desired. In the papers by Bellman and coworkers [156] the authors study the "control aspects of the pioneering efforts in the modeling of chemotherapy. Bischoff and Dedrick [157] described a physiological model for the anticancer drug concentration in anatomical compartments. Swan and Alexandro [158] present the reduction of a physiological model for the chemotherapy of brain tumors.
Problems of the following type remain to be investigated: determining the optimal drug inputs in regional chemotherapy to achieve a cell kill in the tumor while minimizing the toxicity to important anatomical compartments. Gaglio et al. [159] describe two expert systems that are being used for the characterization of optimal adjuvant cancer therapies. Unfortunately, no details of the optimalcontrol application are provided. Maceratini et al. [160] have reviewed expert systems and some of their applications in medicine.
For more details about the role of optimal-control theory, we refer the reader to reference [161], in which four different mathematical models of chemotherapy from the literature are investigated with respect to the optimal control of the drug-treatment schedules. The various models are based on two different sets of ordinary differential equations and contain either chemotherapy, immunotherapy, antiangiogenic therapy or combinations of these. Optimal-control problem formulations based on these models are proposed, discussed and compared. For different parameter sets, scenarios, and objective functions, optimal-control problems are solved numerically with Bock's direct multiple shooting methods.

Role of the Nonlinear Optimization Method
In this paragraph, we classify the optimal-control models with respect to the solution method used.

Analytical solution:
We use the analytical method when the model is composed of linear differential equations; in most optimalcontrol studies, the sufficient and necessary conditions for optimality, the optimality conditions, can be obtained using frameworks such as the Pontryagin Maximum Principle [99]. Zietz and Nicolini [33] use the Pareto version of the Pontryagin Maximum Principle (Yu and Leitmann [162]) to show that a specific type of chemotherapy plan is optimal under certain conditions; their solution depends on the two weighting factors in the Hamiltonian form of the model. Murray [36] uses Pontryagin Maximum Principle to show that the optimal chemotherapy treatment plan is a mixture of an initial bolus dose (an immediate infusion of a single dose) drug application followed by no drug and then continuous infusion.
Another optimal control is solved by analytical solution, for example, by explicit formulation of the treatment plan [30,33,43], or by explicit formulation of the treatment administration period given the dose size [163].

Approximation with analytical intuition:
In some studies, the optimality conditions obtained using the Pontryagin's Maximum Principle are not simple enough to explicitly derive the optimal chemotherapy treatment plan. The optimality conditions provide an idea about the structure of the optimal solutions. However, optimality conditions, which are used in the characterization of the optimal solution, are expressed in terms of adjoint variables. These adjoint variables turn up in the Hamiltonian function that the optimal solution needs to maximize for all instances. These conditions give a large system of differential equations, which is difficult to solve. For this reason, these studies use approximation techniques to determine the optimal chemotherapy treatment plan, such as i) discretization of the continuous decision horizon and solving of the optimal model using control parameterization; and ii) solving the system of equations obtained from the optimality conditions using approximate methods such as Newton's method [164].
Certain optimal controls are solved by approximation with analytical intuition. Several example are by using standard methods with hypothetical data [165], by characterizing the optimal solutions for quadratic and linear controls [145], by using an iterative algorithm that may converge to a local optimum [166], by transforming the model to mixed integer linear programming (MILP) by discretizing the decision horizon and linearizing the nonlinear constraints [15], by characterizing analytically the form of the optimal solution and using discretization and nonlinear programming to solve it [167], by using a linear time-varying approximation technique, or finally, by using an explicit formulation of the optimal treatment plan as a function of variables from the Hamiltonian form of the model that are quantified via approximation [38].

Heuristics:
Although the optimal solution is typically computationally intractable, we develop heuristic algorithms to solve the optimal chemotherapy treatment. Tan et al. [168] use distributed evolutionary computing software that employs the resources of a network of computers to overcome complex chemotherapy treatment planning problems. Floares et al. [147] design an adaptive neural network algorithm to solve this problem. Liang et al. [169] minimize the size of the tumor population using an adaptive elitist-populationbased genetic algorithm, and they design a multimodal optimization genetic algorithm that represents the chemotherapy treatment plans as genes to optimize the chemotherapy plan under various cumulative toxicity functions.
Several optimal controls are solved by heuristics. For example, Iliadis and Barbolosi [118] use simulation to evaluate the performance of possible drug-dose plans for given drug-administration periods. Petrovski and McCall [170] adapted the strength Pareto evolutionary algorithm [170] to evolve a population of treatments spread out along the Pareto efficient frontier and find the Pareto-optimal chemotherapy plans. They also develop a new heuristic algorithm in their further research: a genetic and a particle swarm optimization (PSO) algorithm [171]. The genetic algorithm encodes the multi-drug chemotherapy schedules as binary strings. The PSO algorithm uses a population of candidate-solution particles that swam around the search space. Tse et al. [172], minimize the size of the tumor population using a customized algorithm combining a genetic algorithm and iterative dynamic programming.
All of the methods cited above are used to find a global optimum value, knowing that there are several computational biology research studies [142,173,174] that integrate local search algorithms into a global research algorithm to increase the computing accuracy.

Discussion and Conclusion
The number of cancerous cells, the toxicity, and the drug resistance are the key factors in chemotherapy-treatment planning. Most of the studies on chemotherapy-treatment optimization problem differ from each other according to how they model the interactions among these key factors.
This paper provides a comprehensive review of the relevant literature while considering several other pieces of information not mentioned previously, such as the cancer modeling characteristics and optimization computational methods used to solve these models. We also provide information on the medical relevance, such as real treatment protocols and one-drug/multi-drug involvement.
Our study is oriented toward the optimization of treatment protocols. The optimization methods can still lead to solutions, even if the parameters to identify are numerous. Clinically, however, to identify the parameters, we need many experiences. The number of parameters is large, and therefore, the number of experiences is large; for this reason, it is preferable that the number of parameters is small. It is clear from all of the studies presented that there is a desperate need to create more interaction between mathematicians, clinicians and biologists to correctly identify the model. Of course, these highly interdisciplinary efforts are not easy to achieve. We see the role of mathematical modeling in cancer biology as twofold. On the one hand, mathematical models are able to verify supposed word models suggested by experimentalists.
The field of cancer biology is reaching a stage of maturity at which the next step in the modeling process must be the careful parameterization of a number of the models so that specific experimental predictions can be made and tested in close collaboration between experimentalists and theoreticians. Many problems in this review require the determination of the parameters in the model by fitting the model to data. Clinical tests are performed over short durations, so we require an efficient identification protocol to correctly identify the model parameters. The clinical trials occur over a short duration, whereas the phenomenon is long, occurring over several months.
In the literature, most mathematical models that describe the diffusion of drugs in the body are based on a single drug. Whereas, in real cases, oncologist use a treatment protocol with a combination of drugs. In some cases, these drugs work together in series or parallel mode, or without any interaction between them. Some drugs have physical effects and other have chemical effects. Based on these information, it is necessary to develop these models to consider all these aspects of administered drugs. Although there have been some success stories in the application of mathematical models to cancer biology with chemotherapy treatment, mathematics has much more to offer for the use of mathematical models with real chemotherapy treatment protocols. For example, in Eq. (9), this model describes the use of a single drug. However, by using a real protocol, the treatment of many drugs must be considered together, so we propose to change the form of the equation: This form must be verified by mathematicians and biologists to ensure its reliability. Briefly, one of the advantages of mathematical modeling is its ability to determine a delicate balance when building any model between reliability and realism.
Research in chemotherapy should consider a more realistic objectives function. Currently, most models aim to obtain the chemotherapy treatment plan that minimizes the size of the tumor population while satisfying the constraints on maximum toxicity in a given time interval. However, this is only a part of an optimal treatment plan. In this context, the question of how to define the objective function naturally becomes more important. In this paper we have reviewed all of the objectives functions to include and cover all of the possibilities of cancer phenomena in our optimization method.
It is important to compare the performance of all objective functions to choose the best that can lead to the optimal solutions. Reading various articles shows that all objective functions lead practically to optimal solutions close to each other, without showing the specificity of each objective function. For these reasons, comparative study is important.
The chemotherapy process involves the interactions between the patient's response and the oncologist's scheduling decision. During the chemotherapy session, when drug resistance becomes obvious, the drugs must be changed immediately. During this process, better information about the characteristics of the tumor will be available, which necessitates dynamic and sequential decision making. On the other hand, current models consider only the decision process prior to the initiation of chemotherapy and then calculate the optimal plan of the entire period. For these reasons, we incorporate in our optimization method the evolution of the normal and immune cells to help the oncologist to make a better decision about the administration of the treatment. In real cases, the evolution of normal and cancer cells varies according to the treatment. The relations between the treatment dosage and the number of normal and immune cells are clearly visible in the mathematical models that describe the evolution of these cells.
Additionally, we can improve the optimization problem by adding as constraints the evolution of the white blood cells (WBC) in the body after treatment administration. There are mathematical models that describes the evolution of the number of immune cells under the influence of cancer and drugs. We can add Eq. (66) [118] for mathematical models that do not have an equation that describes the evolution of the number of immune cells. In this case, the oncologist can define if it is necessary to stop the treatment if the number of immune cells is less than a predefined number. (66) where WBC are produced at rate ( ) r t and eliminated by the firstorder process The normal range for the WBC count is 4.3 to 10.8×10 9 cells per liter. A range of 11 to 17×10 9 /L may be considered mild-to-moderate leukocytosis, and a range of 3.0 to 5.0×10 9 /L may be considered mild leucopenia [142].
In the literature, applying the optimization method to the mathematical models that describe the evolution of cancer cells, the authors do not provide any clinical applications. However, in our further research, the usefulness of these models is shown for using real chemotherapy treatment protocols. The major limitation of the existing studies is the lack of clinical realism in the models. Most studies treat chemotherapy-treatment planning as a pure optimization problem. As a result, the practical applications of this area of research are limited. In this paper, we provide several suggestions to improve the research in this area.
Mathematical models must be developed to aid in our understanding of how to implement this work and hopefully to show how to develop new optimal treatment strategies. To accomplish this goal, we can determine a strategy that makes it possible by following it to apply the optimization method using real chemotherapy treatment protocols. This strategy can be described as follows: 1) Define the genre of the cancer treated.
2) Choose the model of cancer evolution with and without treatment.
3) Identify the parameters of this model according to the cancer treated.

4)
Choose a real treatment protocol defined by the oncologist.

5)
Choose the pharmacokinetic model that describes the evolution of the drug concentration after its administration.
6) Define the model of resistance (e.g., white blood cells) added as constraints.

7)
Define the minimal and maximal thresholds for the resistance of treatment and white blood cells.

8)
Resolve the optimization problem using one solution method (we recommend the genetic algorithm).
By using a real treatment protocol, we suggest that the treatment optimization is performed using genetic algorithms (GA) [175] to avoid the problem of oversizing caused by the use of deterministic optimization methods. GA allows us to set an initial population size to a fixed value at the beginning of the search that remains constant throughout the run. However, making it necessary to specify this initial parameter value is problematic in many ways. If it is too small, the GA may not be able to reach high-quality solutions. If it is too large, the GA spends unnecessary computational resources.
We propose the use the genetic algorithm because it is a method that has been adopted in all cases. It allows us to offer a software tool ready for use by oncologists. Additionally, the deterministic models converge according to each case, so it is not a generalized method contrariwise for stochastic models.
Due to the complexity of cancer, chemotherapy-treatment studies should focus on a specific cancer rather than modeling cancers in general. It is important to biologically classify the models, based on the cancer definition and not as a function of the mathematical description (e.g., linear or nonlinear models).
It is evident from all of the studies presented that there is a desperate need to have more interaction between the modelers and either the experimentalists or the clinicians. Of course, such highly interdisciplinary endeavors are not easily accomplished. Another feature of importance is the need to work with models of biological relevance that can be tested on a computer and in a laboratory situation.
At present, there are no commonly used techniques to measure the number of tumor cells in humans when the level is below the minimum level of a tumor that can be diagnosed, (C d ). Some clinicians believe that one should treat the tumor only to the level C d because, as indicated, they cannot measure the size below C d . This means that we should attempt to find a model to describe cancer appearance starting with at least one cell. In the thesis [176], the authors demonstrate, using stochastic models, that the cancer passes through two phases: i) a creation phase, in which they study the probability that a cancer will arise from a healthy population of cells under the influence of mutations, and ii) a phase in which the cancer appears. The first phase is explained by stochastic static models. These models are divided into three types, i) the Moran model [64,65], ii) the Wright-Fisher model [66] and iii) the Moolgavkar, Venzon, and Knudson (MVK) model [69]. Simulations of these models show the uncontrollable evolution of cancer cells. The cancer cell number grows until reaching the population size or to a nonpredictable number.
In this paper, we reviewed i) the mathematical models that describe the evolution of cancer and biological cells; ii) the mathematical models that describe the diffusion of drugs throughout the body that are administered via bolus injection, continuous infusion, and liposomal and thermo liposomal drug delivery; iii) a set of objective functions that affect the results of the optimization method; and iv) a solution method for the optimizations models. These studies imply significant gaps between theoretical research and clinical application in chemotherapy treatments. We note that the application of the optimization method in the real case of chemotherapy is only a tool that helps doctors to determine the best decision for a patient's treatment.
Finally, the question of therapeutic optimization in cancer is vast, and it can be treated in notably different manners that must be adapted for the particular clinical problem at hand. Nevertheless, modeling the target and the means of control while considering the known clinical issues is still a successful way to reduce the cost and time of clinical trials. Room remains for further mathematical developments to meet other challenges, especially for optimizing treatments for more clinical problems. These results will be all the more useful as more links develop between mathematicians and clinicians.