A mathematical justification for metronomic chemotherapy in oncology

We mathematically justify metronomic chemotherapy as the best strategy to apply most cytotoxic drugs in oncology for both curative and palliative approaches, assuming the classical pharmacokinetic model together with the Emax pharmacodynamic and the Norton-Simon hypothesis. From the mathematical point of view, we will consider two mixed-integer nonlinear optimization problems, where the unknowns are the number of the doses and the quantity of each one, adjusting the administration times a posteriori


Introduction
Since the origins of chemotherapy in oncology, the standard and most widely used treatment has been the maximum tolerated dose (MTD) one: "the higher the dose, the better", whose rationale seems very intuitive. As everyone knows, the trade-off lies in the side effects due to drug toxicity, which not only diminishes quality of life for the patient (adding more illness to the existing one), but also make it necessary to impose rest periods between cycles of therapy, which compromises a good resolution of the cancer treatment. For this reason, the design of more effective and less toxic chemotherapeutic regimens (by choosing the timing and fractionation of the doses) is a very active area of research in oncology.
A turning point in cancer chemotherapy treatments can be placed in the year 2000, when a revolutionary alternative approach called "metronomic chemotherapy" (MC) was proposed, see [13] and the references therein. According to National Cancer Institute [7], this is the "treatment in which low doses of anticancer drugs are given on a continuous or frequent, regular schedule (such as daily or weekly), usually over a long time. Metronomic chemotherapy causes less severe side effects than standard chemotherapy. Giving low doses of chemotherapy may stop the growth of new blood vessels that tumors need to grow". In oncology, this last effect is referred to as "anti-angiogenesis".
The historical controversy between efficacy in tumor killing and lack of toxicity has given rise to numerous studies and experimental trials comparing MTD and MC, the main conclusion of which can be summarized by saying that "the body of experimental and clinical evidence, coupled with theoretical considerations [...] point to MC as a preferred course of action", see [5]. Among the theoretical considerations, in addition to the anti-angiogenesis effect, we can mention that MC can activate anti-tumor immunity and minimize therapeutic resistance.
In this work, we will justify mathematically that MC is the best strategy in oncology for most cytotoxic drugs and both curative and palliative treatments, using simple well-established mathematical models and without taking these additional benefits into account. To that end, we will study two different optimization problems that arise when modeling usual chemotherapy treatments for cancer: to minimize the tumor volume at a fixed final time (curative approach) or to minimize the total cumulative dose that a patient must take during a specific period of time for maintaining the tumor volume below a given threshold (palliative approach). We will assume a Gomperztian type tumor growth and the Norton-Simon hypothesis to represent the growth-inhibitory influence due to the cytotoxic chemotherapy effect (see [9,10]) with the classical pharmacokinetics and the Emax pharmacodynamic model. As far as we know, this is the first time that the superiority of MC is mathematically justified for minimally parameterized models. Up to now, the efficiency of MC was related with heterogeneous tumor populations (sensitive/resistant cells) and palliative purposes, see [6] and the references therein.
From the mathematical point of view, the optimization problems under consideration are nonlinear and of mixed-integer type (where we will have the number of treatment doses as an integer variable, N , and the amount of each individual dose as a continuous variable, d i ). In both cases we will give explicit expressions of the approximated optimal solutions which will lead us to choose treatments corresponding to MC. A posteriori, the administration times can be adjusted (in fact, in several ways) completing the description of the treatment. We will illustrate the utility of our approach with numerical experiments concerning the treatment of some brain tumors that exhibit Gompertzian growth (see [14]) using a cytotoxic drug called Temozolomide (TMZ). This drug is the subject of many clinical and biological recent studies (see [2,3,8,11] and their references). The origin of this paper was an academic project [12] submitted by the third author under the supervision of the other authors.

