Modelling the Internalization of Labelled Cells in Tumour Spheroids

Bulletin of Mathematical Biology (1999) 61, 601–623 Article No. bulm.1998.0089 Available online at http://www.idealibrary.com on

Modelling the Internalization of Labelled Cells in Tumour Spheroids K. E. Thompson Department of Mathematics, UMIST, PO Box 88, Manchester M60 1QD, U.K.

H. M. Byrne∗ Division of Theoretical Mechanics, University of Nottingham, Nottingham NG7 2RD, U.K. In this paper a mathematical model is developed to describe the migration of labelled particles within a multicell spheroid. In the model, spatial variations in cell proliferation and death create an internal velocity field which leads to redistribution of the labelled and unlabelled cells. By applying a range of numerical and analytical techniques to the model equations, it is possible to show that, whilst the speed with which the labelled cells migrate through the tumour is independent of the type of cells that are labelled, their limiting distribution depends crucially on whether inert polystyrene microspheres or live tumour cells are labelled. These predictions are shown to be in good qualitative agreement with independent experimental results. c 1999 Society for Mathematical Biology

1.

I NTRODUCTION

In their efforts to understand the complex mechanisms underpinning the development of solid tumours, experimental biologists have designed several assays which they use to ascertain various features of tumour growth. For example, classical assays for studying tumour angiogenesis include the hamster cheek pouch, the rabbit ear chamber and the chick chorioallantoic membrane (Cliff, 1963, 1965; Fan and Polverini, 1997). Due perhaps to its reproducibility, one of the most widely used assays for studying avascular tumour growth is the multicell spheroid system, in which the growth of clusters or aggregates of tumour cells is monitored. Typically, as the cells proliferate and the aggregate expands, cells towards the centre of the cluster become starved of vital nutrients, undergo hypoxia and eventually die. Thus a well-developed spheroid may comprise an outer ring of proliferating cells, an inner ring of viable, but nonproliferating hypoxic cells and a central core of necrotic matter. ∗ Author to whom correspondence should be addressed.

0092-8240/99/040601 + 23

$30.00/0

c 1999 Society for Mathematical Biology

602

K. E. Thompson & H. M. Byrne

Microspheres after 3 days

0.25 0

0

Microspheres after 4 days

0

0

50 100 150 200 250 300 Distance from outer edge νm

0

100

200

300

Labelled cells after 3 days

0.25 0

50 100 150 200 250 300

0.25

Labelled cells after 2 days

0.25 0

50 100 150 200 250 300

Fraction

Fraction

0

Fraction

0.25 0

Fraction

(b) Microspheres after 2 days

Fraction

Fraction

(a)

0

100

200

300

Labelled cells after 4 days

0.25 0

0 100 200 300 Distance from outer edge νm

Figure 1. Typical distribution of (a) polystyrene microspheres and (b) labelled tumour cells 2, 3 and 4 days after internalization [reproduced from McElwain and Pettet (1993) and Dorie et al. (1982; 1986)].

In the early to mid-eighties, Dorie et al. (1982, 1986), used the multicell spheroid assay to test whether parameters relevant to tumour growth kinetics could be estimated by tracking the movement of labelled cells within a growing spheroid. In separate experiments, inert polystyrene microspheres and tumour cells radioactively labelled with tritiated thymidine were allowed to adhere to the outer surface of multicell spheroids of approximately 300–400 µm in radius. A certain fraction of the labelled probes were internalized. Changes in the distributions of the two types of labelled probes over time were studied and typical results are presented in Fig. 1. Referring to Fig. 1(a) we observe that the microspheres appear to migrate towards the centre of the tumour in a wavelike manner, so that by day 4 no microspheres are discernible at the outer tumour boundary, i.e., complete internalization has occurred. Referring to Fig. 1(b) we remark that when labelled tumour cells are used once again a wavelike invasion front is discernible. By contrast with Fig. 1(a), however, complete internalization of the labelled tumour cells does not occur. Indeed, by day 4 the labelled cells appear to be fairly evenly distributed throughout the tumour volume. McElwain and Pettet (1993) developed a mathematical model to explain why the internalization patterns characterizing the microspheres differ from those characterizing the labelled tumour cells. They showed that pressure gradients, caused by nonuniform cell proliferation and death, could explain the wavelike pattern of internalization of the microspheres. In order to explain why the labelled tumour

Cell Internalization in Tumour Spheroids

603

cells were not completely internalized, McElwain and Pettet (1993) postulated that, unlike the inert microspheres, the tumour cells were sensitive to the local oxygen gradient, with cells migrating preferentially, via chemotaxis, towards the source of oxygen. In this paper we provide a different explanation of Dorie et al.’s (1982) results. Rather than invoking chemotaxis, we suggest that the different internalization patterns may be attributed simply to nonuniform cell proliferation and death of the labelled tumour cells. The model that we develop is similar in form to that studied by Ward and King (1997). We assume that the tumour is incompressible and comprises a mixture of unlabelled tumour cells and labelled cells, the latter being either microspheres or tumour cells. As in McElwain and Pettet’s (1993) model, mitosis and cell death create an internal velocity field, which leads to redistribution of the labelled and unlabelled cells and growth of the tumour itself. The key dependent variables involved in the model are therefore the nutrient (or, specifically, oxygen) concentration, the cell velocity, the tumour radius and the densities of the labelled and unlabelled cells. By choosing the cell birth and death terms appropriately we are able to distinguish between cases in which the tumour cells are labelled and those in which the microspheres are labelled. Indeed this is the key difference between our model and that of McElwain and Pettet (1993). They assumed that neither the microspheres nor the labelled tumour cells underwent cell birth and death. In addition, they did not treat the tumour as incompressible. The paper is organized as follows. In Section 2 the mathematical model is developed. Section 3 continues with numerical results. In Section 4 attention focuses on the long-time behaviour of the model equations. In particular, steady, timeindependent solutions of the model equations are constructed. The paper concludes in Section 5 with a summary of the main results and a discussion of the relative merits of the two labelling systems.

2.

T HE M ATHEMATICAL M ODEL

