The Migration of Cells in Multicell Tumor Spheroids

Bulletin of Mathematical Biology (2001) 63, 231–257 doi:10.1006/bulm.2000.0217 Available online at http://www.idealibrary.com on

The Migration of Cells in Multicell Tumor Spheroids G. J. PETTET CiSSaIM, School of Mathematical Sciences, Queensland University of Technology, GPO Box 2434, Brisbane, QLD 4001, Australia C. P. PLEASE, M. J. TINDALL Faculty of Mathematical Studies, The University of Southampton, Southampton SO17 1BJ, U.K. D. L. S. McELWAIN CiSSaIM, School of Mathematical Sciences, Queensland University of Technology, GPO Box 2434, Brisbane, QLD 4001, Australia A mathematical model is proposed to explain the observed internalization of microspheres and 3 H-thymidine labelled cells in steady-state multicellular spheroids. The model uses the conventional ideas of nutrient diffusion and consumption by the cells. In addition, a very simple model of the progress of the cells through the cell cycle is considered. Cells are divided into two classes, those proliferating (being in G 1 , S, G 2 or M phases) and those that are quiescent (being in G 0 ). Furthermore, the two categories are presumed to have different chemotactic responses to the nutrient gradient. The model accounts for the spatial and temporal variations in the cell categories together with mitosis, conversion between categories and cell death. Numerical solutions demonstrate that the model predicts the behavior similar to existing models but has some novel effects. It allows for spheroids to approach a steady-state size in a non-monotonic manner, it predicts self-sorting of the cell classes to produce a thin layer of rapidly proliferating cells near the outer surface and significant numbers of cells within the spheroid stalled in a proliferating state. The model predicts that overall tumor growth is not only determined by proliferation rates but also by the ability of cells to convert readily between the classes. Moreover, the steady-state structure of the spheroid indicates that if the outer layers are removed then the tumor grows quickly by recruiting cells stalled in a proliferating state. Questions are raised about the chemotactic response of cells in differing phases and to the dependency of cell cycle rates to nutrient levels. c 2001 Society for Mathematical Biology

0092-8240/01/020231 + 27

$35.00/0

c 2001 Society for Mathematical Biology

232

G. J. Pettet et al.

1.

I NTRODUCTION

Multicell spheroids have been studied for a number of years as in vitro models of in vivo tumors. Many of their characteristics are well documented (Sutherland, 1988; Knuechel and Sutherland, 1990; Mueller-Klieser, 1997; Kunz-Schugart et al., 1998), and in the past few decades they have become the focus of applied mathematicians who seek to model the growth dynamics of avascular tumors. Cultured from mono-layer cells, the growth of a multicellular spheroid consists of three distinct stages. The first is the proliferation of cancerous cells, which are assisted experimentally into forming spherical aggregations. As the multicell spheroid grows, those cells some distance from the periphery of the spheroid are observed to enter a quiescent phase where they cease to proliferate. Proliferation continues near the outer surface driving the growth of the spheroid. Finally, a region of dead cell matter, the necrotic core, develops at the center of the spheroid. The necrotic core acts as a cell volume loss mechanism such that eventually, the cell loss at the center balances the cell proliferation at the periphery. When such a balance occurs the spheroid obtains a steady state in size. Thus a fully developed tumor spheroid can consist of three distinct regions; a central necrotic core, an inner shell of quiescent cells and an outer rim of proliferating, viable cells. Initial work by Burton (1966) and Greenspan (1972, 1976) focused on reaction–diffusion based models to describe the growth of multicellular spheroids by modeling nutrient supply and cell volume growth and loss. These models involve postulating the existence of a pressure gradient within the tumor mass generated by proliferation and surface tension on the boundary. Cells then exhibit a passive motion driven by the pressure differential that is related to both their position relative to the periphery of the spheroid and the spheroid’s temporal proximity to its steady state. The models successfully reproduce the layered structure observed in multicell spheroids, and have recently been extended and applied to various avascular and vascularized stages of tumor development: [see, for example, Chaplain and Britton (1993), McElwain and Pettet (1993), Byrne and Chaplain (1995), Chaplain (1996), Byrne and Gourley (1997), Please et al. (1998, 1999), Ward and King (1998a,b) and Landman and Please (to appear)]. This recent proliferation of models has contributed to our understanding of certain details of multicell spheroid development. The ability of these models to reproduce the gross structure of multicell spheroids depends primarily on the underlying ideas of Burton (1966) and Greenspan (1972, 1976). There is, however, little experimental data to verify the predicted velocities derived from these models. However, the work of Dorie et al. (1982, 1986) describes experiments to measure differential migration and internalization of labelled cells and inert polystyrene beads. Dorie et al. undertook two separate experiments. One involved monitoring the internalization of polystyrene microspheres adhered to the exterior of EMT6 cell line spheroids. In a second experiment, using similar multicell spheroids, the internalization of cells labelled with 3 H-thymidine was measured. As can be seen in

Migration of Cells

233

Figure 1. Representations of data originally presented by Dorie et al. (1982, 1986) showing the internalization of inert microspheres compared with the internalization of 3 Hthymidine labelled cells. Reproduced with the kind permission of Academic Press Inc.