Associated optimization problems
The Gompertzian law has often been used to simulate the growth of some untreated tumors (see for instance [1]). In mathematical terms this can be described by the following ordinary differential equation (ODE): where Ψ(L) = ξL log θ L .
In the above equations, L(t) represents the tumor volume (or tumor size) at time t, ξ its growth rate and θ the maximum size it can reach (also called carrying capacity for biological systems). In practice, the parameters ξ and θ can also change with time (due to angiogenesis), but here we will focus on the simpler problem in which they are fixed.
As mentioned before, when using a cytotoxic drug to treat the tumor, we will consider the Norton--Simon hypothesis to model its effect; more precisely, we will assume that the growth-inhibiting effect due to the treatment is proportional to the growth rate of the untreated tumor (see [9,10]). Mathematically, this can be written as Here the term Ψ(L(t))ρ(t) represents the growth-inhibitory influence due to the cytotoxic chemotherapy effect and it reflects both, the level of the therapy at time t, ρ(t), and the tumor's sensitivity to therapy. First let us present a general result for the existence and uniqueness of solution for the Cauchy problem associated with (2) (see [4] for the proof). As usual, we will denote by W 1,∞ (0, T ) the Sobolev space of all functions in L ∞ (0, T ) having first order weak derivative (in the distributional sense) also belonging to L ∞ (0, T ). It is well known that W 1,∞ (0, T ) can be identified with C 0,1 [0, T ], the space of Lipschitz continuous functions in [0, T ], after a possible redefinition on a set of zero measure. Theorem 1. Let us assume that ξ, θ, L 0 and T are given positive real numbers, L 0 ∈ (0, θ) and ρ ∈ L ∞ (0, T ). Then, there exists a unique solution L ∈ W 1,∞ (0, T ) of the following Cauchy problem given by Depending on the expression of ρ(t), different pharmacodynamics (PD) can be studied. In this work we consider the classic Emax model: where k 1 and k 2 are fixed positive real numbers and c(t) denotes the concentration of the drug in the tumor at time t. Other choices have been considered in the literature, as the Skipper model where ρ(t) = k 1 c(t). The Emax model seems to be more appropriate from the clinical point of view because it saturates for high concentration values. Another crucial term is related with the pharmacokinetics (PK). We use the following Cauchy problem to describe the PK of the drug: The coefficient λ is the clearance rate and it is related to the half-life of the drug. The second term on the right hand of the equation depends on both the specific drug and the way it is administered. We assume that N doses, The coefficient σ is determined on the basis of drug-, patient-and tumorspecific parameters. For instance, in the case of brain tumors the drug loss during the transport to the brain has to be taking into account, see [2] and the references therein. In other cases, σ = α V D β , where α is the patient's body surface area, V D is the volume of distribution of the drug and β is the patient weight in kg if V D is given in l/kg. Finally, δ(t − t i ) is the Dirac delta distribution concentrated at t i . In principle, it seems reasonable to consider the administration times {t i } N i=1 as variables of the problem, as well as N and the doses {d i } N i=1 . However, this approach would lead us into serious difficulties because the criterion for deciding a priori the characteristics that determine the feasible protocols may be unknown even to the specialist. In practice, this would make the mathematical approach unfeasible. Instead, in a first step, we can consider the administration times as data, fixed according to some pattern. It is clear that this assumption will implicitly affect the solution obtained: for example, a daily regimen will inevitably lead to a solution with small doses, because MTD cannot be administered over a long period due to the serious side effects it produces. Apparently, it looks impossible to get out of this vicious circle, but we will see that this can be done for most cytotoxic drugs by following our approach and readjusting the administration times a posteriori, once N and the doses {d i } N i=1 have been determined. Now, let us present our optimization problems: in this paper we are interested in the solutions of two optimization problems associated with the Cauchy problem (3) using the Emax model (5) and when the following data are given: the final time T , T > t N and the levels of effectiveness and the toxicity of each individual dose, d min and d max , respectively. We assume 0 < d min < d max . To simplify some expressions we will use t N +1 = T . As we have explained before, the administration times t i , i = 1, . . . , N , are fixed a priori following some pattern, but the number N of doses of the treatment is a variable to be determined.
Firstly, we consider a curative approach with the goal of minimizing the tumor volume at the final time T with a fixed cumulative dose, D, that will be divided into N smaller doses, {d i } N i=1 , to be determined. Hence we formulate the following problem with N + 1 variables Using (4) and the fact that L 0 < θ (where L 0 is the initial tumor size), we transform (P 1 ) into the following problem: Now, let's formulate the problem related to a palliative approach in which we want to minimize the total dose (determining N and d i , for i = 1, . . . , N ) while maintaining the tumor volume at T below a given threshold L * > 0 that is not harmful to the patient, allowing him to have an acceptable quality of life. We consider the same bound constraints on the individual doses as in the previous optimization problem: Using again (4) it is easy to check the equivalence: Therefore, we can reformulate the problem (P 2 ) as follows:

