TUMOR GROWTH DYNAMICS WITH NUTRIENT LIMITATION AND CELL PROLIFERATION TIME DELAY

. It is known that avascular spherical solid tumors grow monoton-ically, often tends to a limiting ﬁnal size. This is repeatedly conﬁrmed by various mathematical models consisting of mostly ordinary diﬀerential equations. However, cell growth is limited by nutrient and its proliferation incurs a time delay. In this paper, we formulate a nutrient limited compartmental model of avascular spherical solid tumor growth with cell proliferation time delay and study its limiting dynamics. The nutrient is assumed to enter the tumor proportional to its surface area. This model is a modiﬁcation of a recent model which is built on a two-compartment model of cancer cell growth with transitions between proliferating and quiescent cells. Due to the limitation of resources, it is imperative that the population values or densities of a population model be nonnegative and bounded without any technical conditions. We conﬁrm that our model meets this basic requirement. From an explicit expression of the tumor ﬁnal size we show that the ratio of proliferating cells to the total tumor cells tends to zero as the death rate of quiescent cells tends to zero. We also study the stability of the tumor at steady states even though there is no Jacobian at the trivial steady state. The characteristic equation at the positive steady state is complicated so we made an initial eﬀort to study some special cases in details. We ﬁnd that delay may not destabilize the positive steady state in a very extreme situation. However, in a more general case, we show that suﬃciently long cell proliferation delay can produce oscillatory solutions.


