Continuum versus discrete model: a comparison for multicellular tumour spheroids

Gernot Schaller, Michael Meyer-Hermann

Abstract

We study multicellular tumour spheroids with a continuum model based on partial differential equations (PDEs). The model includes viable and necrotic cell densities, as well as oxygen and glucose concentrations. Viable cells consume nutrients and become necrotic below critical nutrient concentrations. Proliferation of viable cells is contact-inhibited if the total cellular density locally exceeds volume carrying capacity. The model is discussed under the assumption of spherical symmetry. Unknown model parameters are determined by simultaneously fitting the cell number to several experimental growth curves for different nutrient concentrations. The outcome of the PDE model is compared with an analogous off-lattice agent-based model for tumour growth. It turns out that the numerically more efficient PDE model suffices to explain the macroscopic growth data. As in the agent-based model, we find that the experimental growth curves are only reproduced when a necrotic core develops. However, evaluation of morphometric properties yields differences between the models and the experiment.

Keywords:

1. Introduction

In many medical applications, multicellular tumour spheroids constitute a popular experimental model system, since they are closely mimicking avascular tumours within living tissues. This opens the possibility to test the effect of, for example, chemotherapeutic agents under more realistic conditions (Santini & Rainaldi 1999; Ward & King 2003). Tumour cell lines tend to grow indefinitely in two-dimensional culture, which is mainly due to the unlimited supply of nutrients such as oxygen and glucose. Cells situated at a free boundary of the population will always have enough space and nutrients to follow the normal course of the cell cycle. However, if there is not enough space available, cells cease to proliferate—a phenomenon known as contact inhibition. Experimental evidence for the mitosis-inhibiting function of stress has been quantified (Helmlinger et al. 1997). In two dimensions, the total cell number is, in principle, unlimited. This is different, however, if the cells are cultivated in a three-dimensional setup (Folkman & Hochberg 1973). Here, the process of contact inhibition is cooperating with the lack of nutrients emerging in the interior of the cell population (Freyer & Sutherland 1986) and also with the loss of cells at the spheroid surface due to shedding (Landry et al. 1981). Usually, a lack of nutrients causes cells to undergo apoptosis or necrosis—processes that differ in many aspects (Bell et al. 2001), but in terms of population dynamics result in the same outcome: cell death. In contrast to the two-dimensional situation, in three dimensions the emergence of necrosis due to insufficient nutrient supply may lead to a substantially faster transition from the initial exponential growth to polynomial growth, which eventually even settles into saturation (Folkman & Hochberg 1973; Freyer & Sutherland 1986). Multicellular tumour spheroids have a very common structure within their late stages: a core consisting of mainly necrotic cells is surrounded by a viable layer of quiescent cells, which, in turn, is surrounded by a layer of proliferating cells (Greenspan 1972; Mueller-Klieser et al. 2002; see also figure 1). Even such a simple model represents a highly coupled system of many involved factors. This holds true for tissue modelling in general; unfortunately, there is no established general method to treat such systems.

Figure 1

(a) Cross-section through a three-dimensional multicellular tumour spheroid, taken with kind permission from Mueller-Klieser et al. (2002). The inner necrotic core is surrounded by a layer of quiescent cells, which, in turn, is surrounded by a layer of viable cells—the bar represents 250 μm. (b) This qualitative picture can be reproduced in an agent-based model (Schaller & Meyer-Hermann 2005), calculated for 0.28 mM oxygen and 16.5 mM glucose concentrations, and cell diameters of 10 μm.

From the viewpoint of mathematical modelling, the dynamics of the nutrient molecules can conveniently be described by reaction–diffusion equations (Araujo & McElwain 2004). A minor complication is represented by the fact that the diffusion coefficient for larger molecules such as glucose may vary considerably between water and tissue (Casciari et al. 1988). Also, it is not entirely clear on which factors the cellular nutrient uptake rate depends. Experiments point to the direction that these rates may strongly depend on the local nutrient concentration and the local pH value (Casciari et al. 1992).

The situation is more complicated for the cellular dynamics, which is not fully understood at present and, consequently, the used mathematical approaches differ considerably. Within the early approaches, the dynamics of the tumour boundaries is investigated. The model proposed by Greenspan (1972) already sufficed to explain the saturation of growth observed in experiments but was constrained to the previously required morphology (as are many models used for data fitting, see, for example, Groebe & Mueller-Klieser 1996). In contrast, further continuum models treated the cellular movement as a dynamic cellular density (Preziosi 2003; Araujo & McElwain 2004). In most of the theoretical investigations, tumour cells are assumed to move passively, such that the only driving force are stresses resulting from local volume changes due to proliferation, cell growth and cell death (compare Ward & King 1997; Byrne & Preziosi 2003; Schaller & Meyer-Hermann 2005). This is different in other models, where, for example, active chemotactic tumour cell movement has been assumed (Dormann & Deutsch 2002) to explain the experimental observation (Dorie et al. 1982) of viable tumour cells moving towards the necrotic region, or also in Matzavinos et al. (2004), where the movement of tumour-infiltrating cytotoxic lymphocytes is influenced by a chemotactic reaction to a chemokine. Within the continuum approach (Ward & King 1997), two dynamic distributions representing viable and necrotic cells, respectively, have been considered, whereas in an extension, the necrotic cells are rather viewed as a fluid constituted of necrotic debris (Ward & King 1999a). Using the balance equations for mass and momentum, the authors use the local volume changes arising from proliferation and cell death to derive a closed set of equations for the cellular dynamics (Byrne et al. 2003). Interestingly, the resulting equations revert to nonlinear diffusion equations (Breward et al. 2002; Byrne et al. 2003; Byrne & Preziosi 2003)—that had previously been used for animal dispersal (Murray 2002) and tumour growth (Sherratt 2000; Sherratt & Chaplain 2001)—in special cases.

