Cell kinetics in tumour cords studied by a model with variable cell cycle length

Mathematical Biosciences 177&178 (2002) 103–125 www.elsevier.com/locate/mbs

Cell kinetics in tumour cords studied by a model with variable cell cycle length Alessandro Bertuzzi

a,*

, Antonio Fasano b, Alberto Gandolfi a, Doriana Marangi a,1

a

b

Istituto di Analisi dei Sistemi ed Informatica del CNR, Viale Manzoni 30, 00185 Rome, Italy Dipartimento di Matematica ‘U. Dini’, Universita
di Firenze, Viale Morgagni 67/A, 50134 Firenze, Italy Received 31 December 2000; received in revised form 2 August 2001; accepted 21 September 2001

Abstract A mathematical model is developed that describes the proliferative behaviour at the stationary state of the cell population within a tumour cord, i.e. in a cylindrical arrangement of tumour cells growing around a blood vessel and surrounded by necrosis. The model, that represents the tumour cord as a continuum, accounts for the migration of cells from the inner to the outer zone of the cord and describes the cell cycle by a sequence of maturity compartments plus a possible quiescent compartment. Cell-to-cell variability of cycle phase transit times and changes in the cell kinetic parameters within the cord, related to changes of the microenvironment, can be represented in the model. The theoretical predictions are compared against literature data of the time course of the labelling index and of the fraction of labelled mitoses in an experimental tumour after pulse labelling with 3 H-thymidine. It is shown that the presence of cell migration within the cord can lead to a marked underestimation of the actual changes along cord radius of the kinetics of cell cycle progression. Ó 2002 Elsevier Science Inc. All rights reserved. Keywords: Solid tumours; Cell migration; Cell kinetics; Mathematical model

1. Introduction According to a common view on tumour development, solid tumours start their growth as spheroids of proliferating cells that receive oxygen and nutrients necessary to survive and *

Corresponding author. Tel.: +39-06 771 6417; fax: +39-06 771 6461. E-mail address: [email protected] (A. Bertuzzi). 1 Present address: Ekar s.p.a., Via Terenzio 35, 00193 Roma, Italy. 0025-5564/02/$ - see front matter Ó 2002 Elsevier Science Inc. All rights reserved. PII: S 0 0 2 5 - 5 5 6 4 ( 0 1 ) 0 0 1 1 4 - 6

104

A. Bertuzzi et al. / Mathematical Biosciences 177&178 (2002) 103–125

proliferate only through diffusion from the external environment, the neighbouring blood vessels running externally to the tumour. This avascular phase is followed by angiogenesis, stimulated by factors released by tumour hypoxic cells, to develop a vascular network inside the tumour mass itself, that supports its further growth [1]. Experimental models of the avascular phase of tumour growth are the multicellular spheroids grown in vitro [2]. The growth of multicellular spheroids or of solid tumours in their initial avascular phase, up to 1–2 mm in diameter, has been investigated by mathematical models (see the review in [3]). The above view seems to be quite adequate for the initial growth of experimental tumours implanted subcutaneously, because tumour cells are placed in a space that is devoid of vessels, thus requiring angiogenesis to provide vessels to the tumour. For spontaneous tumours or experimental tumours implanted not subcutaneously, there are findings indicating that the tumour can initially grow by coopting existing host vessels. As the tumour grows, however, part of the initial vasculature in the internal region of the tumour regresses causing massive tumour cell death, and subsequently angiogenesis develops [4,5]. The structure of the vascular system that supplies the tumour is complex and irregular in most neoplastic tissues, so it is difficult to study the relationship between the distance from blood vessels and the cell proliferation, and therefore between the extent of vasculature and the tumour growth. In some tumours, however, it is possible to observe cylindrical arrangements of viable tumour cells, surrounded by necrosis, around blood vessels, and this simple geometry makes this study more feasible. These structures are named tumour cords [6–11]. In tumours that show vessel cooption, cuffs of viable tumour cells around the coopted blood vessels have been observed [4]. The mean thickness of the cords (i.e., the distance between vessel wall and necrosis) has been found to be 60–120 lm in different tumours, with a mean radius of central vessel of 10–40 lm [6–11]. The cell proliferation within the tumour cord induces migration of cells towards the periphery: cells are pushed by growing and dividing cells and eventually move into the necrotic zone. A mathematical model that describes the growth of an isolated tumour cord within a normal tissue and the formation of necrotic regions, taking into account the diffusion of a chemical critical for cell viability such as oxygen, was presented in [12]. The model predicts that a constant value for the radius of the cord will be attained in a finite time. Experimental observations have shown that the cell proliferation rate, measured by the fraction of cells in S phase and by the fraction of cells in mitosis, appears to decrease from the vessel wall to the periphery of the cord [6–8], and this phenomenon is likely to be related to the decrease of the concentration of oxygen, nutrients, and other critical chemicals. A mathematical model which describes the cell kinetics in a tumour cord at the stationary state, that is when the cord has attained the maximal radius compatible with the supply of vital substances, and the structure of the cell population is time invariant, was proposed in [13]. In [13], the model was formulated using the age formalism and the variability of the cell cycle phase transit times was disregarded, all proliferating cells being assumed to traverse the cycle in the same time. The variability of the microenvironment within the cord was assumed to affect only the transition of cells to a quiescent postmitotic compartment. The mathematical model for the cell population within a tumour cord at the stationary state, proposed in the present paper, takes into account the cell-to-cell variability of phase transit times and the decrease of the rate of progression across the cell cycle that is likely to accompany the migration of cells from the region close to the central vessel to the periphery of the cord. The time courses of the labelling index (LI) and of the fraction of labelled mitoses (FLM) after pulse la-

A. Bertuzzi et al. / Mathematical Biosciences 177&178 (2002) 103–125

105

belling of the S phase cells, as predicted by the model, are compared with experimental data reported in the literature.

2. The model We will describe a tumour cord around a blood vessel of radius r0 as a cylinder of tumour cells of radius R. At the stationary state, the radius R will represent the cord–necrosis interface. Let r denote the radial distance from the axis of the central blood vessel. As in [13], the changes of parameters and variables along the axial direction of the cord are neglected, and the tumour cord is considered as a continuous medium in which the number of tumour cells per unit volume is assumed constant in space and time. Cell motions in the axial direction are considered to be negligible, and cell migration in the radial direction is described by a velocity field uðr; tÞ, common to all cells, irrespective of cell position within the cycle and of its cycling or quiescent status. The cell cycle of tumour cells is represented, according to the model proposed by Kendall P [14] and Takahashi [15,16], by a sequence of discrete compartments of cell maturity. Let M ¼ i mi be the total number of maturity compartments, mi , i ¼ 1; . . . ; 4, being the number of compartments in G1, S, G2, and M phases, respectively. Let kk , k ¼ 1; . . . ; M, be the exit rate constant from compartment k. The rate constants kk may differ for G1, S, G2, and M phases and we will assume 8 kG1 k ¼ 1; . . . ; m1 ; > > < kS k ¼ m1 þ 1; . . . ; m1 þ m2 ; kk ¼ ð1Þ k ¼ m1 þ m2 þ 1; . . . ; m1 þ m2 þ m3 ; k > > : G2 kM k ¼ m1 þ m2 þ m3 þ 1; . . . ; M: We assume that newborn cells can arrest their progression across the cycle and become quiescent, so that only a fraction hðr; tÞ 2 ½0; 1 of newborn cells enter the G1 phase (see Fig. 1). The dependence of h on r and t reflects the dependence of cell arrest on the concentrations of oxygen, nutrients and other chemicals within the cord. Since in their motion towards the cord periphery cells are likely to experience more severe conditions, we will take the cell arrest as irreversible. Moreover, we assume that the rate of progression across the cycle can depend on these

Fig. 1. Scheme of the cell population model with discrete maturity compartments and quiescent compartment. For the symbols see the text.

106

A. Bertuzzi et al. / Mathematical Biosciences 177&178 (2002) 103–125

