NUTRIENT LIMITATIONS AS AN EXPLANATION OF GOMPERTZIAN TUMOR GROWTH

An intuitive and influential two-compartment model of cancer cell growth proposed by Gyllenberg and Webb in 1989 [18] with transition rates between proliferating and quiescent cells reproduces some important features of the well known Gompertzian growth model. However, other plausible mechanisms may also be capable of producing similar dynamics. In this paper, we formulate a resource limited three-compartment model of avascular spherical solid tumor growth and study its dynamics. The resource, such as oxygen, is assumed to enter the tumor proportional to its surface area and the dead cells form the necrotic core inside the tumor. We show the tumor growth of our model mimics that of Gompertzian model, and the solutions of our model are naturally bounded. We also identify general and explicit expressions of the tumor final sizes and study the stability of the tumor at steady states. In contrast to the Gyllenberg-Webb model, our model confirms that tumor size at the positive steady state is strictly decreasing function of the dead cell removal rate. We also present two intriguing mathematical open questions.

1. Introduction.In [15], Gompertz found that the age distribution l(t) of many communities is given by the now called Gompertz curve l(t) = ke −e b−at , where t denotes age, k > 0, a < 0 and b are constants.This function is the solution of the following simple exponential growth model with exponentially decaying growth rate.
dN (t) dt = r(t)N (t), N (0) > 0, dr(t) dt = −ar(t), r(0) = r > 0. (1) The above system is often referred as the Gompertz growth model.Notice that the solution of the above system takes the form of which resembles a Gompertz curve with b = r a and k = N (0)e r a except here the constant a must be positive.For this reason, the Gompertz curve is also frequently referred as Gompertz growth function.The tumor size is an increasing function that tends to the carrying capacity K = e b/a .Using K, we can show that N (t) = Ke AB , where A = ln(N 0 /K) and B = e −at .

EBRAHEEM O. ALZAHRANI AND YANG KUANG
Many researchers reported that Gompertz curve provided surprisingly good fit to their experimental data on growth of various animals, tumors and organisms.It is observed that the growth of solid tumors seems to follow the Gompertz growth curve very precisely in the early stages (Laird [23], Norton et al., [26]).The key assumption embodied in the Gompertz curve is that due to resource limitation, the cell growth rate decreases exponentially as a function of time.
The Gompertz curve is an empirical finding and hence one wonders why tumors grow according to the Gompertz curve in so many occasions.For this purpose, Burton [6] considered a diffusion based model for spherical tumors which develop necrotic centers.He noticed that the tumor growth limiting resource oxygen must enter the tumor by diffusion.Burton's [6] diffusion-limitation theory also predicts Gompertzian-like growth with the exception that ultimately the tumor grows linearly.Further along this line, Greenspan [17] considered a model that influenced many aspects of early mathematical oncology.Thalhauser et al. [29] reproduced the Gompertzian curve with the intuitive "go or grow" mechanism and explicit nutrient dynamics in a more complex spatiotemporal setting.The model of Frenzen and Murray [13] provides a mechanism explaining Gompertzian curve which seems to be appropriate for tumors for which it has been found that the duration of all the phases of the cell cycle increased with the size of the tumor.However, this is not necessarily the case for many solid tumors since the mean duration of the cell cycle in solid tumors is often relatively constant and that for a given solid tumor, variation in growth rate is probably due to variations in the ratio between proliferating and quiescent cells (Martinez and Griego [24]).
In this paper, we formulate a resource ratio based three-compartment model of avascular spherical solid tumor growth that explicitly models the proliferating, quiescent and dead cell populations that are usually found in an avascular tumor in its whole life span (see Fig. 1).The resource, such as oxygen or glucose, is assumed to enter the tumor proportional to its surface area and the dead cells form the necrotic core inside the tumor with the possibility that dead cells are gradually exuded from the tumor.One of our key questions is if our simple model based on nutrient limitation in the form of resource ratio can generate a family of solutions with Gompertz-like sigmoid (S-shaped) form.We show the tumor growth of our model indeed mimics that of Gompertzian curve.We also identify both mathematical and biological conditions under which the tumor has a finite final size and the proliferating layer width tends to a constant.Moreover, we systematically study the model's nonlinear dynamics and its biological implications.