Fig. 1 two very contrasting results were observed. The essential observation to be made is that the inert polystyrene microspheres are driven towards the center of the spheroid while the labelled cells, although exhibiting a similar tendency for internalization, persisted near the periphery of the spheroid occupied by proliferating cells. Consequently a distinctly bimodal distribution of labelled cells can be seen. These experiments are taken by us to imply that there are both passive and active mechanisms governing cell migration. It is the intention of the authors of this paper to explain these observations. The explanations will highlight the underlying significant mechanisms of internalization and hence will show that interpretation of behavior based on previous models of avascular tumor growth may prove inadequate. McElwain and Pettet (1993) attempted an explanation of this observed behavior in a model incorporating pressure gradients generated by cell proliferation and the chemotactic response to the established gradient in nutrient concentrations as the main mechanisms for providing cell movement. Although the numerical solutions presented ostensibly matched the experimental results of Dorie et al. (1982, 1986), one shortcoming of their model is that the labelled cells do not proliferate. Furthermore, despite the comparable nature of labelled and unlabelled cells, their model does not include the chemotactic response of the unlabelled cells, reserving this capability for the labelled cells alone. More recently, Thompson and Byrne (1999) presented a model including mitosis of labelled cells, but excluding chemotaxis. This model reproduces some of the observed internalization behavior but fails to predict the bimodal distribution of labelled cells as illustrated in Fig. 1. Furthermore, in order to achieve an internal-

234

G. J. Pettet et al.

ization of labelled cells, their model requires the labelled cells to be initially located at some distance from the periphery of the spheroid in contrast to their adhesion on the spheroid surface as described by Dorie et al. In this paper we reproduce results analogous to those presented by Dorie et al. (1982, 1986) by including different chemotactic behavior between proliferating and quiescent cells. In this formulation we model the quiescent and proliferating sub-populations as different cell cycle states with conversion from one state to another dependent upon local nutrient concentration. We can interpret the experiments of Dorie et al. as providing evidence for the self-sorting behavior of two sub-populations of tumor cells, based upon their different chemotactic behavior.

2.

M ODEL F ORMULATION

2.1. Nutrient distribution. We begin by considering a growing multicell spheroid to be surrounded by nutrient-rich serum. For the purposes of modeling we will consider the effect of a single nutrient, usually taken to be oxygen. Assuming the diffusion of oxygen to be very rapid, with diffusion coefficient D, and the consumption of oxygen by living cells to be dependent on the local oxygen concentration C(x, t), the oxygen concentration maintains a quasi-static distribution determined by D∇ 2 C = N σ (C),

(1)

where σ (C) describes the rate of oxygen consumption per cell and N represents the number of living cells per unit volume of the multicell spheroid. 2.2. Cell kinetics. In this paper we make the assumption that the cells that constitute the multicell spheroid may be in one of three possible states; a proliferating state, a non-proliferating (quiescent) state and a dead state. Each of these states may be interpreted as representing components of the cell cycle. It is an intrinsic element of our model that cells will exhibit a set of kinetic and dynamic behaviors characteristic of the part of the cell cycle in which they reside. The proliferative part of the cell cycle is broken down in to a number of phases known as G 1 , S, G 2 and M phases. These individual stages play crucial roles in the growth, synthesis of cellular DNA and subsequent cell division at the end of Mphase. Cells in the proliferative part of the cell cycle grow and divide at intervals dependent upon their size, environment and internal cell cycling mechanisms. Quiescent cells are those in the G 0 phase and persist in a dormant or non-proliferative state. Cells may move to the quiescent phase from the proliferative phase due to a lack of available nutrients although growth inhibitors can also trigger such a conversion. Cells may also return to the proliferative part of the cycle either by an increase in nutrient supply or in response to external stimuli such as growth factors.

Migration of Cells

235

Figure 2. Schematic diagram of the three state model used showing the rate functions for conversion between states. P = proliferating, Q = quiescent, D = dead.

In contrast, lack of nutrient supply to cells in the quiescent state may trigger necrosis and hence cell death. Those cells that undergo apoptotic cell death at the end of G 1 and the onset of S phases may be represented in this model by the conversion of cells in the proliferative state to cells in the death state. Note that the model does not assume all cells in a hypoxic region are quiescent, but rather considers the rate at which cells make the transition between the phases to be oxygen dependent. It is worth noting that equation (1) assumes that all cells use the nutrient at the same rate. For the purposes of the work presented here we assume that the rate at which cells transfer between states is only dependent upon the local extracellular nutrient concentration. The rate functions for each of the possible transfers in our model of cell kinetics are summarized in Fig. 2. More explicitly we have; • cells in the proliferative state proliferate at a rate K B (C) and are converted to the quiescent state at the rate K Q (C) and undergo apoptotic death at a rate K A (C), and • cells in the quiescent state may revert to the proliferating state at a rate K P (C) and undergo necrotic death at a rate K D (C). By considering a given region of the spheroid, at every point the cells can exist in one of three states, namely proliferating, quiescent or dead. We assume that dead cells occupy no space in the tumor. If we further assume that all proliferating and quiescent cells are the same size, that they are incompressible and that the cells are closely packed together, then the number of cells per unit volume, N , remains constant. Note that this assumes that daughter cells quickly grow to their natural size. We introduce the notation that the number of active proliferating cells per

236

G. J. Pettet et al.

Figure 3. Profiles of the concentration within the spheroid at different time steps.

unit volume is P(x, t) while the number of quiescent non-proliferating cells per unit volume is Q(x, t). The assumption that the cells are well packed implies N = P + Q.

(2)

2.3. Cell migration. Assuming cells to be incompressible leads us to model cell movement within a region of well packed cells as being dependent upon the pressure gradients generated by cell proliferation. If cells are in close contact with each other then proliferation will generate forces that create motion in such a way as to make space for the new cells. Although this convective flux of cells is the basis of previous multicell spheroid models we wish to extend these ideas to consider a nutrient dependent chemotactic response that also contributes to the cell velocity. Experiments regarding the chemotactic response of cells as a function of the phases of the cell cycle such as those of Hughes and McCulloch (1991) and Palka et al. (1996) indicate that there is a strong dependence. However, the available evidence does not seem to identify the exact nature of this dependency. In particular, we anticipate that cells undergoing mitosis may be somewhat less motile than those in the quiescent state (Hughes and McCulloch, 1991). It can be argued that quiescent cells will detect the nutrient concentration gradient and seek to move towards areas of higher nutrient concentration. In seeking to do so they will have to compete for space with the proliferating cells as they continually divide and drive cells towards the center of the spheroid.