concentrations, so that the exit rate constants of the maturity compartments may change with r and t, although we exclude that the kk s can vanish. In particular, in view of the larger variability of G1 phase duration with respect to the other cycle phases [17], the cell microenvironment is likely to mainly affect the progression through that phase, so in part of the following treatment only the rate constant kG1 will be considered to be varying. If a stationary state is attained, the concentration of the chemicals along the radial direction will not change with the time, and h and the exit rate constants are functions of the radial distance only. Because of the worsening of the microenvironment when the distance from the central vessel increases, it is reasonable to assume that these functions of r are non-increasing. Although it is likely that cell death can occur within the cord when approaching the periphery, apoptosis within the cord is disregarded, assuming that at the stationary state cell death occurs immediately beyond the radius R. Denoting by nk ðr; tÞ the number of cells in the kth compartment, and by nQ ðr; tÞ the number of quiescent cells, in the unit volume at distance r and time t, from the previous assumptions we can write the following conservation equations: on1 1 o ðrun1 Þ ¼ 2hkM nM  k1 n1 ; þ r or ot onk 1 o ðrunk Þ ¼ kk1 nk1  kk nk ; þ r or ot

ð2Þ k ¼ 2; . . . ; M;

onQ 1 o ðrunQ Þ ¼ 2ð1  hÞkM nM : þ r or ot

ð3Þ ð4Þ

Moreover, denoting by nðr; tÞ the total cell density at r and t, we have on 1 o ð5Þ þ ðrunÞ ¼ kM nM : ot r or By assuming that n is constant and equal to  n, and that uðr0 ; tÞ ¼ 0 (i.e., there is no cell flux across the vessel wall), uðr; tÞ can be derived from Eq. (5) as Z 1 r zkM ðz; tÞnM ðz; tÞ dz: ð6Þ ruðr; tÞ ¼  n r0 In contrast with the present model, we will denote as ‘Takahashi model’ the cell population model with cell cycle structure as in Fig. 1, but without spatial structure and cell migration. In a stochastic framework, according to the Takahashi model cells spend in the maturity compartments random times that are independent and exponentially distributed with parameter kk . Thus, under assumption (1), the time spent in each cell cycle phase has a gamma distribution with, e.g. . The expected values for G1 phase, mean value equal to m1 =kG1 and coefficient of variation m1=2 1 of the number of cells in the maturity compartments of cell cycle and in the quiescent state are given in the Takahashi model by the solution of a set of ODEs whose right-hand side (r.h.s.) is the same as that of Eqs. (2)–(4). We note that, in deriving the distribution of times spent by cells in the cell cycle phases within a tumour cord, it should be taken into account that the maximal lifespan of a cell is equal to the time needed to migrate from the position at birth to the boundary with the necrotic zone (and thus is finite if the cell is born at r > r0 ) and that the rate of exit from a maturity compartment can be varying with the time.

A. Bertuzzi et al. / Mathematical Biosciences 177&178 (2002) 103–125

107

3. Stationary state Experimental data suggest that in most of the observed tumour cords the cord radius is constant with the time. We will assume that in these cords also the maturity distribution is in a stationary state, i.e., it is time invariant, so that nk ðr; tÞ ¼ nk ðrÞ. Let fk ðrÞ ¼ nk ðrÞ=n be the fraction of cells in the kth compartment of maturity at distance r in the stationary state. The conservation equations become 1 d ðruf1 Þ ¼ 2hkM fM  k1 f1 ; ð7Þ r dr 1 d ðrufk Þ ¼ kk1 fk1  kk fk ; r dr

k ¼ 2; . . . ; M;

1 d ðrufQ Þ ¼ 2ð1  hÞkM fM r dr with ruðrÞ ¼

Z

ð8Þ ð9Þ

r

ð10Þ

zkM fM dz: r0

The above equations can be rewritten in the following form:  Z r  1 df1 zkM fM dz ¼ 2hkM fM  k1 f1  kM fM f1 ; r r0 dr  Z r  1 dfk ¼ kk1 fk1  kk fk  kM fM fk ; k ¼ 2; . . . ; M; zkM fM dz r r0 dr  Z r  1 dfQ zkM fM dz ¼ 2ð1  hÞkM fM  kM fM fQ : r r0 dr

ð11Þ

ð12Þ

ð13Þ

Eqs. (7)–(10) are a non-linear system of integro-differential type with degeneracy at r ¼ r0 , since uðrÞ vanishes there. However, we can prove the following. Lemma 1. Suppose that h, k1 ; . . . ; kM are positive constants and 0:5 < h 6 1. Then the system (7)– (10) possesses a unique constant P solution f1 ; . . . ; fM , fQ such that 0 < fk < 1, k ¼ 1; . . . ; M, and    fQ ¼ 2ð1  hÞ 2 ½0; 1Þ. Moreover M k¼1 fk þ fQ ¼ 1. Proof. From (11) and (12), we obtain the following equations: 2hkM fM  k1 f1  kM fM f1 ¼ 0; kk1 fk1  kk fk  kM fM fk ¼ 0;

k ¼ 2; . . . ; M;

ð14Þ ð15Þ

from which we have f1 ¼

2hkM  fM ; k1 þ kM fM

ð16Þ

108

A. Bertuzzi et al. / Mathematical Biosciences 177&178 (2002) 103–125

kk1 ð17Þ fk1 ; k ¼ 2; . . . ; M: kk þ kM fM Supposing fM assigned in the denominators kk þ kM fM , Eqs. (16) and (17) are linear in fk , k ¼ 1; . . . ; M so that a non-zero solution will exist if and only if the determinant of coefficients is zero. The determinant is equal to M Y kk ; 1  2h kk þ kM fM fk ¼

k¼1

so that fM must be a root of the following algebraic equation of order M: M Y ð1 þ ak fM Þ ¼ 2h;

ð18Þ

k¼1

where ak ¼ kM =kk > 0. We note that the polynomial in the left-hand side (l.h.s.) of (18), PM ðfM Þ, has negative roots, PM ð0Þ ¼ 1, it is strictly increasing for fM P 0, and PM ð1Þ > 2 since aM ¼ 1. Having supposed 1 < 2h 6 2, we conclude that (18) has a unique real positive solution that is less than one. At this point, all the other fractions f1 ; . . . ; fM1 can be determined from (16) and (17). It is easy to check that the value of fM provided by the last equation in (17), PMi.e., for k ¼ M, is precisely the positive solution of (18). By summing and (15), we have m¼1 fk ¼ 2h  1 6 1. P (14) k þ fQ ¼ 1.  f From (13), we obtain fQ ¼ 2ð1  hÞ and thus M m¼1 We note that the fractions fk are the same as the fractions obtained from the solution in asynchronous exponential growth of the Takahashi model with the same parameter values. In the general case of h and kk dependent on r, because of the degeneracy at r0 , it is not possible to prescribe data for r ¼ r0 , and the values of the unknown functions fk ðrÞ and fQ ðrÞ at r ¼ r0 (denoted by fk0 , k ¼ 1; . . . ; M, and fQ0 ) must be deduced from suitable conditions. We will impose that the right derivative of the solution at r0 is bounded. Thus the values of fk0 turn out to be precisely the ones given by Eqs. (16)–(18) with h and kk evaluated at r0 , and fQ0 ¼ 2ð1  hðr0 ÞÞ. We will take hðr0 Þ > 0:5, which is necessary and sufficient to have fk0 positive. If we assume that hðrÞ and the coefficients kk ðrÞ have enough regularity, knowing their derivatives up to some order at r ¼ r0 , it is possible to calculate the derivatives at r0 of the possible solutions of (7), (8) and (10) up to the desired order. We show the procedure for the first derivative in Appendix A. We now establish some useful a priori results. Lemma 2. Let fk ðrÞ, k ¼ 1; . . . ; M, and fQ ðrÞ be solutions of the system (7)–(10) with the values specified above for r ¼ r0 . Then we have, for r 2 ½r0 ; R, 0 < fk ðrÞ < 1;

k ¼ 1; . . . ; M;

0 6 fQ ðrÞ < 1;

as long as h, k1 ; . . . ; kM are positive functions, and M X fk ðrÞ þ fQ ðrÞ ¼ 1:

