A MATHEMATICAL MODEL OF CHEMOTHERAPY WITH VARIABLE INFUSION

. A nonautonomous mathematical model of chemotherapy cancer treatment with time-dependent infusion concentration of the chemotherapy agent is developed and studied. In particular, a mutual inhibition type model is adopted to describe the interactions between the chemotherapy agent and cells, in which the chemotherapy agent is modeled as the prey being consumed by both cancer and normal cells, thereby reducing the population of both. Properties of solutions and detailed dynamics of the nonautonomous system are investigated, and conditions under which the treatment is successful or un- successful are established. It can be shown both theoretically and numerically that with the same amount of chemotherapy agent infused during the same pe- riod of time, a treatment with variable infusion may over perform a treatment with constant infusion.


1.
Introduction. Cancer is still one of the leading causes of death worldwide (see, e.g., [3,5,18]). Conventional methods used in cancer treatment include chemotherapy, immunotherapy, radiotherapy and surgery, etc. 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. While being easily applicable, most of the chemotherapy agents attack not only the cancer cells but also other fast-renewing tissues such as skin, bone-marrow, gut, and other digestive epithelia (see, e.g., [1,10,22,23]). This motivates both theoretical and experimental studies to better understand the trade-offs between reducing cancer cells and impacting healthy cells.
In particular, mathematical models have been used extensively to study the effectiveness of chemotherapy treatments, from dynamical point of view, optimization point of view, and compartmental point of view (see, e.g., [1,2,6,11,13,16,19,20,21,22,26]). Since the chemotherapy agents and cells have negative effects on each others' growth rates, their interactions are of the mutual inhibition type [24]. In the context of chemotherapy treatments, such interactions can be understood as a "predator-prey" type relation that the "predator" has a negative growth rate by consuming the "prey", as if the prey is poisoned. Thus, among various mathematical models, the predator-prey type of systems with mutual negative affects have been adopted to describe the interaction between the chemotherapy agent and cells. There are two perspectives in modeling the chemotherapy treatment as a predatorprey system: the first type is to model the chemotherapy agent as the "predator" that kills both normal and cancer cells (see, e.g., [13,22]), and the second type is to model the chemotherapy agent as the "prey" that is consumed by both normal and cancer cells (see, e.g., [1,25]).
For the first type where the chemotherapy agent is considered as the predator, an autonomous model of the interaction between the normal cells and tumor cells with metastasis and time delays was studied in [22], and a nonautonomous model of interactions among tumor cells, normal cells and the chemotherapy agent under time varying environmental conditions was later developed and studied in [13]. It was demonstrated in [13] that treatments with time-dependent infusion of the chemotherapy agent can be more effective than treatments with constant infusion. The goal of this work is to investigate dynamics of a nonautonomous chemotherapy model of the second type, in which the chemotherapy agent is modeled as a prey, and infused into the treatment site with time-dependent infusion. The autonomous counterpart of this model with or without delay was recently studied by Abdulrashi et al. in [2], in which both types of autonomous models were analyzed and compared. Yet no result on the second type nonautonomous system is available up to date. This is the focus of this work.
The paper is organized as follows. In Section 2 we first formulate the underlying nonautonomous dynamical system with parameters of physical meanings and then simplify it to a dimensionless system. In Section 3 we first study the existence and uniqueness of a non-negative and bounded solution of the model and then introduce the basic terminology and theory of nonautonomous dynamical system related to this work. In Section 4 we investigate detailed dynamics of the nonautonomous model, in terms of the existence of a pullback attractor, as well as detailed dynamics within the attractor. Biological interpretations of the resulting mathematical results are also provided. In the end numerical simulations are presented to illustrate the theoretical results.
2. Mathematical model. The model to be developed and studied is based on the idea introduced in [25] and later studied in [1,2], but with time-dependent infusion due to the natural (temporal or random) fluctuation of environments or human control.
2.1. Model formulation. Consider a single site where the cells are treated, e.g., a tumor, with fixed volume V . It is assumed that all cells, as well as the chemotherapy agent, are spatially uniform within the site, i.e., their concentrations do not depend on the location. At any time t denote by N 1 (t), N 2 (t) and C(t) be the concentration of cancer cells, normal cells, and the chemotherapy agent at the treatment site, respectively. Let F in = F out = F be the blood flows brought into and coming out from the tumor site at any time. The novelty and focus of this work is that the chemotherapy is assumed to be infused with blood flow at time-dependent concentration. More precisely, denote by I(t) the concentration of the chemotherapy agent in the blood flowing into the site, where I(t) is a continuous, positive and bounded function that varies with time deterministically or randomly.
Using the idea of [25], the negative effect of the chemotherapy agent on the growth of cells is modeled by a "kill rate" K j (C) (j = 1, 2 for cancer and normal cells, respectively), and the chemotherapy agent is regarded as the "prey" being consumed by both types of cells at rates proportional to the kill rates. In addition, assume that the normal and cancer cells both follow a logistic growth [9,13,17,22] and have Lotka-Volterra type intra-specific competitions between them [8]. These leads to the following nonautonomous system of ODEs describing dynamics of chemotherapy Note that the key difference between the model above and autonomous models in the literature is that the input concentration I is time-dependent. In addition, the difference between the model above and the nonautonomous model studied in [13] lies in that the killing rates K 1 and K 2 are functions of C instead of functions of N 1 and N 2 . More precisely, the functions −r 1 K 1 (C)N 1 and −r 2 K 2 (C)N 2 can be regarded as the interactions that create a positive feedback on both variables in the mutual inhibition relation between the chemotherapy agents and the cells. Meanings and units of parameters b 1 , b 2 , κ 1 , κ 2 , r 1 , r 2 , d 1 and d 2 are listed in Table 1 below. Throughout this paper we adopt the Michaelis-Menten formulation of the killing rates [25]: where K max j is the maximum killing rate of the chemotherapy agent on the cells, and k half j is the concentration of cells corresponding to K j (C) = K max j /2, which is usually referred to as the half saturation rate. Note that K max j is a rate, and has units 1/time and k half j has units of concentration.