Migration of Cells

237

Figure 4. The predicted internalization of proliferating and quiescent cells for the case of χ = −1.0.

As stated previously, it is fundamental to our model of cell migration in the context of multicell spheroids that cells in different states exhibit different dynamics. If we denote the velocity of proliferating cells to be u P (x, t) and the velocity of the quiescent cells to be u Q (x, t) then the two velocities are related by u Q = u P + χ∇C,

(3)

where cells in the quiescent state are taken to be more chemotactically active than their proliferating counterparts, migrating in the direction of increasing nutrient concentration as quantified by the chemotactic coefficient χ. In what follows χ is taken to be constant. Note that it is possible to allow χ to depend on the local oxygen concentration C [for example as discussed in Murray (1993)], in particular to account for reduced chemotactic response in hypoxic regions. However, in hypoxic regions not only is the concentration very small but also different forms of the chemotactic coefficient have little effect upon the overall behavior of the model. 2.4. Cell conservation. Combining the above descriptions we arrive at the following conservation equations for cells in the proliferating state ∂P + ∇ · (u P P) = (K B − K Q − K A )P + K P Q, ∂t

(4)

238

G. J. Pettet et al. 5 4.5 4 3.5

^^ R(t)

3 2.5 2 1.5 1 0.5 0 0

10

20

30

40

50 t^

60

70

80

90

100

ˆ tˆ) from the initial to steady state. Figure 5. The outer cell radius R( and for cells in the quiescent state ∂Q + ∇ · (u Q Q) = K Q P − (K D + K P )Q. ∂t

(5)

Here K B is the birth rate, K Q is the rate to quiescence, K A is the apoptotic rate, K D is the death rate and K P is the rate of conversion to proliferation from quiescence. Using equations (2) and (3) and neglecting the influence of apoptotic cell death by taking K A ≡ 0 we can combine equations (4) and (5) to give ∇ · (u P N + χ (N − P)∇C) = K B P − K D (N − P)

(6)

together with   P ∂P + u P + χ ∇C · ∇ P ∂t N   P P 2 = (N − P) K P − K Q + (K B + K D + χ ∇ C) . N−P N

(7)

It can be seen that the model at this stage consists of a diffusion equation (1) for the nutrient concentration, a quasi-static equation (6) for the velocity and a nonlinear hyperbolic equation (7) for the distribution of cells in the proliferating state.

Migration of Cells

239

Figure 6. The velocity of the outer surface u Q (ˆr = R(tˆ), tˆ) from the initial to steady state.

2.5. Labelled cell migration. In introducing labelled cells to our system we assume that the labelled cells are identical to unlabelled cells except for the labelling and hence can exist in either the proliferating or quiescent states. Thus in any unit volume of the spheroid at any time t the total number of cells is now given by N = P + PL + Q + QL,

(8)

where P L (x, t) is the number of labelled proliferating cells per unit volume and Q L (x, t) is the number of labelled quiescent cells per unit volume. The conservation equations for the labelled cells then become ∂ PL + ∇ · (u P P L ) = (K B − K Q )P L + K P Q L , ∂t

(9)

∂ QL + ∇ · (u Q Q L ) = K Q P L − (K D + K P )Q L . ∂t

(10)

In order to simplify the model we make the assumption that the concentration of labelled cells at any point is small compared to N , and hence we can take the dilute limit such that PL  N, QL  N, which gives N ≈P+Q

(11)

240

G. J. Pettet et al.

and P + P L ≈ P,

Q + Q L ≈ Q.

Given the dilute limit approximation, the velocity of the proliferating cells is still determined by equation (6), so we can rewrite the conservation equation (9) for the labelled proliferating cells as PL ∂ PL + uP · ∇ P L = {(K B + K D + χ∇ 2 C)(N − P) − K Q N − χ∇ P · ∇C} ∂t N +K P Q L .

(12)

For the labelled quiescent cells we use (3) and (10) to obtain   ∂ QL P χ L L 2 + u Q · ∇ Q = −Q (K B + K D + χ∇ C) + K P + ∇ P · ∇C ∂t N N +K Q P L .

(13)

2.6. Microsphere migration. Finally we consider extending the model to account for the presence of inert microspheres. Due to their passive nature the motion can be described by ∂M + ∇ · (u M M) = 0, ∂t

(14)

where M(x, t) is the number of microspheres per unit volume and u M is their associated velocity. Because the microspheres are not actively mobile their velocity is the same as the non-chemotactic proliferating cells so that uM = u P .

(15)

As with the labelled cells, the number of microspheres per unit volume is assumed small compared to the total number of cells and, as was the case with the labelled cells, we therefore utilize the dilute limit approximation. Expanding (14) and substituting for u M using (15) we have M ∂M + u P · ∇ M = − [K B P − K D (N − P) + χ∇ P∇C − χ(N − P)∇ 2 C)]. ∂t N (16)

Migration of Cells

3.

241

O NE - DIMENSIONAL S PHERICAL M ODEL F ORMULATION

We now formulate the complete problem in the case of a spherically symmetric tumor giving the associated boundary and initial conditions. For the distribution of nutrient within the spheroid we have, assuming that the nutrient is consumed by the cells at a rate proportional to the local concentration,   1 ∂ 2 ∂C r = σ0 N C (17) D 2 r ∂r ∂r with σ0 constant and boundary conditions ∂C (r, t) = 0 ∂r C(R, t) = C0

at r = 0, at r = R(t).

Here R(t) is the outer radius of the spheroid at time t and R(0) = R0 is the initially prescribed value. The first of the boundary conditions is a result of the symmetry of the problem, and the second is a consequence of our assumption that the spheroid is supported in a medium that is well supplied with nutrient. Similarly the radial velocity, u P , of the proliferating cells is determined by   1 ∂ 2 KB P K D (N − P) χ ∂ P ∂C χ(N − P) 1 ∂ 2 ∂C (r u P ) = − + − r r 2 ∂r N N N ∂r ∂r N r 2 ∂r ∂r (18) with u P (0, t) = 0