ð19Þ

ð20Þ

k¼1

Proof. First we show that fk > 0. This is true for r ¼ r0 . Suppose that fj is the first function in the set ðf1 ; . . . ; fM Þ that vanishes at the point rj > r0 . For j > 1, using (12) we obtain dfj =drjr¼rj > 0,

A. Bertuzzi et al. / Mathematical Biosciences 177&178 (2002) 103–125

109

unless fj1 ðrj Þ ¼ 0, and for j ¼ 1 we need fM ðrj Þ ¼ 0 to avoid a positive derivative. In other words, all the functions must vanish simultaneously. Now we note that, since ðf1 ; . . . ; fM Þ is by assumption a solution of (11) and (12), we may regard fM as a known function and thus (11) and (12) with k ¼ 1; . . . ; M  1 as a linear system of differential equations, which is certainly not degenerate in a left neighbourhood of rj . Therefore, in such an interval, the only solution which is compatible with fk ðrj Þ ¼ 0 is fk ðrÞ 0, k ¼ 1; . . . ; M  1, which contradicts our assumptions on the initial values. Concerning fQ , it is fQ P 0 at r0 . For r > r0 , from (13) it is dfQ =dr > 0 if fQ ðrÞ PM< 0. Thus fQ cannot become negative. We now prove fk < 1 and fQ < 1 showing that F ¼ k¼1 fk þ fQ ¼ 1. This is true for r ¼ r0 . From Eqs. (11)–(13), the equation satisfied by F is u

dF ¼ ð1  F ÞkM fM : dr

ð21Þ

The latter equation shows that F is increasing in the set {r: F < 1} and decreasing in the set {r: F > 1}. Now, if one of these sets is not empty we see, integrating backward, that the condition F ðr0 Þ ¼ 1 cannot be matched. Then F must be constant and equal to 1.  Remark 1. If hðrÞ vanishes beyond some r , following the previous argument we can still assert that, if rj is the first point where fj vanishes, then all fk with lower index must vanish too and remain zero. However, also in this case we may argue as before, contradicting fk ðr0 Þ > 0 for every k. It is easy to see that the solutions of Eqs. (12) and (13) still respect F ¼ 1 and hence the inequalities fk < 1. The proof of the following existence and uniqueness theorem can be found in Appendix B. Theorem 1. Suppose h, k1 ; . . . ; kM are positive C 2 functions of r in ½r0 ; R, and that 0:5 < hðr0 Þ 6 1. Then the system (7)–(10) has one unique solution, with the functions f1 ; . . . ; fM , fQ taking the values f10 ; . . . ; fM0 , fQ0 for r ¼ r0 . In Appendix B, the concept of solution to the degenerate problem (7)–(10) is also defined. To obtain a numerical solution of Eqs. (7), (8) and (10), we note that these equations can be rewritten as the following set of integral equations: Z r  Z r Z r zkM fM dz f1 ðrÞ ¼ 2 zhkM fM dz  zk1 f1 dz; ð22Þ r0

Z

r

r0

r0

r0

 Z r Z r zkM fM dz fk ðrÞ ¼ zkk1 fk1 dz  zkk fk dz; r0

k ¼ 2; . . . ; M:

ð23Þ

r0

A solution of the above equations can be obtained by an iterative procedure. We define as x the vector of all the values fk ðri Þ, k ¼ 1; . . . ; M, ri being equispaced mesh points in ½r0 ; R, and rewrite Eqs. (22) and (23) as x ¼ GðxÞ, computing the integrals by the trapezoidal rule. This system of algebraic equations was solved using the iteration [18] xnþ1 ¼ xn þ aðGðxn Þ  xn Þ;

ð24Þ

110

A. Bertuzzi et al. / Mathematical Biosciences 177&178 (2002) 103–125

Fig. 2. Cell cycle phase fractions and fraction of quiescent cells in the stationary state as a function of radial distance. Kinetic parameters: m1 ¼ 2, m2 ¼ 9, m3 ¼ m4 ¼ 2, kG1 ¼ 0:25 h1 , kS ¼ 1:0 h1 , kG2 ¼ 1:0 h1 , kM ¼ 2:0 h1 , and hðrÞ ¼ 1=ð1 þ 16ððr  r0 Þ=ðR  r0 ÞÞ4 Þ. Geometric parameters: r0 ¼ 40 lm, R ¼ 100 lm.

where xn denotes the solution at the nth step. The convergence (that resulted very slow indeed) was achieved with small positive a. Fig. 2 represents the steady-state fractions of cells in the G1, S, G2, and M phases, obtained by summing the fractions fk ðrÞ over the appropriate indexes, and the quiescent fraction fQ ðrÞ. The kinetic and geometric parameters of the simulated tumour cord are reported in the figure legend. It can be observed that the ratios of the phase fractions with respect to the fraction of proliferating cells are not constant with r. An experimental technique largely used to study the cell proliferation within tumour cords is based on the measurement of the fraction of cells in S phase at various radial distances, obtained as the fraction of labelled cells immediately after the injection of a radioactive DNA precursor, or of the fraction in M phase, obtained by examination of histological slices of the tumour [6–8]. Thus, the question arises if these data allow us to distinguish whether there is a decelerated progression through G1 phase, as the radial distance increases, or there is an increase in the fraction of postmitotically arrested cells, when all the other cell cycle parameters are constant with r. We prove that this is not possible, see Appendix C, as far as the S and M phase fractions at the stationary state are considered. However, the technique of labelling S-phase cells followed by autoradiography also allows to measure the time course of the fraction of labelled cells as well as the fraction of labelled cells in mitosis. As we will see in the following, the time evolution of these fractions can be different if a decreasing rate of progression in G1 or a true postmitotic cell arrest occurs. 4. Pulse labelling The injection of a labelled DNA precursor will label tumour cells in a way that can be considered equivalent in practice to a true impulsive labelling. If the injection occurs at t ¼ 0, we will assume that at t ¼ 0þ all (and only) S phase cells are labelled. Thereafter, the subpopulation of labelled cells progresses through the cycle and cells undergo division: we will assume that, at each division, cells generated by labelled cells remain labelled.

A. Bertuzzi et al. / Mathematical Biosciences 177&178 (2002) 103–125

111

Let fkl ðr; tÞ be the fraction of labelled cells in the kth compartment at distance r and time t. We will assume here that the exit rate constants kS , kG2 and kM are constant with r. By writing the conservation equations for the labelled subpopulation in all the maturity compartments and in the quiescent compartment, we obtain for the time evolution of the labelled fractions the following equations: of1l 1 o ðruf1l Þ ¼ 2hkM fMl  k1 f1l ; þ r or ot ofkl 1 o l  kk fkl ; þ ðrufkl Þ ¼ kk1 fk1 ot r or

ð25Þ k ¼ 2; . . . ; M;

ofQl 1 o ðrufQl Þ ¼ 2ð1  hÞkM fMl ; þ r or ot where Z r zfM ðzÞ dz; ruðrÞ ¼ kM

ð26Þ ð27Þ

ð28Þ

r0

with the initial conditions  fk ðrÞ k ¼ m1 þ 1; . . . ; m1 þ m2 ; l fk ðr; 0Þ ¼ 0 otherwise;

fQl ðr; 0Þ ¼ 0:

ð29Þ

Since we assume that cell labelling occurs in a population at the stationary state, the functions fk ðrÞ and uðrÞ in the above equations are the fractions of cells and the velocity field in the stationary state, as defined in the previous section. Eqs. (25)–(27) can be rewritten as of1l of l þ u 1 ¼ 2hkM fMl  k1 f1l  kM fM f1l ; ot or ofkl of l l  kk fkl  kM fM fkl ; þ u k ¼ kk1 fk1 ot or

ð30Þ k ¼ 2; . . . ; M;

ð31Þ

