STABILITY ANALYSIS OF A CHEMOTHERAPY MODEL WITH DELAYS

. A chemotherapy model for cancer treatment is studied, where the chemotherapy agent and cells are assumed to follow a predator-prey type rela- tion. The time delays from the instant that the chemotherapy agent is injected to the instant that the treatment is eﬀective are taken into account and dynam- ics of systems with or without delays are compared. First basic properties of solutions including existence and uniqueness, boundedness and positiveness are discussed. Then conditions on model parameters are established for diﬀerent outcomes of the treatment. Numerical simulations are provided to illustrate theoretical results.

1. Introduction. Unarguably, cancer is still one of the deadliest disease to date (see, e.g., [1]). Chemotherapy is a well-known and commonly used method of cancer treatment, that involves the application of a chemotherapy agent to the body of the infected individual thereby attacking the cancerous cells with the aim of curative intent, control, or palliation. However, most of the chemotherapy agents attack not only the infected cells but also the non-infected normal cells and hence may cause severe side effects to the body of the patient (see, e.g., [5,14,8]). For this reason, both theoretical and experimental research on chemotherapy cancer treatments has been carried out to maximize the effectiveness of chemotherapy on tumor cells while minimizing the damage to normal cells [5].
Among theoretical studies on chemotherapy treatments, mathematical modeling, analysis and computation have been used extensively. In particular, chemotherapy treatments have been investigated from dynamical, optimization, and compartmental point of view (see, e.g., [12,13,11,6,2,14,18,8]). For example, an autonomous dynamical system modelling the interaction between the normal cells and tumor cells with metastasis and time delays in the effect of the metastasis was studied by [14]. Based on the idea of [14], a nonautonomous dynamical system was developed and studied in [8], that models the interactions among tumor cells, normal cells and the chemotherapy agent under time varying environmental conditions.
Note that the effect of chemotherapy treatment is not always instantaneous, i.e., there may be a time delay from the moment that the chemotherapy agent is injected to the moment that the cells are killed by the agent. In this work we develop a dynamical model similar to the models in [14] and [8], but with the inclusion of time delays as well as some other modifications. In particular, to focus on the effect of delays, an autonomous dynamical system without the consideration of varying environments will be considered. The aim is to investigate the stability of each meaningful steady state of the chemotherapy model with constant delays, and compare with its analog model without delays.
The paper is organized as follows. In Section 2 we develop the dynamical systems to be studied and introduce the meaning of each parameter. In Section 3 we discuss basic properties of solutions, including the existence and uniqueness, boundedness, and positiveness. In addition, steady state solutions will be calculated and classified. In Section 4 we investigate stability of each steady state solution obtained in Section 3 and provide intuitive biological interpretation. Illustrative numerical simulations are provided in Section 5, and some closing remarks are given in Section 6.
2. Model formulation. The quantities of interests are concentrations of cancer cells, normal (un-infected) cells, and the chemotherapy agent, at a single site for treatment in the body. The tumor site is assumed to be spatially uniform, i.e., the concentrations are the same everywhere over the site. For any time t ≥ 0 let x(t) be the concentration of cancer cells, y(t) be the concentration of normal cells, and z(t) be the concentration of the chemotherapy agent at the treatment site.
The model to be developed is based on those in [14,8], but with key modifications. First, in [14] and [8], the normal and cancer cells are considered as the "predators" that consume the chemotherapy agent and the predation was modeled by consumption functions U (x(t)) and U (y(t)). Here we follow the idea given in [17] and model the cells as the "preys" that are consumed by the chemotherapy agent. In particular, the negative effect of the chemotherapy agent on the growth of cells are described by an uptake function U (z(t)). Similar to [14,8], the normal and cancer cells are both assumed to follow a logistic growth [3,10,14,8] and have Lotka-Volterra type intraspecific competitions between them [4]. The resulting system of ODEs describing dynamics of chemotherapy reads where meanings of the parameters are summarized in Table 1.
In addition to the modifications in the relations between the chemotherapy agent and cells described above, we consider another important factor in the model: time delay. When the negative effect of the chemotherapy agent on the cells are not instantaneous, there exist time delays from the injection of the chemotherapy agent to the loss of the cells. To model such delays, we modify the terms describing the negative effect of the chemotherapy agent on the cells, −αx(t)U (z(t)) and −αy(t)U (z(t)) to −αx(t)U (z(t − τ 1 )) and −αy(t)U (z(t − τ 2 )), respectively, where τ 1 and τ 2 are the time delays in the effect of chemotherapy agent on the cancer cells and normal cells, respectively. The system (1) -(3) then becomes Assume that the above system is subjective to the initial conditions: where τ := max{τ 1 , τ 2 } > 0 and φ is a continuous and nonnegative function on [−τ, 0]. We aim to study the stability of each meaningful steady state of the system (4) -(6) with the initial condition (7), and investigate the effect of delays on the stability. Throughout the paper when explicit calculations are needed it is assumed that the uptake function U takes the Michaelis-Menten or Holling type-II form given by where θ is the half saturation constant (see, e.g., [15,17]). Clearly 0 ≤ U (z) ≤ 1 for every z ≥ 0, and U is an increasing function of z.
3. Basic properties of solutions. In this section we discuss the existence and uniqueness, as well as boundedness and positiveness of solutions to the system (4) -(6) with initial conditions (7). To this end, let and define F : Then the vector fields of system (4) -(6) can be written as and system (4) -(6) along with the initial conditions (7) can be written as Remark 1. The formulation (8) via the definition of F appears to cumbersome, but will facilitate the stability analysis in the sequel.
Throughout the rest of this paper, denote Proof. First, since F (u, v, w) and F u (u, v, w) are both continuous on R 3 and Φ : [−τ, 0] → R 3 + is also continuous, then it follows immediately from standard theory of delay ODEs (e.g., [7,16,9]) that the initial value problem (8) has a local solution To show that the local solution is in fact global, we first verify that u(t; Φ) ∈ R 3 + for all t ∈ [−τ, T ]. In fact, this follows immediately by the continuity of u(t; Φ) and observing that dx dt x=0,y,z≥0 = 0, dy dt y=0,x,z≥0 = 0, dz dt z=0,x,y≥0 = DI > 0.
As a direct consequence, and by a comparison principle Using again the facts that x, y ≥ 0 and U ≥ 0, we have which implies that The inequalities (9) and (10) together shows that the solution u(t; Φ) exists for all t ≥ 0 and is non-negative. Moreover, notice that the upper bounds for x, y, z are independent of t, therefore the solution u(t; Φ) is bounded for every t ≥ 0.
Remark 2. Calculations in the above proof also apply to the system (1) -(3) without delay. Thus the system (1) -(3) also has a unique bounded and nonnegative global solution given any non-negative initial condition.
The steady state solutions of the system (4) -(6) can be classified into four categories: namely the axial state, the preferred state, the failure state, and persistent state, respectively. For clarity of exposition, we will focus on the stability of the axial, the preferred, and the failure state. But note that the same set of techniques can also be applied to analyze the persistent state with more complex computations.
More precisely, a preferred state satisfies and a failure state satisfies Notice that equations (11) are equivalent to with z * satisfying the cubic equation and hence have at least one positive solution since a 3 < 0. Similarly, equations (12) also have at least one positive solution.
4. Stability analysis. In this section we investigate the stability of the axial, preferred and failure states. For comparison purpose, we first discuss the stability of each state when there is no time delay. Throughout this section denote by u * = (x * , y * , z * ) a generic steady state solution of the system (4) - (6), as well as the system (1) (ii) a preferred steady state (0, y * , z * ) is asymptotically stable provided (iii) a failure steady state (x * , 0, z * ) is asymptotically stable provided Proof. (i) When x * = y * = 0 and z * = I, the Jacobian matrix is reduced to with eigenvalues which are all negative under the assumption (13).
(ii) When x * = 0 and y * and z * satisfy equation (11), the Jacobian matrix is reduced to The first eigenvalue of the above matrix is under the assumption (14). The other two eigenvalues of the above matrix coincide with the eigenvalues of the submatrix First notice that by the assumption (15). Therefore all eigenvalues of J are negative, which implies the asymptotic stability of the preferred state (0, y * , z * ).
(iii) When y * = 0 and x * and z * satisfy equation (12), the Jacobian matrix is reduced to The first eigenvalue of the above matrix is under the assumption (16). The other two eigenvalues of the above matrix coincide with the eigenvalues of the submatrix First notice that by the assumption (17). Therefore all eigenvalues of J are negative, which implies the asymptotic stability of the preferred state (x * , 0, z * ).  (14) and (15) are sufficient conditions for the asymptotic stability of the preferred steady state, but the assumption (14) is also a necessary condition for the asymptotic stability of the preferred steady state. Similarly, assumptions (16) and (17) are sufficient conditions for the asymptotic stability of the failure steady state, but the assumption (16) is also a necessary condition for the asymptotic stability of the failure steady state.
Notice that assumption (13) is equivalent to Thus part (i) of Theorem 4.1 indicates that when the killing rate of the chemotherapy agent is higher than a certain threshold (dependent of the largest growth rate of cells, the input concentration of the chemotherapy agent, and the half saturation constant), then both normal and cancer cells will be cleared. The assumptions (14) and (15) together are equivalent to Thus part (ii) of Theorem 4.1 indicates that when the killing rate of the chemotherapy agent is well controlled between β1−δy * U (z * ) and β2 2U (z * ) , the treatment will be successful. This implicitly requires the steady state to satisfy β 2 ≥ 2(β 1 − δy * ), which essentially put a restriction on the growth rates of cancer and normal cells. Similarly part (iii) of Theorem 4.1 indicates that when the killing rate of the chemotherapy agent lies between β2−δx * U (z * ) and β1 2U (z * ) , the treatment will fail. This also implicitly requires the steady state to satisfy β 1 ≥ 2(β 2 − δx * ).
Notice that the conditions in parts (ii) and (iii) depend on the magnitude of each specific steady state. For example, a preferred steady state with smaller quantity of normal cells is easier to achieve than a preferred steady state with larger quantity of normal cells, and a failure steady state with smaller quantity of cancer cells is easier to achieve than a failure steady state with larger quantity of cancer cells.

