Modeling insight into spontaneous regression of tumors

Mathematical Biosciences 155 (1999) 45±60

Modeling insight into spontaneous regression of tumors Andrej Yakovlev *, Kenneth Boucher, James DiSario Huntsman Cancer Institute, Department of Oncological Sciences, University of Utah, 546 Chipeta Way, Suite 1100, Salt Lake City, UT 84108, USA Received 20 January 1998; received in revised form 18 September 1998; accepted 16 October 1998

Abstract The phenomenon of spontaneous regression of benign and malignant tumors is well documented in the literature and is commonly attributed to the induction of apoptosis or activation of the immune system. We attempt at evaluating the role of random e€ects in this phenomenon. To this end, we consider a stochastic model of tumor growth which is descriptive of the fact that tumors are inherently prone to spontaneous regression due to the random nature of their development. The model describes a population of actively proliferating cells which may give rise to di€erentiated cells. The process of cell di€erentiation is irreversible and terminates in cell death. We formulate the model in terms of temporally inhomogeneous Markov branching processes with two types of cells so that the expected total number of neoplastic cells is consistent with the observed mean growth kinetics. Within the framework of this model, the extinction probability for proliferating cells tends to one as time tends to in®nity. Given the event of nonextinction, the distribution of tumor size is asymptotically exponential. The limiting conditional distribution of tumor size is in good agreement with epidemiologic data on advanced lung cancer. Ó 1999 Elsevier Science Inc. All rights reserved. Keywords: Tumor growth; Spontaneous regression; Stochastic modelling; Gompertz birth±death process; Asymptotic results

*

Corresponding author. Tel.: +1-801 585 9543; fax: +1-801 585 5357; e-mail: [email protected]

0025-5564/99/$ ± see front matter Ó 1999 Elsevier Science Inc. All rights reserved. PII: S 0 0 2 5 - 5 5 6 4 ( 9 8 ) 1 0 0 5 2 - 4

46

A. Yakovlev et al. / Mathematical Biosciences 155 (1999) 45±60

1. Introduction It is a well-known fact that cell populations constituting benign tumors, such as colon polyps, not infrequently become extinct [1]. In particular, Knoernschild [2] followed 213 patients with colon polyps for a period of 3±5 years. Some polyps were marked by a small tattoo placed near the base of the polyps. Knoernschild reported 18% of the marked polyps entirely disappeared within the period of observation. By contrast, spontaneous regression of malignant tumors is a rare event because death of the tumor host precludes this event from occurring. Nevertheless malignant tumors do regress occasionally for no apparent reason as evidenced by many clinical observations. The best known example for this type of behavior is neuroblastoma that has the highest rate of spontaneous regression among human malignant tumors [3]. At the same time, neuroblastoma has one of the poorest prognosis when occurring as disseminated neoplasm, accounting for at least 15% of cancer-related deaths in children [4]. Brodeour et al. [5] classi®ed neuroblastomas into three subsets with di€erent biological features and clinical course. One of the types of neuroblastoma is characterized by a lower malignant potential than the other two. Neuroblastomas of this type are prone to spontaneous regression and have a very favorable clinical outcome. The three tumor types appear to be genetically distinct; there is no evidence that one type may evolve into another. A limited therapeutic intervention may sometimes cause spontaneous regression of neuroblastoma. One such case was recently reported by Kullendor€ and Stromblad [6]. In two infants, paresis of the legs due to spinal cord compression by a dumbbell neuroblastoma disappeared within two weeks after surgical excision of only the paraspinal tumor mass. This e€ect was attributed to prompt spontaneous regression of the intraspinal component of dumbbell neuroblastoma. Spontaneous regression of other malignant tumors is also well documented. As one example, a case of biopsy-proven hepatocellular carcinoma, which was considered to be unresectable at initial laparotomy but subsequently regressed without any speci®c treatment, was reported by Grossmann et al. [7]. The authors indicated nine other case reports of spontaneous regression of hepatocellular carcinoma published in the English literature. Similar cases have been reported for intracranial germinoma [8], meningioma [9], pulmonary metastases of renal cell carcinoma [10], primary small cell bronchial carcinoma [11], malignant melanoma and basal cell carcinoma [12], lung metastases from osteosarcoma [13], and hepatic metastases in carcinoid heart disease [14]. This intriguing phenomenon is most commonly explained by the induction of apoptosis or immune responses (see Section 5). To the best of our knowledge, none of the authors have taken account of the stochastic nature of tumor growth in their discussion of possible mechanisms underlying spontaneous regression of tumors. Clearly, one has to rely on mathematical modeling as the

A. Yakovlev et al. / Mathematical Biosciences 155 (1999) 45±60

47