ofQl ofQl ð32Þ þu ¼ 2ð1  hÞkM fMl  kM fM fQl : ot or Eqs. (30)–(32) are a system of first-order and linear equations, whose characteristic lines tend to have vanishing slope as r ! r0 . Thus the initial conditions (29) determine completely the solution for t P 0 and r 2 ðr0 ; R. For r ¼ r0 , by denoting fkl0 ðtÞ ¼ fkl ðr0 ; tÞ and fQl0 ðtÞ ¼ fQl ðr0 ; tÞ, we have df1l0 ð33Þ ¼ 2hðr0 ÞkM fMl0  kG1 ðr0 Þf1l0  kM fM0 f1l0 ; dt dfkl0 l0  kG1 ðr0 Þfkl0  kM fM0 fkl0 ; ¼ kG1 ðr0 Þfk1 dt

k ¼ 2; . . . ; m1 ;

dfml01 þ1 ¼ kG1 ðr0 Þfml01  km1 þ1 fml01 þ1  kM fM0 fml01 þ1 ; dt

ð34Þ ð35Þ

112

A. Bertuzzi et al. / Mathematical Biosciences 177&178 (2002) 103–125

dfkl0 l0 ¼ kk1 fk1  kk fkl0  kM fM0 fkl0 ; dt

k ¼ m1 þ 2; . . . ; M;

dfQl0 ¼ 2ð1  hðr0 ÞÞkM fMl0  kM fM0 fQl0 ; dt with the initial conditions  0 fk k ¼ m1 þ 1; . . . ; m1 þ m2 ; fkl0 ð0Þ ¼ 0 otherwise;

ð36Þ ð37Þ

fQl0 ð0Þ ¼ 0:

ð38Þ

The numerical solution of Eqs. (25)–(29) was computed by a finite-difference scheme over a rectangular grid with fixed mesh size on ½0; T  ½r0 ; R. Time derivatives were approximated by forward differences whereas spatial derivatives by backward differences. The scheme required the knowledge of the solution at r ¼ r0 , and this is given by the solution of Eqs. (33)–(38). The time evolution of the fractions of labelled cells allows to compute the evolution at different radial distances of various observable quantities, such as the total fraction of labelled cells, LIðr; tÞ, and the fraction of labelled cells within the cells in mitosis, FLMðr; tÞ, according to PM l M X k¼Mm þ1 fk ðr; tÞ l l : ð39Þ fk ðr; tÞ þ fQ ðr; tÞ; FLMðr; tÞ ¼ PM 4 LIðr; tÞ ¼ k¼1 k¼Mm4 þ1 fk ðrÞ Numerical simulations have shown that, for large t, LI and FLM converge to a common value which is independent of r. Because Eqs. (25)–(28) admit stationary solutions of the form fkl ðrÞ ¼ cfk ðrÞ and fQl ðrÞ ¼ cfQ ðrÞ with 0 6 c 6 1, the observed behaviour suggests that the solution of Eqs. (25)–(29) tends to a stationary solution with c depending on the initial condition. If it is so, it is easy to see from Eq. (39) that both LIðr; tÞ and FLMðr; tÞ will tend to c independently of r.

5. Application to experimental data We applied the above model to the analysis of experimental data reported by Hirst and Denekamp [7]. These data have been obtained from the KHH mammary carcinoma, a tumour that grows in mice showing a corded structure. As reported in [7], when the tumour reached a diameter of 7.5–10 mm, the animals were injected with 3 H-thymidine and, at different times after label injection, the animals were sacrificed and the tumours excised. Histological sections were scanned and cells categorized in relation to the nearest visible blood vessel. Three classes were distinguished: cells close to vessels, cells close to necrosis and cells in the intermediate zones. The fraction of labelled cells (LI) and the FLM were measured in the three above defined zones up to 50 h after label injection. The LI and FLM data were analyzed under two different hypotheses. First, we assumed that progression across the G1 phase of cell cycle is not affected by the radial distance of cells, whereas h decreases with r. Alternatively, we assumed that kG1 is affected by the radial distance whereas no cell arrest occurs (i.e., hðrÞ ¼ 1). Under the first hypothesis, using a trial-and-error procedure, a reasonable fitting of the data was obtained with the following choice of parameter values: m1 ¼ 2, m2 ¼ 9, m3 ¼ m4 ¼ 2, kG1 ¼ 0:278 h1 , kS ¼ 1:059 h1 , kG2 ¼ 1:333 h1 , kM ¼ 3:333 h1 , and hðrÞ given by

A. Bertuzzi et al. / Mathematical Biosciences 177&178 (2002) 103–125

hðrÞ ¼

1 1 þ 16ððr  r0 Þ=ðR  r0 ÞÞ4

113

ð40Þ

:

The geometric parameters of the cord were r0 ¼ 40 lm and R ¼ 100 lm, according to the average measured diameters. The values chosen for the mi ’s and the exit rate constants are such that the mean values of TG1 , TS , TG2 and TM are 7.2, 8.5, 1.5 and 0.6 h, respectively, that is the phase transit times estimated in [7] for the cells in the inner zone of the cord. Fig. 3 shows the behaviour of LI in the inner and the outer zone of the cord. The predicted time course was calculated at r ¼ 50 lm and at r ¼ 90 lm, respectively. In Fig. 4, the experimental FLM data are shown together with their predictions. We note, in particular, that the theoretical prediction of the FLM curve in the outer zone fails to follow the rising front of the second wave of experimental data. Under the second hypothesis, data were fitted by choosing the following values of the parameters: m1 ¼ 2, m2 ¼ 7, m3 ¼ m4 ¼ 2, kS , kG2 and kM as before, kG1 ðrÞ given by kG1 ðrÞ ¼ 0:1 þ

0:55 1 þ 4:5ððr  r0 Þ=ðR  r0 ÞÞ1:5

h1 ;

ð41Þ

and hðrÞ ¼ 1. The time course of the prediction of the LI data is depicted in Fig. 5, showing a fitting of comparable quality with respect to the fitting obtained with h dependent on r. In the prediction of the FLM curve in the outer zone of the cord, the rising front of the second wave appears now to be delayed as shown by the data (see Fig. 6). We stress that the above sets of fitting parameters are not to be considered as optimal estimates. Since multiple curves have to be fitted (LI and FLM curves at different radial distances), the choice of the index to be minimized in an optimal estimation procedure would require an accurate analysis of data statistics, which is beyond the scope of the present paper. Moreover, the long computation time required to compute the model response was an obstacle to the use of a minimization algorithm. It can be noted that the second wave of the predicted FLM curve in the outer zone (Fig. 6) is slightly delayed although the value of kG1 at 90 lm, as given by Eq. (41), is 0.0243 h1 , which

Fig. 3. Experimental data of LI as a function of time in the inner () and the outer zone (j) of the cord. Data replotted from Hirst and Denekamp [7]. Model predictions with kG1 constant and h dependent on r. The function hðrÞ and the parameter values are reported in the text.

114

A. Bertuzzi et al. / Mathematical Biosciences 177&178 (2002) 103–125

Fig. 4. Experimental data of FLM () as a function of time in the inner (upper panel) and the outer zone (lower panel) of the cord. Data replotted from Hirst and Denekamp [7]. Model predictions with kG1 constant and h dependent on r (––). The function hðrÞ and the parameter values are reported in the text.

Fig. 5. Experimental data of LI as in Fig. 3. Model predictions with h ¼ 1 and kG1 dependent on r. The function kG1 ðrÞ and the parameter values are reported in the text.

A. Bertuzzi et al. / Mathematical Biosciences 177&178 (2002) 103–125

115

Fig. 6. Experimental data of FLM () as in Fig. 4. Model predictions with h ¼ 1 and kG1 dependent on r (––). The function kG1 ðrÞ and the parameter values are reported in the text.