Continuum models neglect the effects of cellular discreteness, even though these may result in qualitatively different population dynamics even in simple Lotka–Volterra systems (Shnerb et al. 2000; Bettelheim et al. 2001). Such effects can be taken into account by cellular automaton (von Neumann 1966) approaches, where cells can either occupy single lattice nodes (Dormann & Deutsch 2002) or several lattice nodes (Turner & Sherratt 2002). In such systems, the interaction rules between the different lattice nodes can be designed intuitively from existing biological observations, especially in the single-node case. Such lattice descriptions have some disadvantages: the predefined lattice geometry violates the isotropy of space. Though the emergence of lattice artefacts can be avoided by introducing stochastic interaction rules, it is still quite difficult to relate the effective parameters of cellular automata to physically relevant quantities. A further class of agent-based models, however, does admit continuous cellular positions—thereby allowing knowledge on microscopic physical interactions to be directly included (e.g. Drasdo 2003; Drasdo et al. 2003; Dallon & Othmer 2004; Schaller & Meyer-Hermann 2004, 2005). However, in order to do this consistently, an intrinsic cellular shape, such as, for example, spherical (Drasdo & Loeffler 2001; Galle et al. 2005; Schaller & Meyer-Hermann 2005), ellipsoid-like (Palsson & Othmer 2000), or polygonal (Honda et al. 2000), must be assigned, which usually makes these models specific to particular cell types. As an advantage, all deviations from the intrinsic shape can be treated as perturbations that give rise to physically motivated interactions, such as repulsive forces due to volume conservation or attractive forces due to surface energy minimization.

All the mentioned models allow the qualitative description of avascular tumour growth. However, this does not translate to a quantitative level and, in particular, the limits of the different methods are not sufficiently understood. In this article, we will as a first step apply a continuum model to multicellular tumour spheroid data in order to explore the advantages and disadvantages of continuum reaction–diffusion models in comparison with an analogously constructed agent-based model (Schaller & Meyer-Hermann 2005). In accordance with the experiment (Freyer & Sutherland 1986), we will consider the dynamics of oxygen and glucose simultaneously.

2. The model

Our model is described by the following partial differential equations (PDEs) of the reaction–diffusion type:Embedded Image(2.1)where Embedded Image represent the concentrations of oxygen and glucose and the densities of viable and necrotic cells, respectively.

The coefficients Embedded Image represent the effective diffusivities of oxygen and glucose. Especially for larger molecules such as glucose in real tissue, the presence of cells and extracellular matrix may significantly change the effective diffusion coefficient for the nutrients (Casciari et al. 1988). This effect is also existent for oxygen (Grote et al. 1977), though to a much smaller extent. Many theoretical studies consider oxygen as limiting nutrient and assume a constant diffusion coefficient (Ward & King 1999b; Breward et al. 2002). Here, we want to treat both nutrients consistently by assuming the simple linearized relationship,Embedded Image(2.2)where Embedded Image and Embedded Image are the measured diffusion constants in water and tissue, respectively, and the cellular threshold density is just the inverse cell volume Embedded Image, where the correction factor 0.74 arises from the closest packing of spheres. The diffusivities for oxygen and glucose become implicitly space and time-dependent, as the populations of necrotic and viable cells evolve. Note also that a nonlinear dependence would be possible, but theoretical investigations (Alvarez-Ramirez et al. 1996) suggest that the assumed linear dependence provides a good approximation.

As assumed, for example, in Landman & Please (2001) and Chen et al. (2001), to enable analytical investigations, we (for simplicity and consistency with Schaller & Meyer-Hermann 2005) assume cells to consume nutrients at constant rates Embedded Image and Embedded Image—regardless of the local nutrient concentration. This contrasts with other models (Casciari et al. 1992; Ward & King 1997, 1999a), and, hence, the values for Embedded Image and Embedded Image should be regarded as values averaged over all configurations occurring within the tumour spheroid.