2.
A resource limitation model for tumor growth.In an effort to better understand the mechanisms that maybe responsible for the Gompertz curve's success in fitting clinical tumor data, Gyllenberg and Webb [18] proposed a twocompartment model of the tumor cells consisting of only proliferating and quiescent cells.They assumed that the dead cells are removed from the tumor instantly.Their key assumption is that in a typical avascular multicellular tumor spheroid, proliferating cells may transit into and out of quiescence.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.The Gyllenberg-Webb model which will be referred as GW model below takes the following form: A center cross section of a typical avascular multicellular tumor spheroid.Q=quiescent layer, P=proliferating layer, D=necrotic core.
with initial conditions, proliferating cells are assumed to proliferate at a constant per capita rate β > 0, and proliferating and quiescent cells are assumed to die at constant rates µ p ≥ 0 and µ q ≥ 0, respectively.While these two assumptions are rather simplistic as both rates are dependent on the limiting resource levels, they are useful to gain an initial understanding of the key assumption which is that the tumor cells transition to and from the quiescent compartment at rates r 0 (N ) and r i (N ), respectively, where both functions are continuous and defined for N ≥ 0. Since stress often increase with tumor size and stressed cells tend to enter and stay at a quiescent state, Gyllenberg and Webb made the following qualitative assumptions regarding the properties of the transition functions r 0 (N ) and r i (N ).
(A1): r 0 (N ) ≥ 0 and r 0 (N ) ≥ 0 for all N > 0 and, (A2): r i (N ) ≥ 0, but r i (N ) ≤ 0, and They found that their model assumptions imply diminishing of the growth fraction (i.e.proportion of proliferating cells), a phenomenon found in most tumors.They performed extensive biologically well motivated mathematical analysis of this model and presented a host of insightful biological observations based on their mathematical findings for various interesting special cases.Naturally, all these model-based findings are likely model specific and may not hold for other plausible models.We would like to formulate an alternative model that include possibly more realistic features related to resource limitation effects on various key rates in our model.We will compare and contrast our model dynamics with that of the GW-model.More specifically, we will perform both preliminary analysis and global qualitative analysis of our model below.The qualitative analysis will be guided and complemented by extensive numerical and bifurcation results.
In our tumor model below, we let N (t) = P (t) + Q(t) + D(t), where P and Q are as in the GW-model and D represents the dead cells inside the tumor.The resource, such as oxygen, is assumed to enter the tumor proportional to its surface area that assumed to take the form of S(t) = aN (t) θ where θ is a scaling parameter with values in [2/3, 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 its effect on cell state switching between proliferating state and quiescent state may shape the growth of an avascular tumor growth.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 [14]) on R(t) which yields 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. [29]), then θ = 1.
Our model below makes an additional natural modification to the Gyllenberg-Webb model (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))), and quiescent cells transition to proliferating cells at the rate of h(u(R(t))).For simplicity, in the rest of this paper, we assume a = 1 and rename u(R(t)) as r(t).In other words, We will also refer r(t) as resource ration for the proliferating cells.Cells transit from quiescent state to proliferating state when they restart the proliferation process.
For simplicity, we assume that the proliferation rate of such newly reemerged cells from quiescent state is proportional to that of the proliferation rate of the current proliferating cells, Since proliferating cells are most likely enter quiescent state before death, we shall simply assume proliferating cells do not die (i.e., µ p = 0 in model ( 3)).All parameters below are nonnegative constants.This may result in the following threecompartment model.
where initial conditions are, Throughout the rest of this paper, we assume that (F): f (x) is continuously differentiable.f (0) = 0, f (x) > 0 for all x > 0 and, We present here some examples of functions f (x) and g(x) We would like to point out that these are ad hoc functions.Mechanistically derived forms of these functions are highly desirable for future more realistic tumor growth models.