would correspond in the Takahashi model to a mean TG1 of 82.3 h (whereas kG1 at 50 lm is 0.321 h1 , corresponding to a mean TG1 of 6.2 h). This behaviour is due to the cell migration within the cord. To show the effect of cell migration on the time course of FLM, we compared the response of our model at r ¼ 50 and 90 lm to the FLM curves predicted by the Takahashi model (i.e., a model without cell migration) having kG1 equal to kG1 ðrÞ at r ¼ 50 and 90 lm, respectively, and the other parameters unchanged. Fig. 7 shows this comparison: whereas the FLM curve in the inner zone slightly differs from that predicted by the Takahashi model, the difference is very marked when the FLM curve in the outer zone of the cord is considered, where the FLM predicted by the present model shows a moderate delay of the second wave. This marked difference can be explained considering that the labelled mitoses observed in the outer zone of the cord are the mitoses of cells that have spent the G1 phase at smaller radial distances where the progression across G1 phase is faster. Because of cell migration, the mean cell cycle duration of cells dividing at the periphery will thus be shorter than that of non-moving cells characterized by the cell kinetic parameters of the cord periphery. Due to their movement, cells progress through the cycle according to time-varying cell kinetic parameters. An analysis of the FLM curves for spatially homogeneous cell populations under non-stationary conditions can be found in [19].

116

A. Bertuzzi et al. / Mathematical Biosciences 177&178 (2002) 103–125

Fig. 7. Comparison between the FLM curves predicted by the model of the tumour cord (- - - -) and the FLM curves predicted by the Takahashi model (––): inner zone (upper panel), outer zone (lower panel). For the values of the parameters see the text.

Without taking into account cell migration, previous analyses of the FLM curves, based on cell population models similar to the Takahashi model, had found no difference [6] or only slight differences [7] in the kinetic parameters of cell cycle in regions proximal to the central vessel and in regions at the periphery of the cord. 6. Concluding remarks The analysis of labelling data here proposed suggests that, because of the decay of oxygen and nutrient concentration in tumour cords, the progression across cell cycle is slowed down as the distance from the central vessel increases. We have found that the retarded progression suggested by the experimental data can be achieved in the model by suitably lowering only the rate of progression in G1. In the presence of decreased rates of progression through cell cycle phases at the periphery of tumour cords, the kinetic differences between inner and outer zones can be largely masked by cell migration when data from pulse labelling are considered. Thus, more in general, our results in-

A. Bertuzzi et al. / Mathematical Biosciences 177&178 (2002) 103–125

117

dicate that a correct analysis of data from tumour microregions requires that the possible cell migration through these regions has to be taken into account. Finally, we note that the continuous-medium approach here adopted was able to describe substantially the experimental data considered, although the size of the particles (the cells) is comparable in the tumour cords to the spatial scale. An approach based on the behaviour of individual cells could, however, improve the description of cell motion, in particular near the wall of the central vessel. Acknowledgements The support of the Progetto Strategico del CNR ‘Metodi e Modelli Matematici per lo Studio dei Fenomeni Biologici’ is gratefully acknowledged. Appendix A We show here how to compute the first right derivatives at r ¼ r0 , denoted as fk00 , of the solutions fk of Eqs. (7), (8) and (10), assuming that the functions hðrÞ and kk ðrÞ have enough regularity. Let us insert in (7), (8) and (10) the first-order approximations hðrÞ ’ h0 þ h00 ðr  r0 Þ; 0

00

kk ðrÞ ’ k0k þ k00 k ðr  r0 Þ;

0

where h ¼ hðr0 Þ, h ¼ h ðr0 Þ, and mation fk ðrÞ ’ fk0 þ fk00 ðr  r0 Þ;

k0k

¼ kk ðr0 Þ,

k00 k

¼

k ¼ 1; . . . ; M;

k0k ðr0 Þ,

and look for the first-order approxi-

k ¼ 1; . . . ; M:

ðA:1Þ

From div u ¼ kM fM , we deduce 0 div u ’ k0M fM0 þ ðk0M fM00 þ k00 M fM Þðr  r0 Þ;

ðA:2Þ

and, from Eq. (10), 1 r2  r02 ’ k0M fM0 ðr  r0 Þ: u ’ k0M fM0 ðA:3Þ r 2 Using all the above approximations in (7) and (8), the zero-order terms cancel because of (14) and (15) written with h and kk evaluated at r0 . Equating the first-order terms we arrive at the system 00 00 0 0 00 0 0 0 0 00 0 0 f100 ð2k0M fM0 þ k01 Þ ¼ f10 ðk00 M fM þ k1 Þ þ 2fM ðh kM þ h kM Þ þ fM ð2h kM  f1 kM Þ;

ðA:4Þ

00 00 0 0 0 00 0 0 00 k ¼ 2; . . . ; M; fk00 ð2k0M fM0 þ k0k Þ ¼ fk0 ðk00 M fM þ kk Þ þ kk1 fk1  fM fk kM  kk1 fk1 ; 00 00 00 which defines fk , k ¼ 2; . . . ; M  1, recursively in terms of f1 and fM by means of ! ! k1 k k1 k X Y X Y k0j k0j 0 00 0 00 0 0 bki fki  kM fM bki fki kk fk ¼  a a j¼ki j j¼ki j i¼0 i¼0 ! 0 k2 Y k k X Y kj 0 0 f100 ; k ¼ 2; . . . ; M; þ k00 f þ k ki1 ki1 1 a j i¼0 j¼ki j¼2

ðA:5Þ

ðA:6Þ

118

A. Bertuzzi et al. / Mathematical Biosciences 177&178 (2002) 103–125

00 0 00 where aj ¼ 2k0M fM0 þ k0j and bk ¼ k00 M fM þ kk . Using (A.6) for k ¼ M and replacing f1 by the expression deduced from (A.4), we obtain the equality: ! # " M1 M X Y k0j 0 00 0 kM fM 1 þ fMi a j¼Mi j i¼0 ! M1 M M2 Y M X Y X k0j 0 0 bMi fMi þ k00 ¼ Mi1 fMi1 a j¼Mi j i¼0 i¼0 j¼Mi ! M Y k0j 00 00 0 0 00 0 00 0 0 0 0 ½f10 ðk00 þ ðA:7Þ M fM þ k1 Þ þ 2fM ðh kM þ h kM Þ þ kM fM ð2h  f1 Þ: a j j¼1

At this point we note that ! M Y k0j f0 ð2h0  f10 Þ < 1  1 0 < 1; a 2h j¼1 j owing to (18) and taking into account that aj > k0M fM0 þ k0j , and that 1

f10 k01 ¼ <1 2h0 k01 þ k0M fM0

because of (16). Thus, by rearranging (A.7), it can be seen that the coefficient of fM00 is positive. Therefore (A.7) can be solved for fM00 and all the other unknowns fk00 can be calculated. We remark 00 that, if h00 ¼ k00 k ¼ 0, k ¼ 1; . . . ; M, then fk ¼ 0 for all k. Following the above procedure higher order derivatives can be calculated, of course their expressions will be much more complicated. Appendix B. Proof of Theorem 1 The proof goes through the following steps: (i) we consider a family of regularized Cauchy problems obtained from (11)–(13) by adding e > 0 to the coefficients of the derivatives:  Z r  1 df1 ¼ 2hkM fM  k1 f1  kM fM f1 ; zkM fM dz þ e ðB:1Þ r r0 dr 

 dfk zkM fM dz þ e ¼ kk1 fk1  kk fk  kM fM fk ; dr r0  Z r  1 dfQ zkM fM dz þ e ¼ 2ð1  hÞkM fM  kM fM fQ ; r r0 dr 1 r

Z

r

k ¼ 2; . . . ; M;

ðB:2Þ ðB:3Þ

with initial data f10 ; . . . ; fM0 , fQ0 ; (ii) we show that letting e # 0 we obtain a unique limit, which satisfies (7)–(10) and that we define to be the solution of the degenerate system (7)–(10) with the values f10 ; . . . ; fM0 , fQ0 at r ¼ r0 . ðeÞ ðeÞ ðeÞ Step (i): First of all, we establish the existence and uniqueness of a solution f1 ; . . . ; fM , fQ to the regularized non-linear Cauchy problem (B.1)–(B.3). This can be done by means of the following fixed point argument. For fixed e, consider the set of functions:

A. Bertuzzi et al. / Mathematical Biosciences 177&178 (2002) 103–125

  Uðr Þ ¼ u 2 C½r0 ; r  : uðr0 Þ ¼ fM0 ; 0 6 u 6 1

119

ðB:4Þ