Starting from mass balance equations for some medium C of the form Embedded Image, where Embedded Image denotes the velocity of the medium and Embedded Image the local net production rate (compare Landman & Please 2001), one recovers the law of diffusion by assuming that the velocity is directly proportional to the gradient of concentration Embedded Image. In other theoretical studies (Breward et al. 2002; Byrne et al. 2003; Byrne & Preziosi 2003), this direct proportionality is effectively replaced by cellular interactions, but also from these equations a nonlinear diffusion equation arises as a special case. Thus, the cellular diffusion coefficient Embedded Image can be regarded as a rough measure for describing at least passive cellular mobility. It may be influenced by the total cellular density Embedded Image as well. In analogy to a model for animal dispersion as presented in Murray (2002), we chooseEmbedded Image(2.3)where Embedded Image is assumed, compare also Sherratt (2000) and Sherratt & Chaplain (2001) for a similar approach. It should be noted that the above choice does not directly incorporate effects of adhesion, as the diffusion coefficient will always be positive. In the case m>0, the above choice ensures that the net diffusion coefficient will increase if the cellular density exceeds the dense packing (thus effectively modelling intercellular repulsion), and that it will decrease if the total cellular concentration is very small (indirectly mimicking intercellular adhesion). Naturally, for m=0, the normal reaction–diffusion equation is recovered.

The coefficient Embedded Image denotes the cellular proliferation rate. In the simplest approach, it should reflect at least the phenomenon of contact inhibition (Helmlinger et al. 1997). Via a substantial decrease in cell size (Roose et al. 2003), the external stress leads to decreased proliferation rates. This has already been analysed using agent-based modelling approaches (Galle et al. 2005). As a simple choice for the functional dependence, we consider a continuous dependence of the proliferation rate on the available volumeEmbedded Image(2.4)where Embedded Image is the largest cellular proliferation rate that can be observed in the experiment under consideration (Freyer & Sutherland 1986). The constant Embedded Image can be related to a different agent-based model (Galle et al. 2005) viaEmbedded Image(2.5)where the parameter Embedded Image denotes the critical compression, above which cell proliferation is contact-inhibited. The thus assumed direct relation of cell compression and proliferation differs from other approaches (Landman & Please 2001; Chen et al. 2001; Byrne & Preziosi 2003), where the net proliferation rate solely depends on the local nutrient supply by becoming negative for low nutrient concentrations to model apoptosis and/or necrosis. In addition, there the necrotic regions are not defined by the existence of necrotic cells, but rather defined using a balance of cell and necrotic fluid pressure. The function of the compression is then to transport cells from regions with a net positive growth rate to regions with a net negative growth rate, which facilitates global growth saturation for some parameters.

We have assumed necrotic and viable cells to move in the same way. This is due to our view of necrotic cells as inanimate balloons still having similar mechanical properties as viable cells, but neither exhibiting proliferation nor consumption of nutrients. The transformation of viable cells into necrotic ones is here solely determined by the local nutrient concentration. In order to be consistent with our previous agent-based attempt (Schaller & Meyer-Hermann 2005), we chose the local product of both concentrations to be the critical parameter inducing necrosis:Embedded Image(2.6)where Embedded Image is a transition rate and Embedded Image denotes the minimum nutrient concentration product to sustain cellular life functions without triggering necrosis or apoptosis. Note that in the analogous agent-based model approach, the parameter Embedded Image has no direct correspondent, since, there, cells turned necrotic immediately (i.e. after one time-step of, on average, 10 s) and also deterministically. Under the precondition that nutrient uptake rates are independent of the local nutrient concentration, equation (2.6) was the simplest assumption that allowed us to fit four growth curves for different nutrient supply (Freyer & Sutherland 1986) simultaneously in the agent-based model (Schaller & Meyer-Hermann 2005). Note that necrosis is a faster process than cell proliferation, which implies that Embedded Image should hold; assuming a necrosis duration of roughly 3 h (Sengelaub & Finlay 1982) and a minimum cell doubling time of 15 h (Freyer & Sutherland 1986), one would estimate Embedded Image versus Embedded Image. Also, if necrosis shall effectively reduce the local cell number, the inequality must hold strictly, since for low cell concentrations and poor nutrient support, the net proliferation rate is given by Embedded Image. Necrotic cells decay, thereby expelling their content. In living organisms, this leads to processes such as inflammation, whereas in in vitro experiments their content is passed into the fluid phase. The parameter γ represents the necrosis removal rate that has been fixed to Embedded Image in order to compare with the agent-based model (Schaller & Meyer-Hermann 2005). There, this value had been chosen in order to obtain a reasonable cell density within the necrotic core (compare figure 1). Within many other models, necrotic cells are not considered a separate transient state but are immediately assigned to the fluid phase (compare Ward & King 1999a; Landman & Please 2001; Breward et al. 2002; Byrne et al. 2003). Here as an approximation, the pressure of this fluid phase is assumed to be equilibrated at all times, therefore its existence has no effect other than the sink term Embedded Image.

3. Discussion of the model equations

Since the coefficients α and β depend on the concentrations themselves, it is obvious that even for m=0, the model equations (2.1) constitute a nonlinear system that may give rise to non-trivial dynamics. Some model parameters can be determined from independent experiments and others have to be estimated by fitting the model to experimental data (compare table 1). This is especially interesting for the nutrient uptake rates Embedded Image and Embedded Image and the critical compression Embedded Image, since it is difficult to measure these parameters directly. In addition, other hypotheses concerning the functional form of α and β and the value of m can, in principle, be tested to see whether they provide better agreement with the experimental data.

View this table:
Table 1