Optimal treatments
In view of our optimization problems it is convenient to have an explicit formulation for the concentration of the drug. (6) is given by

Proposition 1. The solution of the Cauchy problem
Proof. Applying Laplace transform to problem (6) we get and consequently, This implies that which is equivalent to (12).
We continue by obtaining the expression for the integral term appearing in (8) and (11). In this respect, we present the following result: Lemma 1. Let us assume that c(t) is given by (12) and ρ(t) by (5). Then, it is verified that Proof. Using (12) and the definition of ρ, we derive (13) from the following equalities:

Curative approach
We start by considering the problem (P 1 ) defined in Section 1 (in which we want to minimize the tumor volume at time T with a fixed cumulative dose D). In this case, using (8), (13), the positivity of k 1 /λ and the monotonicity of the logarithm function, we reformulate the problem (P 1 ) as follows: Note that it is a mixed-integer nonlinear optimization problem with an integer variable, N , and N continuous variables, d i .
In order to simplify the objective function f 1 , along the rest of the work, we will assume the main hypothesis that we are considering cytotoxic drugs verifying where s denotes the minimum of the elapsed time periods; this is Let us stress that this condition holds (at least) when λ is sufficiently large (or equivalently, the half-life of the drug is sufficiently small) and that most of the usual cytotoxic drugs satisfy this condition, even when s is small. Assuming (15), the following approximation for the objective function holds: that leads to the mixed-integer nonlinear programming problem: At this point, it is very important to underline that the administration times {t i } N i=1 do not appear in the formulation of (P 1 ). This is a crucial property that will allow us to break the vicious circle mentioned in the previous section: for most cytotoxic drugs we will first determine N and the doses as the solution of (P 1 ) and then assign appropriate administration times, checking that (15) holds.
Using the bound constraints and the total dose value, in order to have a non-empty set of feasible points, we must assume In particular, hypothesis (18) implies that only a finite number of values for variable N have to be taken into account to solve the optimization problem. Namely, where x denotes the greatest integer less than or equal to x and x the least integer greater than or equal to x. Now, for each fixed feasible value of N , we can formulate the following continuous nonlinear programming problem: It is not difficult to prove that equal-dosage is the optimal solution of the previous problem: Proof. For any feasible point d, using the relation between the geometric and arithmetic means, we get and therefore,f is a feasible point for (P N 1 ) such that the corresponding objective function value reaches the above upper bound, so the conclusion of the theorem follows straightforwardly.
As a consequence, the optimal solution of the approximate mixed-integer optimization problem (P 1 ) leads to the longest feasible treatment: Corollary 1. Let us assume that d min > 0 and (18). Then the optimal solution of (P 1 ) is given by (N ,d), withN = D/d min andd i = D/N , for i = 1, . . . ,N .
Proof. This can be deduced combining Theorem 2 with the fact that the auxiliary function x is strictly increasing in (0, +∞) and noticing that Furthermore, we can study the dependence of the objective function optimal value with respect to d min . Using previous expressions, we derivê On the one hand, this implies that to increasef 1 (N ,d) (and therefore to decrease the tumor size at the fixed final time) is interesting to take d min as small as possible, in line with the principles of MC. But, on the other hand, from the clinical point of view, it is known that for any drug there is a minimum effective dose (MED) (the smallest dose with a discernible useful effect) that must be taken into account. In general, the estimation of the MED (and even its definition) remains controversial and this carries over to the determination of the optimal d min in the practice.