Boundedness of tumor.
A basic requirement for any plausible population model is that the population values shall be nonnegative and bounded under appropriate conditions.We shall show below that this is readily satisfied by model (10).
From the quiescent cell growth equation, we see that if the tumor initially has only proliferating cells, then in an infinitesimal time, some proliferating cells will enter the quiescent state, and some proliferating and quiescent cells will die in the same period.In other words, for very small positive value t, P (t), Q(t) and D(t) are all positive if P 0 > 0. Hence, without loss of generality, we can assume initial population values P 0 , Q 0 and D 0 are all positive.The following proposition confirms that with positive initial values, the solutions of model (10) stay positive when exist.
Proposition 1. Solutions of model (10) with positive initial values will stay positive.
Proof.If the proposition is false, then at some time the solution turn negative, then there is a first time t 1 > 0 such that one of the three solution components becomes zero.Without loss of generality, assume first that P (t 1 ) = 0 and min{P (t), Q(t), D(t)} > 0 for t ∈ (0, t 1 ).We have which implies that a contradiction to the assumption that P (t 1 ) = 0. Similar contradictions can be derived if Q(t 1 ) = 0 or D(t 1 ) = 0 for some time t 1 > 0. This completes the proof of the proposition.
In the Gyllenberg-Webb model (3), if l 0 < β, then dP dt ≥ (β − l 0 )P, which implies that P (t) ≥ P (0)e (β−l0)t .Hence, mathematical conditions must be placed on l 0 to ensure that the tumor size is bounded.However, the sizes of most avascular tumors are limited due to the resource limitation.Naturally, we would like to see no or minimum requirements to ensure the boundedness of the population densities in model (10).
(15) Therefore, Similarly, we can see that D(t) ≤ D(0) + k D P m (t).The rest of this lemma is obvious.
The above lemma implies that if P (t) is bounded, so are Q(t), D(t) and hence N (t).If P (t) is unbounded, then obviously N (t) > P (t) is also unbounded.We are now ready to state and prove the following boundedness result for the solutions of model (10).Theorem 3.2.Solutions of model (10) are bounded provided that θ ∈ (0, 1).In which case, there is a constant L > 0 such that if N (0) < L, then N (t) ≤ L.
Proof.Since f (x) is continuously differentiable and f (0) = 0.For any ε > 0, there is a constant By the previous lemma and the assumption that θ ∈ (0, 1), we see that there is a P ε > 0 such that if N (0) < P ε and P (t) ≥ P ε , then We claim that if N (0) < 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 (17), we see that r(t 1 ) < δ.We have a contradiction, proving the claim.Let L = (2 + k Q + k D )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.We shall also consider the existence and locations of nonnegative steady states of the model (10) which may help us to better understand how the model parameters and functions change the steady state tumor size and structure.
We consider first the existence and locations of nonnegative steady states of model (10).Clearly (0, 0, 0) is the trivial steady state, which affirms the reality that if no proliferating cells present at the beginning, then the tumor size stay at zero.However, the model system does not admit a Jacobian at (0, 0, 0) which will prevent one to study the local stability of model ( 10) at (0, 0, 0) via its linearized system at origin.Let E * = (P * , Q * , D * ) be a nontrivial steady state of model (10).It is easy to see that if such a steady state E * exists, all components of E * must be positive.
It is easy to see that dN dt = f (r)P − dD.
At steady state, we must have f (r * )P * − dD * = 0. From the D equation, we have µQ * − dD * = 0. Hence, f (r * )P * = dD * = µQ * and From the Q equation, we have which implies that It is easy to see that F (x) is a strictly decreasing function in [0, ∞), F (0) = g(0) > 0, and lim x→∞ F (x) = −f M (αf M + µ) < 0. Hence, F (x) = 0 has a unique solution r * .Observe that r * does not depend on the dead cell removal rate d.It is easy to see that lim Observe that Equation ( 23) is equivalent to which implies that Therefore, It is easy to see that Clearly, g(r * ) − f (r * ) > 0. An implicit differentiation of Equation ( 21) with respect to µ at the value r * yields Equation ( 21) and r * = (N * ) θ /P * imply there is a unique positive value P * which is given by Therefore, E * = (P * , Q * , D * ) = (P * , f (r * ) µ P * , f (r * ) d P * ) exists and is unique.Observe that all components of E * and hence the tumor size N * at the positive steady state are strictly decreasing function of the dead cell removal rate d.This removes an unrealistic result from the Gyllenberg-Webb model (3) which asserts that tumor size at the positive steady state does not depend on the dead cell removal rate d (see [2]).
Equation ( 21) also implies that at the positive steady state, the percentage of proliferating cells in the tumor is The elementary mathematical findings of Equations ( 28), ( 29) and ( 31) can be summarized into the following proposition on the steady state tumor size and structure.
Proposition 2. If an avascular tumor grows following the model (10), then the percentage of proliferating cells in a tumor at the positive steady state increases when the rate α of quiescent cells returns to proliferating state or the dead cell removal rate d increases.In addition, as the dead cell removal rate d increases to infinity, the percentage of proliferating cells in a tumor at the positive steady state approaches a constant; while as the dead cell removal rate d or the death rate of quiescent cells µ decreases to zero, the percentage of proliferating cells in a tumor at the positive steady state approaches to zero.
The last statement of the above proposition agrees to the Gyllenberg-Webb model ( 3), for which one can readily show that as the death rate of quiescent cells decreases to zero, the percentage of proliferating cells approaches to zero.However, the Gyllenberg-Webb model (3) also implies a slow dying quiescent population will increase the size of the tumor at the steady state level ( [2]).Bifurcation diagram (see Figure 4 below) clearly shows that this is also true for our model (10).However, this is yet to be confirmed mathematically.