Following Adam (1987), Adam and Maggelakis (1990), Byrne and Chaplain (1998), and Greenspan (1972, 1976), we consider an avascular solid tumour whose growth is regulated by an externally-supplied nutrient. In order to mimic Dorie et al.’s experimental assay, we assume that the tumour comprises a dynamic mixture of unlabelled and labelled cells. The unlabelled cells are always live, proliferating tumour cells whereas the labelled cells are either inert polystyrene microspheres or proliferating tumour cells. Mitosis and cell death create an internal velocity field which causes the labelled and unlabelled cells to redistribute themselves within the growing tumour. Our aim here is to develop a model that reproduces Dorie et al.’s experimental results regarding the patterns of internalization of inert microspheres and tumour cells within an avascular tumour. At the same time the model should explain why the microspheres and the tumour cells give rise to different migration patterns.

604

K. E. Thompson & H. M. Byrne

The mathematical model that we study is stated below in dimensionless form and is similar in derivation to that developed by Ward and King (1997). For simplicity, a one-dimensional cartesian geometry is employed. This corresponds to the growth of a slab of tumour cells, rather than a spheroid. As we show, this does not affect the qualitative form of our results or predictions. We also assume that sufficient numbers of labelled and unlabelled cells are present within the tumour to justify the adoption of a continuum model, in which all physical variables depend on distance from the centre of the tumour x and time t. We denote the density of unlabelled tumour cells by n(x, t), the density of labelled cells by m(x, t), the cell velocity by v(x, t) and the nutrient concentration by c(x, t). These dependent variables are defined within the growing tumour whose outer boundary is located at x = R(t) at time t. Equations governing the evolution of n, m, v, c and R are presented below. They derive from the principle of conservation of mass. 2.1. Unlabelled tumour cells, n(x, t). We assume that the dominant factors governing the evolution of the unlabelled cells are random motion, convective transport (due to the underlying velocity field), cell proliferation and cell death. Introducing µn and Sn (c, n) to denote, respectively, the assumed constant random motility coefficient and the net proliferation rate of the unlabelled cells (cell proliferation rate minus cell death rate), we deduce that n(x, t) satisfies the following partial differential equation: ∂ 2n ∂n ∂ + (vn) = µn 2 + Sn (c, n). ∂t ∂x ∂x

(1)

2.2. Labelled tumour cells, m(x, t). Following (1), we assume that the factors governing the evolution of the labelled cells are random motion, convective transport and, where appropriate, cell proliferation and cell death. Introducing µm and Sm (c, m) to denote, respectively, the assumed constant random motility coefficient and the net proliferation rate of the labelled cells, we deduce that m(x, t) satisfies the following partial differential equation: ∂m ∂ 2m ∂ + (vm) = µm 2 + Sm (c, m). ∂t ∂x ∂x

(2)

As we explain below, the functional form used to describe Sm (c, m) depends on whether microspheres or tumour cells have been labelled. 2.3. Cell velocity, v(x, t). By assuming that the tumour, which comprises a mixture of labelled and unlabelled cells, is incompressible we deduce that n + m = constant at all points within the tumour. Adding equations (1) and (2) leads to the following partial differential equation which describes the evolution of the cell velocity v(x, t): (n + m)

∂v ∂ 2n ∂ 2m = µn 2 + µm 2 + Sn (c, n) + Sm (c, m). ∂x ∂x ∂x

(3)

Cell Internalization in Tumour Spheroids

605

2.4. Tumour radius, R(t). For the simple, 1D cartesian geometry under consideration, the outer tumour boundary grows at a rate which is equal to the velocity there. Thus the equation governing R(t) can be written: dR = v(R, t). dt

(4)

2.5. Nutrient concentration, c(x, t). As it is diffusing through the tumour, we assume that the nutrient is consumed at the rates 0n (n) and 0m (m) by the unlabelled and labelled cells respectively. Following Adam (1987), Byrne and Chaplain (1998), Greenspan (1972), McElwain and Pettet (1993) and Ward and King (1997), we make a quasi-steady approximation and assume that the nutrient satisfies the following reaction-diffusion equation: 0=

∂ 2c − 0n (n) − 0m (m). ∂x2

(5)

The neglection of time-derivatives in (5) reflects the fact that, because they are typically much smaller than the tumour cells, the nutrient molecules rapidly redistribute themselves to take account of any changes in the size of the tumour. 2.6. Reaction terms. Throughout this paper we assume that cell mitosis and apoptosis (or programmed cell death) contribute to the net proliferation rate of the unlabelled cells (Kerr et al., 1972; Folkman and Hochberg, 1973). For simplicity we assume further that the tumour does not contain a necrotic core. Introducing positive constants s and d, we assume that, for the unlabelled tumour cells, mitosis and apoptosis occur at the rates scn and dn, respectively. Thus the rate of cell proliferation is proportional to the local nutrient concentration whereas apoptosis occurs at the rate dn. In addition, we assume that the unlabelled cells consume nutrient at a rate 0n (n) = 0n where 0 is a positive constant. Combining these ideas, in equations (1), (3) and (5) we fix: Sn = (sc − d)n = (sc − d)(1 − m)

and

0n = 0n.

(6)

Our choices of Sm and 0m depend on whether the labelled cells are polystyrene microspheres or tumour cells. Being inert, the microspheres do not consume nutrient, proliferate or die. Hence if the polystyrene microspheres are labelled in equations (2), (3) and (5) we fix: Inert microspheres labelled Sm = 0 = 0m .

(7)

In contrast, when live tumour cells are labelled, it is natural to assume that their proliferation and nutrient consumption rates will be identical to those used for the unlabelled tumour cells. Hence, following (6), in this case we fix:

606

K. E. Thompson & H. M. Byrne

Tumour cells labelled Sm = (sc − d)m

and

0m = 0m.

(8)

2.7. Model simplification. In practice, it is possible to exploit the assumed incompressibility of the tumour to eliminate one of the dependent variables, n or m, say. Assuming that the model equations have been nondimensionalized so that n + m = 1, and for comparison with Dorie et al.’s experiments, it is appropriate to eliminate n = 1 − m. It is also convenient to reorder the governing equations. In this way the following simplified model can be obtained: 0=

∂ 2c − 0(1 − m) − 0m , ∂x2

∂ 2m ∂v = (µm − µn ) 2 + (sc − d)(1 − m) + Sm , ∂x ∂x dR = v(R, t), dt

(9) (10) (11)