Remark 1.
i) Choosing the Gompertzian law for the untreated tumor growth is not a crucial aspect for deriving our results. Other typical ODEs with sigmoid-type solutions may be reasonable choices, such as (1) with with γ > 0. For γ = 1, we recover the well-known logistic ODE, meanwhile for small γ > 0 , this is an approximation of the Gompertz ODE, see [1]. All these options lead to the same results.
ii) Instead, it can be shown that the Norton-Simon hypothesis is crucial for our argument. There are other classical hypotheses such as the log-kill one. It states that a given dose of chemotherapy kills the same fraction of tumor cells regardless of the size of the tumor at the time of treatment, which leads to the following problem instead of (3). Although the log-kill hypothesis continues to appear frequently in the literature, it was shown in [9] to be inconsistent with some clinical experiences.
iii) We have studied a related problem for continuous (non-discrete) drug infusion in [4] (labeled (OP 2 ) with G = G 2 ). There, it was proved the appearance of a constant maintenance infusion rate (during a fairly long time interval) in the expression of the optimal control. This can be seen as a piecewise continuous version of our previous Corollary 1.
From a theoretical point of view, we have seen that the optimal curative treatment (i.e. aiming to minimize the tumor size at a given time with a fixed cumulative dose using drugs for which (16) holds) is the longest feasible treatment with equal individual doses. This is the case for drugs with a short half-life such as Temozolomide (TMZ), an oral alkylating agent with antitumor activity in high and low grade gliomas (HGG and LGG, respectively). In the monotherapy phase for adult patients with newly-diagnosed glioblastoma multiforme, TMZ is usually administered consecutively for 5 days of every 28 days (this is denoted by 5/28d) for a maximum of 6 cycles with individual doses of 150 mg/m 2 in first cycle and 200 mg/m 2 in the others. So, taking into account this usual treatment (UT), for our numerical experiments we have fixed T = 210 days for the final time and D = 5,750 mg/m 2 for the cumulative dose.
The parameter values used for our numerical computations are showed in Table 1. Following the reference cited in [2] we know that TMZ has a mean elimination half-life of 1.8 hours (i.e. t 1/2 = 0.075 days). The clearance rate λ is estimated using that λ = log(2)/t 1/2 = 9.242 days −1 . The growth rate ξ depends on the type of tumor; the value ξ = 5.51e−3 could correspond to an HGG with a volume-doubling time about 126 days (see [14] for growth dynamics of glioblastomas). Following [2], we have used σ = 1.6(2.5e−3) ≈ 4e−3 that corresponds to a woman with a body surface of 1.6 m 2 and it takes into account the fraction of TMZ getting to her brain interstitium. On the other hand, it is well known that k 1 represents the maximum effect of the drug on the tumor and k 2 is the effective concentration, i.e concentration producing 50% of the maximum effect. We have chosen k 2 = 0.36 mg/l, a value that is consistent with predicted peak concentrations of TMZ from [2,11]. Moreover, we have taken k 1 = 60 for illustrative purposes. Let us note that the parameter k 1 does not play any role either in the optimization process of this subsection nor in the main hypothesis (15), but it does influence the evaluation of the tumor size. Let us stress that (15) holds in this case, because d max e −λs ≈ 0.01938 k 2 = 90, with s = 1. Finally, let us underline that the value of θ is unnecessary here, because the initial datum and the results are given rescaled in the form L 0 /θ and L(t)/θ, resp. (see Table 1 and Figure 1). In Table 2 we present some numerical results to illustrate and compare the optimal treatments corresponding to (P 1 ) without using the approximation (16) (i.e. (14)) with those associated to (P 1 ) (given by Theorem 2 and Corollary 1) in the case of individual doses between d min = 100 mg/m 2 and d max = 200 mg/m 2 . To solve the problems associated with (P 1 ) when the value for N is fixed, we have used FMINCON (a nonlinear programming solver provided in MATLAB's Optimization Toolbox) with the SQP algorithm and the value 1.e−8 for the optimality tolerance. Here, we present results only for the feasible values of N taking into account the length of the treatment interval and the number of dosage times for the 5/28d schedule: N ∈ [29, 40]. This number N appears in the first column; then we present two columns with the individual doses of the treatment described in Theorem 2 and the ratio of the final tumor volume to the initial volume, L(T )/L 0 , for an initial medium size tumor (25% of carrying capacity); the next three columns report results given by FMINCON: the minimum and maximum values reached by the individual doses of each treatment,d min andd max , and its corresponding ratio L(T )/L 0 . Let us note at this point that, for each value of N , the individual doses obtained with FMINCON are almost identical and the optimal solution is reached administering 40 doses (the longest 5/28d treatment, as Theorem 2 and Corollary 1 state, using the approximation (16)). The last line of the table is associated to the usual treatment (UT) described above using N = 30 doses. It is noteworthy that with it the final tumor volume is 91% of the initial volume while our optimal treatment provides 73%. At the top of Figure 1 we see the evolution of the tumor volume, L(t)/θ, for a selection of treatments of Table 2.
The most relevant results of this section are given in Table 3 and at the top of Figure 2, where we illustrate the influence of the level d min and compare several optimal treatments corresponding to (P 1 ) requiring the same cumulative dose (D = 5,750 mg/m 2 ). The rest of the parameters are given in Table 1. As we have mentioned previously, first we have fixed d min and solved (P 1 ); once we have determined the corresponding solution (N and {d i }N i=1 ) we have selected the most appropriate schedules for TMZ used in some studies, see [8]: the 7/14d with doses of 150 mg/m 2 , the 21/28d with doses of 100 mg/m 2 and 75 mg/m 2 and the daily regimen with 50 mg/m 2 doses. The results of the penultimate column of Table 3 have been calculated with the times given by these schedules. Let us underline that the specific values of the administration times do not affect the optimal solution and the only relevant factor is the maximum number of doses that can be administered with the chosen d min . We observe that the greatest tumor shrinkage is achieved with the longest treatment which also has the lowest doses, in line with MC. Moreover, three of the treatments of Table 3 are more effective than the optimal treatment of Table 2, all of them improving the usual treatment UT.
As predicted by the Norton-Simon hypothesis [10], it is apparent that   Tables 2 and 4). the optimal solution is more "dose-dense" than other feasible protocols, see the top of Figure 2 (the period between cycles is shorter). Nevertheless, concerning the "dose-intensity" (defined as the total dose administered during a treatment, divided by its duration) it can be seen in the last column of Table 3 that the optimal protocol has the lowest one.   Tables 3 and 5).