Assumptions. By the physical meanings of parameters listed in
3. Properties of solutions. In this section, we first investigate basic properties of solutions to the system (4) -(6) including existence, uniqueness, boundedness and non-negativeness of the solution. We then provide a basic introduction on concept and theory of nonautonomous dynamical systems required in the sequel.
3.1. Basic properties of solutions. In this subsection we prove that system (4) -(6) has a unique global solution under the initial condition Moreover, we will prove that the solution is non-negative and bounded for all time t ≥ t 0 . For convenience, write u(t) := (x(t), y(t), z(t)) and u 0 = (x 0 , y 0 , z 0 ).
Lemma 3.1. The ODE system (4) -(6) with initial condition (9) has a unique bounded solution u(t; Proof. First it is straightforward to rewrite (4) -(6) as the following ODE on R 3 , Since µ(t) is both continuous and bounded, function g is continuous in t and locally Lipschitz in u. It then follows immediately from the classical theory of ODEs (see, e.g., [12]), that equation (10) has a unique local solution u(t; i.e., the positive quadrant R 3 + is positively invariant for u. Therefore by continuity of solutions, any solution trajectory that starts from u 0 ∈ R 3 + at t 0 will stay nonnegative for all It then follows immediately that Moreover, by using Assumption (A1) we have which implies that The inequalities (12) and (13) and the existence of local solutions, together imply that given any initial condition u 0 = (x 0 , y 0 , z 0 ) ∈ R 3 + the equation (10) has a unique solution defined for all t ≥ t 0 and remains in the bounded region The proof is complete.

Preliminaries on nonautonomous dynamical systems.
In this subsection we provide introductory material of nonautonomous dynamical systems (see, e.g., [4,7,14,15]) required in the sequel. In particular, we will introduce the process formulation of nonatuonomous dynamical systems and concepts and theory on pullback and forward attractors. Denote by (t, t 0 ) ∈ R 2 ≥ , which satisfies (i) initial value property: ϕ(t 0 , t 0 , u) = u for all u ∈ R d and any t 0 ∈ R; (ii) two-parameter semigroup property: for all x ∈ R d and (t 2 , t 1 ), The existence of a pullback attractor follows from that of a pullback absorbing family, which is usually more easily determined.
The proof of the following proposition is well known, see e.g., [15]. Proposition 1. Suppose that a process ϕ on R d has a ϕ-positively invariant pullback absorbing family Λ = {Λ(t) : t ∈ R} of nonempty compact subsets of R d . Then ϕ has a unique global pullback attractor A = {A(t) : t ∈ R} with its component sets determined by If Λ is not ϕ-positively invariant, then 4. Dynamics of the nonautonomous chemotherapy model. First of all, due to the existence and uniqueness of a global solution to the system (4) -(6), we can define a process {ϕ(t, t 0 ) where u(t; t 0 , u 0 ) is the solution of (4) -(6) with the initial condition u(t 0 ) = u 0 . Moreover, it is straightforward to check that the process defined above is continuous and hence all concepts and theory introduced in the subsection 3.2 can be applied.
In what follows, we first establish the existence of a pullback attractor, and then investigate detailed structures of the attractor and provide their biological insights.

Existence of pullback attractors.
In this subsection we first construct a positive invariant absorbing set for the process {ϕ(t, t 0 )} (t,t0)∈R 2 ≥ defined in (14), stated in the Lemma below.
The following theorem follows directly from Proposition 1.

Remark 4.3.
Notice that the estimations (18) -(25) hold both forwardly and pullback, i.e., for t 0 fixed with t → ∞, as well as for t fixed with t 0 → −∞. The set Λ is both a pullback absorbing set and a forward absorbing set. Although this does not necessarily ensures the existence of a forward attractor (see, e.g., [7]), it can still be used to investigate forward dynamics of the system.

4.2.
Detailed dynamics within the attractor. Theorem 4.2 provides the existence of a pullback attractor for the process {ϕ(t, t 0 )} (t,t0)∈R 2 ≥ defined by the solution of system (4) - (6). In fact, since Λ is ϕ-positively invariant, the component subsets of the attractor A are defined by In this subsection we investigate detailed structure of A , with both mathematical and biological interpretations.
Then the pullback attractor A has a singleton component subset Proof. First note that dx dt x=0 = 0. Then for any x > 0, using the lower bound of z in (25), we have z It then follows immediately from the assumption (26) that is negative definite. Thus the x component of all trajectories in the nonnegative quadrant R 3 + approaches 0 asymptotically. Similarly, the y component of all trajectories in the nonnegative quadrant R 3 + all approaches 0 asymptotically provided α 2 z m > β 2 , which is equivalent to the assumption (27).
With x(t) = 0 and y(t) = 0, the equation (6) becomes which can be solved to get The proof is complete.
Remark 4.5. The singleton trajectory z * (t) is obtained by fixing z 0 and letting t 0 approach −∞. Notice that the chemotherapy agent does not exists until the treatment starts, thus z 0 = 0 for t 0 < 0 and µ(t) = 0 for t < t 0 . While it seems that z * (t) then depends on the starting time t 0 , it is in fact a function of t dependent on the definition of µ(t) which is given.
Theorem 4.7. Assume that (27) holds and Then the pullback attractor contains points inside the strictly positive subspace {(x, y, z) ∈ R 3 + : x > 0, y = 0, z > 0}. Proof. We look at the derivative of x(t) at any ε < 1/γ 1 . Using Lemma 4.1 In particular picking ε ≤ α1θ1 β1γ1(θ1+2µ M ) , then under the assumption (30), which implies that dx(t) dt x=ε > 0. Thus x(t) ∈ [ε, 2/γ 1 ] for all t ≥ t 0 and the attractor contains points inside {(x, y, z) ∈ R 3 + : x ≥ ε, y = 0, z > 0}. The proof is complete. (26) is satisfied and all normal cells will die out if the assumption (27) is satisfied. The assumption (26) is equivalent to b 1 < K max 1 r 1 , and the assumption (27) is equivalent to b 2 < K max 2 r 2 where r 1 and r 2 can be thought of as a portion of the maximal killing rate on the cancer and normal cells, respectively. More importantly r 1 and r 2 depend on the minimum infusion concentration min t≥t0 I(t).

Biological interpretations. Theorem 4.4 says that all cancer cells will die out if the assumption
Theorem 4.6 provides sufficient conditions for a successful treatment, i.e., all cancer cells are killed but normal cells still remain. The assumption (29) is equivalent to b 2 > K max 2 + 2d 2 κ 1 , which means that the per capital birth rate of normal cells has to be large enough to cover the maximal killing rate of the chemotherapy agent on the normal cells and twice the intra-specific competition created by all cancer cells carried by the environment. In the special case where d 2 = 0, this reduces to b 2 > K max 2 only. Theorem 4.7 provides sufficient condition for a failed treatment, i.e., all normal cells are killed but cancer cells are remaining. The assumption (30) is equivalent to b 1 > K max 1 + 2d 1 κ 2 , which means that the per capital birth rate of cancer cells is even larger than the maximal killing rate of the chemotherapy agent on the cancer cells and twice the intra-specific competition created by all normal cells carried by the environment. In the special case where d 1 = 0, this reduces to b 1 > K max 1 only. It is implied by the theoretical results above that the success or failure of a chemotherapy treatment is mostly determined by the relations between the per capita growth rate of cells, the maximum killing rate, i.e., effectiveness of the chemotherapy agent on cells. The carrying capacity of cells also affect the results, but according to the strength of intra-specific competitions.
However, it is worth mentioning that after a closer look at the computations in the proof of Theorem 4.6, the assumption (29) can be weakened to which is equivalent to This means that for the normal cells to remain while all cancer cells are cleared, the per capita growth rate of normal cells does not really need to be much larger than the maximum killing rate of the agent on the normal cells. In fact, it only needs to be faster than a percentage R 2 of the maximum killing rate on the normal cells, which is determined by the relation between the maximum input concentration of the chemotherapy agent and the half saturation concentration of the consumption function of normal cells.
Similarly, for the cancer cells to remain while all normal cells die, the per capita growth rate of cancer cells does not really need to be much larger than the maximum killing rate of the agent on the cancer cells. In fact, it only needs to be faster than a percentage R 1 of the maximum killing rate on the cancer cells, which is determined by the relation between the maximum input concentration of the chemotherapy agent and the half saturation concentration of the consumption function of cancer cells.
These bring in the effect of control on the input concentration I(t), as well as a major difference between nonautonomous and autonomous models.

4.4.
Comparison to the autonomous counterpart. For comparison purpose, we analyze the autonomous counterpart of the system (4) -(6), in which µ(t) ≡μ. In particular, we exam the sufficient conditions for a successful treatment and a failure treatment and compare to the nonautonomous results. For reader's convenience, we state the autonomous system below.
Note that all computations in Lemma (4.1) still hold for the above system, and hence we can also focus our attention on the positive invariant set (34) Recall that a major difference between autonomous and nonautonomous systems is that solutions of autonomous systems depend only on the time elapsed, t − t 0 , while solutions of nonautonomous sytems depend on both t 0 and t. In general, nonautonomous systems do not possess constant equilibria as autonomous systems do. But there may exist entire trajectories of nonautonomous systems which can be regarded as the time-dependent counterpart of equilibria for autonomous systems. For example, (0, 0,μ) is one equilibrium for the autonomous system (31) -(33) that is asymptotically stable under the assumptions (26) and (27), while (0, 0, z * (t)) is an entire trajectory of the nonautonomous system (4) -(6) that attracts all other solutions under the assumptions (26) and (27).