∂m ∂m ∂ 2m +v = [µm − (µm − µn )m] 2 + (1 − m)Sm − (sc − d)m(1 − m), ∂t ∂x ∂x (12) where 0m and Sm are defined by (7) if microspheres are labelled and by (8) if tumour cells are labelled. 2.8. Boundary and initial conditions. Equations (9)–(12) are closed by imposing the following boundary and initial conditions: c = c∞ ∂c =0 ∂x v=0 µm

at x = R(t),

(13)

at x = 0,

(14)

at x = 0,

(15)

∂m − vm = 0 at x = 0, R(t), ∂x m(x, 0) = m in (x) and R(0) = R0 , prescribed.

(16) (17)

In equation (13) c∞ is the external nutrient concentration. Thus (13) states that the nutrient concentration is continuous across x = R(t). Equations (14), (15) and the first equation of (16) reflect the assumed symmetry of the tumour about x = 0. The second equation of (16) states that the flux of labelled cells across the outer tumour boundary is zero. Equations (17) define the initial distribution of labelled cells within the tumour and also its initial radius.

Cell Internalization in Tumour Spheroids

3.

607

N UMERICAL R ESULTS

Equations (9)–(17) define a moving boundary problem which may be solved numerically by mapping the equations onto a fixed spatial domain and using finite differences to construct approximate solutions to the transformed equations (Crank, 1988). In practice, the tranformation is effected by introducing 0 ≤ ρ = x/R(t) ≤ 1 and τ = t and rewriting equations (9)–(17) in terms of ρ and τ . Using the chain rule, it is possible to show that the partial derivatives transform as follows: ∂ ∂ρ ∂ ∂τ ∂ 1 ∂ = + = , ∂x ∂ x ∂ρ ∂ x ∂τ R ∂ρ ∂τ ∂ ∂ρ ∂ ∂ ρ dR ∂ ∂ ρv(1, τ ) ∂ ∂ = + = − = − , ∂t ∂t ∂τ ∂t ∂ρ ∂τ R dt ∂ρ ∂τ R ∂ρ where we have used (11) to substitute for ddtR . In this way we deduce, for example, that the transformed versions of equations (9) and (12) are given by: 0=

1 ∂ 2c − 0(1 − m) − 0m , R 2 ∂ρ 2 ∂m ∂ 2m 1 ∂m 1 + [v(ρ, τ ) − ρv(1, τ )] = 2 [µm − (µm − µn )m] 2 ∂τ R ∂ρ R ∂x +(1 − m)Sm − (sc − d)m(1 − m).

In the absence of accurate parameter estimates, O(1)-estimates were adopted for the dimensionless parameters needed to perform the numerical simulations, and qualitative agreement with Dorie et al.’s experimental results sought. Repeated simulations, obtained using a range of parameter values, suggest that the qualitative form of the solutions depends crucially on whether the labelled particles are inert microspheres or live tumour cells. Typical simulations demonstrating the different types of behaviour observed are presented in Figs 2 and 3. These were obtained by assuming that initially the labelled cells had penetrated a finite depth into the tumour. Thus in (17) we fix   2  2  x x  3 1− 1− x ∈ [R L (0), RT (0)], m(x, 0) = m in (x) = R L (0) RT (0)  0 otherwise. R1 In this expression 3 is chosen so that 0 m(x, 0)d x = 1 and R L (0) = 0.7R(0) and RT (0) = 0.9R(0) denote the positions of the leading and trailing edges of the labelled cells at t = 0. Figure 2 shows how the nutrient concentration, cell velocity and labelled cell density evolve when inert polystyrene microspheres are labelled. From Fig. 2(a)

608

K. E. Thompson & H. M. Byrne (a) 1

1 t=1

c(x, 1)

c(x, 0)

t=0

0.8 0.6 0.4

0.8 0.6

0

0.2 0.4 0.6 0.8 Scaled distance, x/R(t)

1

0.4

1

0

1

0.6 0.4 0

0.2 0.4 0.6 0.8 Scaled distance, x/R(t)

1

1

0.2

1

t=5

c(x, 5)

c(x, 4)

1

0.6 0.4 0

1

0.8 0.6 0

0.2 0.4 0.6 0.8 Scaled distance, x/R(t)

0.8

t=4

0.4

1

t=3

0.8

c(x, 3)

c(x, 2)

t=2

0.2 0.4 0.6 0.8 Scaled distance, x/R(t)

0.2

0.4

0.6

0.8

Scaled distance, x/R(t)

1

0.8 0.6 0.4

0

0.4

0.6

0.8

Scaled distance, x/R(t)

Figure 2. Numerical results describing the internalization of labelled polystyrene microspheres within a growing tumour. Evolution of (a) nutrient concentration, (b) cell velocity and (c) labelled microspheres at times t = 0, 1, 2, 3, 4, 5 (in dimensionless units) after internalization. From the simulations we observe that the nutrient concentration and cell velocity rapidly approach steady state profiles whereas the microspheres gradually accumulate at the centre of the tumour. Parameter values: s = 10, c∞ = 1.0, d = 6.0, cnec = 0.1, 0 = 0.5, µn = 0.01, µm = 0.008, R(0) = 1.4.

and (b) we note that c(x, t) and v(x, t) rapidly adopt time-independent forms, with the nutrient concentration increasing monotonically with x and the velocity field having an internal minimum near x = 0.6R(t). We note also that as v(x, t) ≤ 0 convective transport will drive cells towards the centre of the tumour. Figure 2(c) shows how the microspheres are internalized within the tumour. Initially they spread out (t = 1) and then they are transported by the velocity field, in a wavelike manner, towards the centre of the tumour (t = 2) where they eventually aggregate (t = 3, 4, 5). We remark that the numerical results presented in Fig. 2 show how the dependent variables c, v and m vary over time and the scaled distance 0 ≤ ρ = x/R(t) ≤ 1. Replotting these solutions in terms of 0 ≤ x ≤ R(t) does not alter the qualitative form of the solutions: it simply alters the width of the region containing the microspheres. Further, as we show below (see Fig. 4), the tumour radius R rapidly attains a time-independent equilibrium value. For all subsequent times the spatial

Cell Internalization in Tumour Spheroids

609

(b) 1

0.5 t=1

v(x, 1)

v(x, 0)

t=0

0.5 0 –0.5

0

0.2 0.4 0.6 0.8 Scaled distance, x/R(t)