with r > r0 . For a given u in this set, solve the non-linear non-degenerate Cauchy problem for r 2 ½r0 ; r   Z r  1 dw1 ¼ 2hkM wM  k1 w1  kM wM w1 ; zkM u dz þ e ðB:5Þ r r0 dr  Z r  1 dwk zkM u dz þ e ðB:6Þ ¼ kk1 wk1  kk wk  kM wM wk ; k ¼ 2; . . . ; M; r r0 dr wk ðr0 Þ ¼ fk0 ;

k ¼ 1; . . . ; M:

ðB:7Þ

It is useful to remark that if we divide (B.5) and (B.6) by the coefficient of the derivatives, the resulting system is such that all the coefficients on the r.h.s. can be estimated by constants independent of r and u 2 U (and of course depending on e, which however is fixed). From well-known classical theorems (see, e.g., [20]) we can say that the system (B.5)–(B.7) has one unique solution ðw1 ; . . . ; wM Þ in an interval ½r0 ; r , bounded independently of the choice of u. We can even say that 0 < wk < 1 for k ¼ 1; . . . ; M, on the basis of the same a priori argument we have used in Lemma 2. This implies two main consequences: (a) The solution ðw1 ; . . . ; wM Þ can be continued up to r , as it follows again from the classical theory, since the solution can terminate only at a point where one of the wk has no limit (but they are all Lipschitz continuous and in addition we have uniform bounds for jdwk =drj) or blows up (but they are all bounded). (b) For any r > r0 , if we put ~ ¼ wM ¼: Mu; u

ðB:8Þ

we have defined a mapping M from Uðr Þ into itself. If we now take two functions u1 , u2 2 Uðr Þ and we compute the corresponding solutions ðiÞ ðiÞ ð1Þ ð2Þ ðw1 ; . . . ; wM Þ, i ¼ 1; 2, we can easily see that the differences dk ¼ wk  wk satisfy a linear differential system with zero Cauchy data, in which the free terms are such that their sup–norm is estimated in terms of ku1  u2 kC½r0 ;r  . Therefore we immediately conclude that: ðaÞ the mapping M is completely continuous with respect to the sup–norm (remember that the image of U is a set of uniformly Lipschitz continuous functions); ðbÞ there is some constant K, independent of the ~1  u ~ 2 kC½r0 ;r  6 Kðr  r0 Þ ku1  u2 kC½r0 ;r  . Thanks to ðaÞ we can choice of u1 , u2 , such that ku apply Schauder’s theorem (the set U is bounded, closed and convex) for any r , implying the existence of at least one fixed point, which necessarily gives a solution to the original regularized problem. Thanks to ðbÞ we can say that, if r  r0 is sufficiently small, the mapping M is contractive and the solution is unique. This same argument can be used to exclude the possibility that the solution can bifurcate, say at r, for some larger value of r , since we may repeat the procedure starting from r < r and sufficiently close to r, and construct a contraction mapping in a sufficiently small interval. Therefore, we have proved that the regularized system is uniquely solvable over an ðeÞ fk < 1, k ¼ 1; . . . ; M, are satisfied (as a property of arbitrary interval ½r0 ; r . The inequalities PM0 < ðeÞ ðeÞ all the wk ’s irrespective of u 2 U), and k¼1 fk ðrÞ þ fQ ðrÞ ¼ 1 (arguing as in Lemma 2). Remark 1 following Lemma 2 also applies to the regularized system. ðeÞ Step (ii): For any e > 0 all the derivatives dfk =dr for r ¼ r0 are zero. Therefore, we may not expect that they converge uniformly as e tends to zero in the vicinity of r0 , since the derivatives of the possible solution will be generally different from zero (see Appendix A). However, we show that the

120

A. Bertuzzi et al. / Mathematical Biosciences 177&178 (2002) 103–125 ðeÞ

derivatives of fk remain bounded uniformly in e. Writing the regularized system isolating the first derivatives on the l.h.s and integrating over ðr0 ; rÞ, r0 6 r0 < r, we can derive a system of integral ðeÞ ðeÞ ðeÞ equations for the ratios .k ðr; r0 Þ ¼ ½fk ðrÞ  fk ðr0 Þ=ðr  r0 Þ, in which linear approximations of h and kk ’s are taken around r0 (indeed high-order terms in r  r0 do not involve e in a critical way). For instance, for r0 ¼ r0 (actually the only critical point), and with the notations of Appendix A, we have Z rn o ðeÞ ðeÞ ðeÞ ðeÞ ðeÞ ðeÞ .1 ðr; r0 Þ ¼  K1 ðg; rÞ.1 ðg; r0 Þ þ H1 ðg; rÞ.M ðg; r0 Þ dg þ K1 ðrÞ; ðB:9Þ r0

Z r ðeÞ ðeÞ ðeÞ ðeÞ ðeÞ Kk ðg; rÞ.k ðg; r0 Þ  Gk1 ðg; rÞ.k1 ðg; r0 Þ .k ðr; r0 Þ ¼  r0  ðeÞ ðeÞ ðeÞ þ Hk ðg; rÞ.M ðg; r0 Þ dg þ Kk ðrÞ; k ¼ 2; . . . ; M;

ðB:10Þ

where ðeÞ

Kk ðg; rÞ ¼

ðeÞ

H1 ðg; rÞ ¼

ðeÞ

Hk ðg; rÞ ¼

ðeÞ

k0k þ k0M fM0 h i; e ðr  r0 Þ k0M fM0 þ gr 0 k0M ðf10  2h0 Þ h i; e ðr  r0 Þ k0M fM0 þ gr 0 k0 f 0 hM k i; e ðr  r0 Þ k0M fM0 þ gr 0

k ¼ 2; . . . ; M;

k0 h k1 i; e ðr  r0 Þ k0M fM0 þ gr 0

k ¼ 2; . . . ; M;

Gk1 ðg; rÞ ¼

ðeÞ K1 ðrÞ

ðeÞ Kk ðrÞ

k ¼ 1; . . . ; M;

1 ¼ r  r0 1 ¼ r  r0

Z

r r0

Z

r r0

00 0 0 00 0 00 0 0 0 2h0 k00 M fM þ 2h kM fM  k1 f1  kM f1 fM ðg  r0 Þ dg; e þ k0M fM0 ðg  r0 Þ 00 0 00 0 0 0 k00 k1 fk1  kk fk  kM fk fM ðg  r0 Þ dg; e þ k0M fM0 ðg  r0 Þ

k ¼ 2; . . . ; M:

It can be seen that all the free terms are uniformly bounded and all the kernels are uniformly integrable. For instance, we have   Z r k0k þ k0M fM0 1 ðeÞ Kk ðg; rÞ dg ¼ 1  logð1 þ Y Þ ; Y k0M fM0 r0 with Y ¼ k0M fM0

r  r0 ; e

A. Bertuzzi et al. / Mathematical Biosciences 177&178 (2002) 103–125

121

the function in brackets being increasing from 0 to 1 when Y varies from 0 to þ1. In this way, we ðeÞ can establish that around r0 the ratios .k are bounded uniformly in e, so that the derivatives are also uniformly bounded. Hence, the family of regularized solutions is compact and there are uniformly convergent sequences. Now we show that there is in fact one unique limit, thus concluding the proof of the theorem. It P ðeÞ ðeÞ ðeÞ f is enough to show the result for fk , k ¼ 1; . . . ; M, because M k¼1 k ðrÞ þ fQ ðrÞ ¼ 1. Take any ðe1 Þ ðe2 Þ pair ðe1 ; e2 Þ and consider the differences dk ¼ fk  fk , k ¼ 1; . . . ; M. Our aim is to prove that for any r > 0 we have kdk k < r (k  k will denote the sup–norm on ½r0 ; R), provided that je1  e2 j is sufficiently small. The functions dk ðrÞ satisfy the system    Z r  ðe2 Þ Z 1 dd1 1 r df1 ðe1 Þ zkM fM dz þ e1 zkM dM ðzÞ dz þ e1  e2 þ r r0 r r0 dr dr ðe1 Þ ðe2 Þ ¼ 2hkM dM  k1 d1  kM dM f1  kM fM d1 ; ðB:11Þ 

   ðe2 Þ Z Z 1 r ddk 1 r dfk ðe1 Þ zkM fM dz þ e1 zkM dM ðzÞ dz þ e1  e2 þ r r0 r r0 dr dr ðe1 Þ ðe2 Þ ¼ kk1 dk1  kk dk  kM dM fk  kM fM dk ; k ¼ 2; . . . ; M:

ðB:12Þ

For k ¼ M we may perform a formal integration, starting from a point r1 > r0 and obtaining the expression    Z r  Z r Z 1 z 0  bM;2 ðzÞ e1  e2 þ bM;1 ðzÞ dz þ z kM dM ðz0 Þ dz0 dM ðrÞ ¼ dM ðr1 Þ exp  z r0 r1    Z r r1 þ bM;3 ðzÞdM1 ðzÞ exp  ðB:13Þ bM;1 ðz0 Þ dz0 dz; z

where ðe Þ

ðe Þ

kM ð1 þ fM 1 þ fM 2 Þ > 0; bM;1 ðrÞ ¼ R r ðe1 Þ 1 zk f dz þ e M 1 M r r0 ðe Þ

bM;2 ðrÞ ¼ R r 1 r

r0

bM;3 ðrÞ ¼ R r 1 r

r0

dfM 2 =dr ðe Þ

zkM fM 1 dz þ e1 kM1 ðe Þ

zkM fM 1 dz þ e1

;

> 0:

We can derive from (B.13) a Gronwall inequality for dM ðrÞ of the form Z r Z r jdM ðrÞj 6 jdM ðr1 Þj þ AM ðr; zÞjdM1 ðzÞj dz þ BM ðr; zÞje1  e2 j dz r1 r1 Z r Z z þ BM ðr; zÞ kM ðz0 ÞjdM ðz0 Þj dz0 dz r1

with

r0

ðB:14Þ

122

A. Bertuzzi et al. / Mathematical Biosciences 177&178 (2002) 103–125

 AM ðr; zÞ ¼ bM;3 ðzÞ exp 

Z

r

 bM;1 ðz Þ dz ; 0

0

ðB:15Þ

z ðe Þ

BM ðr; zÞ ¼ R z 1 z

r0

kdfM 2 =drk ðe Þ

z0 kM fM 1 dz0 þ e1

 exp 

Z

r 0

0



bM;1 ðz Þ dz :

ðB:16Þ

z

Interchanging the integration order in the last term of (B.14) we can rewrite it as # Z r "Z r BM ðr; zÞ dz kM ðz0 ÞjdM ðz0 Þj dz0 ; r0

cðz0 Þ

where cðz0 Þ ¼ max½z0 ; r1 . Therefore, the kernel of the Gronwall inequality is Z r 0 0 KM ðr; z Þ ¼ kM ðz Þ BM ðr; zÞ dz:

ðB:17Þ

cðz0 Þ

To guarantee that the kernel KM ðr; z0 Þ is bounded irrespective of the choice of r1 , e1 , e2 , we have to ðeÞ see that it remains bounded as r1 tends to r0 and e1 tends to zero. Since the derivative of fM is ðeÞ uniformly bounded, two positive constants h and l1 will exist such that kM ðrÞfM ðrÞ > l1 for r 2 ½r0 ; r0 þ h and irrespective of e. Then, for z < r0 þ h, we have  Z r  ðe Þ kdfM 2 =drk l2 0 BM ðr; zÞ 6 r0 exp  dz 0 l ðz  r0 Þ þ e1 z l3 ðz  r0 Þ þ e1 r 1 " #l2 =l3 ðe Þ kdfM 2 =drk l3 ðz  r0 Þ þ e1 ¼ r0 ; ðB:18Þ l ðz  r0 Þ þ e1 l3 ðr  r0 Þ þ e1 r 1 where l2 ¼ min kM and l3 ¼ max kM over the interval ½r0 ; R. Formula (B.18) shows clearly that when r1 tends to r0 and e1 tends to zero, the integral of BM with respect to z remains bounded since the integral from r1 to r0 þ h is bounded. Therefore, the kernel KM ðr; z0 Þ is bounded irrespective of the choice of r1 , e1 , e2 . We can make a similar calculation for AM and we arrive at the following conclusion: kdM kr 6 CðrÞðjdM ðr1 Þj þ je1  e2 j þ kdM1 kr Þ;

ðB:19Þ

where k  kr denotes the sup up to r, and CðrÞ is an increasing function independent of r1 , e1 , e2 . Using (B.19) in (B.12) for k ¼ M  1 we arrive at a similar result for dM1 and, proceeding recursively to the last step, we obtain from (B.11) a Gronwall inequality for d1 with the free terms containing only e1  e2 and dk ðr1 Þ, k ¼ 1; . . . ; M. The resulting estimate is ! M X kd1 kr 6 CðrÞ je1  e2 j þ jdk ðr1 Þj ; ðB:20Þ k¼1

which propagates to d2 ; . . . ; dM , going once more through (B.12) in the reverse order. ðeÞ ðe Þ ðe Þ ðe Þ Owing to the uniform Lipschitz continuity of the fk , since jfk 1 ðr1 Þ  fk 2 ðr1 Þj 6 jfk 1 ðr1 Þ  ðe2 Þ fk ðr0 Þj þ jfk ðr1 Þ  fk ðr0 Þj, we can make the dk ðr1 Þ as small as desired independently of e1 , e2 , by

A. Bertuzzi et al. / Mathematical Biosciences 177&178 (2002) 103–125

123

taking r1 sufficiently close to r0 . Thus, for any r > 0, we can choose je1  e2 j in such a way that kdk k < r, k ¼ 1; . . . ; M.  Appendix C The following proposition shows that, at the stationary state, the same fractions in S, G2 and M phases as a function of the radial distance can be obtained with different pairs of the functions hðrÞ and kG1 ðrÞ, when kS , kG2 and kM are constant with r. Proposition 1. Given kG1 ðrÞ, hðrÞ, and kS , kG2 and kM constant, let fk ðrÞ, k ¼ 1; . . . ; M, be the solution of Eqs. (7), (8) and (10). For every constant h1 , 0:5 < h1 6 1, a function k~G1 ðrÞ can be found such that Eqs. (7), (8) and (10), with ðC:1Þ h~ðrÞ ¼ h1 ; k~S ¼ kS ; k~G2 ¼ kG2 ; k~M ¼ kM ; have the solution f~k ðrÞ, k ¼ 1; . . . ; M, with the following property: f~k ðrÞ ¼ fk ðrÞ ;

k ¼ m1 þ 1; . . . ; M:

Proof. In view of (C.1), the fractions f~k ðrÞ will satisfy the equations 1 d ðr~ uf~1 Þ ¼ 2h1 kM f~M  k~G1 f~1 ; r dr 1 d ðr~ uf~k Þ ¼ k~G1 f~k1  k~G1 f~k ; r dr

k ¼ 2; . . . ; m1 ;

ðC:2Þ

ðC:3Þ ðC:4Þ

1 d ðr~ uf~m1 þ1 Þ ¼ k~G1 f~m1  km1 þ1 f~m1 þ1 ; r dr

ðC:5Þ

1 d ðr~ uf~k Þ ¼ kk1 f~k1  kk f~k ; r dr

ðC:6Þ

with r~ uðrÞ ¼ kM

Z

k ¼ m1 þ 2; . . . ; M;

r

zf~M dz:

ðC:7Þ

r0

For r ¼ r0 , the values f~k0 of the solution f~k of Eqs. (C.3)–(C.7) will be given by f~10 ¼ f~k0 ¼

2h1 kM k~0G1 þ kM f~M0 k~0G1 k~0G1 þ kM f~M0

f~m01 þ1 ¼

f~M0 ; 0 ; f~k1

ðC:8Þ

k ¼ 2; . . . ; m1 ;

k~0G1 f~m0 ; km1 þ1 þ kM f~M0 1

ðC:9Þ

ðC:10Þ

124

A. Bertuzzi et al. / Mathematical Biosciences 177&178 (2002) 103–125