only way to evaluate the role of random e€ects in this process. In this paper, a tumor is modeled as a stochastic system, and we explore the properties of such a system, which are coherent with the phenomenon of spontaneous tumor regression.

2. A stochastic model of tumor growth Based on many experimental and clinical observations the common belief is that tumors initially grow rapidly (exponentially or superexponentially) but increasingly slow their growth as they become larger, so that the mean tumor size tends to a constant level as time increases. Consistent with this pattern is the Gompertz growth kinetics of the expected number of tumor cells. Luebeck and Moolgavkar [15] developed a stochastic model of clonal expansion in carcinogenesis proceeding from the Gompertz birth-death process proposed by Tan [16]. The authors applied the model to study the clonal growth of enzyme altered foci induced by a chemical carcinogen (N-nitrosomorpholine) in the liver of rats. These lesions, although premalignant in character, show a growth pattern which is much better described by the Gompertz mean growth model than by the exponential one. Thus, a retardation of growth rate can be seen even in early stages of carcinogenesis. The analysis by Luebeck and Moolgavkar [15] also suggests that the clonal growth of initiated cells is signi®cantly accelerated (superexponential) rather than exponential for nontoxic doses of N-nitrosomorpholine. In this paper, we explore similar ideas as applied to the growth of malignant tumors. There are other, probably more natural, ways to model the general pattern of tumor growth, but we decide on the Gompertztype stochastic process in order to retain model tractability (see Ref. [15] for discussion). Our model is based on the following assumptions: (1) The growth of a tumor begins with a single clonogenic malignant cell at time s. Generally speaking, this moment of time is a random variable and its distribution can be speci®ed, for example, by the two-stage model of carcinogenesis, known as the Moolgavkar±Venzon±Knudson model [17,18], but this will not be necessary for purposes of this paper. In what follows s will be thought of as a nonrandom moment of time. (2) In the small interval of time …t; t ‡ DtŠ the initiator malignant cell divides into two cells of the same type with probability a…t†Dt ‡ o…Dt†, and it transforms (di€erentiates) into one di€erentiated cell with probability b…t†Dt ‡ o…Dt†. The di€erentiation pathway is irreversible and results in cell death. A di€erentiated cell dies with probability k…t†Dt ‡ o…Dt†. (3) The usual independence hypotheses for the birth-and-death process are adopted.

48

A. Yakovlev et al. / Mathematical Biosciences 155 (1999) 45±60

The stochastic process thus de®ned is a particular case of the Markov branching process with two types of cells. For this process, the probability generating function for the cell progeny is given by h1 …s1 ; s2 ; t† ˆ ps21 ‡ …1 ÿ p†s2 ;

h2 …s1 ; s2 ; t† ˆ 1;

…1†

where p ˆ p…t† ˆ

a…t† : a…t† ‡ b…t†

We also introduce the notation h1 …t† ˆ a…t† ‡ b…t†;

h2 …t† ˆ k…t†:

…2†