–0.5

0.2 0.4 0.6 0.8 Scaled distance, x/R(t)

1

0.2 0.4 0.6 0.8 Scaled distance, x/R(t)

1

0.2

1

t=3

0

–1

–0.5 –1

0

0.2 0.4 0.6 0.8 Scaled distance, x/R(t)

0.5

–1.5

1

t=5

v(x, 5)

0

–0.5 –1 0

0

0.5

t=4

0 v(x, 4)

0

0.5 v(x, 3)

v(x, 2)

–1.5

1

t=2

0

–1.5

–0.5 –1

0.5

–1.5

0

0.2

0.4

0.6

0.8

1

–0.5 –1 –1.5

0

Scaled distance, x/R(t)

0.4

0.6

0.8

Scaled distance, x/R(t)

Figure 2. Continued

transformation is the same. For these reasons, the numerical results were presented in terms of the scaled variables. Figure 3 shows how the migration of the labelled particles changes when live tumour cells are used in place of the microspheres. The evolution of the nutrient concentration and cell velocity are qualitatively similar to those presented in Fig. 2(a) and (b) and therefore are not included. Comparing Figs 2(c) and 3(a), we note that initially (0 < t < 2) the speed of migration of the wave of invading cells is the same for the microspheres and tumour cells. Once they have penetrated to the centre of the tumour, unlike the microspheres, the labelled tumour cells redistribute themselves evenly throughout the tumour. Recalling that the motivation for Dorie et al.’s experiments was to relate the migration of labelled cells to kinetic parameters associated with the unlabelled tumour cells, a series of numerical experiments were performed in which the proliferation rate constant s was varied. For each choice of s the loci of the outer tumour boundary R(t) and the leading and trailing edges of the wave of migrating cells [R L (t) and RT (t)] were determined. Results that pertain to labelled microspheres are presented in Fig. 4. As expected, the limiting size of the tumour radius increases with s. Further, as s increases so does the proliferation-driven velocity and, hence, the time

610

K. E. Thompson & H. M. Byrne

(c) 0.2

0.2 t=1

0.15

m(x, 1)

m(x, 0)

t=0

0.1 0.05 0

0

0.2 0.4 0.6 0.8 Scaled distance, x/R(t)

0.15

0

0.2

t=2

m(x, 3)

m(x, 2)

0

1

0.1 0.05

0.2 0.4 0.6 0.8 Scaled distance, x/R(t)

1

0.2 0.4 0.6 0.8 Scaled distance, x/R(t)

1

0.2

1

t=3

0.15 0.1 0.05

00

0.2 0.4 0.6 0.8 Scaled distance, x/R(t)

0.2

00

1

0.2

t=4

0.1 0.05 0

t=5

0.15 m(x, 5)

0.15 m(x, 4)

0.1 0.05

0.2

0

0.15

0.2

0.4

0.6

0.8

1

0.1 0.05 0

0

Scaled distance, x/R(t)

0.4

0.6

0.8

Scaled distance, x/R(t)

Figure 2. Continued

taken for the labelled cells to reach the centre of the tumour decreases (from t ∼ 2.0 dimensionless units when s = 10 to t ∼ 1.5 when s = 14). As, in each case, the tumours are growing when the labelled cells are introduced (at t = 0), there is an initial transient period, whose duration increases with s, during which both R L (t) and RT (t) move outwards, towards the tumour boundary and away from the centre of the tumour. In all cases there is, however, a significant period (0.5 < t < 1.0) during which R L (t) and RT (t) both migrate towards the origin with approximately the same, constant speed. In order to determine whether these migration speeds can provide information about the proliferation rates of the unlabelled tumour cells, we estimated, for a range of values of s, the average speeds of R L (t) and RT (t) over the period 0.5 < t < 1.0. The results, which are presented in Fig. 5, show that, for all values of the proliferation rate s, the speed with which the leading edge of the front of labelled cells moves is greater than that of the trailing edge. This is consistent with the numerical results presented in Fig. 2(c) where the width of the pulse of labelled cells initially increases slightly. For small values of s (6 < s < 8) we notice that the migration speeds decrease to a minimum value near s = 8. The reduction in speed as s increases towards s ∼ 8 is due to the initial value of the tumour radius, R(0). When

Cell Internalization in Tumour Spheroids 0.1

0.1 t=1

m(x, 1)

m(x, 0)

t=0

0.05

0

0

0.2 0.4 0.6 0.8 Scaled distance, x/R(t)

0.1

0.05

0

1

0

0.1

0.05

0

0.2 0.4 0.6 0.8 Scaled distance, x/R(t)

0.1

0

1

0

0.1 m(x, 5)

m(x, 4)

0.2 0.4 0.6 0.8 Scaled distance, x/R(t)

1

0.2

1

t=5

0.05

0

1

0.05

t=4

0

0.2 0.4 0.6 0.8 Scaled distance, x/R(t) t=3

m(x, 3)

m(x, 2)

t=2

0

611

0.2

0.4

0.6

0.8

Scaled distance, x/R(t)

1

0.05

0

0

0.4

0.6

0.8

Scaled distance, x/R(t)

Figure 3. Numerical results describing the internalization of labelled tumour cells within a growing tumour. The evolution of the nutrient concentration, the cell velocity and the tumour radius (not shown) are qualitatively similar to those presented in Fig. 2. The labelled tumour cells gradually distribute themselves uniformly throughout the tumour volume. Parameter values: s = 10, c∞ = 1.0, d = 6.0, cnec = 0.1, 0 = 0.5, µn = 0.01, µm = 0.008, R(0) = 1.4.

s < 8 the steady state tumour radius is smaller than the initial tumour radius and, hence, the tumour initially shrinks. As s increases the steady state radius increases in magnitude, exceeding R(0) when s > 8. Provided that s is sufficiently large (s > 8 in Fig. 5) the average speeds of R L (t) and RT (t) (and, hence, also their average) increase approximately linearly with s. Hence we deduce that measurements of the speeds of migration of the labelled cells can provide information about the growth kinetics of the underlying multicell spheroid, provided that the tumour increases in size during the experiment. In practice, by comparing the speeds of migration of microspheres through multicell spheroids from different cell lines, all subject to the same external nutrient concentration, it should be possible to predict which cell line has the most rapid turnover rate. Similar simulations, performed with labelled tumour cells, yielded profiles for R L (t) which were indistinguishable from those obtained using the microspheres. However, as we observe by comparing Figs 2(c) and 3(a), because the labelled tumour cells become rapidly distributed throughout the tumour, RT (t) ∼ R(t).