(Communicated by Yuan Lou)
Abstract. It is known that avascular spherical solid tumors grow monotonically, often tends to a limiting final size. This is repeatedly confirmed by various mathematical models consisting of mostly ordinary differential equations. However, cell growth is limited by nutrient and its proliferation incurs a time delay. In this paper, we formulate a nutrient limited compartmental model of avascular spherical solid tumor growth with cell proliferation time delay and study its limiting dynamics. The nutrient is assumed to enter the tumor proportional to its surface area. This model is a modification of a recent model which is built on a two-compartment model of cancer cell growth with transitions between proliferating and quiescent cells. Due to the limitation of resources, it is imperative that the population values or densities of a population model be nonnegative and bounded without any technical conditions. We confirm that our model meets this basic requirement. From an explicit expression of the tumor final size we show that the ratio of proliferating cells to the total tumor cells tends to zero as the death rate of quiescent cells tends to zero. We also study the stability of the tumor at steady states even though there is no Jacobian at the trivial steady state. The characteristic equation at the positive steady state is complicated so we made an initial effort to study some special cases in details. We find that delay may not destabilize the positive steady state in a very extreme situation. However, in a more general case, we show that sufficiently long cell proliferation delay can produce oscillatory solutions.
1. Introduction. In order to understand the possible mechanisms for the Gompertz curve's surprising close fit to clinical tumor data sets, Gyllenberg and Webb [13,17] proposed a two-compartment model of the tumor cells consisting of proliferating and quiescent cells. When the tumor is small, its cell proliferation process takes much longer than that of removing dead cells. For simplicity, they assumed that the dead cells are removed from the tumor instantly. Their key assumption is that in an avascular multicellular tumor spheroid, proliferating cells may become quiescent when a limiting nutrient level is low and get out of quiescence when that limiting nutrient level rebounds.
Let P (t) and Q(t) be the densities of proliferating and quiescent cells, respectively. Define N (t) = P (t) + Q(t) which represents the total amount of cells in a tumor. Gyllenberg and Webb [13] proposed the following plausible model (GWmodel ) to describe the interaction dynamics of proliferating and quiescent cells: with initial conditions, For convenience, the authors made a few general but simplistic assumptions. They assume proliferating cells grow at a constant rate β > 0, and proliferating and quiescent cells die at constant rates µ p ≥ 0 and µ q ≥ 0, respectively. In reality, these rates are most likely depend on the dynamic limiting resource levels. They assume that proliferating tumor cells transit to and from the quiescent compartment at rates in the form of continuous functions r 0 (N ) and r i (N ), respectively. Since starvation due to resource limitation often increase with tumor size and starved cells tend to enter and stay at a quiescent state, the authors made the following assumptions.
(A1): r 0 (N ) ≥ 0 and r 0 (N ) ≥ 0 for all N > 0 and, A key finding of their model assumptions is the diminishing of the growth fraction (the ratio of proliferating cells to the total tumor cells), conforming a phenomenon observed in most tumors. They also performed extensive mathematical analysis of this model and presented some interesting biological implications stemming from their mathematical findings. However, the assumption that proliferating cells maintain a constant proliferating rate of β is a big contrast to well known models such as Gompertz model where the cell growth rate is assumed to be exponentially decreasing due to nutrient limitation [12]. Many well established cell growth models place some constraints on the cell proliferating rate by assuming the rate is a decreasing function of the tumor density. In general, model-based findings are likely to be model specific. To ensure modelbased findings represent robust phenomena, we need to see similar properties appear in closely related plausible models. This objective was pursued in several recent papers. In [25], the author included a discrete time delay in the proliferation term in the GW-model (1.1) and performed stability analysis of the resulting model. In [2], the dead cells are assumed to form a necrotic core and hence the model takes the form of a system of three ordinary differential equations and the motivation was to contrast the dynamics to that of the GW-model (1.1). In addition to including the dead cell population, the authors of [3] also assumed that the nutrient enters the tumor by diffusion through the surface and their objective was to study the nutrient dependent tumor growth dynamics. The main findings of [3] include that the modified model generates naturally bounded solutions and an explicit expression of the tumor final sizes can be established.
In this paper, we would like to formulate yet another alternative model that includes both nutrient limitation and a discrete time delay in cell proliferation. We assume that cell division takes an average of τ units of time.
2. Model formulation. As in [3], we assume that the resource, such as oxygen, that enters the tumor is proportional to its surface area. Therefore, we can assume it takes the form of S(t) = αN (t) θ where θ is a scaling parameter with values in [1/2, 1]. Let R(t) be the limiting resource and u(R(t)) be the proliferating cells' limiting nutrient uptake rate function, then we have Our main objective is to understand how resource limitation and the cell division time delay affect the growth of an avascular tumor. In general, nutrient uptake is a process that takes place much faster than proliferation, state switching and death. Hence, we may apply the usual quasi-steady-state argument (for example, see [11]) on R(t) which yields where we assumed α = 1 for simplicity. If a tumor takes the form of a sphere, then θ = 2/3 (see Bertalanffy [5]) while if it takes the form of a tumor cord (the from of a tumor which grows by accumulating layers of cells enveloping a blood vessel, see Thalhauser et al. [23]), then θ = 1. Following [3], we assume that proliferating cells proliferate at a rate of f (u(R(t))), proliferating cells transition to quiescent cells at the rate of g(u(R(t))). However, due to the persistent nutrient limitation, we assume that the rate of quiescent cells transition to proliferating cells is negligible.
For simplicity, we let One can think r(t) represents the resource ration for proliferating cells.
Since proliferating cells are assumed to enter quiescent state before death, we assume proliferating cells do not die (i.e., µ p = 0 in model (1.1)). All parameters below are nonnegative constants. This may result in the following two-compartment model.
The initial conditions are, Observe that if the tumor initially has only proliferating cells, then in an infinitesimal time, some proliferating cells will enter the quiescent state. In other words, for a very small positive value t, P (t) and Q(t) are both positive if P 0 > 0 but Q(0) = 0. Hence, without loss of generality, we assume Q 0 > 0. We assume that These conditions are slightly more general than those for f (x) and g(x) in [3].
Here are some possible examples of functions f (x) and g(x) Carefully derived forms of f (x) and g(x) are highly desirable for studying the growth dynamics of more realistic tumor growth models.
In the rest of this paper, we assume that θ ∈ (0, 1).
3. Boundedness of tumor size. As pointed out in [3], if l 0 < β in the Gyllenberg-Webb model (1.1), then one can show that P (t) ≥ P (0)e (β−l0)t . Hence, artificial conditions must be placed on l 0 or other parameters to ensure that the tumor size is bounded. However, the sizes of most avascular tumors are naturally limited due to the resource limitation. It is thus desirable to see no or minimum requirements to ensure the boundedness of the population densities in model (2.4). In general, boundedness of a solution is a global and important issue which can have important consequence in establishing global existence of solutions, perform global bifurcation analysis and global stability studies. Even without time delay, the proof of boundedness of solution in [3] is nontrivial, requiring an approach mimicking a Razumikhin type argument only seen in theory of delay differential equations [16]. Our proof of the boundedness of solution for model (2.4) is essentially a modification of that presented in [3]. First, we would like to point out that with positive initial values, the solutions of model (2.4) stay positive. Proof. If at some time the solution turn negative, then there is a first time t 1 > 0 such that one of the solution components becomes zero. Without loss of generality, assume first that P (t 1 ) = 0, a contradiction to the assumption that P (t 1 ) = 0. Similar contradictions can be derived if Q(t 1 ) = 0 for some time t 1 > 0. This completes the proof of the proposition.
Observe that if θ ∈ (0, 1), then the function F (x) = is strictly decreasing for large values of P (t) and tends to zero as P (t) tends to infinity. By the previous lemma and the assumption that θ ∈ (0, 1), and g(0) > ε > 0, we see that there is a P ε > 0 such that if N (0) < P ε and P (t) ≥ P ε , then If the claim is false, then there is a first time t 1 > 0 such that P (t 1 ) = P ε and P (t) < P ε for t ∈ [0, t 1 ). Therefore, we must have P (t 1 ) ≥ 0. By (3.4), we see that r(t 1 ) < δ. We have a contradiction, proving the claim. Let L = Q(0) + (1 + k Q )P ε , we see from Lemma 3.1 that N (t) ≤ L. This completes the proof of the theorem.