Now we have to specify the rates a and b as functions of time. In order to describe the Gompertz shape of the mean growth curve, Luebeck and Moolgavkar [15] speci®ed the `net' proliferation rate, d ˆ a ÿ b, as d…t ÿ s† ˆ beÿa…tÿs† ;

t P s P 0;

…3†

while the ratio c ˆ b=a …c < 1† was kept constant. Under these conditions both a and b, given by d cd ; bˆ ; 1ÿc 1ÿc tend to zero as t ! 1. The model will become more ¯exible and realistic if we just require that a and b tend to some ®nite limit, not necessarily equalling zero, as time tends to in®nity. Speci®cally, it suces to assume that a…t† and b…t† are bounded measurable functions satisfying the following conditions: aˆ

Zt lim a…t† ˆ lim b…t† ˆ c;

t!1

t!1

lim

‰a…u† ÿ b…u†Š du ˆ A…s†;

t!1

…4†

s

where c is a non-negative constant and A…s† is a bounded function of s. It is natural to assume that A…s† takes on only positive values. A pertinent example is given by the following form of the rates a and b: a…t ÿ s† ˆ c ‡ b1 eÿa…tÿs† ;

b…t ÿ s† ˆ c ‡ b2 eÿa…tÿs† ;

…5†

where b1 , b2 and a are positive constants and b1 P b2 . This speci®cation of a and b satis®es condition (3), of course. In Section 4, it will become clear that, under the model considered here, the condition (4) are sucient for the extinction probability for tumor cells to tend to 1 as time increases, providing c > 0. Consider the two-dimensional stochastic process Z…t† ˆ …Z1 …t†; Z2 …t††, where Z1 …t† is the number of proliferating tumor cells and Z2 …t† is the number of di€erentiated cells at time t. Introduce the probability generating function U…s; s; t† ˆ …U1 …s; s; t†; U2 …s; s; t††;

…6†

A. Yakovlev et al. / Mathematical Biosciences 155 (1999) 45±60

with the components X Ui …s; s; t† ˆ PrfZ…t† ˆ kjZ…s† ˆ ei gsk ;

49

i ˆ 1; 2; jsj 6 1;

…7†

k

where sk ˆ sk11 sk22 , e1 ˆ …1; 0†, e2 ˆ …0; 1†, and the summation in Eq. (7) is over the set of all points in R2 with non-negative integer coordinates. Suppose that a cell of type i, i ˆ 1; 2, existing at time t undergoes transformation within the interval …t; t ‡ DtŠ with probability hi …t†Dt ‡ o…Dt†, and let hi …s; t† be the generating function of the numbers of cells of the two di€erent types born at transformation. Given in Ref. [19], pp. 113±114, is the following system of backward equations for the generating functions (7): oUi …s; s; t† ˆ ÿhi …s†‰hi …U…s; s; t†; s† ÿ Ui …s; s; t†Š; i ˆ 1; 2; …8† os with the initial condition Ui …s; t ÿ 0; t† ˆ si . With hi and hi given by Eqs. (1) and (2), respectively, the Eq. (8) assume the form oU1 …s; s; t† ˆ ÿa…s†U21 …s; s; t† ÿ b…s†U2 …s; s; t† ‡ ‰a…s† ‡ b…s†ŠU1 …s; s; t†; os oU2 …s; s; t† ˆ ÿk…s†‰1 ÿ U2 …s; s; t†Š: os

…9†

These equations provide the basis for all considerations that follow. 3. The moments We need explicit formulas for the ®rst and the second moments of the processes Z1 (t), Z2 (t) and their covariance to be employed in the proof of the main theorem in Section 4. The existence of the moments follows from general results of the theory of branching processes [20]. The marginal probability generating function of the process Z1 (t) can be obtained from Eq. (9) by setting s1 ˆ s, s2 ˆ 1. Introduce the notation: /1 …s; s; t† ˆ U1 …s; 1; s; t† and /2 …s; s; t† ˆ U2 …s; 1; s; t†. Then it follows from Eq. (9) that /1 ˆ /1 …s; s; t† and /2 ˆ /2 …s; s; t† satisfy the following equation: o/1 os

ˆ ÿa/21 ÿ b ‡ …a ‡ b†/1 ;

/2 ˆ 1:

…10†

In a similar manner, we set s1 ˆ 1, s2 ˆ s to obtain ow1 ow2 ˆ ÿaw21 ÿ bw2 ‡ …a ‡ b†w1 ; ˆ ÿk‰1 ÿ w2 Š: …11† os os where w1 ˆ w1 …s; s; t† ˆ U1 …1; s; s; t† and w2 ˆ w2 …s; s; t† ˆ U2 …1; s; s; t†. To ®nd the expected number of proliferating (type-1) cells we di€erentiate Eq. (10) with respect to s and then set s ˆ 1. The resultant di€erential equation for M11 …s; t† ˆ EfZ1 …t†jZ1 …s† ˆ 1g is solved easily yielding the expression

50

A. Yakovlev et al. / Mathematical Biosciences 155 (1999) 45±60

M11 …s; t† ˆ exp

8 t
s

9 =

‰a…u† ÿ b…u†Š du ; ;

…12†

which is well-known in the theory of nonhomogeneous birth±death processes [19]. Di€erentiating Eq. (11) with respect to s and setting s ˆ 1 we have the equation for the expected number of di€erentiated (type-2) cells M12 …s; t† ˆ EfZ2 …t†jZ1 …s† ˆ 1g: 8s 9
M12 …t; t† ˆ 0;

…13†

which has the solution 8x 9 Zt Zt
s

…14†

x

From Eqs. (4) and (12) it follows that for any ®nite s, lim M11 …s; t† ˆ eA…s† :

…15†

t!1

Suppose that a ®nite limit k0 ˆ limt!1 k…t† exists and k0 > 0, then recalling Eq. (4) again and applying l'H opital's rule to Eq. (14) we get c lim M12 …s; t† ˆ eA…s† : …16† t!1 k0 Thus we can state that M11 …s; t† k0 ! ; t ! 1; …17† k0 ‡ c M11 …s; t† ‡ M12 …s; t† for any ®nite s. We will further explore this convergence in the next section. Representative shapes of the functions M11 …0; t† and M12 …0; t† are shown in Fig. 1. This numerical example uses the rates a…t† and b…t† as speci®ed in formulas (5). Now we are in a position to derive formulas for the variances of the processes Z1 and Z2 . It is a matter of direct veri®cation to prove that the second equation in Eq. (11) has the following solution: Zt w2 …s; s; t† ˆ 1 ‡ …s ÿ 1† exp ÿ

k…x† dx; s

so that Eq. (11) can be rewritten as ow1 ˆ ÿaw21 ÿ bw2 ‡ …a ‡ b†w1 ; os

Zt w2 ˆ 1 ‡ …s ÿ 1† exp ÿ

k…x† dx: …18† s

A. Yakovlev et al. / Mathematical Biosciences 155 (1999) 45±60

51

Fig. 1. Typical behavior of some characteristics of the model. (A) The function M11 …t†  103 computed for a…t† and b…t† speci®ed by formula (5) with b1 ˆ 20, b2 ˆ 2, c ˆ 100; k ˆ 1; (B) M12 …t†  105 computed for the same set of parameters; (C) The extinction probability p0 …t† computed using formula (31) and the same set of parameters.

First we need to ®nd the factorial moments 2

D11 …s; t† ˆ Ef‰Z1 …t†Š jZ1 …s† ˆ 1g ÿ M11 …s; t†; 2

D12 …s; t† ˆ Ef‰Z2 …t†Š jZ1 …s† ˆ 1g ÿ M12 …s; t†:

52

A. Yakovlev et al. / Mathematical Biosciences 155 (1999) 45±60

Di€erentiating Eqs. (10) and (18) twice with respect to s and setting s ˆ 1, we have the di€erential equations

and

oD11 …s; t† 2 ˆ ÿ2a…s†M11 …s; t† ÿ ‰a…s† ÿ b…s†ŠD11 …s; t† os

…19†

oD12 …s; t† 2 ˆ ÿ2a…s†M12 …s; t† ÿ ‰a…s† ÿ b…s†ŠD12 …s; t†; os

…20†

with the initial conditions: D11 …t; t† ˆ D12 …t; t† ˆ 0. These equations can be solved using the integrating factor exp fq…x; t†g with Zt ‰b…u† ÿ a…u†Š du:

q…x; t† ˆ

…21†

x

The resulting solutions are of the form D11 …s; t† ˆ 2eÿq…s;t†

Zt

2 a…x†M11 …x; t†eq…x;t† dx;

…22†

2 a…x†M12 …x; t†eq…x;t† dx;

…23†

s

D12 …s; t† ˆ 2e

ÿq…s;t†

Zt s

while the corresponding variances are given by 2 …s; t†; D11 …s; t† ˆ VarfZ1 …t†jZ1 …s† ˆ 1g ˆ D11 …s; t† ‡ M11 …s; t† ÿ M11

…24†

2 …s; t†: D12 …s; t† ˆ VarfZ2 …t†jZ1 …s† ˆ 1g ˆ D12 …s; t† ‡ M12 …s; t† ÿ M12

…25†

In like manner we can ®nd K…s; t† ˆ EfZ1 …t†Z2 …t†jZ1 …s† ˆ 1; Z2 …s† ˆ 0g from the equations for the joint probability generating functions oU1 …s1 ; s2 ; s; t† ˆ ÿa…s†U21 …s1 ; s2 ; s; t† ÿ b…s†U2 …s1 ; s2 ; s; t† os ‡ ‰a…s† ‡ b…s†ŠU1 …s1 ; s2 ; s; t†; Zt U2 …s1 ; s2 ; s; t† ˆ 1 ‡ …s2 ÿ 1† exp ÿ k…x† dx;

…26†

s

see formulas (9) and (18). Di€erentiating Eq. (26) once with respect to s1 and once with respect to s2 , and setting s1 ˆ s2 ˆ 1, we derive the equation for K…s; t†: oK…s; t† ˆ ÿ2a…s†M11 …s; t†M12 …s; t† ÿ ‰a…s† ÿ b…s†ŠK…s; t†: os

…27†

A. Yakovlev et al. / Mathematical Biosciences 155 (1999) 45±60

53

Using the same integrating factor exp fq…x; t†g and taking the condition: K…t; t† ˆ 0 into account, we obtain from Eq. (27) K…s; t† ˆ 2eÿq…s;t†

Zt

a…x†M11 …x; t†M12 …x; t†eq…x;t† dx;

…28†

s

and the covariance results from the formula: CovfZ1 …t†; Z2 …t†jZ1 …s† ˆ 1; Z2 …s† ˆ 0g ˆ K…s; t† ÿ M11 …s; t†M12 …s; t†: Eqs. (19), (20) and (27) can also be solved by the variation of constant method. 4. Asymptotic properties In this section, we con®ne our consideration to the case s ˆ 0 without loss of generality. In reference to the function A…s†, introduced in Eq. (4), we will use the notation: A ˆ A…0†. For simplicity we also assume that the death rate for di€erentiated cells, k…t†, is constant in time, i.e. k…t† ˆ k, although the results that follow are readily extendable to the case where k…t† tends to a ®nite (strictly positive) limit as time tends to in®nity. The equation for /1 …s; t† ˆ EfsZ1 …t† jZ1 …0† ˆ 1g in Eq. (10) was solved by Kendall [21]. The solution is /1 …s; t† ˆ 1 ÿ

sÿ1 ; …s ÿ 1†G…t† ÿ g…t†

…29†

where Zt g…t† ˆ exp fq…t†g;

G…t† ˆ

a…u†eq…u† du;

…30†

0

and q…t† ˆ q…0; t† is de®ned in Eq. (21). Unfortunately, it has not proved feasible to obtain an analytical solution of the system of Eq. (11). It immediately follows from Eq. (29) that the extinction probability p0 …t† ˆ Pr fZ1 …t† ˆ 0g is given by p0 …t† ˆ 1 ÿ

1 1 ˆ1ÿ : Rt g…t† ‡ G…t† eq…t† ‡ 0 a…u†eq…u† du

…31†

Recall the conditions (4) imposed on the rates a and b. Suppose c > 0, then it is clear from Eq. (31) that p0 …t† ! 1 as t ! 1, i.e., the population of tumor cells eventually becomes extinct with probability 1 within the framework of our model. This is not particularly surprising since by virtue of conditions (4) the process Z1 …t† is asymptotically a critical branching stochastic process. The rate of convergence of p0 …t† is quite slow as can be seen from

54

A. Yakovlev et al. / Mathematical Biosciences 155 (1999) 45±60

eA ; t ! 1: …32† ct A typical behavior of the extinction probability as a function of time is shown in Fig. 1(C). Suppose Z1 …0† ˆ 1 and introduce the processes 1 ÿ p0 …t† 

W1 …t† ˆ ‰1 ÿ p0 …t†ŠZ1 …t†

W2 …t† ˆ ‰1 ÿ p0 …t†ŠZ2 …t†:

…33†

Note that EfW1 …t†g ! 0 but EfW1 …t†jZ1 …t† > 0g ! exp …A† as t ! 1. Given the event fZ1 …t† > 0g, the process W1 …t† converges weakly to an exponentially distributed random variable W, while the process W2 …t† converges in the same sense to …c=k†W . To prove this assertion we need the following lemma by Hanin et al. [22]. Lemma 1. Let R…t; s; u† be a bounded measurable function de®ned for t P s P 0; u P 0 and such that R…t; s; u† ! a as t; s; u ! 1. Also, let q be a positive t!1 q…t ÿ u†=q…t† ˆ 1 for every u P 0, R 1 continuous function such that Rlim t and 0 q…u† du ˆ 1. Denote Q…t† ˆ 0 q…u†du. Then 1 Q…t†

Zt R…t; t ÿ u; u†q…u† du ! a;

t ! 1:

0

Now we can formulate the main result of this section. Theorem 1. If conditions (4) with c > 0 are met, then lim PrfW1 …t† > ujZ1 …t† > 0g ˆ eÿe

ÿA u

t!1

;

k ÿA u

lim PrfW2 …t† > ujZ1 …t† > 0g ˆ eÿce

t!1

;

u P 0;

…34†

u P 0:

…35†

Proof. Introduce the conditional generating function W…s; t† ˆ EfsZ1 …t† jZ1 …t† > 0g ˆ

/1 …s; t† ÿ p0 …t† : 1 ÿ p0 …t†

Using Eq. (29) we have sg…t† g…t† ÿ …s ÿ 1†G…t† and the Laplace±Stieltjes transform of the conditional distribution of Z1 …t† is given by W…s; t† ˆ

^ t† ˆ W…s;

eÿs g…t† ; g…t† ÿ …eÿs ÿ 1†G…t†

Rs P 0:

…36†

Recalling Eq. (33), we have the Laplace±Stieltjes transform of the distribution of W1 …t†,

A. Yakovlev et al. / Mathematical Biosciences 155 (1999) 45±60

^ t† ˆ X…s;

55

eÿ‰1ÿp0 …t†Šs g…t† : g…t† ÿ …eÿ‰1ÿp0 …t†Šs ÿ 1†G…t†

It is clear from Eqs. (30) and (31) that lim ‰1 ÿ p0 …t†ŠG…t† ˆ 1

t!1

and, therefore, ^ t† ˆ lim X…s;

t!1

eÿA ; eÿA ‡ s

t ! 1:

…37†

The right-hand side of Eq. (37) represents the Laplace transform of the exponential distribution density with parameter eÿA . This proves the asymptotics of W1 …t† given by Eq. (34). Let us prove Eq. (35). It follows from Eqs. (21) and (22) that the function D11 …t† ˆ D11 …0; t† can be represented as Zt ÿ2q…t† a…x†eq…x† dx; D11 …t† ˆ 2e 0

where q…t† ˆ q…0; t†. Therefore, it is easy to see that D11 …t†  2ceA t; Recalling Eq. (32) we get the limit ‰1 ÿ p0 …t†ŠEfZ12 …t†g ! 2e2A ;

t ! 1:

t ! 1:

It can be veri®ed that M11 …x; t† ! 1, M12 …x; t† ! c=k and exp fq…x; t†g ! 1 as x; t ! 1. Applying the lemma with q  1 to formulas (23) and (28) with s ˆ 0, we have D12 …t†  2…c3 =k2 †eA t; K…t†  2…c2 =k†eA t as t ! 1. Therefore, the following limits hold: c2 2A c e ; ‰1 ÿ p0 …t†ŠEfZ1 …t†Z2 …t†g ! 2 e2A ; 2 k k as t ! 1. Proceeding from Eq. (33) we have h  i2 c W1 …t† ÿ W2 …t† Z1 …t† > 0 06E k h  i2 c 2 Z1 …t† ÿ Z2 …t† Z1 …t† > 0 ˆ ‰1 ÿ p0 …t†Š E k h i2  c 6 ‰1 ÿ p0 …t†ŠE Z1 …t† ÿ Z2 …t† k c2 ˆ 2 ‰1 ÿ p0 …t†ŠEfZ12 …t†g ‡ ‰1 ÿ p0 …t†ŠEfZ22 …t†g k c c2 c2 c2 ÿ 2 ‰1 ÿ p0 …t†ŠEfZ1 …t†Z2 …t†g ! 2 2 e2A ‡ 2 2 e2A ÿ 4 2 e2A k k k k ˆ 0; ‰1 ÿ p0 …t†ŠEfZ22 …t†g ! 2

which completes the proof of the theorem.

56

A. Yakovlev et al. / Mathematical Biosciences 155 (1999) 45±60

Corollary 1. The sum W1 …t† ‡ W2 …t† converges in distribution (conditionally on Z1 …t† > 0) to U ˆ …c ‡ k=k†W . Indeed, h  c i2 W1 …t† ‡ W1 …t† Z1 …t† > 0 E W1 …t† ‡ W2 …t† ÿ k h  i2 c ˆ E W2 …t† ÿ W1 …t† Z1 …t† > 0 ! 0; t ! 1: k Corollary 2. Given Z1 …t† > 0, the process Y …t† ˆ ˆ

Z1 …t† ‡ Z2 …t† EfZ1 …t†jZ1 …t† > 0g ‡ EfZ2 …t†jZ1 …t† > 0g W1 …t† ‡ W2 …t† ‰1 ÿ p0 …t†ŠEfZ1 …t†jZ1 …t† > 0g ‡ ‰1 ÿ p0 …t†ŠEfZ2 …t†jZ1 …t† > 0g

converges weakly to a random variable V as t ! 1. The distribution of V is exponential with parameter 1. Obviously, we have ‰1 ÿ p0 …t†ŠEfZ1 …t†jZ1 …t† > 0g ! eA ;

t!1

…38†

and c ‰1 ÿ p0 …t†ŠEfZ2 …t†jZ1 …t† > 0g ! eA ; t ! 1: k Now we need to simplify our notation as follows:

…39†

X …t† ˆ W1 …t† ‡ W2 …t†; m…t† ˆ

1 : ‰1 ÿ p0 …t†Š‰EfZ1 …t†jZ1 …t† > 0g ‡ EfZ2 …t†jZ1 …t† > 0gŠ

By Corollary 1, the process X …t† tends in distribution to the random variable U as t tends to in®nity, and the distribution function FU …u† ˆ Pr fU 6 ug is continuous. It follows from a result of P olya (see Ref. [29], p. 86) that supjFX …t† …y† ÿ FU …y†j ! 0; y2R

t ! 1:

…40†

It also follows from Eqs. (38) and (39) that m(t) tends to m ˆ …k=c ‡ k†eÿA . It remains to show that for every x 2 R, FY …t† …x† ˆ Fm…t†X …t† …x† ! FmU …x† ˆ FV …x† as t ! 1. Indeed, we have the chain of inequalities

A. Yakovlev et al. / Mathematical Biosciences 155 (1999) 45±60

   x  x ÿ FU jFm…t†X …t† …x† ÿ FmU …x†j ˆ FX …t† m…t† m     x x ÿ FU 6 FX …t† m…t† m…t†    x  x ÿ FU ‡ FU m…t† m    x  x ÿ FU 6 supjFX …t† …y† ÿ FU …y†j ‡ FU m…t† m y2R

57

ˆ I1 …t† ‡ I2 …t†: But I1 …t† ! 0 as t ! 1 by virtue of Eq. (40), and I2 …t† ! 0 by continuity of the function FU …t†; this completes the proof of Corollary 2. Remark 1. We refer the reader to Sevast'yanov's paper [23] for asymptotic properties of one-dimensional critical Markov branching processes and to the paper by Jagers [24] for relevant results regarding an age-dependent, but temporally homogeneous, branching process with two types of cells. 5. Discussion The phenomenon of spontaneous regression of malignant tumors is most commonly ascribed to apoptosis or immune attacks induced by unknown causes. Some clinical reports are consistent with the idea that apoptosis may be a possible cause of spontaneous regression and remission among cancer patients. For example, Ozeki et al. [25] reported a female patient with a part of liver mass which, when surgically removed, was believed to be spontaneous complete necrosis of hepatocellular carcinoma. No viable cells were histologically detectable in the resected hepatic tissue. However, a recent study [26] of DNA fragmentation in the course of neuroblastoma progression provided evidence against this hypothesis. In those authors' opinion, it is unlikely that apoptosis is a basic mechanism of neuroblastoma regression. There is evidence in favor of the hypothesis that spontaneous regression of tumors may be immunologically mediated. Halliday et al. [12] observed a larger number of CD4+ T lymphocytes in®ltrated regressing primary melanomas and basal cell carcinomas than their non-regressing counterparts. The number of T cells positive for interleukin 2 receptor, but not for transferrin receptor, was increased as well, indicating that the in®ltrating T lymphocytes were in the early stage of their activation. It remains unknown why only some tumors are strongly attacked by the immune system and what factor is responsible for this selective e€ect.

58

A. Yakovlev et al. / Mathematical Biosciences 155 (1999) 45±60

The results of our work indicate that the role of random e€ects must not be underestimated when one seeks to explain why tumors regress spontaneously. The stochastic model considered here, quite simplistic as it is, provides a certain basis for exploration of such e€ects. The model is descriptive of the fact that tumors are inherently prone to spontaneous regression due to the random nature of their development. To evaluate this model adequacy we analyzed data on tumor sizes in 392 patients with advanced lung cancer. The patients were identi®ed through the Utah Cancer Registry. The estimated mean tumor size was 50.2 (in relative units). Since the exponential distribution suggested by Corollary 2 is conditional on the event fZ1 …t† > 0g, it appears to be more appropriate for describing the distribution of large tumor sizes. For this reason, the distribution and the corresponding data on tumor volume were truncated from the left at 40. Another reason to do so is that the process of tumor detection may have an impact on the shape of tumor size distribution [27,28]. Since the detection rate increases monotonically with tumor size, truncating the distribution of tumor size at detection would make this impact smaller. The truncated exponential distribution provides a good ®t to the observed tumor sizes (Fig. 2) as evidenced by the Chi-square goodness of ®t test yielding a signi®cance level of greater than 0.3. It is encouraging that our basic assumptions are not in con¯ict with real data of this type.

Fig. 2. Fitting the exponential distribution to data on tumor sizes in patients diagnosed with advanced lung cancer.

A. Yakovlev et al. / Mathematical Biosciences 155 (1999) 45±60

59

Acknowledgements We are very grateful to Dr J. Simone (University of Utah) for stimulating and fruitful discussions of biomedical implications of this work. We would like to thank Dr L. Hanin (Idaho State University) for some helpful suggestions. We are also indebted to the reviewers whose comments have led to substantial improvements in the manuscript. The research of A.Y.Y. and K.B. is supported in part by NCI Cancer Center Support Grant 5P30 CA 42014, and by a sub-contract of NCI grant 1 P01 Ca 76466 held by the Fred Hutchinson Cancer Research Center (FHCRC). The contents of this paper are solely the responsibility of the authors and do not necessarily represent the ocial views of the FHCRC or the Awarding Agency. A part of this research was carried out while A.Y.Y. was visiting the Institute of Applied Mathematics and Statistics, University of W urzburg, as a recipient of the Alexander von Humboldt Research Award.

References [1] B. Hofstad, M. Vatn, Growth rate of colon polyps and cancer, Gastrointest. Endosc. Clini. North Am. 7 (1997) 345. [2] H.E. Knoernschild, Growth rate and malignant potential of colonic polyps: Early results, Surg. Forum 14 (1963) 137. [3] R.P. Castleberry, Biology and treatment of neuroblastoma, Pediatr. Clin. North. Am. 44 (1997) 919. [4] K.K. Matthay, Neuroblastoma: A clinical challenge and biological puzzle, CA Cancer J. Clin. 45 (1995) 179. [5] G.M. Brodeur, J.M. Maris, D.J. Yamashiro, M.D. Hogarty, P.S. White, Biology and genetics of human neuroblastomas, J. Pediatr. Hematol. Oncol. 19 (1997) 93. [6] C.M. Kullendor€, L.G. Stromblad, Spontaneous regression of intraspinal neuroblastoma, Pediatr. Surg. Int. 12 (1997) 305. [7] M. Grossman, R. Hoermann, M. Weiss, K.W. Jauch, H. Oertel, A. Staebler, K. Mann, D. Engelhardt, Sponateneous regression of hepatocellular carcinoma, Am. J. Gastroenterol. 90 (1995) 1500. [8] M. Ide, M. Jimbo, M. Yamamoto, S. Hagiwara, M. Aiba, O. Kubo, Spontaneous regression of primary intracranial germinoma. A case report, Cancer 79 (1997) 558. [9] P.H. Phillips, N.J. Newman, Here today ... gone tomorrow, Surv. Ophthalmol. 41 (1997) 354. [10] S.D. Bos, H.J. Mensink, Spontaneous caval tumor thrombus necrosis and regression of pulmonary lesions in renal cell cancer, Scand. J. Urol. Nephrol. 30 (1996) 489. [11] F.X. Roithinger, M. Wallner, C. Baldinger, R. Puchner, L. Leitner, W. Hobling, F.C. Prischl, Spontaneous regression of a small cell bronchial carcinoma, Pneumologie (in German) 49 (1995) 509. [12] G.M. Halliday, A. Patel, M.J. Hunt, F.J. Tefany, R.S. Barnetson, Spontaneous regression of human melanoma/nonmelanoma skin cancer: Association with in®ltrating CD4+T cells, World J. Surg. 19 (1995) 352. [13] Y. Ogihara, K. Takeda, T. Yanagawa, Y. Hirasawa, Spontaneous regression of lung metastases from osteosarcoma, Cancer 74 (1994) 2798.

60

A. Yakovlev et al. / Mathematical Biosciences 155 (1999) 45±60

[14] D. Rayson, H.C. Pitot, L.K. Kvols, Regression of metastatic carcinoid tumor after valvular surgery for carcinoid heart disease, Cancer 79 (1997) 605. [15] G. Luebeck, S.H. Moolgavkar, Stochastic analysis of intermediate lesions in carcinogenesis experiments, Risk Analy. 11 (1991) 149. [16] W.Y. Tan, A stochastic Gompertz birth-death process, Statist. Prob. Lett. 4 (1986) 25. [17] S.H. Moolgavkar, A.G. Knudson, Mutation and cancer: a model for human carcinogenesis, J. Nat. Cancer Inst. 66 (1981) 1037. [18] S.H. Moolgavkar, D.J. Venzon, Two-event models for carcinogenesis: Incidence curves for childhood and adult tumors, Math. Biosci. 47 (1979) 55. [19] T.E. Harris, The Theory of Branching Processes, Springer, Berlin, 1963. [20] P. Jagers, Branching Processes with Biological Applications, Wiley, New York, 1975. [21] D.G. Kendall, On the generalized birth-and-death process, Ann. Math. Statist. 19 (1948) 1. [22] L.G. Hanin, S.T. Rachev, A.D. Tsodikov, A.Y. Yakovlev, A stochastic model of carcinogenesis and tumor size at detection, Adv. Appl. Prob. 29 (1997) 607. [23] B.A. Sevast'yanov, The theory of branching random processes, Uspehi Matemat. Nauk. (in Russian) 6 (1951) 47. [24] P. Jagers, The proportions of individuals of di€erent kinds in two-type populations. A branching process problem arising in biology, J. Appl. Prob. 6 (1969) 249. [25] Y. Ozeki, T. Matsubara, K. Tateyama, M. Kokubo, H. Shimoji, M. Katayama, Spontaneous complete necrosis of hepatocellular carcinoma, Am. J. Gastroenterol. 91 (1996) 391. [26] G.P. Tonini, K. Mazzocco, A. di Vinci, E. Geido, B. de Bernardi, W. Giaretti, Evidence of apoptosis in neuroblastoma at onset and relapse. An analysis of a large series of tumors, J. Neurooncol. 31 (1997) 209. [27] R. Bartoszy nski, A modeling approach to metastatic progression of cancer, in: J.R. Thompson, B.W. Brown (Eds.), Cancer Modeling, Marcel-Dekker, New York, 1987, p. 237. [28] J.P. Klein, R. Bartoszy nski, Estimation of growth and metastatic rates of primary breast cancer, in: O. Arino, D.E. Axelrod, M. Kimmel (Eds.), Mathematical Population Dynamics, Marcel-Dekker, New York, 1991, p. 397. [29] K. L. Chung, A Course in Probability Theory, 2nd ed., Academic Press, New York, 1974.