kk1 0 ; k ¼ m1 þ 2; . . . ; M; ðC:11Þ f~k1 kk þ kM f~M0 where k~0G1 ¼ k~G1 ðr0 Þ. We will search the value k~0G1 that yields f~M0 ¼ fM0 . Expressing f~M0 by means of (C.8)–(C.11), we obtain the equation similar to (18) ! !m1 M Y kk þ kM f~M0 k~0G1 þ kM f~M0 ¼ 2h1 ðC:12Þ 0 ~ kk k f~k0 ¼

k¼m1 þ1

G1

and, for fM0 , we have the equation ! !m1 M Y kk þ kM fM0 k0G1 þ kM fM0 ¼ 2h0 ; 0 k k k k¼m1 þ1 G1

ðC:13Þ

where k0G1 ¼ kG1 ðr0 Þ and h0 ¼ hðr0 Þ. From Eqs. (C.12) and (C.13), by setting f~M0 ¼ fM0 , we obtain  1 1=m1 0 h kG1 þ kM fM0 k~0G1 þ kM fM0 ¼ ; ðC:14Þ h0 k0G1 k~0G1 with h1 =h0 2 ð0:5; 2Þ, from which k~0G1 can be found, resulting k~0G1 < k0G1 for h1 =h0 > 1. Since f~M0 ¼ fM0 , from Eq. (C.11) we have f~k0 ¼ fk0 , k ¼ m1 þ 1; . . . ; M  1, so property (C.2) holds for r ¼ r0 , and from (C.10) we have k~0G1 f~m01 ¼ k0G1 fm01 . For r > r0 , we note that if k~G1 ðrÞf~m1 ðrÞ ¼ kG1 ðrÞfm1 ðrÞ;

ðC:15Þ

~ uðrÞ ¼ ruðrÞ ¼ then R rEqs. (C.5)–(C.7) have the solution fk ðrÞ ¼ fk ðrÞ, k ¼ m1 þ 1; . . . ; M and thus r~ ~ ~ kM r0 zfM dz. Setting kG1 ðrÞ ¼ kG1 ðrÞfm1 ðrÞ=fm1 ðrÞ and substituting uðrÞ to u~ðrÞ into Eqs. (C.3) and (C.4), we obtain the ODE system (with degeneracy at r ¼ r0 ) df~ f~ ðC:16Þ u 1 ¼ 2h1 kM fM  kG1 fm1 1  kM fM f~1 ; dr f~m 1

u

df~k f~ f~ ¼ kG1 fm1 k1  kG1 fm1 k  kM fM f~k ; dr f~m1 f~m1

k ¼ 2; . . . ; m1 :

ðC:17Þ

The solution of the above system, whose existence and uniqueness can be shown following the line of Theorem 1, will give f~m1 ðrÞ from which k~G1 ðrÞ is found through (C.15). It can be easily seen that the values f~k0 , k ¼ 1; . . . ; m1 , obtained from (C.8)–(C.11) when k~0G1 is found from (C.14), are the values at r0 that can be deduced from (C.16) and (C.17). Thus the value k~G1 ðr0 Þ obtained from the solution of (C.16) and (C.17) is the value found from (C.14).  Note that, taking h1 ¼ 1, the above result states that cell arrest into the quiescent compartment can be fully mimicked, as far as the phase fractions in S, G2 and M phases are concerned, by a suitably varying rate of progression through G1 in the absence of any cell arrest. In a symmetric way, it seems possible to construct a function h~ðrÞ that substitutes for a decreasing kG1 ðrÞ. Thus we may advance the following conjecture:

A. Bertuzzi et al. / Mathematical Biosciences 177&178 (2002) 103–125

125

Conjecture 1. Given kG1 ðrÞ, hðrÞ regular enough and decreasing, and kS , kG2 and kM constant, let fk ðrÞ, k ¼ 1; . . . ; M, be the solution of Eqs. (7), (8) and (10). A function h~ðrÞ 2 ð0; 1, 0:5 < h~ðr0 Þ 6 1, and a positive constant k1G1 can be found such that Eqs. (7), (8) and (10), with k~S ¼ kS ; k~G2 ¼ kG2 ; k~M ¼ kM ; have the solution f~k ðrÞ, k ¼ 1; . . . ; M, with the following property k~G1 ¼ k1G1 ;

f~k ðrÞ ¼ fk ðrÞ;

k ¼ m1 þ 1; . . . ; M:

References [1] J. Folkman, Tumor angiogenesis, Adv. Cancer Res. 43 (1985) 175. [2] W. Mueller-Klieser, Multicellular spheroids. A review on cellular aggregates in cancer research, J. Cancer Res. Clin. Oncol. 113 (1987) 101. [3] J.A. Adam, General aspects of modeling tumor growth and immune response, in: J.A. Adam, N. Bellomo (Eds.), A Survey of Models for Tumor-Immune System Dynamics, Birkh€ auser, Boston, MA, 1997, p. 187. [4] J. Holash, P.C. Maisonpierre, D. Compton, P. Boland, C.R. Alexander, D. Zagzag, G.D. Yancopoulos, S.J. Wiegand, Vessel cooption, regression, and growth in tumors mediated by angiopoietins and VEGF, Science 284 (1999) 1994. [5] G.D. Yancopoulos, S. Davis, N.W. Gale, J.S. Rudge, S.J. Wiegand, J. Holash, Vascular-specific growth factors and blood vessel formation, Nature 407 (2000) 242. [6] I.F. Tannock, The relation between cell proliferation and the vascular system in a transplanted mouse mammary tumour, Br. J. Cancer 22 (1968) 258. [7] D.G. Hirst, J. Denekamp, Tumour cell proliferation in relation to the vasculature, Cell Tissue Kinet. 12 (1979) 31. [8] J.V. Moore, H.A. Hopkins, W.B. Looney, Tumour-cord parameters in two rat hepatomas that differ in their radiobiological oxygenation status, Radiat. Environ. Biophys. 23 (1984) 213. [9] J.V. Moore, P.S. Hasleton, C.H. Buckley, Tumour cords in 52 human bronchial and cervical squamous cell carcinomas: inferences for their cellular kinetics and radiobiology, Br. J. Cancer 51 (1985) 407. [10] K.H. Falkvoll, The relationship between changes in tumour volume, tumour cell density and parenchymal cord radius in a human melanoma xenograft exposed to single dose irradiation, Acta Oncol. 29 (1990) 935. [11] D.G. Hirst, V.K. Hirst, B. Joiner, V. Prise, K.M. Shaffi, Changes in tumour morphology with alterations in oxygen availability: further evidence for oxygen as a limiting substrate, Br. J. Cancer 64 (1991) 54. [12] A. Bertuzzi, A. Fasano, A. Gandolfi, A mathematical model for the growth of tumor cords incorporating the dynamics of a nutrient, in: N. Kenmochi (Ed.), Free Boundary Problems: Theory and Applications II, Gakkotosho, Tokyo, 2000, p. 31. [13] A. Bertuzzi, A. Gandolfi, Cell kinetics in a tumour cord, J. Theor. Biol. 204 (2000) 587. [14] D.G. Kendall, On the role of variable cell generation time in the development of a stochastic birth process, Biometrika 35 (1948) 316. [15] M. Takahashi, Theoretical basis for cell cycle analysis. I. Labelled mitosis wave method, J. Theor. Biol. 13 (1966) 202. [16] M. Takahashi, Theoretical basis for cell cycle analysis. II. Further studies on labelled mitosis wave method, J. Theor. Biol. 18 (1968) 195. [17] G.G. Steel, Growth Kinetics of Tumours: Cell Population Kinetics in Relation to the Growth and Treatment of Cancer, Clarendon, Oxford, 1977. [18] J.M. Ortega, W.C. Rheinboldt, Iterative Solution of Nonlinear Equations in Several Variables, Academic Press, New York, 1970. [19] A.Yu. Yakovlev, N.M. Yanev, Transient Processes in Cell Proliferation Kinetics, Lecture Notes in Biomathematics, vol. 82, Springer, Berlin, 1989. [20] P. Hartman, Ordinary Differential Equations, Wiley, New York, 1964.