612

K. E. Thompson & H. M. Byrne 2

2

1.5

1.5

R(t) R(t)

RT (t)

RT (t)

1

1

RL(t)

0.5

0.5

0

0

RL(t)

(s = 10)

(s = 8)

0.5

1

1.5

2

0

0

0.5

Time, t 2

2

RL(t)

(s =14)

0.5

2

0.5

RL(t) (s =12)

0

1.5

RT (t)

1

1

0

2

1.5

RT (t)

0.5

1.5

R(t)

R(t)

1.5

1 Time, t

1 (Time, t )

1.5

2

0

0

0.5

1 (Time, t )

Figure 4. Sequence of figures showing how the evolution of the outer tumour boundary, and the leading and trailing edges of the invading front of microspheres depend upon the proliferation rate, s = 8.0, 10.0, 12.0, 14.0. As s increases the limiting value of R(t) increases and the time taken for the microspheres to reach the centre of the tumour decreases. The qualitative form of the profiles for the leading and trailing edges of the invading front also changes with s. Parameter values: c∞ = 1.0, d = 6.0, cnec = 0.1, 0 = 0.5, µn = 0.01, µm = 0.008, R(0) = 1.4.

Thus experiments involving labelled tumour cells will yield only estimates of the average speed of the leading edge R L (t). Consequently, the resulting estimates of the spheroid’s proliferation rate will be less robust than estimates obtained from experiments involving microspheres. In this respect, experiments involving labelled microspheres may be preferred to those involving labelled tumour cells.

4.

A NALYTICAL R ESULTS

In this section attention focuses on the initial wavelike migration of the labelled cells and their limiting, or long-time, behaviour. Guided by the numerical results of Section 3, and in order to simplify the analysis, we will assume that the nutrient concentration, cell velocity and tumour radius have evolved to time-independent, steady-state profiles. In addition, and to ensure that the analysis is tractable, we focus attention on a special case for which the total number of labelled particles is

Cell Internalization in Tumour Spheroids

613

1.4

Estimate of migration speed

1.2 VL

1 0.8 0.6

VAV

0.4 VT

0.2 0

6

7

9

8

10

11

12

13

14

15

Proliferation rate, s Figure 5. Figures showing how, over the period 0.5 < t < 1.0 the average migration speeds of the leading and trailing edges of the pulse of labelled microspheres (VL and VT , respectively) depend upon the proliferation rate, s. Also shown is the average speed. Vav = (VL + VT )/2. The results show that the various speeds (VL , VT and Vav ) are all approximately linearly dependent on the external nutrient concentration. Parameter values: c∞ = 1.0, d = 6.0, cnec = 0.1, 0 = 0.5, µn = 0.01, µm = 0.008, R(0) = 1.4

small and may be characterized by the small parameter 0 <   1. Specifically, we seek solutions to equations (9)–(17) of the form: c = c0 (x) + O(), R = R0 + O(),

v = v0 (x) + O(), m = m 0 (x, t) + O( 2 ).

(18)

[Note: the assumption of incompressibility enables us to establish that n = 1 − m 0 (x).] By equating to zero coefficients of O(1) in equations (9)–(11) we obtain equations for c0 , v0 and R0 . As, by assumption, m ∼ O(), these equations are independent of m 0 and, hence, their form will be unaffected by whether microspheres or tumour cells are labelled. The simplicity of the kinetic terms employed in the model ensures that equations (9) and (10) are integrable at leading order. In this way we deduce that c0 and v0 are defined in terms of the equilibrium tumour radius R0 as follows: 0 c0 (x) = c∞ − (R02 − x 2 ) 2



and

0 d v0 (x) = sx c∞ − − s 2

 R02

x2 − 3

 . (19)

614

K. E. Thompson & H. M. Byrne

Using (11), with

d R0 dt

= 0, we deduce that  0 = R0

 d 0 R02 , c∞ − − s 3

(20)

and hence that for a nontrivial solution s R0 =

  d 3 c∞ − . 0 s

Using this expression it is clear that as s increases (or d decreases) the limiting, steady-state tumour radius increases—this agrees with the numerical results presented in Fig. 4. Now, for non-necrotic solutions we require that min0≤x≤R0 c0 (x) = c0 (0) ≥ cnec . Following Byrne and Chaplain (1998) and substituting with our expression for R0 , we deduce, from (19), that a necessary condition for obtaining non-necrotic solutions is c∞ − cnec

  3 d c∞ − . ≥ 2 s

With R0 determined by (20), the expression for v0 can be rewritten v0 (x) = −

s0x 2 (R0 − x 2 ). 6

Using this result it is possible to show that v0 (x) has a minimum turning point √ at x = R0 / 3 ∼ 0.58R0 . This is in good agreement with the numerical results presented in Fig. 2(b). As equation (12) is trivially satisfied at O(1) when m = m 0 (x), we must equate to zero coefficients of O() in order to determine m 0 (x). In this way we deduce that ∂m 0 ∂m 0 ∂ 2m0 + S0 (c0 , m 0 ), + v0 = µm ∂t ∂x ∂x2

(21)

where the form of the source term S0 (c0 , m 0 ) depends on the type of particles that are labelled:  −(sc − d)m = dv0 m 0 0 0 S0 (c0 , m 0 ) = dx 0

labelled microspheres, labelled tumour cells.

(22)

Cell Internalization in Tumour Spheroids

615

4.1. Wavelike migration. Now, regardless of whether the labelled particles are microspheres or tumour cells, they are always much larger than nutrient molecules. Hence their random motility coefficients will be much smaller than for the nutrient, that is in (21) 0 < µm  1. When considering the initial, wavelike internalization of the labelled particles we exploit this fact by neglecting random motion in equation (21). Then, to leading order, the initial migration of the labelled cells can be described approximately by the following nonlinear wave equation: ∂m 0 ∂m 0 + v0 = S0 (c0 , m 0 ). ∂t ∂x

(23)