Palliative approach
In this section we are interested in the optimization problem (P 2 ) i.e. we want to minimize the total administered dose with a constraint on the tumor size at the final time: L(T ) ≤ L * , where L * is a given level beforehand. Using (10) and Lemma 1, that constraint is equivalent to the following inequality: being f 1 the objective function of problem (14) and Let us remember that T R depends on the data ξ, θ, L 0 and L * (see (10)). Therefore, using (11), the optimization problem (P 2 ) can be reformulated as follows: with d min > 0 andT R ∈ R defined in (22). Using the approximation (16) in the nonlinear constraint we arrive to the following nonlinear programming problem withf 1 defined in (17). As one can see, (P 2 ) is also a mixed-integer nonlinear optimization problem with an integer variable N and N continuous variables, d i .
In this case, let us remark that the set of feasible values for N is infinite: using the general inequality constraint and the upper bound constraints, the following inequality should be satisfied In fact, only a finite subset has practical interest, i.e. those values that also verify because at least the pair (N max , d min , . . . , d min Nmax ) is feasible and any other feasible pair with larger value of N will produce a greater value of the objective function f 2 . Therefore, to have a non-empty set of feasible points we will assume Now, as in the previous subsection, we will focus on the associated continuous problems generated when the value of the integer variable N is fixed: with the functionf N 1 defined in (20). In the following result we show that the optimal treatment for (P N 2 ) has also equal doses and moreover, we obtain the explicit formula to determine them: Proof. First of all, it is immediate to check that the proposed dose vector is feasible for the problem (P N 2 ), if N verifies (27). On the other hand, using once more the relation between the geometric and arithmetic means, the doses of any feasible treatment for (P N 2 ) verify the following inequalities Taking into account that the lower bound of (29) is equal to f N 2 (d), the conclusion follows.
Finally, using (16), we determine the optimal solution of the approximate mixed-integer optimization problem of the palliative approach.
Corollary 2. Let us assume that d min > 0 and (25). Then a global optimal solution of (P 2 ) problem is given by (N ,d) with Proof. This can be deduced combining Theorem 3 with the fact that the auxiliary function ϕ 2 (x) = x(e 1/x − 1) is strictly decreasing in (0, +∞) and for N ≤ N max −1 and using that (d min , . . . , d min

Nmax
) is the solution of (P Nmax 2 ).

