HOW DO CELL CROWDING AND STARVATION AFFECT AVASCULAR TUMOR GROWTH OF THE EMT6/RO TUMOR?

. Here we re-examine experimental in vitro data for the EMT6/Ro tumor volume and viable rim thickness. This shows that the growth speed is constant in time, and independent of nutrient concentration at large concentrations, but that the viable rim thickness increases in this high concentration regime. We then present a simple mechanistic reaction-diﬀusion equation that includes crowding and starvation eﬀects, and show that the model qualitatively captures the experimental observations. More-over, the model predicts that the cancer cell concentration is characterised by a wave pulse (soliton), and the pulse shape is explored through zero’th order perturbation analysis corresponding to large wave speeds. It is shown that this zero’th order term is dominant for the experimental condition, further indicating that the non-linear reaction governs the pulse characteristic shape. Finally, at low nutrient concentrations we ﬁnd that the front is a pulled-front, that is, the growth speed is determined by the cell multiplication at the wave front edge in accordance with earlier modelling approaches.


Introduction
Tumor growth can be separated into three growth stages: the avascular, vascular and metastatic stages [23]. In the first avascular stage the tumor is divided into three cell layers; the outermost layer consists of proliferating cells and the tumor growth is envisioned to be due to cell multiplication in this layer by consumption of nutrients [3]. Behind this proliferating cell layer there is a layer of living quiescent cells that do not multiply, but are in a hibernating state. At the tumor centre there are dead cells forming the so-called necrotic core [3]. This three-layer picture is commonly referred to as the multi-cellular spheroid (MCS) model. It is important to note that proliferating and quiescent cancer cells are not clearly distinguishable, and these two layers should be considered as a single layer of living cells wherein the fraction of proliferating cells is varying [15,23]; we will refer to this single layer as the viable rim.
The common understanding of avascular tumor growth is that the tumor volume initially increases exponentially with time, and after this initial phase the growth rate decreases, and the volume eventually reaches a constant plateau [2,3,9,12,15,19,22]. The first exponential growth is a consequence of the fact that all proliferating cancer cells are exposed to approximately constant nutrient resources, however, as the tumor  [15] and punctured line a super-imposed line with slope of 3. The tumor is EMT6/Ro mouse mammary tumor cell grown in vitro [8,9]. (b) The corresponding viable rim thickness versus glucose concentration. The glucose concentration varies by a factor of 32. Data are taken from Mueller-Klieser et al. [17]. The tumor specifications are the same as in (a).
volume increases the proliferating cells in the viable rim consume the nutrients and starve the cells closer to the tumor center, hence, the growth rate decreases. For long times the growth eventually stops and the entire growth dynamics may be described empirically, for example, by the Gompertz function [3].
In Figure 1a we re-plot data from Jiang et al. [15] showing the avascular solid tumor volume, V = V (t), for the EMT6/Ro tumor. The tumor is grown in vitro in a homogeneous gel and do not attach to any surface before growth, hence, it can grow in all three directions. Notice that both axes are in logarithmic scales and that the superimposed line has slope of 3. From this, one observes that V ∝ t 3 in the time span of at least 1 to 22 days. Equivalently, as the tumor is to a good approximation spherical [3,15] the tumor radius, R, is proportional with time. If the tumor radius is R 0 as the growth reaches the power law regime, the tumor radius can then be written as R = R(t) = st + R 0 , where s is a characteristic growth speed. This interpretation of the data results in the growth speed (or rate) of s = 19.6 ± 0.6 µm per day; the error estimator is the standard error assuming the slope is normal distributed.
The speed is, within statistical uncertainty, independent of the oxygen and glucose concentrations for the concentrations investigated by Jiang et al. A natural first approach is to assume that the growth rate is proportional to the tumor surface area, i.e., proportional to V 2/3 which agrees with behavior V ∝ t 3 [10,16]. However, this coarse-grained picture does not provide details about, for example, the nutrient penetration depth and viable rim thickness. In fact, to our knowledge only few other authors [6,10,14,16,21,25,26] actually discuss this simple power law.
Jiang et al. [15] also showed that the viable rim thickness, L, is constant as a function of tumor diameter, that is, L is constant with respect to time. This contradicts earlier findings by Freyer and Sutherland [9] who conclude that the viable rim thickness decreases as the tumor grows, however, from their data it is questionable whether this claim is supported within statistical uncertainties. We therefore follow Jiang et al. and conjecture that the viable rim thickness remains constant after some short time at specific glucose and oxygen concentrations. In Figure 1b we plot data for the viable rim thickness of the EMT6/Ro tumor as function of the nutrient concentration; data are from [17], see also comment [20]. As shown, the rim thickness increases as the nutrient concentration increases. Comparing the glucose concentrations in Figure 1a and b one sees that the concentrations used in (a) corresponds to the high concentration regime in (b), where the viable rim thickness is increasing. In this work we seek to qualitatively model the behaviour shown in Figure 1.
From the description above, we envision the viable cancer cell concentration as a traveling wave-pulse propagating through the nutrient enriched gel at constant speed. Ahead of the pulse the system is homogeneous, and behind the pulse the necrotic core is located. Such systems can be modelled through a reaction-diffusion equation [13,18], see e.g. Jiang et al. [15] for a specific example, and this approach is also sought here. There exists an overwhelming large body of previous modelling efforts on the growth of avascular tumors, see for example [2,4,12,19,[23][24][25]. With time the models have become highly detailed involving many mechanisms and can often only be analysed numerically due to the non-linear nature of the problem and the many dynamical degrees of freedom. Also, the parameter space in these models can be very large, thus, complicating the interpretation and insight into the underlying mechanisms that are responsible for the growth phenomena. In this paper we will move in the opposite direction and propose a simple mechanistic model for the MCS growth that emphasizes the main relevant underlying mechanism(s) behind the data. Such a simplified model can also act as an integral part of a larger modelling effort for the multi-stage tumor growth development.

The model
As an initial approximation we assume that the tumor is perfectly spherical and that the concentrations are independent of the angles. The model dynamical variables are then defined by the nutrient (glucose) concentration N = N (r, t), living cancer cell concentration C = C(r, t), and dead cell concentration D = D(r, t). r is the radial distance to the sphere centre and t is the time. We will not include the oxygen dependency here. There is no bulk flow present in the system, and also neglecting effects from cell volume changes we can disregard advection, hence the general balance equation [5] reads where φ represents one of the concentrations, f φ the reaction term, and j φ the flux of φ. We first address the reaction terms, f φ . We stress again that there is no clear distinction between proliferating and quiescent cells [15], thus, we simply let both these cell states make up the pool of living cells C; this also follows Ward and King [25]. The cell concentration must be bounded, i.e., there exists a finite cell carrying capacity, C 0 . Moreover, we introduce a metabolic saturation state wherein the cell division rate is bounded by an upper limit with respect to nutrient concentration. The cell division rate, r 1 , then follows the functional form where k 1 is a second order rate constant, and K 1 > 0, hence, it acts as a Michaelis-type constant. We let the cell death rate, r 2 , be dependent on the nutrients available. At zero concentration the cells are in a starvation state, where the death rate is largest. As the nutrient concentration increases the death rate decreases and approaches an inherent minimum death rate. Again, we describe this using the functional form [25] where k 2 is the corresponding first order rate constant, and K 2 > 0 the Michaelis constant. N 2 is a positive dimensionless constant fulfilling 0 < N 2 < 1.
We can now write the reaction terms for the nutrient concentration, f N , living cell concentration, f C , and dead cell concentration, f D . If ν is the stochiometric relationship between nutrient consumption and cell production we have f C is also referred to as the mitotic rate. While the algebraic expressions for the cell division and cell death rates are of the Michaelis-Menten type, they are not related to enzyme kinetics, nor are they unique. However, they do offer intuitive mathematical expressions for the cell dynamics and, also, ease the analysis below compared to non-algebraic expressions.
Next we address the cell mobility. The cells do not attach to any surface and are thus mobile, and we will assume that only the living cells feature mobility as the dead cells are caged in the necrotic core. The living cell mobility is caused by a crowding effect meaning that the cells migrate towards lower concentration regions. In the most general case the cells move down gradients of both living and dead cell concentrations, and importantly, in the presence of a dead cell concentration gradient a flux of cells can naturally be realised if C > 0 at that point. The cell flux j C must therefore be a function of C as well as the gradients ∇C and ∇D. Applying simple linear constitutive relations we write the flux j C = j C (r, t) of C as where D C is the mobility or diffusion coefficient for the living cells, and α is a weighting function describing the migration due to gradients in dead cells. We will denote this the cross-mobility or cross-diffusion effect, where the weighting function fulfills α ≥ 0 and α = 0 for C = 0. Notice that cell crowding enters both the reaction term, equation (2.2), via the carrying capacity and the migration term. The nutrients will also diffuse in presence of gradients, and we here ignore any cross-diffusion effects. The flux, j N = j N (r, t), of N is therefore simply Fickian where D N is the corresponding diffusion coefficient. Substitution of f N , f C , f D , j C , and j N into the balance equation for N, C, and D in three dimensional spherical geometry gives assuming that the diffusion coefficients are constants. Introducing the variable transformation we get the dimensionless and final form It is worth highlighting that for D > 1 the systems has high nutrient diffusivity compared to cell migration, and large values for β indicates a faster cell multiplication rate compared to the death rate. The tumor growth can be initialised by a non-zero initial concentration of cancer cells around x = 0, c(x, 0) = c 0 (x) and a uniform nutrient distribution n(x, 0) = n ∞ , where n ∞ is the constant concentration ahead of the front. We assume that the number of dead cells is initially zero, d(x, 0) = 0, and we let the boundary conditions for n and c be given by Equations (2.9) together with initial and boundary conditions, equation (2.10), form the basic tumor growth model.

Zero order perturbation analysis
In this section we focus on the regime where the cross diffusion effect can be ignored, i.e., we set α = 0 such that the dynamics decouples from the dead cell concentration d. We will comment on this choice later in Section 3. After a transient time the tumor radius is sufficiently large while the gradient in concentrations are bounded such that in the wave we can ignore the differential operator (1/x)∂/∂x that enters in the diffusion terms. The validity of the simplification is addressed numerically in Section 3, where the different biophysical effects are studied. In this regime the model, equation (2.9), reads Assume that a wave exists and that this wave has, after the transient, a constant propagation speed s. We can now introduce the moving frame coordinate z = x − sτ , and we have from the chain rule such that c = c(z) and n = n(z). Thus, following the wave in the moving frame, the system is described through two coupled ordinary differential equations. Notice that z → −∞ is the moving frame coordinate behind the wave and z → ∞ is the moving frame coordinate ahead of the front, hence, we have lim z→∞ n(z) = n ∞ and lim z→∞ c(z) = 0. It is possible to introduce a perturbation parameter, ε = 1/s 2 . By a second change of variable, see e.g. [18], and we arrive at a singular perturbation problem. Notice that ε → 0 as s → ∞, that is, the perturbation analysis valid at small ε-parameter values correspond to large wave speeds. From equation (2.13) we see that in this large speed limit, diffusion has small effect on the wave characteristics which is controlled by the reaction terms. However, diffusion enters the problem implicitly as equations (2.13) are formed under the fundamental assumption that a wave exists. Now, let n and c be given by the power series where n k and c k are functions of ξ. Substituting into Equations (2.13), and ordering in powers of ε we have for zeroth order Note, only the reaction terms remain on the right-hand side. For ν = 0 equations (2.15) allow for fixed points fulfilling (1 − N 2 n 0 /(n 0 + κ 2 ))c 0 = 0. Since we restrict the investigation to positive nutrient concentrations we have the constraint n 0 ≥ 0, hence, c 0 has fixed point at 0. Due to the boundary conditions n 0 → n ∞ as ξ → ∞. The solution to equation (2.15) is not known to the authors, and it is therefore solved numerically by using backwards integration with n 0 = n ∞ and c 0 = 10 −4 as initial condition. Notice that equations (2.15) is an initial value problem which is now independent of the boundary values at r = 0, or equivalently x = 0.
For the numerical treatment one needs to consider the model's characteristic concentration scales. Assuming a cell volume of 10 2 −10 3 µm 3 one has, with ideal packing, on order of 10 13 −10 14 cells per liter. This corresponds to a molarity, C 0 , on order of 10 −9 −10 −10 M. The molarity of the nutrient is on order of N = 10 −2 M, meaning that the dimensionless nutrient concentration, n ∞ = N/C 0 , ahead of the wave is on order of 10 7 −10 8 . We will ease ideal packing and use the concentration 10 5 ≤ n ∞ ≤ 3 × 10 6 unless otherwise stated. This particular choice does not affect the overall results, for example, letting n ∞ = 10 7 and varying the model parameters appropriately gives approximately the same solution curves. Figure 2a shows the (n 0 , c 0 )-phase plane and four solutions. Figure 2b shows the corresponding profiles; notice that the origin of the ξ-coordinate is arbitrary. One clearly sees the existence of a wave-pulse, i.e., a maximum value for c 0 . The persistent pulse shape is thus due to the non-linear reaction terms, that is, the phenomena is a soliton.
In the regime κ 2 n κ 1 and 0 < c 1 it is possible to find an analytical expression for the maximum in terms of the model parameters. First, note in this regime we have and where ∆N = 1 − N 2 , that is, we are well below the system cell carrying capacity and starvation point. The zeroth order perturbation term is then Notice, that this perturbation equation does not bound the cell concentration to c ≤ 1, that is, c can exceed the carrying capacity. The solution for c 0 at relatively low nutrient concentration is shown in Figure 2b as punctured line. Compared to equation (2.15) the tail (decreasing ξ) decays slowly and the pulse thickness increases. This is due to the simple linear decay for the death rate in equation (2.18b). For n 0 > 0 we have from Equations (2.18) This equation is separable, and integrating with lower bounds c 0 = 0 and n 0 = n ∞ , that is, applying the values ahead of the front we get This function has a maximum at n 0 = ∆N κ 1 /β with value Figure 2b by the filled circle for n ∞ = 3 × 10 5 . This expression for the maximum concentration is convex always fulfilling max(c 0 ) ≥ 0 with max(c 0 ) = 0 whenever n ∞ β = κ 1 ∆N . For n ∞ β > κ 1 ∆N the first term on the right-hand side dominates and max(c 0 ) ∝ n ∞ . This result is, of course, only valid for low cell and nutrient concentrations.

Wave pulse speed
We now leave the perturbation analysis, i.e., we are not limited to large wave-pulse speeds and continue the investigation of the non-perturbed system. Specifically, we study the regime α = 0, κ 2 n κ 1 , and 0 < c 1. Furthermore, assume that in the wave front-edge the change in nutrient concentration is dominated by the reaction, such that for C cells produced the nutrient amount νC is consumed, i.e. n ≈ n ∞ − νc in the front edge. Substitution into equation (2.9b) leads to a single dynamical equation for the cell dynamics in the front-edge after the initial transient. In the moving frame where c = c(z) this becomes This equation is on the form of a standard auto-catalysis process with decay [11,13], and can be studied using straight forward linear stability analysis. To this end we identify two fixed points, namely, c min = 0 and c max = βn ∞ − ∆N κ 1 βν .

(2.24)
This implies an imposed constraint βn ∞ ≥ ∆N κ 1 since c ≥ 0. This is the region of increasing maximum cell concentration with respect to nutrient as shown above. The former fixed point c min has eigenvalues and, again, as c ≥ 0 this fixed point cannot be a spiral point which is fulfilled if From this we define the minimum allowed speed by This linear analysis is equivalent to the assumption that at the front edge c decays proportional with e −λz , λ being the eigenvalue at c min . The concentration constraint is illustrated in Figure 3 showing the solution curves for equation (2.23). Also note that for βn ∞ > ∆N κ 1 the fixed point is a stable node. The latter fixed point c max has eigenvalues which is then a saddle point if βn ∞ > ∆N κ 1 . Thus, we expect that for α = 0, κ 2 n κ 1 , a wave front exists with a lower bounded speed given approximately by equation (2.27) if (2.29) For this relation we also have from equation (2.21) that max(c 0 ) < n ∞ /ν. Importantly, this analysis is based on an unbounded reaction term which is only valid in the limit of sufficiently small concentrations. This in turn leads to an unbounded wave speed with respect to nutrient concentration. Even if s min is not bounded with respect to n ∞ , we conjecture that the wave-pulse speed resembles a pulled front mechanism, such that the speed increases with nutrient concentration, at least, in the sufficient low nutrient concentration regime.

Numerical solution
The model equation system, equation (2.9), poses a relatively simple numerical problem. The system is dissipative and an explicit first order time integrator with a second order spatial finite difference scheme (FTCS) [7] suffice for the numerical integration. In order to control the integration error an adaptive time step algorithm is implemented. The error estimator is given by the maximum of the spatial difference in the cell concentration, i.e., err(τ ) = max(|c(x, τ ) ∆τ − c(x, τ ) ∆τ /2 |), where c(x, τ ) ∆τ is the cell concentration at position x and time τ using a time step ∆τ and c(x, τ ) ∆τ /2 is the corresponding value using half the time step. The relative acceptance tolerance is set to 0.1 per cent. The spatial resolution of the dimensionless spatial coordinate is set to ∆x = 0.25. The system length L sys is around 20 to 40 times the viable rim thickness. Finally, in accordance with the boundary conditions, equations (2.10), we apply homogeneous Neumann boundary conditions at x = 0. For x = L sys we use n(L sys , τ ) = n ∞ and c(L sys , τ ) = d(L sys , τ ) = 0. The validity of the scheme settings was tested by increasing the system size L sys , decreasing the error estimator, and changing the spatial resolution showing no measurable effects.
In the first simulation series we investigate the full model, equation (2.9), by assuming where k is sufficiently large; here we use k = 50. This corresponds to a smoothed step function modelling no cross-migration at zero cell concentration to cross-migration α = α 0 for non-zero concentrations. In particular, α goes from zero to α 0 over the cell concentration range 0 to 0.1, resulting in a full coupling effect even in the front's edge. The system is initialized using three different initial conditions: (i) A Gaussian c(x, 0) = e −x 2 /10 , (ii) a step function c(x, 0) = 1 for 0 ≤ x ≤ 2 otherwise c(x, 0) = 0, and (iii) a linear function c(x, 0) = −x/2 + 1 for 0 ≤ x ≤ 2 otherwise c(x, 0) = 0. The different initial conditions may not all be experimentally relevant, but are used here to investigate whether the wave features are sensitive to the initial conditions. We let n(x, 0) = n ∞ and d(x, 0) = 0. Figure 4(a) shows the living cell concentration profile for five different times in the case of initial condition (i); the first profile, corresponding to τ = 0, is scaled with a factor 1/2 for clarity. A propagating wave-pulse is clearly seen to develop in agreement with the zero order perturbation analysis above, Figure 2. The tumor radius is defined to be the point, x = R, where the wave-pulse has its maximum; this is illustrated for time τ = 25 in Figure 4a. Figure 4b plots the corresponding radius R as a function of time. After approximately τ = 50 the speed is within 5 per cent of the final speed (within numerical accuracy). Using the experimental data, Figure 1a, where the speed is constant within statistical uncertainty after approximately 24 hours, we can from equation (2.8) define the characteristic time scale 1/k 2 = 1/2 h. From the time scale we can estimate the characteristic length scale D C /k 2 in term of the dimensionless diffusion coefficient D. Jiang et al. [15] estimate D N to be on order of 10 −7 cm 2 h −1 , hence, we obtain D C /k 2 = D N /(Dk 2 ) ≈ 2µm/ √ D. Importantly, the final wave-pulse profile and speed is observed to be independent of the three different initial conditions, indicating that the model features a stable wave-pulse solution (solition).
Recall that in the analysis in Section 2.1 we ignore the curvature, i.e., the results are only valid for times where the operator (1/x)∂/∂x is sufficiently small. At τ = 25, approximately 50 hours, the terms (1/x)∂c/∂x and ∂ 2 c/∂x 2 were evaluated numerically; the former contribution (i.e. the tumor curvature) to the total diffusion is less than 4 per cent. Thus, the effect of curvature is neglectable after a relative short transient time compared to the experimental time span, Figure 1a, and we will simply ignore this term in the treatment from here on. We then proceed and study the wave-pulse in the moving frame coordinate z. This is done by monitoring the pulse maximum and shift it such that the wave maximum is fixed at z = 2/3L sys . The total shift length is stored and from this the speed is calculated. Figure 5 shows the wave-pulse for α 0 = 0, 1, and 2 in the moving frame after convergence to constant profile and speed (within numerical accuracy). It is seen that, for these parameters, the pulse profile does not change dramatically when the cross-migration term is included. Especially, the wave front-edge is very similar. Also, the front speeds for α 0 = 0 and α 0 = 2 are evaluated to 1.25 and 1.56, respectively, that is, the speed does not feature a significant dependency on the exact value of α 0 . In line with Section 2.2, we conjecture that this low dependency of the cross diffusion is because the pulse speed is dominated by the reactions occurring at the edge (denoted a pulled front), i.e., away from the region of non-zero dead cell gradient. We therefore conclude that the pulse features are not significantly dependent on the crowding due to the dead cells, and we continue without the cross-migration term, reducing the problem to that of equation (2.11).
We return to the experimental data in Figure 1, and recall that while the wave propagation speed is constant with respect to nutrient concentration in the interval from 5.5 to 16.5 mM glucose, the rim thickness increases. Figure 6 shows the wave-pulse speed and rim thickness as function of nutrient concentration. The rim thickness is defined as the pulse width, where c is 10 per cent of the maximum value (see inset in Fig. 6b); the qualitative behavior is independent of the exact definition. At relative low concentrations the speed increases with increasing nutrient concentration, but convergences to constant in the large nutrient concentration regime. On the other hand, the wave-pulse width increases linearly in the entire nutrient concentration interval studied here.
The regime of constant speed and increasing rim thickness is characterized by a large cell concentration in the pulse, C ≈ C 0 . This is also seen in the inset in Figure 6b. We therefore attribute the maximum allowed cell concentration (or crowding) to be the limiting factor for the speed. Moreover, the model assumes that the cell death rate decrease with increasing nutrient concentration, i.e., for large concentrations the nutrients penetrates deeper into the tumor resulting in a cell death rate decrease in these regions. This picture is supported by the mitotic rate which is plotted in Figure 7 for different nutrient concentrations. Positive mitosis (cell multiplication) is observed to be located only at the wave front-edge, and negative mitosis (cell death) is located the pulse back-end. The negative mitosis is shifted backwards as the nutrient concentration increases, thus, increasing the rim thickness. Figure 8 compares the perturbation equations, equations (2.15) and (2.18), with the numerical solution. In (a) where n ∞ = 10 5 we are in the low cell concentration limit and we expect equation (2.18) to show fair agreement with simulation results. This is also the case at the pulse front-edge, however, the pulse tail decays slowly compared to simulation data. This is due to the linear cell death rate of equation (2.18). This linearity is not assumed in equation (2.15) and the agreement with the simulation result is much better. In (b) we have  n ∞ = 10 6 which is outside the limits of equation (2.18); this is evident from the maximum cell concentration that goes above the bound C 0 , or equivalently c = 1. On the other hand, equation (2.15) agrees perfectly with the output from the numerical simulation. The agreement is better for the high nutrient concentration regime which is consistent with larger speed s and thus smaller perturbation parameter ε. This strongly indicates that the pulse characteristics is dominated by the reaction term in the domain of large wave speeds.
Equation (2.27) predicts that the minimum allowed wave-pulse speed increases as the square root of the nutrient concentration. Recall, that strictly the prediction is (i) limited to the regime κ 2 n κ 1 and 0 < c 1, and (ii) based on the assumption that the nutrient concentration at the wave front-edge is given by the cell consumption alone. Figure 9 compares s min and s found in the numerical simulations. First, it can be seen that s min , in fact, overestimates the speed, which is likely due to assumption (ii) since the large gradients  in n is ignored. As expected, however, the difference goes to zero as the nutrient concentration goes to zero. Interestingly, the relative difference is approximately linear indicating different functional forms for s and s min ; this is also evident in Figure 6 where s is seen to be bounded and converge to a constant value.

Conclusion
In this paper we revisited the experimental data for the growth of EMT6/Ro avascular solid tumor from Mueller-Klieser et al. [17], Freyer and Sutherland [9], and Jiang et al. [15], see Figure 1. It was shown that the tumor radius grows linearly in time, and that the viable rim thickness increases with increasing glucose concentration. From this we proposed a reaction-diffusion model, which includes the cross-migration effects between dead and living cells, a cell carrying capacity, as well as metabolic saturation and starvation.
Through numerical explorations it was shown that the model features a propagating wave-pulse with a constant speed. It was also shown that in the high nutrient concentration regime the wave speed is approximately constant, but the pulse thickness, i.e., the tumor viable rim thickness, increase. This is also in agreement with experimental data. The model features are present for zero cross-migration and also when ignoring the tumor curvature. Therefore, we conjecture that the experimental observations are due to (i) maximum cell carrying capacity which constrains the pulse speed, and (ii) a starvation phenomenon, where the cell death rate increases as the nutrient concentration decreases.
In the limit of large wave-pulse speed, it is possible to make a perturbation analysis. It was shown from this analysis that the pulse shape is primarily determined by the reaction terms in the model. This is in good agreement with the conjecture above.
Assuming that the pulse speed is determined by the cell multiplication at the very front-edge it is possible to derive a minimum allowed wave speed. The predicted minimum speed was observed in fair agreement with simulation in a limited parameter regime and for low nutrient concentrations indicating that the pulse is in this case a so-called pulled-front.
It is important to point out that recently Brüningk et al. [1] showed that the HCT116 tumor (human colon cancer) also features a linear growth regime for up to 25 days, hence, the present model is not limited to describing the mechanisms for the RMT6/Ro tumor growth, but is more general.
This work focuses on the affects from varying nutrient concentration. A natural continuation of this would be inclusion of effects from, for example, oxygen concentrations and vascular development. In line with the underlying philosophy used here, insight into important mechanisms in these extended and more complicated systems may still be achieved by low dimensional models, which allow for more than numerical analysis and, thus, a deeper understanding of the tumor development and treatments.