From (23) it is clear that the propagation of the labelled cells is driven by the velocity field. This, in turn, is created by differential cell proliferation and death rates. Now, for the asymptotic limit under consideration, v0 is independent of the choice of cell labelling. Hence we deduce that the speed of migration of the labelled cells is also independent of whether microspheres or tumour cells are used. This prediction is in good agreement with the numerical results of Section 3 (see Fig. 4). Given that (23) is independent of both µm and µn , the random motility coefficients for the labelled and unlabelled cells respectively, we deduce that the initial migration of the labelled cells will be unaffected by the size of the labelled particles. This is consistent with Dorie et al.’s experimental results (1982, 1986). It is possible to construct explicit solutions to equation (23) using the method of characteristics (Williams, 1980). The details of these solutions depend on the choice of S0 . It is appropriate to construct these solutions by introducing new independent variables p, q to parametrize the characteristic curves and the initial data respectively. With m 0 (x, t) = m char ( p, q), when the microspheres are labelled 0 equation (23) can be written in terms of these variables as dt = 1, dp

s0x 2 dx = v0 = − (R0 − x 2 ), dp 6

dm char s0 2 0 , = (R − 3x 2 )m char 0 dp 6 0

with the initial data now prescribed at p = 0 and parametrized by q in the following way: t ( p = 0, q) = 0,

x( p = 0, q) = q,

m char ( p = 0, q) = m in (q). 0

Integrating along the characteristic curves subject to the initial data, we obtain the following parametrized representation of the solution:  t ( p, q) = p,  2 2 −s0 R02 p/3 2 2 −s0 R02 p/3 2 2 x ( p, q) = q e /{(1 − q /R0 ) + q /R0 e }, 2 2 −s0 R02 p/3 3/2  s0 R02 p/3 2 2 m char ( p, q) = m (q)e {(1 − q /R ) + q /R e } . in 0 0 0

(24)

616

K. E. Thompson & H. M. Byrne

Eliminating p, q from (24) our solution can be rewritten in terms of the original independent variables p 2 m 0 (x, t) = m in (x/ Q) e−s0 R0 t/3 Q −3/2 ,

(25)

where 

 x2 x2 2 Q(x, t) = 1 − 2 e−s0 R0 t/3 + 2 . R0 R0 In a similar manner it is possible to show that when the tumour cells are labelled p m 0 (x, t) = m in (x/ Q).

(26)

In Fig. 6 we use (25) and (26) to present sketches of these approximate solutions. From Fig. 6(a) we observe how the pulse of labelled cells decreases in width as it travels towards the centre of the tumour. In order to preserve the total number of labelled cells the height of the pulse increases. Eventually the height of this peak becomes so large that our asymptotic expansion ceases to be valid. Given that the initial data has compact support (m 0 (x, 0) > 0 ∀ x ∈ (x L , x T )), an alternative way to see this is by following the characteristics through the endpoints of the domain 2 of compact support. From (24) we deduce that x L (t), x T (t) ∼ O(e−s0 R0 t/6 ) as t → ∞, i.e., the width of the domain of compact support shrinks exponentially. At the same time the height of the peak must increase in order to preserve the total 2 number of microspheres within the tumour: m 0 ∼ O(es0 R0 t/6 ). Eventually the spatial gradients become so large that our approximate solution ceases to be valid because second order effects such as random motility can no longer be neglected. Figure 6(b) depicts the corresponding situation when the tumour cells are labelled. As in Fig. 6(a) the width of the pulse of labelled cells decreases as it migrates towards the centre of the tumour. However, as the total number of cells is no longer preserved, our approximation predicts that the total number of cells actually falls. As in Fig. 6(a), the approximation breaks down when the width of the peak becomes so small that terms involving second spatial derivatives, such as random motility, can no longer be neglected. 4.2. Limiting distributions of labelled cells. In the previous subsection approximate solutions were constructed to describe the initial migration of the labelled cells within the tumour and the limitations of these approximate solutions were discussed briefly. From the numerical simulations of Section 3 (see Figs 2 and 3) it is clear that for both labelling systems the particles eventually adopt a well-defined structure. In this subsection we focus on the behaviour of equation (21) in the limit as t → ∞ for which ∂t∂ = 0. As our analysis will show, it is important to retain the

Cell Internalization in Tumour Spheroids

617

(a) t = 2.0

Labelled tumour cells, m(x,t)

0.25

0.2

0.15 t = 1.5

0.1

0.05

t = 1.0 t = 0.0

t = 0.5

0

0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.4

1.6

Distance from centre of tumour, x (b) 0.025 t = 1.0

t = 0.5

t = 0.0

Labelled tumour cells, m(x, t)

0.02

0.015

0.01

0.005

0

0

0.2

0.4

0.6

0.8

1

1.2

Distance from centre of tumour, x

Figure 6. Sketches depicting solutions of approximate equations that describe initial migration of labelled particles within the tumour. (a) When microspheres are labelled, the pulse of labelled cells decreases in width and increases in height as it migrates towards the centre of the tumour. (b) When the tumour cells are labelled, the pulse decreases in width whilst remaining of constant height as it migrates towards the centre of the tumour. Parameter values: s = 10.0, c∞ = 1.0, c˜ = 6.0, cnec = 0.1, 0 = 0.5.

618

K. E. Thompson & H. M. Byrne

random motility term in this asymptotic limit. Thus we now consider equation (21) with ∂t∂ = 0 in order to investigate the long-time behaviour of the labelled particles: 0 = µm

d 2m0 dm 0 + S0 . − v0 2 dx dx

(27)

When the labelled cells are inert microspheres 

dv0 S0 = −(sc0 − d)m 0 = − m0 dx



equation (27) can be rewritten 0=

d dx

 µm

 dm 0 − v0 m 0 . dx

Integrating with respect to x and imposing (16) yields dm 0 v0 m 0 s0x 2 = =− (R − x 2 )m 0 . dx µm 6µm 0

(28)

As, for physically realistic solutions, m 0 ≥ 0 we deduce that at equilibrium m 0 (x) is monotonically decreasing with x. Integrating again with respect to x yields   s0 m 0 (x) = M ∗ exp − x 2 (2R02 − x 2 ) , 24µm

(29)

where the constant M ∗ is chosen so that the total number of labelled cells within the tumour remains fixed, that is: Z  0

1

Z m 0 (x)d x = 0

1

Z m(x, 0)d x =

1

m in (x)d x.

0