4.
Size and structure of a tumor at the positive steady state. The boundedness of solution suggests that the size of a nutrient limited avascular tumor may eventually tend to or oscillate around a positive steady state value. Hence, we would like to consider the existence and locations of nonnegative steady states of the model (2.4). We would also like to understand how the model parameters and functions change the steady state tumor size and structure. Since we do not consider cell mutation in our model, it is natural to expect that (0, 0) is a tumor free steady state, which amounts to say that if there are no proliferating cells initially, then the tumor size stay at zero. In order to ensure that (0, 0) indeed is a steady state for model (2.4), mathematically we need to define f (r(t)) = f (N θ (t)/P (t)) and g(r(t)) = g(N θ (t)/P (t)) at (P (t), Q(t)) = (0, 0). For tiny tumors, N (t) is approximately P (t). Hence Therefore, we define f (r(t)) = f M when (P (t), Q(t)) = (0, 0) and g(r(t)) = 0 when (P (t), Q(t)) = (0, 0). Notice that the model system does not admit a Jacobian at (0, 0) which will prevent us from studying the local stability of model (2.4) at (0, 0) via its linearized system at origin.
Let E * = (P * , Q * ) be a nontrivial steady state of model (2.4). It is easy to see that if such a steady state E * exists, all components of E * must be positive.
Adding together the two equations of the model (2.4) yields At steady state, we must have f (r * )P * − µQ * = 0, and From the Q equation, we have g(r * )P * = µQ * = f (r * )P * , which implies that F (r * ) ≡ f (r * ) − g(r * ) = 0. Therefore, E * = (P * , Q * ) = (P * , f (r * ) µ P * ) exists and is unique. Equation (4.1) also implies that at the positive steady state, the percentage of proliferating cells in the tumor is We summarize our mathematical findings into the following proposition on the steady state tumor size and structure. We wonder what will happen when the time delay presents. The routine approach of studying the local stability of a steady state involves the computation of the Jacobian at the steady state. Unfortunately, r(t) = (P + Q) θ (t)/P (t) is not differentiable at the trivial steady state E 0 = (0, 0). In other words, we do not have a standard Jacobian obtained by partial derivatives of the model rate functions at the origin. Fortunately, from Lemma 3.1, we see that if Q(0) = 0 and P (0) is very small, then for small t > 0, which tends to infinity as P (t) tends to zero. This and the fact that lim r→∞ g(r) = 0 imply that near the trivial steady state E 0 = (0, 0), r(t) is very large. Therefore, near the origin, the model (2.4) dynamics can be approximated by the following decoupled system Hence the stability of the trivial steady state E 0 = (0, 0) is determined by that of The characteristic equation of Equation (5.3) takes the simple form of Observe that for any given τ ≥ 0, G(0) = −f M < 0 and lim λ→∞ G(λ) = ∞. This implies that G(λ) = 0 always has at least one positive root, indicating E 0 = (0, 0) is unstable. We summarize this as the following proposition.
In the following, we consider the local stability of the positive steady state E * = (P * , Q * ). Recall that f (r * ) = g(r * ). We will use the following notations.
f * = f (r * ) = g(r * ), Since r = N θ /P , we see that System (5.6) can be rewritten as The characteristic equation of system (5.7) takes the form of To proceed further, we would like to take advantage of the following results presented in [16].
3. If c 2 > d 2 , and that one of the following two conditions is not true: ; then (5.9) has no imaginary root. In this case, the real part of any characteristic root of (5.9) stays negative for any value τ > 0 if the real parts of the two roots of the characteristic equation (5.9) are negative when τ = 0.
By comparing the characteristic equations (5.8) and (5.9), we can establish the following correspondences . These rather involved expressions suggest that a complete analysis of the number of characteristic roots with positive real part for characteristic equation (5.8) will be complicated. Indeed, our extensive simulation work confirms this. Figure 1 below presents some bifurcation diagrams illustrating the intriguing dynamics exhibited by model (2.4) with f (r) = kr ar+1 , g(r) = c r+m and the cell proliferation time delay τ as the bifurcation parameter. The parameter values are k = 2, a = 1, m = 2, µ = 0.5, c = 1, θ = 1/2. The positive steady state appears to be globally attractive for short time delay but loses its stability for larger values of τ . As cell proliferation time delay increases, tumor size oscillates more noticeably and at a lower level. Notice that the percentage of the average amount of proliferating cells decreases as τ increases.
To avoid the complexity involved in the stability analysis, we consider below two special cases that are mathematically tractable and interesting.
Consider first an extremely special case when f (r * ) = g (r * ) = 0. In this case, we have