Model parameters. (Model parameters that have reproduced the best fit to experimental growth curves—compare figure 5—are given in the first and second columns. For the first number in brackets, the uncertainties represent estimates for upper parameter bounds, within which the calculated Χ2 should not differ by more than 5% from the minimum value given. The corresponding Hessian matrix at the minimum has been estimated by varying the parameters over a range of 25%, compare appendix B. For the other values in parentheses in the second column, the exponent m has been fixed to m=0. If applicable, the third column contains the analogous model parameters of the agent-based approach (Schaller & Meyer-Hermann 2005). Note that in this respect the cellular radius of 5 μm as assumed in the agent-based model directly corresponds to the value used for Cthresh if maximum cellular packing density of 74% is assumed. Within the last column, further literature containing information about the parameters is gathered and the type of the parameter is specified further; for more detailed explanations, see the text.)

Since most tumour spheroids are approximately spherical (compare Santini & Rainaldi 1999), we assume, in accordance with many other authors (e.g. Greenspan 1972; Ward & King 1997; Breward et al. 2002; Byrne & Preziosi 2003), spherical symmetry in all equations, which simplifies the numerical solution of (2.1) considerably, for details see appendix A. It is possible to extract simple characteristics of the model from direct observation of (2.1) under radial symmetry. By assuming that the solutions for large radii will travel outwards without changing shape, one can transform the system (2.1) to travelling wave coordinates by replacing the dependence on r and t by Embedded Image, where v is the wavespeed that shall be determined as part of the solution (compare Ward & King 1997; Sherratt 2000; Murray 2002). In particular, the equation for the viable cells then takes the formEmbedded Image(3.1)where Embedded Image denotes differentiation with respect to z and Embedded Image is obtained from equation (2.3). Usually, one will start with initial conditions, where there is an identically vanishing population of necrotic cells, a localized distribution of viable cells and nutrients in abundance. In analogy to Ward & King (1997), one can then for the special case m=0 (implying Embedded Image) and within regions with sufficiently large nutrient concentrations (such that Embedded Image and sufficiently small cell concentrations (such that Embedded Image) find a solution for the above equation that propagates outward with the wavespeed (compare also Murray 2002 and references therein) Embedded Image, which is also confirmed in the numerical solution (see figure 2).

Figure 2

The concentration of viable cells resembles travelling waves—here calculated for Embedded Image, Embedded Image, Embedded Image and Embedded Image. The curves advancing from left to right represent densities of viable cells after 4.63, 9.26, 13.89, 18.52 and 23.15 days. For m=0 (solid dark lines), the asymptotic maximum wave velocity (obtained for large radii and times) is 19.6 μm d−1, which is also found in the numerical solution. By a stability analysis of (3.1), one can show that the concentration must reach the critical density Embedded Image, which has been marked together with the threshold density Embedded Image by dashed horizontal lines. For m=2 (dashed grey lines), the apparent wave velocity of the numerical solution decreases to Embedded Image. In addition, the cells tend to aggregate leading to more localized distributions, reflected by the steeper decaying cell front.

In the case of m>0, the situation is more difficult. The general structure of equation (3.1), however, suggests that in the vicinity of m=0 no singularities occur, and therefore one can expect that the solutions will exhibit a continuous transition. This is well reflected in figure 2: for m=0, the normal wave velocity is recovered, whereas for m>0, a travelling wave is still found but the wave front is steeper and the apparent wave velocity is decreased considerably. Note that the existence of a travelling wave with lower bounds on amplitude and width will reproduce the experimentally observed transition from initial exponential to polynomial growth correctly (see also Ward & King 1997). However, a saturation in terms of the total viable cell number,Embedded Image(3.2)can in such a case not be observed, since the radial volume element would still lead to polynomial growth. This inability to predict growth saturation is a limitation of the present continuum model (compare also Sherratt 2000, where a minimum wavespeed is derived). More involved models—where cell movement is not described by positive definite diffusion but where the mass balance equations include specific interactions and the cellular material created by necrosis may either leak outwards into the external matrix or can be reused by the viable cells—predict growth saturation for some regions in parameter space (Ward & King 1999a; Landman & Please 2001; Breward et al. 2002; Byrne & Preziosi 2003).

4. Results

As for the computational domain, a sphere with radius of 500 μm—divided into 100 concentric shells—and Dirichlet boundary conditions for all populations have been considered. The choice of Dirichlet boundary conditions for the nutrient concentrations is motivated by the experiment (Freyer & Sutherland 1986), where the measured nutrient concentrations outside the spheroids remained approximately constant between the periodic refills of nutrient. Dirichlet boundary conditions on the cell population also emerge from multiple-phase modelling approaches (e.g. Breward et al. 2002; Byrne et al. 2003). In addition, we have verified that the assumption of vanishing Dirichlet boundary conditions for the cell populations did not make a difference to no-flux von Neumann boundary conditions in the observed time range of 25 days. Since the diffusion coefficients of oxygen and glucose are of orders of magnitude larger than typical diffusivities of cellular movement, the pseudo-steady-state approximation has been applied for the reaction–diffusion equations of the nutrients. This approximation relies on the fact that the corresponding diffusion length Embedded Image was always larger than the largest spheroid radius encountered, and is also adopted by many other authors (e.g. Ward & King 1997; Chen et al. 2001; Breward et al. 2002; Byrne & Preziosi 2003). The volume integral (3.2) of all densities yields the total cell numbers for viable and necrotic cell populations or the total number of glucose and oxygen molecules in the computational domain, respectively. The total number of viable cells has then been compared with experimental growth curves on EMT6/Ro tumour cells (Freyer & Sutherland 1986) by varying the model parameters. Following the assumption that whole spheroid populations in the experiment were grown from single tumour cells, all simulations have been started by setting the initial concentration of viable cells in the innermost (spherical) volume compartment, such that the total initial cell content was 1 (compare also Ward & King 1999b). It should, however, be noted that there was no experimental control in Freyer & Sutherland (1986) that the harvested spheroids were monoclonal. In contrast, the described culturing methods make it highly probable that cells aggregate before forming a spheroid. We have also started our model with different initial cell numbers and found that the growth curves were mainly shifted in time. From the mathematical point of view, a model describing Embedded Image cells by continuum densities has certainly been extrapolated beyond the region of its validity. However, the associated error in the number of cells is limited, as one can see by directly comparing with discrete exponential growth. In a first run, with m=0 being fixed, the deviation from experimental data has been minimized by varying the model parameters using Powell's method (Press et al. 1994), see appendix B. In a second run, m has been varied as a fit parameter as well. Powell's method has been started from different parameter sets and the best minimum has been kept. Figures 3–5 relate to these two fits and the corresponding parameters can be found in table 1.

Figure 3

Cell densities for different nutrient concentrations after 23 days of simulation time. The dark bold lines denote simulations, where the exponent m has been determined as a fit parameter, whereas the grey thinner lines represent simulations with m=0. The cellular threshold density has been marked by a horizontal light grey dashed line. In each panel, the two dashed vertical lines correspond to the largest distance (cell centre) of a necrotic and a viable cell from the spheroid centre within the agent-based model (Schaller & Meyer-Hermann 2005). Thus, they constitute an upper bound on the radii of the necrotic core and the tumour spheroid, respectively, in the agent-based approach. (a) Cellular distributions for 0.07 mM oxygen and 0.8 mM glucose concentrations. As the concentrations never reach Embedded Image, contact inhibition does not play a role within this scenario. (b) Cellular distributions for 0.07 mM oxygen and 16.5 mM glucose concentrations. A thin layer of partially contact-inhibited cells emerges. (c) This is similar for 0.28 mM oxygen and 0.8 mM glucose concentrations. (d) Cellular distributions for 0.28 mM oxygen and 16.5 mM glucose concentrations. As one would expect, the spheroid exhibits a thick layer of contact-inhibited cells.

Figure 4

Glucose concentration calculated for 0.28 mM oxygen and 16.5 mM glucose boundary concentrations and m being determined as fit parameter. During tumour growth, the concentration of glucose drops rapidly as soon as the density of viable cells reaches significant concentrations. The region of steepest descent coincides with the travelling localization of viable cells (rescaled population density shown after 23 days, thin full line).

Figure 5

Temporal dynamics of the number of viable cells for different nutrient concentrations. All simulations have started with a single cell. Symbols represent data read off from Freyer & Sutherland (1986). Solid lines represent simulations with m being used as fit parameter, whereas for the dashed lines (only contrasting for nutrient starvation) m=0 has been fixed.

(a) Qualitative analysis

As can be expected for the full system (2.1), in the case of m being fixed to 0, the cellular concentration resembles a travelling wave, with amplitude and width depending on the nutrient concentrations. In addition, since now cell death is allowed, the total cell density at a given distance first rises to a maximum value depending on the nutrient concentrations and then decays due to the removal of necrotic cells. Generally, the wave of viable cells is followed by a wave of necrotic cells as soon as the local nutrient concentrations fall below the critical threshold Embedded Image (see figure 3). Since in the model (2.1) necrotic cells decay, the population of necrotic cells soon assumes a wave-like shape as well. For m=0, the velocity of this wave solely depends on the cellular diffusion coefficient and the maximum proliferation rate (compare figure 2), and therefore it is equal between test runs with different nutrient concentrations. Consequently, the calculated spheroids will in this case have the same macroscopic size for all nutrient concentrations, as is confirmed in figure 3. Half of the maximum spheroid cell density is reached at a radius of 350 μm—regardless of the nutrient concentration. In this model ansatz, the differences in the overall cell number will then result from the local density amplitude only. The assumption m=0 is therefore in direct contradiction with experimental data and the agent-based model (Schaller & Meyer-Hermann 2005), where considerable size differences are found for spheroids grown in different nutrient concentrations.

At a first glance one might expect from figure 2 that if one allows for a variable diffusion exponent m>0, the effective velocity of the travelling waves should decrease considerably. However, in figure 3 this is not the case, except for poor nutrient support. For m>0, the spheroid grows to a radius of 300 μm, whereas in the case of nutrient abundance the spheroid front can be found at 380 μm. In cases where only one nutrient is limited, the spheroids are negligibly smaller at roughly 370 μm. One should keep in mind that both the exponent m and the diffusion constant Embedded Image have been varied as fit parameters in this scenario. With the resulting increased exponent m, the apparent diffusion coefficient will be smaller than the diffusion constant Embedded Image for cell densities below the threshold density Embedded Image, which is especially important for the advancing tumour front. Therefore, the net effect of an increased exponent m will at the tumour boundary be a decreased apparent diffusion coefficient. This tendency, however, is balanced by an increased diffusion constant Embedded Image. In order to fit the overall cell numbers to experimental data, the apparent diffusion coefficient at the tumour front is similar for both scenarios (m fixed and m varying). The local cellular density is bound by Embedded Image and the only way to harbour enough cells within the tumour spheroid is to increase the velocity of the propagating wave. Therefore, it does not come as a surprise that the macroscopic spheroid sizes in the scenarios where the threshold density Embedded Image is reached (see, for example, figure 3b–d) do not differ very much. The difference between m=0 and m variable becomes manifest only in the case of nutrient starvation, where contact inhibition does not play a role at all. The resulting size differences of the spheroids reproduce the agent-based model (Schaller & Meyer-Hermann 2005) and the experimental data (Freyer & Sutherland 1986) much better, even though the qualitative agreement is still not satisfactory. For example, the spheroid size differences in the agent-based model are also hardly reproduced in nutrient scenarios of the continuum approach, where the threshold density is reached. A possible improvement resulting from the choice of m>0 is the more sharply pronounced transition at the spheroid boundary, since this agrees much better with experimental observations and the agent-based model (compare figure 1). However, it should also be kept in mind that the experimental procedure of obtaining the cross-sections of spheroids might remove cells that are only loosely bound to the spheroid aggregate, and thereby may also lead to a flawed picture of realistic tumour spheroids.

As in the agent-based model (Schaller & Meyer-Hermann 2005), different mechanisms dominate the population dynamics for different nutrient concentrations (figure 3). In the case of nutrient starvation, the cellular concentrations never reach densities that might induce contact inhibition, and the only mechanism inhibiting cellular growth is necrosis. In all other test configurations, the cellular density locally exceeds the threshold density Embedded Image. Cells residing within regions where the cellular concentration lies above Embedded Image have a decreased proliferation rate, and are therefore interpreted as quiescent within this model. Consequently, both contact inhibition and necrosis play an important role in all other scenarios. The first initially dominates if nutrients are provided in abundance, whereas the latter dominates if nutrients are rare.

The model will always predict the emergence of a necrotic core as soon as a sufficient number of cells has been developed. As quiescent cells also consume nutrients here, it is to be anticipated that the number of quiescent cells will eventually be outstripped by the number of necrotic cells. Also in this respect, the present continuum model poorly differentiates between the scenarios where only one of the nutrients is provided in smaller quantities. For both scenarios, a thin contact-inhibited layer followed by a large necrotic core is predicted. In contrast, in the agent-based model (Schaller & Meyer-Hermann 2005), contact inhibition is much more pronounced when glucose is provided in abundance and oxygen is restricted in comparison to the case where more oxygen than glucose is provided. To our knowledge, however, the quantitative experimental signature is too poor to differentiate between the models in this aspect.

Generally, in order to correctly fit the global cellular growth curves, a necrotic core always emerges. Therefore, like in the agent-based model, this fact can be considered a qualitative prediction, since the emergence of a necrotic core is not a stringent ingredient of the model. The sensitivity of the fit routine is also underlined by the fact that when the cellular density approaches the maximum carrying capacity, the minimization procedure leads to cellular diffusivities that produce similar apparent wave velocities—regardless of the value of m.

Note that for all curves in figure 3, the continuum approach overestimates the radii of the spheroids. This may be caused by the positive-definite diffusion coefficient (2.3), which does not describe adhesion correctly. In more realistic models (Landman & Please 2001; Byrne & Preziosi 2003), the relation between the apparent diffusion coefficient and adhesive and repulsive parameters becomes explicit. However, even in these multiphase approaches, the parameters for the cellular interactions remain effective continuum parameters. It would be interesting to see how individual single-cell-based parameters relate to these effective continuum values.

When considering only the total growth curves, the visual difference between m=0 and m>0 is negligible (see also figure 5). Therefore, a more refined analysis makes sense when improved quantitative experimental signatures become available.

Since the evolution of the nutrients is coupled to the cellular distribution, the concentration of nutrients will fall, depending on the concentration of viable cells (see figure 4 for glucose and variable m). The region of steepest descent coincides with the localization of viable cells, which is also found in other models (compare Ward & King 1997; Sherratt & Chaplain 2001). A similar relation is found for the oxygen concentration. However, due to the larger diffusion constants in both water and tissue, the gradients tend to be smaller (data not shown). Also for m=0, the qualitative appearance hardly changes. Owing to the slightly increased glucose uptake rates (see table 1), the baselines are a bit lower than in figure 4 (data not shown).

(b) Quantitative analysis

The total cell numbers in figure 5 demonstrate that the implicit assumption of the nutrient concentration product being a critical candidate for necrosis assumption suffices to explain experimental data within the level of their own uncertainty—independent of the value of m. The difference in the choice of m reflects itself only in the curve for 0.07 mM oxygen and 0.8 mM glucose, but is much more pronounced in the corresponding cell density distribution in figure 3.

The model is not equally sensitive to changes of the parameters. Some parameters (such as diffusion coefficients of oxygen and glucose in tissue) have been determined in independent experiments (Grote et al. 1977; Casciari et al. 1988). It should be kept in mind that diffusion coefficients generally depend on temperature and on the (likewise temperature-dependent) viscosity of the used medium. Since oxygen is a small molecule, its diffusion coefficient is hardly affected by the presence of additional boundary conditions, such as cell membranes. We have also tried how a uniform oxygen diffusion constant would change the results by starting a fit with Embedded Image. This mainly resulted in a decreased oxygen uptake rate (Embedded Image) and a decreased minimum nutrient product (Embedded Image). This was accompanied with an increase of the minimum Embedded Image-value by 2%, which hardly made a visual difference. Thus, the ansatz used in most of the literature to assume the diffusive properties of oxygen as being hardly changed by the presence of tissue is well justified. This is different if glucose is used as the limiting nutrient. However, when moving boundary conditions are used at the spheroid surface, this difference will not manifest itself within spheroids where the necrotic core and the viable layers are well separated. In order to lead to a different outcome, the support of the glucose consumption terms must overlap with regions of non-vanishing glucose gradients.

Generally, the proliferation rate (for the determination of Embedded Image) of tumour cells may vary considerably even within a well-known cell line. It may be speculated that mutated clones of previously well-defined cell lines may have different proliferation rates. Also, in addition to the previously discussed problem of poorly controlled initial conditions, in many experiments the culturing conditions may also not be under good control. Therefore, it does not come as a surprise that the cell doubling times obtained from the literature cited in table 1 range between 13 and 17 h (Casciari et al. 1992).

The necrosis removal rate γ has been fixed in order to be consistent with the agent-based model (Schaller & Meyer-Hermann 2005). There, the parameter value had been chosen in order to obtain reasonable cell densities of the necrotic core (compare also figure 1b). In the agent-based model with moderate cell adhesion, the number of viable cells was hardly sensitive on this parameter. In the present model (where adhesion only decreases but does not revert the diffusion), this is similar. Since the fit is done to the number of viable cells and since these are only weakly coupled to the density of necrotic cells, the quality of the fit to the number of viable cells is hardly sensitive to γ over several orders of magnitude. Therefore, for the detailed fixation of the parameter γ, quantified data on the densities of the necrotic cores would be necessary.

Independent of a varying exponent m, the obtained diffusion constants Embedded Image are considerably larger than those in the agent-based models (Galle et al. 2005; Schaller & Meyer-Hermann 2005), where values of about Embedded Image have been used. This is related to the fact that in the present continuum model, the diffusion constant also has to include the effects of intercellular repulsion, whereas in the agent-based models the diffusive movement was represented by additional stochastic forces that were added to the (deterministic) cell–cell interactions. The obtained value for the apparent diffusion coefficient is well within the range proposed in Matzavinos et al. (2004) for continuum approaches, where for tumour cells with radii between 4 and 15 μm and proliferation rates between 0.1 and 1.0 d−1, the estimated apparent diffusion coefficient (resulting from movement and proliferation) ranges within Embedded Image. The quality of the fit to the overall cell number can be increased by leaving the exponent m as a fit parameter. However, as has been argued before, the model is very sensitive to the cellular diffusion constant (compare also the small variation width in the table) via the cell density distribution: Embedded Image is automatically tuned to obtain the necessary wave velocity, which manifests itself in a threefold increase in Embedded Image when m increases from m=0 to 0.73 (see table 1). Since the movement-related parameters m and Embedded Image are used in (2.3) to mimic elastic and adhesive cellular properties, one can state that the agent-based models have an advantage in this respect, since they inherently incorporate experimentally accessible parameters on a single-cell basis. However, the existence of a continuum model with similar results proves that on this macroscopic scale (by comparing the total cell numbers), a discrimination is possible neither between the PDE models of constant and varying diffusivity nor between continuum and agent-based models. In the present example of multicellular tumour spheroids, in addition the morphology must be compared. This is only possible if it is given in quantitative form, i.e. reliable measurements including error estimates of relative sizes of the quiescent layer, the necrotic core and the proliferating rim under well-defined initial and boundary conditions.

The necrosis transition rate Embedded Image mainly determines the width of the region where viable and necrotic cells coexist, and has little impact on the produced cell number over several orders of magnitude (compare also table 1) if the intrinsic requirement Embedded Image is respected. Therefore, the fact that the derived necrosis transition rate Embedded Image—implying transition durations of roughly 2 h—is well within the expected range (Sengelaub & Finlay 1982) should be interpreted as a direct consequence of choosing this value as a starting position for the fitting procedure rather than a prediction of the model.

The critical compression Embedded Image primarily influences the growth curve when nutrients are supplied in abundance. In Schaller & Meyer-Hermann (2005), the criterion used for initiating contact inhibition is based on normal tensions instead of cellular compression, and can therefore not directly be compared. The present model predicts contact inhibition to set in at a much stronger compression as other agent-based models (Galle et al. 2005), where a value of Embedded Image is reported. This may be due to two reasons: first, the dense packing—as assumed by taking the correction factor 0.74 into the definition of Embedded Image—might neither be realized in realistic multicellular tumour spheroids (Freyer & Sutherland 1986) nor in the agent-based model. Then, it is clear that the crude ansatz Embedded Image cannot be correct. Second, one should keep in mind that the pressure produced by the proliferating rim is not correctly transmitted with a diffusion approach. For example, at least for strongly compressed tissue, the pressure distribution will obey a wave equation rather than a diffusion equation that can even be solved in equilibrium (compare Chen et al. 2001; Roose et al. 2003). The approach used here is only valid for highly dissipative systems. The microscopic interactions inherent to the agent-based model lead—in combination with an adaptive time-step—to a more realistic stress relaxation.

The fit is very sensitive to the last three parameters in table 1. Though all these parameters have the same order of magnitude as in Schaller & Meyer-Hermann (2005), the glucose uptake rate Embedded Image is by a factor of 2–3 smaller in the present model. It is not surprising that a smaller glucose uptake rate than in the agent-based model must come with an increased nutrient threshold Embedded Image to fit the experimental data. However, the continuum model—started with the parameters of the agent-based model—does not produce a fit of similar quality (data not shown). The difference may be a consequence of the different cell density distribution in both models, which is due to the replacement of cellular adhesion and repulsion by an effective and purely repulsive diffusion coefficient. The absolute value of the cellular glucose consumption rate Embedded Image compares with experimental data on EMT6/Ro spheroids (Wehrle et al. 2000)—where Embedded Image (n=6) have been measured—and also with other tumour cell lines (Walenta et al. 2000). Wehrle et al. (2000) report the original experimental uptake rates in Freyer & Sutherland (1985) to be in the range of Embedded Image. These values support the glucose uptake rate found in the agent-based model. A value for the consumption of oxygen is given in Wehrle et al. (2000), where Embedded Image is reported. Depending on the local glucose concentration, Casciari et al. report oxygen uptake rates ranging within Embedded Image and glucose uptake rates within Embedded Image for EMT6/Ro cells (Casciari et al. 1992).

If one neglects possible usage of glucose other than for reaction with oxygen to produce energy, a complete combustion of glucose would require the ratio of oxygen and glucose uptake rates to be 6 : 1 in the steady state. However, in the agent-based model, this is actually contradicted with ratios of about 1 : 5. In the present model, this surprising effect is still existent, but less pronounced with a ratio of 1 : 2. This should not be too disturbing, since it is well known that tumour cells do not completely combust glucose, thereby leading to acidification of the tumour environment. Experimental estimates concerning a different cell line point to a ratio of 1 : 1 (Kunz-Schughart et al. 2000), whereas Wehrle et al. (2000) report a ratio of 1 : 4 for EMT6/Ro cells. Usually, uptake rates of both oxygen and glucose will depend on the local nutrient concentration, and thereby indirectly on the EMT6/Ro spheroid size under consideration, which makes experiments necessary under better defined conditions.

5. Conclusion

It has been demonstrated that when considering macroscopic data, continuum models of the reaction–diffusion type completely suffice to explain global properties such as growth curves if saturation is not sufficiently quantified.

In most qualitative properties, the continuum model reproduces the analogous agent-based model. For example, the emergence of a necrotic core was necessary to fit the data correctly for all nutrient configurations. However, it could be shown that the necessity to fit experimental data considerably changed the cellular diffusion constant—depending whether linear or nonlinear diffusion was assumed. The assumption of a constant cellular diffusion coefficient (m=0) does not correctly reproduce the experimental observation that the spheroid size changes with the nutrient supply (Freyer & Sutherland 1986). A concentration-dependent diffusivity (m>0) does not display this drawback, but still cannot sufficiently reproduce varying spheroid sizes of the agent-based approach, where with sufficiently large adhesion even saturation could be reproduced. The used positive-definite diffusion coefficient therefore represents a limitation of the model. Instead, cellular interactions should be used to determine the concentration-dependence of the apparent diffusion coefficient (or flux) in the cellular mass balance equations.

Generally, one is unable to discriminate between microscopic models of cellular interactions on the macroscopic scale as long as quantitative morphological data of the tissue are not available. Since continuum models have the advantage of being easy to solve computationally, the first step towards more refined tissue models should be to analyse the phase space with a simple model. The agent-based models have then their advantage in the fact that microscopic cellular properties that have been determined from independent experiments (such as adhesion and repulsion) can be inherently included to yield improved predictions. In addition, as agent-based models can be built on physically realistic assumptions, they bear the potential to test the relevance of microscopic cell properties for tissue development.

Editors' note

Please see also related communications in this focussed issue by Byrne et al. (2006) and McDougall et al. (2006).

Acknowledgments

G.S. is indebted to T. Beyer for many valuable discussions. The authors would like to thank anonymous reviewers for valuable suggestions. M.M.-H. was supported by a Marie Curie Intra-European Fellowship within the Sixth EU Framework Programme.

Footnotes

  • One contribution of 15 to a Theme Issue ‘Biomathematical modelling II’.

References

View Abstract