From (29) we deduce that at equilibrium the microspheres are concentrated at the centre of the tumour (x = 0). This agrees with Dorie et al.’s experimental results (1982, 1986), and the numerical results of Section 3 (see Fig. 2). In addition, we note that as µm → 0 this equilibrium profile becomes steeper and more localized, approximating a delta-function in the limit as µm → 0. We remark further that if random motion is neglected (µm = 0) then our model does not admit a classical steady solution. As discussed above, in this case the cells accumulate at the origin in a region (or boundary layer) whose width decreases exponentially with time. Eventually m 0 (0, t) ∼ O( −1 ) and the asymptotic expansion breaks down. In practice, and on the basis of Dorie et al.’s experimental results, it seems reasonable to assume that the cells undergo a minimal amount of random motion, and that this 0 physical effect regularizes any sharp increases in ∂m within such a boundary layer. ∂x

Cell Internalization in Tumour Spheroids

619

In (29), the multicell spheroid’s growth parameters are combined in a single parameter grouping, s0/µm . Consequently, fitting (18) to the experimentally observed limiting distribution of labelled microspheres will yield an estimate of s0/µm . Now, the numerical simulations and the analytical solutions presented above indicate that the average migration speed of the microspheres during the early stages of the experiment yield estimates of the parameter grouping s0 which is independent of µm . Combining these two results should enable us to derive estimates of s0 and µm . When the labelled cells are tumour cells (S0 = 0) equation (27) can be written: 0 = µm

 Z d 2m0 dm 0 dm 0 ∗ − v exp − ⇒ = C 0 dx2 dx dx

x

 v0 (x)d ¯ x¯ ,

where the constant of integration C ∗ is determined by imposing (15). This yields C ∗ = 0 ⇒ m 0 = constant. Hence we deduce that at equilibrium the labelled tumour cells are uniformly distributed throughout the tumour volume. This agrees with Dorie et al.’s experimental results (1982, 1986), and the numerical simulations of Section 3 (see Fig. 3). Given that the limiting distribution of the labelled tumour cells is spatially uniform it yields no additional information about the growth kinetics of the underlying multicell spheroid. This contrasts with the situation when the microspheres are labelled.

5.

C ONCLUSIONS

In this paper we have developed a mathematical model that describes the evolution of an avascular tumour, or multicell spheroid, that contains a mixture of labelled and unlabelled cells, growing in response to an externally-supplied nutrient. As such, the model builds upon existing work (Greenspan, 1972, 1976; Adam, 1987; Adam and Maggelakis, 1990; Byrne and Chaplain, 1998; Ward and King, 1997). The main aim of the work was to provide a simple mathematical description of experimental results obtained by Dorie and co-workers (1982, 1986). In a series of experiments they showed that the patterns of internalization of labelled cells within a growing multicell spheroid depends upon the type of particles that are labelled. In particular, when inert microspheres are labelled they accumulate at the centre of the tumour, whereas labelled tumour cells distribute themselves uniformly throughout the tumour (see Fig. 1). McElwain and Pettet (1993) developed a mathematical model to explain Dorie et al.’s results. Their model provided a good description of the internalization of the inert microspheres. In order to explain the internalization pattern of the tumour cells, McElwain and Pettet invoked chemotaxis, postulating that the tumour cells migrate preferentially up nutrient concentration gradients. In order to validate

620

K. E. Thompson & H. M. Byrne

their assertion, one could use a Boyden chamber assay to track the paths of tumour cells through gels pretreated with spatially-varying nutrient concentrations (Boyden, 1962; Zigmond, 1978). Whilst such experiments may prove McElwain and Pettet correct and show that tumour cells migrate preferentially towards sources of nutrients such as oxygen, to our knowledge no such experiments have yet been published. There is, however, experimental evidence showing that certain cell lines migrate towards certain chemicals, such as thrombospondin (Taraboletti et al., 1987). For this reason, chemotaxis was not included in our mathematical model. Rather, we distinguished between the inert microspheres and the tumour cells by choosing the cell birth and death rates appropriately. (McElwain and Pettet neglected cell birth and death for both types of labelled cells.) Numerical simulations of the full model show that this distinction can explain the different internalization patterns observed by Dorie et al. (1982, 1986) (see Figs 2 and 3). Indeed, comparing Figs 1–3 we note that our model yields results which are qualitatively similar to those obtained by Dorie et al. (1982, 1986). In order to explain why the migration patterns and limiting distributions of the microspheres and tumour cells differ, a simplified caricature model was developed. Guided by numerical simulations, which showed the nutrient concentration, the velocity field and the tumour radius rapidly relaxing to steady, time-independent profiles (see Figs 2 and 4), we assumed steady profiles for c, v and R. We also assumed that the total number of labelled cells within the tumour was small. The resulting caricature model defined the distribution of the labelled cells within the tumour in terms of known quantities. Approximate solutions of the caricature model indicated that the speed with which the labelled cells migrated was independent of the labelling system. This agrees with Dorie et al.’s findings that the way in which inert particles moved towards the centre of the spheroid is largely independent of probe size (1982, 1986). Combining our analytical and numerical results, we conclude that it is possible to extract information about the growth parameters of the host multicell spheroid by following the migration of labelled cells through it. In particular, as we show in Fig. 5, the average migration speed of the labelled cells is proportional to the proliferation rate. By constructing time-independent solutions of the caricature model we were also able to show how the limiting distribution of the labelled cells depends upon the labelling system. In particular, we predict that the microspheres accumulate at the centre of the tumour and that the tumour cells distribute themselves evenly throughout the tumour. As they are inert, the microspheres are simply convected by the proliferation-driven velocity field towards the centre of the tumour where, being unaffected by the local nutrient concentration, they accumulate. The labelled tumour cells are also transported by convection towards the centre of the tumour. However as they migrate the labelled tumour cells undergo proliferation and apoptosis at rates which depend on the local nutrient concentration. Consequently they leave a trail of labelled daughter cells behind them as they migrate away from the nutrient-rich outer tumour boundary.

Cell Internalization in Tumour Spheroids

621