The characteristic equation (5.8) becomes
which can be rewritten as This is a critical case 2(i) on page 72 in [16] which concludes that except λ = 0, all other roots of the characteristic equation (5.10) have negative real parts. In other words, in this special case, proliferation delay may not induce instability, or generate periodic solution for model (2.4). Next, we consider a more general but still a rather special case when only g (r * ) = 0. In this case, we have . We claim that in this case, d > c. To prove this claim, we need only to show that r * P µ + f * r * Q > 0. To this end, we evaluate this expression in terms of f * , N * and P * . and g(r) = c r+m using the cell proliferation time delay τ as the bifurcation parameter. The parameter values are k = 2, a = 1, m = 2, µ = 0.1, c = 1, θ = 2/3. The positive steady state appears to be globally attractive for short time delay but lost its stability for larger values of τ . As cell proliferation time delay increases, tumor size oscillates more noticeably and at a lower lever. Notice that the percentage of the average amount of proliferating cells decreases as τ increases.
6. Discussion. A given avascular tumor is likely to have a final size due to the limitation of nutrient supply. A typical model of avascular tumor growth takes the form of a system of nonlinear autonomous ordinary differential equations (ODE).
Mathematically, the final tumor size, if exist, is the steady state tumor size. ODE models tend to confirm the notion that as time goes by, tumor size increases and stabilizes near the positive steady state. However, cell growth is limited by nutrient and its proliferation may incur a significant time delay. Motivated by the recent work of [3], we formulate a nutrient limited model of tumor growth with a discrete cell proliferation time delay and study its dynamics. We seek to answer the natural questions including if and how the cell proliferation delay may cause tumor size to oscillate. Our numerical results as demonstrated by the bifurcation diagram Figure  1 indicate under certain conditions, cell proliferation delay may cause tumor size to oscillate. Moreover, the detailed analysis of a more tractable special case also confirms that cell proliferation delay induced destabilization can happen. However, a complete and general stability investigation of the positive steady state is nontrivial. Due to the many possible scenarios, such analysis is difficult to study without specific expressions of the general functions of f (x) and g(x). We leave this as an open mathematical question.
The need of a better understanding of the dynamics of avascular tumor growth models stems from their potential applications in monitoring tumor growth with or without treatments. Many efforts are made to study the suitability of the Gompertz curve and other simple models in fitting lab data sets [10,15,20]. The Gyllenberg-Webb two dimensional nonlinear ordinary differential equation model (1.1) is regarded as an excellent choice for applications due to its balanced simplicity, reality and richness in growth dynamics ( [13,14,19]). We hope to test model (2.4) with clinical data in the near future.