at r = 0,

and when χ < 0, when χ > 0,

u P (R, t) = d R(t) dt u P (R, t) + χ∇C =

d R(t) dt

at r = R(t), at r = R(t).

For the first case above (χ < 0) the boundary condition corresponds to a surface dominated by proliferating cells driven by chemotaxis to the surface and the surface moving with the speed of the cells there (we assume no cells are lost into the surrounding medium). In the second case the boundary condition corresponds to the reverse scenario of the outer surface dominated by chemotactically driven quiescent cells. The case of χ = 0 is degenerate but easily dealt with. To solve for the unlabelled cells we have   ∂P χ P ∂C ∂ P + uP + ∂t N ∂r ∂r     P 1 ∂ P 2 ∂C = (N − P) K P − K Q + KB + KD + χ 2 r (19) N−P N r ∂r ∂r

242

with

G. J. Pettet et al.

when χ < 0, when χ > 0,

P(r, 0) = 1, P(r, 0) = 0,

0 ≤ r < R(t) 0 ≤ r < R(t)

for given σ0 , N , C0 , R0 , K B , K D , K Q , K P and χ . The second initial condition given here dictates that the spheroid is initially made of cells all in the quiescent state, however, other conditions have been considered. We note that this initial condition may be reasonable as, before the onset of aggregation in to a spheroid, cells are essentially non-proliferative. Once P is known we can solve for Q and its associated velocity u Q (r, t) using Q=N−P uQ = uP + χ

(20) ∂C . ∂r

(21)

The labelled cell distribution is found by solving ∂ PL ∂ PL + uP = KP QL ∂t ∂r     PL 1 ∂ ∂C ∂ P ∂C + KB + KD + χ 2 r2 (N − P) − K Q N − χ N r ∂r ∂r ∂r ∂r (22) ∂ QL ∂ QL + uQ = KQ PL ∂t ∂r     QL 1 ∂ ∂ P ∂C 2 ∂C − KB + KD + χ 2 r P + KP N + χ . N r ∂r ∂r ∂r ∂r

(23)