5.
Stability of a tumor at steady states.We now study the stability of the trivial steady state and the positive steady state.
We consider first the local stability of E 0 = (0, 0, 0).Unfortunately, this steady state does not admit Jacobian since the model is not differentiable at the origin due to the resource ratio expression r(t) = N θ /P.Fortunately, Lemma 3.1 comes to our rescue.From Lemma 3.1, we see that if Q(0) = D(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 origin, r(t) is very large and hence f (r(t)) > f M /2 > 2g(r(t)).Therefore, near the origin, P (t) > f M P (t)/4, (33) implying P (t) grows at least exponentially near the origin and hence (0, 0, 0) is unstable.
We now consider the local stability of E * = (P * , Q * , D * ).The mathematical question of establishing a good set of sufficient conditions for its global stability appears to be formidable and remains open.In fact, even obtaining a set of biologically meaningful conditions for its local stability is a highly involved exercise.
For convenience, we will use the following notations.
Linearizing the model at this steady state yields the following Jacobian, At steady state, we have (37) From Routh-Hurwitz criterion, we know that all eigenvalues of J have negative real parts if all the following three conditions hold: (i) trA < 0, (ii) detA < 0, (iii) ∆ ≡ detA − (trA) where A kk is the determinant of the 2 × 2 matrix obtained by removing the k-th row and k-th column from A (see page 234 in [8]).
Computing the determinant of A by expanding the last row yields Routine computation also shows that (40) A 33 is more complicated and we separate its expression into two parts. where and (45) Therefore, we have established the following local stability result for the positive steady state E * .
The above theorem is rather difficult to use.In the following we consider a probably very plausible scenario: quiescent cells reenter the proliferating states only when there is abundant nutrient supply that typical happens when the tumor has its own blood supply.Since we are considering only the growth dynamics of an avascular tumor, we can simply assume that α = 0.With this reasonable assumption, we shall show below that the positive steady state is always locally asymptotically stable.To this end, we need the following lemma.Lemma 5.2.Assume α = 0, then the determinant of A in Equation (35) is negative.
Proof.Observe that when α = 0, at the steady state, we have A direct computation of the determinant A 33 yields, Since µQ * = g * P * and 0 < θ < 1, we see that It is easy to see that det Simple algebraic derivation yields Replacing g * P * by dD * in the above equation yields This proves the lemma.
We are now ready to state and prove the following theorem.
Proof.By now we know that both the trace and the determinant of the Jacobian matrix A are negative.To establish the theorem, we need only to show that ∆ ≡ detA − (trA) It is easy to see that A kk > dA 33 .

6.
Discussion.A solid tumor grows in two stages.It begins with a slow growing avascular growth period due to the nutrient limitation, followed by an accelerated vascular growth phase that signals an advanced stage for the tumor.In the avascular growth period, the tumor's nutrient supply is provided by the slow diffusion process and hence the tumor's size is often confined into a region with a few millimeters in diameter.Its growth in this phase often is well captured by a Gompertz growth curve.
Over the last several decades, mathematical oncologists were intrigued by the suitability of the Gompertz curve in fitting many lab data sets and made efforts to mechanistically model this phase of tumor growth with the hope to understand why a simple curve like Gompertz curve is often better than some of the more sophisticated models (Frenzen and Murray [13], Kozusko and Bajzer [20], Marusic and Vuk-Pavlovic [25], Norton [26]).The Gyllenberg-Webb two dimensional nonlinear ordinary differential equation model (3) stands out due to its balanced simplicity, reality and richness in growth dynamics (Gyllenberg and Webb [18] and [19]).
However, a closer look with the mathematical analysis lens, as presented by the authors in Gyllenberg and Webb [18] indicates there are still some missing biological links.They assumed that dead cells are removed from tumor instantly by undisclosed process and the model implies that avascular tumor final size does not depend on the dead cell removal rate even if one includes it in a more realistic three dimensional version of their model (Alzahrani et al. [2]).Moreover, the Gyllenberg-Webb model (3) model requires some artificial conditions to ensure the solutions are bounded.Most importantly, while it implicitly cites the nutrient limitation inspired the model's selection of the state transition functions (r 0 (N ) and r i (N )), the resource availability was not explicitly accounted for.This was addressed in a recent paper by Wallace and Guo [30] where the author presented interesting numerical findings of three specific and preliminary models (without model formulation details) with explicit nutrient limitation due to surface diffusion.A further examination of their work may yield biologically more realistic and mathematically tractable models.Observe that all the components as well as the tumor size grow like a sigmoid curve and tending to constants, which is the hallmark of Gompertz tumor growth.
In this paper, we modified the Gyllenberg-Webb model (3) by explicitly expressing the limiting resource as a function (R(t) = N θ (t)) proportional to the tumor surface area and the proliferating cells grow according to the per cell resource availability (r(t) = N θ (t)/P (t)).This explicit nutrient dynamics incorporation ensures solutions are naturally bounded.However, the model remains largely a phenomenological one.It is highly desirable to identify specific biological mechanisms that may lead to the formulation of state switching functions f and g.Despite its simplicity, our model does present some interesting and tractable mathematical challenges largely due to the nonlinearity embedded in the r(t) function.For example, our numerical results suggest that the positive steady state is globally stable (Figs.2,3 and 4), however, we still do not even have a complete result on its local stability and have had trouble to obtain any realistic global stability result for the positive steady state.We are also unable to show mathematically that, like in the Gyllenberg-Webb model (3), a slow dying quiescent population will increase the size of the tumor at the steady state level.We thus present the following two open question.
Open question 2: Assume that a tumor growth can be modeled by (10), then a slow dying quiescent population (with smaller value of µ) will increase the size of the tumor at the steady state level (see Fig. 2).Akin to the common sense that energy often serve as the driving force for many physical processes, one can think a limiting resource drives cells growth, movement, transition, or evolution.It is thus important for tumor models to implicitly or explicitly include such limiting resources.Some recent tumor modeling efforts along this line are presented in Eikenberry et al. [10] for a glioblastoma model, in Eikenberry et al. [11] for a melanoma model, Eikenberry et al. [9] for a evolutionary prostate cancer model, and Thalhauser et al. [29] for a tumor cord model.More recently, some resource-explicit tumor models are also preliminarily validated by clinical data (Portz et al. [27]; Everett et al. [12]).For more realistic and mathematically challenging avascular tumor models, one can incorporate additional biological processes and many such models are discussed in the classical book of Adam and Bellomo [1] and the comprehensive review article of Araujo and McElwain [3].A natural extension of our model is to include time delays to account for the time delays in cell transition and cell death processes and the resulting dynamics can be studied using recent developed mathematical theories (Kuang [21], Kuang [22], Beretta and Kuang [4], Gourley and Kuang [16]).In reality, tumor growth is often limited by space; and in clinical applications, tumor spatial distribution must be modeled.In these situations, one shall consider formulating partial differential equation tumor growth models (Sherratt and Chaplain [28], Byrne et al. [7]).
dA 33 + µr D Ag * − (trA) 3 k=1 A kk .(55)It is easy to see that all the four terms in trA = Ar P + Br Q − µ − d are negative and A ii are all positive for i = 1, 2, 3. Hence, −

Figure 2 .
Figure 2. Two sets of solutions of model (10) with f (r) = kr ar+1 and g(r) = c r+m .Except for the dead cell removal rate d, the parameters for both panels are the same.They are α = 0.1, k = 2, a = 1, m = 2, µ = 0.5, c = 1, θ = 2/3.For panel (a), d = 1 and for panel (b), d = 0.1.The initial conditions are P (0) = 0.1, Q(0) = D(0) = 0 for both panels.Observe that all the components as well as the tumor size grow like a sigmoid curve and tending to constants, which is the hallmark of Gompertz tumor growth.

Figure 3 .
Figure 3. Bifurcation diagrams of model (10) with f (r) = kr ar+1 and g(r) = c r+m using the dead cell removal rate d as the bifurcation parameter.The parameter values are α = 1, k = 2, a = 1, m = 2, µ = 0.5, c = 1, θ = 2/3.The positive steady state appears to be globally attractive and the components of the steady state are decreasing with respect to d.

Figure 4 .
Figure 4. Bifurcation diagrams of model (10) with f (r) = kr ar+1 and g(r) = c r+m using the quiescent cell death rate µ as the bifurcation parameter.The parameter values are the same as in the previous figure except d = 1 and µ varies.The positive steady state appears to be globally attractive and the components of the steady state are decreasing with respect to µ.