For completeness, our key model predictions are summarized below: (1) Differences in the rates of cell proliferation and cell death between the microspheres and tumour cells are sufficient to explain the different internalization patterns observed by Dorie et al. (see Figs 2 and 3). (2) The average speed of migration of the invading pulse of labelled cells can provide information about the kinetic parameters of the underlying multicell spheroid. (3) The microspheres migrate as a pulse with marked leading and trailing edges, whereas the tumour cells yield only a marked leading edge. Consequently, parameter estimates relating to the multicell spheroid’s kinetic parameters obtained by labelling microspheres should be more robust than estimates obtained by labelling tumour cells. (4) At equilibrium the microspheres accumulate at the centre of the tumour, whereas the tumour cells are uniformly distributed throughout the tumour. (5) Being spatially uniform, the limiting distribution of the labelled tumour cells provides no information about the growth kinetics of the multicell spheroid. By contrast, the limiting distribution of the labelled microspheres can be used to estimate parameters associated with the tumour’s growth kinetics. Whilst both labelling systems can provide information about the underlying multicell spheroid, the above results suggest that experiments involving microspheres may provide more useful information and, hence, may be preferable to those involving labelled tumour cells. There are many ways in which our mathematical model could be extended or modified. Natural extensions which would bring the model more in line with Dorie et al.’s experiments include the following: reformulating the model to describe the growth of a three-dimensional, radially-symmetric multicell spheroid; considering more developed tumour structures, involving necrotic and hypoxic regions; studying how the internalization patterns are affected when the unlabelled cells are irradiated or pretreated with a putative anticancer drug. It would also be interesting to see how our model predictions are affected by incorporating effects to describe heterotypic interactions between the different types of cells. More generally, by reinterpreting the definitions of the labelled and unlabelled cells, it should be possible to use the model to study the evolving spatial structure of avascular tumours that comprise two or more subpopulations. The different subpopulations may refer to cells that possess normal and mutant copies of a particular gene, such as p53 (Royds et al., 1998). Given that almost 50% of solid tumours in vivo possess tumour cells whose p53 status is abnormal, by extending our model to describe such situations we hope to be able to provide clinicians with useful insight into some of the complex mechanisms underpinning cancer growth.

622

K. E. Thompson & H. M. Byrne

A CKNOWLEDGEMENTS The authors are grateful to Drs Claire Lewis and Janice Royds at the University of Sheffield for helpful discussions regarding this work. KET acknowledges financial support from the EPSRC (Grant Number: GR/L41967).

R EFERENCES

} Adam, J. A. (1987). A mathematical model of tumour growth. II Effects of geometry and spatial uniformity on stability. Math. Biosci. 86, 183–211. } Adam, J. A. and S. A. Maggelakis (1990). Diffusion regulated growth characteristics of a spherical prevascular carcinoma. Bull. Math. Biol. 52, 549–582. } Boyden, S. V. (1962). The chemotactic effect of mixtures of antibody and antigen on polymorphonuclear leukocytes. J. Exp. Med. 115, 453–466. } Byrne, H. M. and M. A. J. Chaplain (1997). The importance of inter-cellular adhesion in the development of carcinomas. IMA J. Math. Med. Biol. 14, 305–323. } Byrne, H. M. and M. A. J. Chaplain, (1998). Necrosis and apoptosis: distinct cell loss mechanisms in a mathematical model of avascular tumour growth. J. Theor. Med. 1, 223–236. } Cliff, W. J. (1963). Observations on healing tissue: a combined light and electron microscope investigation. Trans. R. Soc. Lond. B246, 305–325. } Cliff, W. J. (1965). Kinetics of wound healing in rabbit ear chambers: a time lapse cinemicroscopic study. Q. J. Exp. Physiol. Cog. Med. Sci. 50, 79–89. } Crank, J. (1988). Free and Moving Boundary Problems, Oxford: Oxford University Press. } Dorie, M. J., R .F. Kallman, D. F. Rapacchietta, D. Van Antwerp and Y. R. Huang (1982). Migration and internalization of cells and polystyrene microspheres in tumour cell spheroids. Exp. Cell. Res. 141, 201–209. } Dorie, M. J., R. F. Kallman and M. A. Coyne (1986). Effect of Cytochalasin B Nocodazole on migration and internalization of cells and microspheres in tumour cells. Exp. Cell. Res. 166, 370–378. } Fan, T. D. and P. J. Polverini (1997). In vivo models of angiogenesis, in Tumour Angiogenesis, R. Bicknell, C. E. Lewis and N. Ferrara (Eds), Oxford: Oxford University Press. } Folkman, J. and M. Hochberg (1973). Self-regulation of growth in three dimensions. J. Exp. Med. 138, 745–753. } Graeber, T. G., C. Osmanian, D. E. Hausman, C. J. Koch. S. W. Lowe and A. J. Giaccia (1996). Hypoxia-mediated selection of cells with diminished apoptotic potential in solid tumours. Nature 379, 88–91. } Greenspan, H. P. (1972). Models for the growth of a solid tumour by diffusion. Stud. Appl. Math. 52, 317–340. } Greenspan, H. P. (1976). On the growth and stability of cell cultures and solid tumours. J. Theor. Biol. 56, 229–242. } Kerr, J. F. R., A. H. Wyllie and A. R. Currie (1972). Apoptosis: a basic biological phenomenon with wide-ranging implications in tissue kinetics. Br. J. Cancer 26, 239–257. } McElwain, D. L. S. and G. J. Pettet (1993). Cell migration in multicell spheroids: swimming against the tide. Bull. Math. Biol. 55, 655–674.

Cell Internalization in Tumour Spheroids

623

} Royds, J. A., S. K. Dower, E. E. Qwarnstrom and C. E. Lewis (1998). Responses of tumour cells to hypoxia: role of p53 and NFkB. J. Clin. Pathol.: Mol. Pathol. 51, 55–61. } Taraboletti, G., D. D. Roberts and L. A. Liotta (1987). Thrombospondin-induced cell migration: haptotaxis and chemotaxis are mediated by different molecular domains. J. Cell. Biol. 105, 2409–2415. } Ward, J. P. and J. R. King (1997). Mathematical modelling of avascular-tumour growth. IMA J. Math. Appl. Med. Biol. 14, 39–69. } Williams, W. (1980). Partial Differential Equations, Oxford: Oxford University Press. } Zigmond, S. H. (1978). Chemotaxis by polymorphonuclear leukocytes. J. Cell. Biol. 77, 269–289. Received 31 August 1998 and accepted 11 November 1998