Constant time delay. First we linearize the system (4) -(6) about
where the coefficient matrices The characteristic equation then reads (i) the axial steady state (0, 0, I) is asymptotically stable provided (ii) a preferred steady state (0, y * , z * ) is asymptotically stable provided (iii) a failure steady state (x * , 0, z * ) is asymptotically stable provided Proof. (i) When x * = y * = 0 and z * = I, the equation (19) becomes The three simple solutions to the characteristic equation (25) are They are all negative when (20) holds, which implies immediately that the axial steady state (0, 0, I) is asymptotically stable.
(ii) When x * = 0 the equation (19) becomes where y * and z * satisfy the equations (11). Equation (26) is equivalent to where First notice that when the assumption (21) holds, by (11) the solution to the equation (27) satisfies Then notice that the polynomial p(λ) has two roots: By using (11) again we have In addition by the assumption (22) and (11), and consequently It then follows from Corollary 4.10 of [16] that Reλ < 0 for every root λ of the equation (28), and together with (29) we conclude that the steady state (0, y * , z * ) satisfying (11) is asymptotically stable.
(iii) When y * = 0 the equation (19) becomes where x * and z * satisfy equations (12). Equation (30) is equivalent to where Under the assumption (23), using (12) we have the solution of equation (31) satisfies Notice that the polynomial p(λ) has two roots, where by using (12) again we have By the assumption (24) and (12), and consequently It then follows from Corollary 4.10 of [16] that Reλ < 0 for every root λ of the equation (32), and together with (33) we conclude that the steady state (x * , 0, z * ) satisfying (12) is asymptotically stable.  (23) and (17) are sufficient conditions for the asymptotic stability of the failure steady state, but the assumption (24) is also a necessary condition for the asymptotic stability of the failure steady state.
Part (i) of Theorem 4.2 has the same interpretation as part (i) of Theorem 4.1. Since the conditions imposed are necessary and sufficient conditions, we obtain an important information that the stability of the axial steady state is not affected by the time delay of the treatment.
To further interpret these assumptions, we rewrite (21) as and rewrite (23) as Then parts (ii) and (iii) of Theorem 4.2 provide two interpretations depending on the strength of inter-specific competitions. When the inter-specific competition is strong, in the sense that δ > κ 2 /β 2 and δ > κ 1 /β 1 , then assumptions (34) and (22) can be combined to be and assumptions (35) and (24) can be combined to be One special case of (36) happens when cancer cells grow slower than normal cells, i.e., β 1 < β 2 , then the treatment will be successful as long as y * is more than half of the environmental carrying capacity of normal cells. Similarly one special case of (37) when cancer cells grow faster than normal cells, i.e., β 1 > β 2 , then the treatment will be a failure as long as x * is more than half of the environmental carrying capacity of cancer cells.
On the other hand, when the inter-specific competition is weak, in the sense that δ < κ 2 /β 2 and δ < κ 1 /β 1 , then assumptions (34) and (22) can be combined to be which may only happen when β 2 > β 1 and β2−β1 β2−κ2δ > 1/2. When the inter-specific competition is negligible, this implies that the normal cells need to grow at least twice as fast as the cancer cells do for a possible successful treatment. Similarly assumptions (35) and (24) can be combined to be which may only happen when β 1 > β 2 and β1−β2 β1−κ1δ > 1/2. When the inter-specific competition is negligible, this implies that the cancer cells need to grow at least twice as fast as the normal cells do for a possible failed treatment.
5. Numerical simulations. In this section we include some numerical results with parameters satisfying the conditions constructed in the previous section. 5.1. Axial steady state. The first set of parameters simulated satisfy the sufficient and necessary conditions for the axial steady state (0, 0, I). In particular, the parameters are chosen as that satisfy the assumptions in (13). The initial conditions are set to be x 0 = 5, y 0 = 3, φ(t) = 4 for t ∈ [−3, 0]. Evolution of the concentrations of normal and cancer cells, and the chemotherapy agent is illustrated in Figure 1. It can be clearly seen that concentration of both normal and cancer cells tend to 0 as time goes on. To closer exam the effect of delays, we compare the concentration of cancer and normal cells of the above example with the special case when τ 1 = τ 2 = 0, shown in Fig. 2. Though concentrations of both cancer and normal approach zero, with or without delays, differences in cell concentrations between the cases with delays and the cases without delays can be clearly seen. In this particular example, the chemotherapy treatment approaches the axial steady state faster with the presence of delays.  The second set of parameters simulated satisfy the sufficient and necessary conditions for a preferred steady state (0, y * , z * ). In particular, the parameters are chosen as There is only one real positive solution to the equations (11), y * = 5.2422 and z * = 0.7168, that satisfy the assumptions (14) and (15). Evolution of the concentrations of normal and cancer cells, and the chemotherapy agent is illustrated in Figure 3, where it can be seen that the concentration of cancer cells tends to 0 while the concentration of normal cells tends to y * = 5.2422. The treatment is hence successful. It is also interesting to see that although the cancer cells become vanish soon, it takes a much longer time for the normal cells to recover to their steady state.  Similarly we compare the concentration of cancer and normal cells of the above example with the special case when τ 1 = τ 2 = 0 in Fig. 4, where differences in cell concentrations for the cases with delays and the cases without delays can be clearly seen. In this particular example, the cancer cells are removed slightly faster when there are delays, whereas the normal cells recover much faster when there are no delays.
Evolution of the concentrations of normal and cancer cells, and the chemotherapy agent is illustrated in Figure 5, where it can be seen that the concentration of normal cells tends to 0 while the concentration of cancer cells tends to x * = 2.728.
A comparison between the concentration of cancer and normal cells of the above example and the special case when τ 1 = τ 2 = 0 is shown in Fig. 6. The difference in cancer cells for the cases with or without delays is clearly seen, whereas the difference in normal cells for the cases with or without delays is not detectible.  6. Closing remarks. We constructed sufficient and necessary conditions for the stability of the axial steady state (0, 0, I), which are indifferent for the system with or without delay. We also constructed sufficient conditions for the stability of the preferred and failure steady states, respectively, which turned out to coincide. The indifference of stability conditions is mainly due to the special bounded structure of the consumption function U , that mitigate the effect of delays. However, we cannot conclude that the time delay does not affect the stability of non-axial steady states. In fact, the numerical experiments presented above clearly show the difference in dynamics of the chemotherapy model with or without delays. Further numerical experiments indicated that stability of preferred or failure steady states that do  In other words, a stable preferred state with y * < κ 2 /2 or a stable failure state with x * < κ 1 /2 when there are no delays, can become unstable with delays, and vice versa. In a subsequent work, conditions on the magnitudes of delays will be established for stability of various states by constructing appropriate Lyapunov functions, which will provide more insights in the effect of delays for chemotherapy treatments.