In this model we take the labelled cells to be applied at time t ∗ in the narrow region (1 − )R(t ∗ ) ≤ r ≤ R(t ∗ ) to represent their placement onto the surface. Hence the conditions are ( P0 (1 − )R(t ∗ ) ≤ r ≤ R(t ∗ ) L ∗ P (r, t ) = 0 0 ≤ r < (1 − )R(t ∗ ) ( Q 0 (1 − )R(t ∗ ) ≤ r ≤ R(t ∗ ) L ∗ Q (r, t ) = 0 0 ≤ r < (1 − )R(t ∗ ). Finally the microsphere distribution is determined by solving  ∂M ∂M M + uP =− K B P − K D (N − P) ∂t ∂r N   1 ∂ ∂ P ∂C 2 ∂C +χ − χ (N − P) 2 r ∂r ∂r r ∂r ∂r

(24)

Migration of Cells

243

with M(r, t ) = ∗

(

M0 0

(1 − )R(t ∗ ) ≤ r ≤ R(t ∗ ) 0 ≤ r < (1 − )R(t ∗ ).

The model for the unlabelled cells comprises a diffusion equation describing the nutrient concentration (17), a quasi-steady equation describing the velocity (18) and equation (19) for the cell phases. The cell phase can have a wave velocity given, from (19), by (u P + (χ /N )P∇C) which, due to its dependence on P, may create shocks within the system. Such shocks will mark the sharp boundary between proliferating and quiescent regions in the tumor. The wave speed for proliferating and quiescent labelled cells is u P and u P +χ ∇C, respectively, and, as these do not depend on the labelled concentrations due to our dilute approximation, the labelled cell concentrations cannot develop any shock behavior. The behavior of the microspheres is similar to the labelled cells.

4.

S OLUTION OF THE M ODEL

To demonstrate the typical behavior of the model the results presented here relate to a particularly simple model of the dependence of the rates in the problem. For the remainder of this paper we have taken K B = k B C,

K D = k D (C0 − C)

with k B , k D constant

so that the proliferation rate is maximum on the outer boundary and the death rate increases towards the center of the spheroid, and K P = kPC

K Q = k Q (C0 − C),

with k Q , k P constant

so that the cells tend to convert to a proliferating state when nutrient levels are high and become quiescent when nutrient levels are low. For this situation the model can be put into a simple non-dimensional form by taking a typical length scale associated with the nutrient diffusion and a typical time scale associated with the proliferation rate at the outer surface

L=

s

D , σ0 N

T =

1 . k B C0

The problem is then solved using variables rˆ = r

1 , L

1 Rˆ = R , L

tˆ = t

1 , T

uˆ P = u P

T , L

244

G. J. Pettet et al.

C Cˆ = , C0

1 Pˆ = P , N

1 PˆL = P L , P0

1 Qˆ L = Q L , P0

1 Mˆ = M . M0

The model in this form depends on just seven non-dimensional parameters, namely, kD =

kD , kB

kQ =

kQ , kB

kP =

kP , kB

which represent the relative rates of the various processes, and s D R0 R0 = = R0 , L σ0 N

χ σ0 N χ χ= = , 2 kB L kB D

,

τ=

t∗ = t ∗ k B C0 , T

which describe the size of the chemotactic response, the initial size of the spheroid, the thickness of the layer of labelled cells and microspheres put on the outer surface of the spheroid and the time at which these are placed there. The full problem is  ˆ 1 ∂ 2 ∂C r ˆ = Cˆ rˆ 2 ∂ rˆ ∂ rˆ

(25)

with boundary conditions ∂ Cˆ (ˆr , tˆ) = 0 ∂ rˆ

at rˆ = 0,

ˆ R, ˆ tˆ) = 1 C(

ˆ tˆ) at rˆ = R(

ˆ and initial data R(0) = R0 . 1 ∂ 2 ˆ ˆ (ˆr uˆ P ) = Cˆ Pˆ − k D (1 − C)(1 − P) rˆ 2 ∂ rˆ  ˆ ∂ Pˆ ∂ Cˆ 1 ∂ ∂ C 2 ˆ +χ − χ (1 − P) rˆ ∂ rˆ ∂ rˆ rˆ 2 ∂ rˆ ∂ rˆ

(26)

with uˆ P (0, tˆ) = 0

at rˆ = 0

and when χ < 0, when χ > 0,

ˆ tˆ) = d R(t ) uˆ P ( R, d tˆ ˆ tˆ) + χ ∂ Cˆ = uˆ P ( R, ∂ rˆ ˆ ˆ

ˆ tˆ) d R( d tˆ

ˆ tˆ), at rˆ = R( ˆ tˆ). at rˆ = R(

(27)

Migration of Cells

245

   ∂ Pˆ ∂ Cˆ ∂ Pˆ ˆ k P Cˆ − k Q (1 − C) ˆ Pˆ + + uˆ P + χ Pˆ = (1 − P) 1− Pˆ ˆ ∂ r ˆ ∂ r ˆ ∂t    ˆ 1 ∂ ∂ C 2 ˆ +χ 2 Pˆ Cˆ + k D (1 − C) rˆ ∂ rˆ rˆ ∂ rˆ

(28)

and when χ > 0 when χ < 0

ˆ r , 0) = 0, P(ˆ ˆ r , 0) = 1, P(ˆ

ˆ tˆ) 0 ≤ rˆ < R( ˆ tˆ) 0 ≤ rˆ < R(

  ˆ  ∂ PˆL ∂ PˆL ˆ + χ 1 ∂ rˆ 2 ∂ C (1 − P) ˆ + uˆ P = k P Cˆ Qˆ L + PˆL Cˆ + k D (1 − C) ∂ rˆ rˆ 2 ∂ rˆ ∂ rˆ ∂ tˆ  ∂ Pˆ ∂ Cˆ ˆ −k Q (1 − C) − χ (29) ∂ rˆ ∂ rˆ

  ˆ  1 ∂ ∂ Qˆ L ∂ Qˆ L 2 ∂C ˆ ˆ L L ˆ ˆ ˆ + uˆ Q = k Q (1 − C) P − Q C + k D (1 − C) + χ 2 rˆ Pˆ ∂ rˆ rˆ ∂ rˆ ∂ rˆ ∂ tˆ  ∂ Pˆ ∂ Cˆ +k P C + χ (30) ∂ rˆ ∂ rˆ with the initial conditions of ( 1 PˆL (ˆr , τ ) = 0 ( 1 Qˆ L (ˆr , τ ) = 0

ˆ ) ≤ rˆ ≤ R(τ ˆ ) (1 − ) R(τ ˆ ) 0 ≤ rˆ < (1 − ) R(τ ˆ ) ≤ rˆ ≤ R(τ ˆ ) (1 − ) R(τ ˆ ) 0 ≤ rˆ < (1 − ) R(τ

and finally  ˆ ˆ ∂ Mˆ ∂ Mˆ ˆ ˆ + χ ∂ P ∂C + uˆ P = − Mˆ k B Cˆ Pˆ − k D (1 − C)(1 − P) ∂ rˆ ∂ rˆ ∂ rˆ ∂ tˆ  ˆ  ˆ 1 ∂ rˆ 2 ∂ C −χ (1 − P) rˆ 2 ∂ rˆ ∂ rˆ with ( ˆ ˆ ˆ r , τ ) = 1 (1 − ) R(τ ) ≤ rˆ ≤ R(τ ) M(ˆ ˆ ). 0 0 ≤ rˆ < (1 − ) R(τ

(31)

246

G. J. Pettet et al.

4.1. Numerical methods. After mapping the moving boundary problem to a fixed grid the finite volume method was used to calculate the numerical solutions to the steady-state diffusion equation (25) and to integrate equation (26) to obtain the active cell velocity. The nonlinear hyperbolic equations describing the unlabelled cell distributions exhibit shock behavior and the source terms in these equations make obtaining accurate solutions with hyperbolic solvers difficult. In order to generate reasonable solutions a central difference scheme was employed with an additional small diffusive term included. In all solutions the diffusive coefficient in this added term has the value D = 5.0 × 10−3 . Accurate solutions to the labelled cell and microsphere problems were obtained using the weighted average flux method (WAF) by Toro (1989) with the ‘Superbee’ flux limiter as defined by Roe (1985). These higher order hyperbolic equation methods give more accurate solutions than simple upwind schemes but such schemes give very similar results. The outer boundary position was updated at each time step by a simple explicit Euler method applied to the appropriate outer boundary condition (27). A total of 150 spatial grid points were used with a fixed time step of 1 × 10−3 . 4.2. Numerical solutions. To illustrate the typical behavior of the model in one dimension, we consider the solutions for the three cases of the chemotactic response: χ < 0, χ = 0 and χ > 0. Firstly, the value of a number of important parameters in the model need to be considered. The time scale used corresponds approximately to the cell cycle time of those cells in the proliferating state, known to vary from between 8 to 24 hours (Tubiana, 1971). Incorporated in our model are a number of rate parameters whose value or dependence upon the extracellular environment is unknown. However, a few observations of some likely relationships between these parameters assists us in determining reasonable values to be adopted for these solutions. Cells may be expected to convert from the proliferating state to the G 0 state more quickly than they return, hence we take k P < k Q . As a multicell tumor spheroid will grow to a radius typically several times the diffusion length of oxygen in tissue, we expect the death rate k D to be small relative to the proliferation rate. Finally, the results reported by Dorie et al. (1982, 1986) indicate that inert microspheres are internalized at a rate of approximately half a spheroid diameter in 4 days. In line with these observations we adopt the following parameter values; R0 = 3.0,

k D = 0.1,

k Q = 0.9

and

k P = 0.05,

where the initial radius of the spheroid is taken to be somewhere near its expected steady-state value. The model has been tested for a range of parameter values (for both χ < 0 and χ > 0) with respect to the cell kinetic constants. In the cases described below the parameter being investigated was varied and the remaining cell kinetic parameters set at the nominal values. The varying of k P in the range 0 ≤ k P ≤ 10 resulted in little change to the overall form of the solution for both χ < 0 and χ > 0. As k P increased the

Migration of Cells

247

Figure 7. The evolution of the proliferating cell distribution from the initial to steady state.

time to steady state was longer and the final tumor size was larger. However, the final size became insensitive for very large values of k P . These observations can be interpreted physically as follows. For small values of k P , i.e., k P < 0.5, few quiescent cells return to the proliferating state and the effects of quiescence and cell death dominate. However, for larger value of k P , i.e., k P > 1, the number of proliferating cells is far greater than quiescent cells due to k Q < k P . The effect of varying k Q was found to be very similar to that of varying 1/k P . The effect of increasing the death rate k D was to reach the steady state in a considerably shorter time. Indeed it would probably be biologically infeasible that the death rate would be near that of the birth rate of the cells, i.e., k D  k B . For each simulation the nutrient profile at a number of different time steps as predicted by equation (25) is shown in Fig. 3. It is clear from this plot that there is a significant drop in nutrient supply to the inner regions of the tumor compared to the outer boundary and this is consistent with experimental findings (MuellerKlieser and Sutherland, 1982; Mueller-Klieser and Freyer, 1986). The case of χ = −1.0. This result leads to the ‘classical’ description of a multicellular spheroid: an outer periphery of proliferating cells surrounding an inner region of quiescent cells. Whilst the model produces this result it is unable to reproduce the internalization characteristics observed by Dorie et al. in their experiments. This is because the velocity profiles of both the proliferating and quiescent cells are consistently negative. This means that all labelled cells and inert microspheres will be convected towards the center of the spheroid, although at marginally

248

G. J. Pettet et al.

Figure 8. The change in the quiescent cell distribution with time, from the initial to steady state.

different velocities for varying values of χ. After a period of time no labelled cells will be found on the outer periphery as demonstrated in Fig. 4, a result which was consistent for a range of negative χ values. The case of χ = 0.0. In this case u P = u Q and the model is therefore not unlike that proposed by Ward and King (1998a), whereby a system of two different cell types with the same cellular velocity is considered. Following on from the intrinsic wave speed of the problem given by equation (7), this system will not exhibit any shock-like discontinuities as is expected for χ 6= 0. Furthermore, the resultant internalization characteristics for labelled cells, as predicted by our model, will be consistent with those for the case of χ < 0. The case of χ = 0.7. Figures 5 and 6 illustrate the development of the multicell tumor spheroid modeled by equations (25) to (31) as it grows in diameter from its initial state with outer radius R0 = 3.0 towards its steady-state value of approximately 4.5. Figures 5 and 6 show an initial transient where the tumor shrinks due to the assumption that cells are initially only in the quiescent phase. Figures 7 and 8 illustrate the distribution of proliferating and quiescent cell densities respectively. Figures 9 and 10 show the calculated cell velocities at corresponding times. From these results we observe that at steady state the cells near the spheroid outer boundary are dominated by those in the quiescent phase, and that a noticeable band of proliferating cells has developed a small distance into the spheroid body. This phenomenon is further illustrated by considering a plot

Migration of Cells

249

Figure 9. The proliferating cell velocity as a function of radial distance at a number of time steps.

ˆ in Fig. 11. Figure 10 also shows that of the mitotic index (determined by K B P) at steady state there is a division of the quiescent cell population into two subpopulations. These populations consist firstly of those cells which actively migrate towards the outer boundary (by chemotactic response) faster than the inward velocity induced by the proliferation and secondly of those whose ability to actively migrate has been debilitated, by lack of nutrient, and are swept towards the center of the spheroid. Whilst initially not intuitive this result has been observed to some degree by Mansbridge et al. (1992) who considered the effect of EGF-stimulation on A431 spheroids. Labelled cells (both proliferating and quiescent) and microspheres were introduced into the spheroid at time tˆ = τ with τ = 200.0. Figures 12–15 are, respectively, the frequency distributions for proliferating labelled cells, quiescent labelled cells, the total distribution of labelled cells and the distribution of microspheres at various times. It can be noted that there is a significant difference between the distributions of labelled cells and microspheres predicted by this model. Figure 14 shows the slow development of an internalized population of labelled cells with a significant sub-population of labelled cells remaining adjacent to the spheroid boundary. The details of this model suggest that the sub-population that is adjacent to the spheroid boundary is dominated by quiescent cells. Figure 15 clearly shows the internalization of the inert microspheres as a coherent pulse. It should be noted that Dorie’s original data show the frequency of cells plotted against the spheroid radius. We infer that this frequency relates to the fraction of

250

G. J. Pettet et al.

Figure 10. The quiescent cell velocity as a function of radial distance at a number of time steps.

labelled cells or microspheres at a specific radius from the center of the spheroid. To account for the geometry it is necessary to measure frequency in the sense of F ∝ RR 0

rˆ 2 φ(ˆr , tˆ) rˆ 2 φ(ˆr , tˆ)d rˆ

.

(32)

Note that as a result of this definition, when considering the labelled cell distribution it is not the simple sum of the frequency of each cell type.

5.

R ESULTS AND C ONCLUSIONS

We now discuss the behavior predicted by the model and its physical implications. Numerical solutions of the model indicate the following general behavior to be inherent from the effects we have included. • The tumor grows to a fixed size which it reaches after a long time. However, the changes in the tumor need not be monotonic towards this size. Instead the dynamics created by spatial variations of the distribution of the two cell states allows the tumor to both undershoot and overshoot the final size depending on its initial size. Figure 5 shows the size of a tumor shrinking briefly before growing towards a steady state due, in this case, to the assumption that cells are initially all in the quiescent phase of the cell cycle.

Migration of Cells

251

Figure 11. The mitotic index through the spheroid. • The model predicts a narrow layer of cells with a significant mitotic index near the outer surface of the tumor. Such mitosis requires that there are cells in the proliferating state and that there is sufficient nutrient for mitosis to occur. The model predicts that most cells on the periphery of the tumor are in the quiescent phase, driven there by chemotaxis, while the center of the tumor is also predominantly filled with quiescent cells. The distribution of quiescent cells is given in Fig. 8 along with the corresponding proliferating cell distribution in Fig. 7. In Fig. 11 the resulting mitotic index is plotted. • The variation in chemotactic response of the two states results in a selfsorting of the cells within the tumor. There is a movement of quiescent cells out and the resulting inward motion of the proliferating cells. The velocity of cells in the two states is shown in Figs 9 and 10. An average individual cell can be considered to start near the outer surface in a quiescent state. The high nutrient concentration there makes the cell transform to a proliferating state and it then proceeds to pass through the various stages of the cycle the whole time being swept into the tumor by the surrounding, more motile, quiescent cells. The cell proliferates, but as the nutrient level decreases, it either drops into the quiescent state to once again climb out to the outer surface or stalls in the cell cycle and moves further into the spheroid. The overall growth of the tumor is therefore tightly controlled by the rate of conversion between the phases. • If the outer layer of quiescent cells of a tumor is removed, perhaps by radiotherapy, then there will be a sudden increase of nutrient levels within the

252

G. J. Pettet et al.

Figure 12. The distribution of labelled proliferating cells.  = 0.1, τ = 200.

tumor. Cells stalled in the proliferating state will start to cycle making the tumor grow rapidly before shrinking back to its steady state. Much of the motivation for this research has been the experimental results of Dorie et al. (1982, 1986) regarding motion of microspheres and labelled cells within tumors. Many of the qualitative properties of the results cannot be explained by any of the existing theories. In addition we know of no similar experimental studies done to elucidate the internal motion of cells within tumors. Notwithstanding these caveats the model presented here has been used to predict the motion of cells and microspheres within a tumor near steady state. The model shows that the speed at which microspheres move into the tumor is dominated by the chemotactic coefficient while the spread of their distribution is governed by the width of the mitotic region. Figure 15 shows the predicted inward motion of the microspheres from the outer surface. When considering the motion of labelled cells we have assumed they are all initially in the quiescent state. The subsequent development of the distribution of both labelled quiescent cells and labelled proliferating cells is shown in Figs 12 and 13. This shows the inward motion of both distributions with the quiescent cells having two parts, one being related primarily to the initial distribution and the other being cells further into the tumor whose motion is delayed relative to the microspheres. This delay is because this part of the distribution arises only due to quiescent labelled cells becoming proliferating labelled cells, which are then carried into the tumor and subsequently become quiescent. The experimental data only measured the overall distribution of labelled cells and so the sum of proliferating and quiescent densities is given in Fig. 14. These plots show the qualitative

Migration of Cells

253

Figure 13. The distribution of labelled quiescent cells.  = 0.1, τ = 200.

properties of Dorie’s observations including sweeping in of all microspheres. In addition, the same model shows a slightly delayed movement of labelled cells into the interior with evidence of a significant number of labelled cells remaining near the surface. Numerical experiments have been carried out using this model to examine the behavior as parameter values are varied over physically reasonable ranges. The behavior described above continues to occur with variations in the size of the resulting spheroid, the thickness of the mitotic layer and the preponderance of proliferating cells in the interior. Extensions to account for more detailed understanding of the various rate constants and in particular their dependence on the nutrient concentration, of which little is known quantitatively, should therefore allow improved agreement with observations. Moreover, the initial state of the labelled cells used in the experiments is unknown. One of the underlying principles of this model has been that the kinetics of the labelled cells is identical to that of the unlabelled cells in the same state. However, we note from Dorie’s results that indeed a number of experiments were conducted with cells of different types from those in the spheroids. Furthermore, Dorie introduced drugs whose influence on the rate parameters and chemotactic response we have not attempted to quantify. These aspects can only be accounted for by further experimental input. In considering the behavior of this model we have neglected the effects due to creation of a necrotic region. The main emphasis of the model is to consider the behavior within the viable rim and it is anticipated that the results presented here

254

G. J. Pettet et al.

Figure 14. The distribution of labelled cells.  = 0.1, τ = 200. can be transferred to models where necrosis is explicitly included. Extensions to include necrosis are currently under consideration. We note that for all the examples considered in this paper we have taken cells in the quiescent state to be more motile than those in the proliferating state. The experimental evidence although supporting this assumption is not, however, sufficient to exclude the alternative case of (χ < 0), which may be true for some cell lines. In this case the model predicts that the microspheres are still internalized in a similar manner. However the internal structure is very different. The quiescent cells are pushed into the center of the spheroid and the proliferating cells are almost all restricted to a layer near the outer surface. In almost all cases there is a sudden jump, characterized by a shock created by the nonlinear hyperbolic equations, between the outer rim and the inner region. The model presented here has incorporated a very simple representation of the cell cycle, accounting for quiescence, proliferation and motility, into the classical models of tumor dynamics. Specifically, cells in the proliferating state direct their energy towards the activities required for proliferation while cells in the quiescent state more actively undergo chemotactic motion. The resulting predictions give an explanation for observed motion of cells in spheroids involving a novel mix of conventional ideas. The model indicates that overall growth of the tumor can be controlled not only by direct manipulation of the mitotic rate but, more importantly, by modification, and specifically reduction, of the conversion rates between quiescent and proliferating states. In addition it shows that removal of an outer layer of cells allows stalled proliferating cells to then be recruited into rapidly increasing the size of a tumor.

Migration of Cells

255

Figure 15. The distribution of passive microspheres. ACKNOWLEDGEMENTS CPP and MJT wish to thank Dr Margaret Peddie of the Department of Cell Sciences, the University of Southampton, for informative and useful discussions regarding this work. CPP appreciates the support provided by a QUT Visiting Fellowship during this period of research. MJT acknowledges the support of an Overseas Research Scholarship (ORS) and University of Southampton Faculty Studentship in providing support for him while undertaking this work. GJP thanks QUT for support, in the form of sabbatical leave, to complete this work.

R EFERENCES Burton, A. (1966). Rate of growth of solid tumors as a problem of diffusion. Growth 30, 157–176. Byrne, H. and M. Chaplain (1995). Growth of nonnecrotic tumors in the presence and absence of inhibitors. Math. Biosci. 130, 151–181. Byrne, H. and S. Gourley (1997). The role of growth factors in avascular tumour growth. Math. Comput. Modelling 26, 35–55. Chaplain, M. (1996). Avascular growth, angiogenesis and vascular growth in solid tumours: The mathematical modelling of the stages of tumour development. Math. Comput. Modelling 23, 47–87. Chaplain, M. and N. Britton (1993). On the concentration profile of a growth inhibitory factor in multicell spheroids. Math. Biosci. 115, 233–243.

256

G. J. Pettet et al.

Dorie, M., R. Kallman and M. Coyne (1986). Effect of cytochalasin b, nocodazole and irradiation on migration and internalization of cells and microspheres in tumor cell spheroids. Exp. Cell Res. 166, 370–378. Dorie, M., R. Kallman, D. Rapacchietta, D. van Antwerp and Y. Huang (1982). Migration and internalization of cells and polystrene microspheres in tumor cell spheroids. Exp. Cell Res. 141, 201–209. Greenspan, H. (1972). Models for the growth of a solid tumor by diffusion. Stud. Appl. Math. 51, 317–340. Greenspan, H. (1976). On the growth and stability of cell cultures and solid tumors. J. Theor. Biol. 56, 229–242. Hughes, F. and C. McCulloch (1991). Quantification of chemotactic response of quiescent and proliferating fibroblasts in boyden chambers by computer-assisted image analysis. J. Histochem. Cytochem. 39, 243–246. Knuechel, R. and R. Sutherland (1990). Recent developments in research with human tumor spheroids. Cancer J. 3, 234–243. Kunz-Schugart, L., M. Kreutz and R. Knuechel (1998). Multicellular spheroids: A threedimensional in vitro culture system to study tumour biology. Int. J. Exp. Pathol. 79, 1–23. Landman, K. A. and C. Please. Tumour dynamics and necrosis: Surface tension and stability. IMAJ. Math. Appl. Med. Biol. (submitted). Mansbridge, J., R. Knuechel, A. Knapp and R. Sutherland (1992). Importance of tyrosine phosphates in the effects of cell-cell contact and microenvironments on egfstimulated tyrosine phosphorylation. J. Cell. Physiol. 1, 433–442. McElwain, D. and G. Pettet (1993). Cell migration in multicell spheroids: Swimming against the tide. Bull. Math. Biol. 55, 655–674. Mueller-Klieser, W. (1997). Three-dimensional cell cultures: From molecular mechanisms to clinical applications. Am. J. Physiol. 273:4, 1109–1123. Mueller-Klieser, W. and J. Freyer (1986). Influence of glucose and oxygen supply conditions on the oxygenation of multicellular spheroids. Br. J. Cancer 53, 345–353. Mueller-Klieser, W. and R. Sutherland (1982). Oxygen tensions in multicell spheroids of two cell lines. Br. J. Cancer 45, 256–264. Murray, J. (1993). Mathematical Biology, 2nd edn, New York: Springer Verlag. Palka, J., B. Adelmann-Grill, P. Francz and K. Bayreuther (1996). Differentiation stage and cell cycle position determine the chemotactic response of fibroblasts. Folia Histochem. Cytobiol. 34, 121–127. Please, C., G. Pettet and D. McElwain (1998). A new approach to modelling the formation of necrotic regions in tumours. Appl. Math. Lett. 11, 89–94. Please, C., G. Pettet and D. McElwain (1999). Avascular tumour dynamics and necrosis. Math. Models Methods Appl. Sci. 9, 569–579. Roe, P. (1985). Some contributions to the modelling of discontinuous flows, in Large-scale Computations in Fluid Mechanics, Lectures in Applied Mathematics 22, B. Engquist, S. Osher and R. Somerville (Eds), Providence, RI: American Mathematical Society, pp. 163–193.

Migration of Cells

257

Sutherland, R. (1988). Cell and environment interactions in tumor microregions: The multicell spheroid model. Science 240, 177–184. Thompson, K. and H. Byrne (1999). Modelling the internalisation of labelled cells in tumour spheroids. Bull. Math. Biol. 61, 601–623. Toro, E. (1989). A weighted average flux method for hyperbolic conservation laws. Proc. R. Soc. Lond. A 423, 401–418. Tubiana, M. (1971). The kinetics of tumour cell proliferation and radiotherapy. Br. J. Radiol. 44, 325–347. Ward, J. and J. King (1998a). Mathematical modelling of avascular-tumour growth ii: Modelling growth saturation. IMA J. Math. Appl. Med. Biol. 15, 1–42. Ward, J. and J. King (1998b). Mathematical modelling of the effects of mitotic inhibitors on avascular tumour growth. J. Theor. Med. 1–25. Received 30 December 1999 and accepted 19 November 2000