Remark 2.
i) Let us note that both possibilities presented in Corollary 2 can occur in practical situations: for example, it can be observed in Table 5 that option a) is true for d min ∈ {50, 100}, while the option b) holds for d min ∈ {75, 150}.
ii) For the trivial cases we recover the expected solutions. For instance, when L * → θ − , we derive that eT R → 0 and the nonlinear constraint is always verified and therefore the solution of (P 2 ) is given by (N ,d) = (1, d min ), as one can expect from the beginning.
iii) We can also observe that N max is a decreasing function of d min ; in fact, N max → +∞ when d min → 0. Moreover, in case a), we have and d * ≈ d min for N max sufficiently large. From here we deduce that case a) is the prevalent one when d min is small enough.
iv) Once more, we can now study the dependence of the optimal cumulative dose with respect to d min . Taking (30) into account and previous remark, we conclude that and to decrease f 2 (N ,d) we have to take d min as small as possible. Again, this is in line with the principles of MC. Table 4 and the plot at the bottom of Figure 1 are devoted to the palliative approach with TMZ. We have considered the threshold L = 0.1813θ which is associated to the size of the tumor at the final time T with the best treatment of the curative approach (see Table 2). Let us remember that in this case we are considering treatments with a maximum of 40 doses associated to the 5/28d schedule, d min = 100 mg/m 2 and the parameter values of Table 1. Moreover, Table 4 presents two kind of optimal solutions. For each value of N (number of doses of the treatment) the following three columns are devoted to problems (P N 2 ) (that use (16)) and the last four columns correspond to optimal solutions obtained with FMINCON for problems (P N 2 ), setting the value for N and without using that approximation. Columns with labelsd min andd max show the minimum and maximum doses of the treatment given by FMINCON, respectively. Let us note that both values in each row are almost equal. Moreover, for each value of N , the results for the problems (P N 2 ) are good approximations of the results obtained with FMIN-CON. Third and seventh columns (with labels f N 2 and (f N 2 ) opt ) present the total dose associated to each optimal solution. We check that the longest treatment is the best from the point of view used in this subsection. Finally, let us note that the total dose of the approximate solution (5,749.24 mg/m 2 ) is less than the FMINCON optimal value (5,749.95 mg/m 2 ), but the constraint on the tumor volume is not satisfied (because the tumor size level is slightly exceeded). The plot at the bottom of Figure 1 shows the Table 4: Results for (P 2 ) and TMZ with L = 0.1813θ and d min = 100 mg/m 2 with schedule 5/28d. Last but not least, in Table 5 and at the bottom of Figure 2 we illustrate the influence of the level d min . Now, we compare the optimal treatments corresponding to the palliative approach. Again, first we have fixed d min and solved (P 1 ); once we have determined the corresponding solution (N and {d i }N i=1 ) we have selected the more appropriated schedules used for TMZ. All the treatments in the table keep the tumor below the given L level, and the longest treatment (that also has the lowest doses) is the best because it has the lowest cumulative dose. These results are also consistent with MC. The corresponding tumor volumes, L(t)/θ , for a medium-sized initial tumor (25% of carrying capacity) with the treatments in Table 5 are shown in Figure 2.
The computations were performed on a 3.9 GHz Core i5-8265U machine, with 8 GB RAM running under the 64-bit version of Windows 10 and MAT-LAB R2019b.

Conclusions
Throughout this work we have given a mathematical justification supporting metronomic chemotherapy (MC) as the best option for most cytotoxic drugs (i.e. when the main hypothesis (15) holds) and for curative and palliative approaches in oncology. First, we have obtained explicit expressions of global optimal solutions of (P 1 ) and (P 2 ) problems (given by Corollary 1 and Corollary 2) that are approximations of (P 1 ) and (P 2 ), respectively. Once more, we emphasize that the administration times do not appear in the formulation of (P 1 ) and (P 2 ) and for both cases the optimal solutions lead to the longest feasible treatments with equal doses, which are in line with MC. We further deduce mathematically that it is convenient to take the effectiveness level d min as small as possible, but above MED (i.e. ensuring a minimum effect!). Numerical experiments with TMZ support the high quality of the solutions of (P 1 ) and (P 2 ) with respect to those of (P 1 ) and (P 2 ), respectively, and illustrate the theoretical results. Moreover, the optimal treatments calculated for both curative and palliative approaches are very similar. Finally, with our method the administration times are adjusted a posteriori (once the number of doses and the individual doses of the treatment are determined). In our work we have used a simple mathematical model being crucial the Norton-Simon hypothesis, the Emax model, the PK of the drug given by (6) and the main hypothesis (15). In the future, we would like to analyze several generalizations and extensions from a mathematical point of view, including some relevant factors such as drug resistance and anti-angiogenesis effect or the case when (15) fails. We think that both theoretical and numerical results would be very useful for clinical purposes.