Iterated birth and death process as a model of radiation cell survival

Mathematical Biosciences 169 (2001) 89±107

www.elsevier.com/locate/mbs

Iterated birth and death process as a model of radiation cell survival Leonid G. Hanin * Department of Mathematics, Idaho State University, Pocatello, ID 83209-8085, USA Received 6 April 2000; received in revised form 8 September 2000; accepted 18 September 2000

Abstract The iterated birth and death process is de®ned as an n-fold iteration of a stochastic process consisting of the combination of instantaneous random killing of individuals in a certain population with a given survival probability s with a Markov birth and death process describing subsequent population dynamics. A long standing problem of computing the distribution of the number of clonogenic tumor cells surviving a fractionated radiation schedule consisting of n equal doses separated by equal time intervals s is solved within the framework of iterated birth and death processes. For any initial tumor size i, an explicit formula for the distribution of the number M of surviving clonogens at moment s after the end of treatment is found. It is shown that if i ! 1 and s ! 0 so that isn tends to a ®nite positive limit, the distribution of random variable M converges to a probability distribution, and a formula for the latter is obtained. This result generalizes the classical theorem about the Poisson limit of a sequence of binomial distributions. The exact and limiting distributions are also found for the number of surviving clonogens immediately after the nth exposure. In this case, the limiting distribution turns out to be a Poisson distribution. Ó 2001 Elsevier Science Inc. All rights reserved. Keywords: Birth and death process; Branching process; Clonogenic tumor cell; Fractionated cancer radiotherapy; Limiting distribution; Probability distribution; Probability generating function; Tumor recurrence

1. Introduction The purpose of this work is to solve, under natural biological assumptions, the following problem of major theoretical and practical importance in radiation oncology that was posed in [1] as early as in 1990: What is the distribution of the number of clonogenic tumor cells surviving

*

Tel.: +1-208 236 3293; fax: +1-208 236 2636. E-mail address: [email protected] (L.G. Hanin).

0025-5564/01/$ - see front matter Ó 2001 Elsevier Science Inc. All rights reserved. PII: S 0 0 2 5 - 5 5 6 4 ( 0 0 ) 0 0 0 5 4 - 7

90

L.G. Hanin / Mathematical Biosciences 169 (2001) 89±107

fractionated radiation? Lack of theoretical results leading to a closed-form analytic expression for this distribution, even under very simplistic models of cell kinetics between radiation exposures, was the most critical deterrent to evaluating and increasing the ecacy of existing clinical practices of cancer treatment and to developing relevant methods of statistical data analysis. We proceed from the following model of tumor population kinetics, see e.g. [2]. A tumor initially comprising a non-random number i of clonogenic cells is exposed to a fractionated radiation schedule consisting of n instantaneously delivered equal doses D separated by equal time intervals s. It is assumed that every clonogen survives each exposure to the dose D with the same probability s ˆ s…D†, given that it survived the previous exposures, and independently of other clonogens. Observe that it is not supposed that all cells forming a tumor are necessarily clonogenic. We also assume that irradiated tumor cells are killed instantaneously. According to a common radiobiological convention, a tumor cell is killed if it is incapable of producing a viable clone. Between the exposures, clonogens proliferate or die spontaneously independently of each other with constant rates k > 0 and m P 0 (de®ned in the sense of Markov processes), respectively. Finally, it is assumed that cell death between radiation exposures is instantaneous, and that descendants of a clonogenic cell are clonogenic. Similar models arise when a tumor cell population is exposed to chemical or immunological agents, as well as in problems of population dynamics in the presence of emigration and external or internal killing factors. Suppose that radiation exposures occur at time moments 0; s; . . . ; …n ÿ 1†s. Let M be the random number of clonogens (including the original clonogenic cells and their descendants) that are alive at time ns, that is, at time s after the nth exposure. Random variable (r.v.) M can be viewed as the ®nal state of a stochastic process de®ned as the n-fold iteration of the combination of exposure to the dose D with the homogeneous birth and death Markov process with the birth rate k and death rate m. This n-stage stochastic process will be referred to in what follows as the iterated birth and death process. Also, denote by L the number of surviving clonogens immediately after delivery of the nth fraction of radiation. Although r.v. L is biologically more natural than M, the latter is more tractable mathematically. Alternatively, r.v.s M and L can be characterized in terms of the Galton±Watson branching process (see e.g. [3, Chapter 1], and [4, Chapter 2]). Denote by Zs the state at time s of the homogeneous birth and death Markov process with the birth rate k, death rate m, and initial population size i ˆ 1 (that can be also viewed as a continuous time branching process with the probabilities of no and two o€spring equal to m=…k ‡ m† and k=…k ‡ m†, respectively). Then the distribution of r.v. M coincides with the distribution of the nth generation size in the Galton± Watson process with the initial population size i and the o€spring distribution de®ned as the distribution of a r.v. that is equal to Zs with probability s and to 0 with probability 1 ÿ s. Similarly, r.v. L is equidistributed with the size of the …n ÿ 1†st generation in the Galton±Watson process where the initial population size follows the binomial distribution B…i; s† and the o€spring distribution is B…Zs ; s†. In what follows, however, we will use the iterated birth and death process approach because it is more straightforward and technically simple. Most of the previous theoretical work was related to the case k ˆ m ˆ 0, where no cell proliferation and cell death occur between radiation exposures, and was primarily focused on the extinction probability p~0 :ˆ Pr…L ˆ 0† (which is also referred to in radiation oncology as the tumor control probability). In this case, the probability that a given clonogenic cell will survive the n-dose

L.G. Hanin / Mathematical Biosciences 169 (2001) 89±107

91

irradiation schedule is equal to sn . Hence, M ˆ L, and these r.v.s follow the binomial distribution B…i; sn †. The initial number i of clonogenic cells in a tumor is typically very large (1 cm3 of a solid tumor contains about 109 cells [5]; also, a clinically detectable tumor is estimated to contain at least 105 clonogenic cells probably ranging up to 109 cells or even more [6]). Furthermore, the survival probability sn is small because the dose D, usually 1±2 Gy, and the number of fractions n, normally 20±40, are chosen to achieve just that. Therefore, the distribution B…i; sn † is close to the Poisson distribution with parameter h ˆ isn . It is this Poisson approximation that is used for planning fractionated cancer radiotherapy in current clinical practice, see e.g. [2,7,8]. The validity of this Poisson approximation was ®rst questioned in [1]. The idea of this work was that due to cell proliferation the distribution of the number of surviving clonogens must deviate from the binomial distribution B…i; sn † which makes the approximating Poisson distribution P …h† not well justi®ed. Extensive numerical simulations conducted in [1] con®rmed that in the absence of cell death between radiation exposures (m ˆ 0) the mean and variance of r.v. L are distinct from those predicted by the binomial distribution B…i; sn †. On this basis, the authors made a conclusion that the distribution of r.v. L cannot be approximated by a Poisson distribution. In response to this work, it was suggested in [9] to use branching stochastic processes as an adequate tool for studying the distribution of r.v. L (for an extensive discussion of this idea, the reader is referred to [10, pp. 31±40]) but the model proposed in this work was too simplistic. Numerical algorithms for computing the extinction probability p~0 in the case m ˆ 0 are described in [8,11]. The work [11] also contains a more general algorithm for numerical evaluation of the probability generating function (p.g.f.) of r.v. L. A numerical procedure for computing the extinction probability p~0 within a more elaborate model accounting not only for cell division but also for cell loss, cell di€erentiation, and cell cycle e€ects is presented in [6]. The paper [2] represents the ®rst attempt at theoretical approach to describing the distribution of the number of surviving clonogens. It contains explicit formulas for p0 :ˆ Pr…M ˆ 0† and the p.g.f. of r.v. M obtained within the framework of homogeneous birth process. The following problems are therefore of theoretical and practical interest. Problem 1. Find explicit computationally feasible formulas for the distribution of the ®nal state M of the n-fold iterated birth and death process and for its counterpart L. Problem 2. Compute the limiting distributions of r.v.s M and L when i ! 1; s ! 0 and isn ! h; 0 < h < 1: In the present communication, we solve Problems 1 and 2. For an at length discussion of biomedical and statistical implications of the results of this work, the reader is referred to [12]. A model based on our ®ndings was successfully applied to statistical analysis of data on prostate cancer post-treatment survival [13]. The structure of the paper is as follows. In Section 2, we compute p.g.f. of the iterated birth and death process and derive some of the immediate consequences for the r.v.s M and L. Their distribution is identi®ed in Section 3, which leads to a solution of Problem 1. Section 4 addresses limit theorems for the iterated birth and death process that solve Problem 2. A numerical example corroborating our theoretical ®ndings is presented in Section 5. Finally, in Section 6, we brie¯y discuss our results from biomedical and statistical data analysis perspectives.

92

L.G. Hanin / Mathematical Biosciences 169 (2001) 89±107

2. Probability generating function of the iterated birth and death process In this section, we will compute p.g.f. of the ®nal state M of the n-fold iterated birth and death process with i initial cells. As a consequence, we will ®nd p.g.f. of the number L of surviving clonogens immediately after the nth exposure. Since p.g.f. is an indispensable tool for analysis of Markov processes with the state space Z‡ (see e.g. [14, pp. 184±192]), we ®rst recall some wellknown facts about them. For a non-negative integer valued r.v. X, denote by uX its p.g.f. de®ned by 1 X X Pr…X ˆ m†zm ; jzj 6 1; uX …z† :ˆ E…z † ˆ mˆ0

where E stands for the expectation. The probabilities Pr…X ˆ m† are related to the p.g.f. of X by Pr…X ˆ m† ˆ

Dm uX …0† ; m!

m P 0:

…1†

In particular, uX …0† ˆ Pr…X ˆ 0†. Also observe that uX …1† ˆ 1; EX ˆ u0X …1†;

…2†

Var X ˆ u00X …1† ‡ u0X …1† ÿ ‰u0X …1†Š2 ;

…3†

and where Var denotes the variance. The expectation and variance of X can also be expressed through the logarithmic derivatives of uX as follows: EX ˆ …ln uX †0 …1†

and Var X ˆ …ln uX †00 …1† ‡ …ln uX †0 …1†:

…4†

For a function u whose range is contained in the domain, we denote by u…n† its n-fold composition with itself. We start our calculation of p.g.f. uM with the case i ˆ 1; n ˆ 1: Since each clonogenic cell survives irradiation in dose D with the probability s, the number of survivors for a single clonogen exposed to the dose D is a Bernoulli r.v. with p.g.f. b…z† ˆ 1 ÿ s ‡ sz:

…5†

According to our assumptions, expansion of a surviving clonogenic cell is governed by a homogeneous birth and death process with the birth rate k and death rate m. It was shown in [15] (see also [16, pp. 165±166]) that the state of such process at time s, that is, the size S at time s of the clone emerging from a single cell, has a generalized geometric distribution of the form Pr…S ˆ 0† ˆ r;

Pr…S ˆ m† ˆ …1 ÿ r†…1 ÿ q†qmÿ1 ;

m P 1:

Parameters r and q, 0 6 r; q < 1; of this distribution are related to the birth and death rates through the formulas rˆ

m…1 ÿ a† k ÿ am

and q ˆ

k…1 ÿ a† ; k ÿ am

…6†

L.G. Hanin / Mathematical Biosciences 169 (2001) 89±107

93

where it is denoted hereafter a :ˆ e…mÿk†s : P.g.f. of r.v. S is given by 1 X r ÿ …r ‡ q ÿ 1†z : uS …z† ˆ r ‡ …1 ÿ r†…1 ÿ q†qmÿ1 zm ˆ 1 ÿ qz mˆ1

…7†

Observe that k ÿ am ˆ 0 only in the critical case k ˆ m, that is, when a ˆ 1. In this case, rˆqˆ

ks 1 ‡ ks

and uS …z† ˆ

ks ‡ …1 ÿ ks†z ; 1 ‡ ks ÿ ksz

which coincides with the corresponding limits in (6) and (7). It follows from (2) and (7) that ES ˆ u0S …1† ˆ

1ÿr ˆ aÿ1 : 1ÿq

…8†

The number X of o€spring of a single irradiated clonogenic cell at time s after exposure has p.g.f. u ˆ 1 ÿ s ‡ suS ˆ b  uS , that is, in view of (5) and (7), u…z† ˆ 1 ÿ s ‡ s

r ÿ …r ‡ q ÿ 1†z 1 ÿ s…1 ÿ r† ÿ ‰q ÿ s…1 ÿ r†Šz ˆ : 1 ÿ qz 1 ÿ qz

…9†

Therefore, r.v. X follows a generalized geometric distribution. By (2) and (8) we have s l :ˆ EX ˆ u0 …1† ˆ su0S …1† ˆ ˆ se…kÿm†s : …10† a Observe that p.g.f. u of r.v. X given in (9) is a fractional linear function. Conversely, every fractional linear p.g.f. of the form …a ÿ bz†=…c ÿ dz† with c; d 6ˆ 0; c 6ˆ d; emerges from a generalized geometric distribution. A standard argument about p.g.f. of a branching process (see e.g. [14, pp. 190±192]) implies that  i …11† uM ˆ u…n† : To compute the function u…n† explicitly, we need a more convenient formula for the p.g.f. u expressed in (9) through parameters s…1 ÿ r† and q. In the generic case a 6ˆ 1 we set 1 ÿ s…1 ÿ r† k ÿ ma ÿ s…k ÿ m† ˆ ; …12† x :ˆ q k…1 ÿ a† while in the critical case a ˆ 1 we de®ne by continuity x :ˆ …1 ÿ s ‡ ks†=…ks†: It is easy to see that x is a ®xed point of the function u: Since 1ÿxˆ

…s ÿ a†…k ÿ m† ; k…1 ÿ a†

…13†

we have 0 < x < 1 for l > 1; x ˆ 1 for l ˆ 1; and x > 1 for 0 6 l < 1: Observe that lÿ1 l…1 ÿ x† qˆ and s…1 ÿ r† ˆ : lÿx lÿx This allows us to rewrite (9) in terms of parameters l and x as follows: x…l ÿ 1† ÿ …xl ÿ 1†z : u…z† ˆ l ÿ x ÿ …l ÿ 1†z

…14†

94

L.G. Hanin / Mathematical Biosciences 169 (2001) 89±107

An explicit expression for the function u…n† is given in formula (15) below. Its remarkable feature is that the action of the semigroup Z‡ on the function u written in form (14) by compositions is simply n 7! ln ; n 2 Z‡ : An equivalent but more cumbersome formula can be found in [9, p. 9], and (with a proof) in [17, pp. 296±298], [18, p. 7], and [19, p. 57]. In the case m ˆ 0; formula (15) can be found in a slightly di€erent form in [2]. Proposition 1. Suppose that l 6ˆ 1. Then, for every natural number n, u…n† …z† ˆ

x…ln ÿ 1† ÿ x…ln ÿ 1†z : ln ÿ x ÿ …ln ÿ 1†z

…15†

Proof. We will use induction in n. For n ˆ 1, Eq. (15) coincides with (14). Suppose it holds for some n P 1: Then we have by (14) and (15) u…n‡1† …z† ˆ …u…n†  u†…z† ˆ

x…ln ÿ 1† ÿ …xln ÿ 1†u…z† An ÿ Bn z ˆ ; ln ÿ x ÿ …ln ÿ 1†u…z† Cn ÿ Dn z

where after an elementary calculation we ®nd that An ˆ x…l ÿ x†…ln ÿ 1† ÿ x…l ÿ 1†…xln ÿ 1† ˆ x…1 ÿ x†…ln‡1 ÿ 1†; Bn ˆ x…l ÿ 1†…ln ÿ 1† ÿ …xl ÿ 1†…xln ÿ 1† ˆ …1 ÿ x†…xln‡1 ÿ 1†; Cn ˆ …l ÿ x†…ln ÿ x† ÿ x…l ÿ 1†…ln ÿ 1† ˆ …1 ÿ x†…ln‡1 ÿ x†; Dn ˆ …l ÿ 1†…ln ÿ x† ÿ …xl ÿ 1†…ln ÿ 1† ˆ …1 ÿ x†…ln‡1 ÿ 1†: Since l 6ˆ 1; we also have x 6ˆ 1. Therefore, the above formulas for An ; Bn ; Cn and Dn lead to the expression for u…n‡1† required in (15). This completes the proof of Proposition 1. Remark 1. By passage to limit in (15) as l ! 1, we ®nd that for l ˆ 1 u…n† …z† ˆ

nk…1 ÿ a† ÿ ‰nk…1 ÿ a† ÿ a…k ÿ m†Šz : nk…1 ÿ a† ‡ a…k ÿ m† ÿ nk…1 ÿ a†z

According to (11) and (15), p.g.f. of r.v. M is given by  i x…ln ÿ 1† ÿ …xln ÿ 1†z : uM …z† ˆ ln ÿ x ÿ …ln ÿ 1†z Thinking of the number of iterations n as being ®xed, we will also write uM in the form  i a ÿ bz ; uM …z† ˆ c ÿ dz

…16†

…17†

where a :ˆ x…ln ÿ 1†;

b :ˆ xln ÿ 1;

c :ˆ ln ÿ x;

d :ˆ ln ÿ 1:

…18†

L.G. Hanin / Mathematical Biosciences 169 (2001) 89±107

95

Using (17) and (18) we evaluate the extinction probability  a i  x…ln ÿ 1† i p0 ˆ Pr…M ˆ 0† ˆ uM …0† ˆ ˆ : c ln ÿ x In particular, in the pure birth setting …m ˆ 0†; we have x ˆ …1 ÿ s†=…1 ÿ a† which leads to the following result obtained in [2]:  i …1 ÿ s†…ln ÿ 1† : p0 ˆ …1 ÿ a†ln ‡ s ÿ 1 Formula (16) allows us to ®nd p.g.f. of r.v. L: In view of uL ˆ ‰u…nÿ1† Ši  b and (5), !i  i ~ 1 ÿ x ÿ s ‡ axln ÿ …axln ÿ s†z a~ ÿ bz uL …z† ˆ ˆ ; ~ 1 ÿ x ÿ s ‡ aln ÿ …aln ÿ s†z c~ ÿ dz

…19†

where a~ ˆ 1 ÿ x ÿ s ‡ axln ;

b~ ˆ axln ÿ s;

c~ ˆ 1 ÿ x ÿ s ‡ aln ;

d~ ˆ aln ÿ s:

In particular, the tumor control probability p~0 ˆ Pr…L ˆ 0† is given by !i  i a~ 1 ÿ x ÿ s ‡ axln ˆ : p~0 ˆ c~ 1 ÿ x ÿ s ‡ aln

…20†

…21†

From (16) and (19) we conclude that in the case i ˆ 1 r.v.s M and L follow generalized geometric distribution. In the case of continuous irradiation with an arbitrary dose rate, p.g.f. for the number of surviving clonogens at the end of treatment and the corresponding tumor control probability were computed within the framework of continuous time birth and death process in [20]. Let U be any non-negative integer valued r.v. with a p.g.f. of the form a ÿ bz : uU …z† ˆ c ÿ dz Then D 2dD and u00U …1† ˆ ; EU ˆ u0U …1† ˆ 2 …c ÿ d† …c ÿ d†3 where D :ˆ ad ÿ bc: Therefore, by (3) Var U ˆ

D D2 c2 ÿ d 2 ÿ D ‡ ÿ ˆ D : …c ÿ d†3 …c ÿ d†2 …c ÿ d†4 …c ÿ d†4 2dD

…22†

Let also V :ˆ U1 ‡


‡ Ui ; where Uk ; 1 6 k 6 i; are independent identically distributed (i.i.d.) r.v.s having the same distribution as U. Then EV ˆ i EU ; Var V ˆ iVar U; and  i a ÿ bz uV …z† ˆ : …23† c ÿ dz For V ˆ M; a straightforward calculation based on (18) yields c ÿ d ˆ 1 ÿ x; c ‡ d ˆ 2ln ÿ x ÿ 1; and

96

L.G. Hanin / Mathematical Biosciences 169 (2001) 89±107

DM :ˆ ad ÿ bc ˆ ln …1 ÿ x†2 :

…24†

Therefore, EM ˆ

iDM …c ÿ d†

2

ˆ iln :

…25†

Also, invoking (22) and (13) we ®nd that, for l 6ˆ 1 and a 6ˆ 1; Var M ˆ iln …ln ÿ 1†

1‡x k…2 ÿ s ÿ a† ‡ m…s ÿ a† ˆ iln …ln ÿ 1† : 1ÿx …k ÿ m†…s ÿ a†

…26†

In a similar manner, for V ˆ L; we derive from (20) that c~ ÿ d~ ˆ 1 ÿ x; c~ ‡ d~ ˆ 1 ÿ x ÿ 2s ‡ 2aln ; and ~c ˆ aDM ˆ aln …1 ÿ x†2 : DL :ˆ a~d~ ÿ b~

…27†

Hence, EL ˆ aEM ˆ ialn

…28†

and, for l 6ˆ 1 and a 6ˆ 1; Var L ˆ ialn

1 ÿ x ÿ 2s ‡ aln …1 ‡ x† : 1ÿx

…29†

The values of Var M and Var L in the critical cases l ˆ 1 or/and a ˆ 1 can be easily computed by passage to appropriate limits in (26) and (29). 3. Distribution of the ®nal state of the iterated birth and death process In this section, we are going to compute the distribution pm ; m P 0; of any non-negative integer valued r.v.V with a p.g.f. of form (23). In particular, our argument applies to r.v.s M and L. We write the fractional linear function involved in (23) as follows:   a ÿ bz b D 1 ˆ 1ÿ
; …30† c ÿ dz d bd z ÿ c=d where D ˆ ad ÿ bc: Then  i X   k i b D i k uV …z† ˆ …ÿ1† …z ÿ c=d†ÿk : k d bd kˆ0 Di€erentiating this equality m times and setting z ˆ 0 we obtain on the basis of (1) that  i  m X  k i   b d D i k‡mÿ1 ; m P 0: pm ˆ k m d c bd kˆ1

…31†

Since the number i is very large, this formula can barely be used for computing the distribution of r.v. V. Fortunately, the probabilities pm can be represented in a much more ecient way (see

L.G. Hanin / Mathematical Biosciences 169 (2001) 89±107

97

Proposition 4) that will allow us to ®nd in Section 4 limiting distributions of r.v.s M and L when i ! 1 and s ! 0: Introduce the polynomials  i   X i k‡mÿ1 k Pm …x† :ˆ …32† x ; m P 0: k m kˆ1 Then it follows from (31) that  i  m   b d D Pm ; pm ˆ d c bc

m P 0:

…33†

From (32) we see readily that i

P0 …x† ˆ …1 ‡ x† : Also,   i X i k k x ˆ x‰…1 ‡ x†i Š0 ˆ ix…1 ‡ x†iÿ1 : P1 …x† ˆ k kˆ1

…34† …35†

Eqs. (34) and (35) invite us to de®ne functions Qm ; m P 0; by Pm …x† ˆ …1 ‡ x†iÿm Qm …x†;

m P 0:

…36†

Then Q0 …x† ˆ 1

and Q1 …x† ˆ ix:

…37†

Indeed, functions Pm and Qm depend on i: However, since in this section i is assumed to be ®xed, this dependence will not be shown notationally. In the next statement we ®nd a simple recurrence relation for polynomials Pm : Proposition 2. Pm‡1 …x† ˆ

1 ‰mPm …x† ‡ xPm0 …x†Š; m‡1

m P 0:

Proof. According to (34) and (35), relation (38) holds for m ˆ 0. In view of (32)  i   i   X i k ‡ m k X i k…k ‡ 1†


…k ‡ m† k Pm‡1 …x† ˆ x ˆ x ; m P 0: k m‡1 k …m ‡ 1†! kˆ1 kˆ1 Hence for m P 1 Z xmÿ1 Pm‡1 …x† dx ˆ

 i   xm X xm i k‡mÿ1 k x ‡C ˆ Pm …x† ‡ C: m m ‡ 1 kˆ1 k m‡1

Di€erentiating this equality we obtain xmÿ1 Pm‡1 …x† ˆ

1 ‰mxmÿ1 Pm …x† ‡ xm Pm0 …x†Š; m‡1

which implies (38). Proposition 2 is proved.

…38†

98

L.G. Hanin / Mathematical Biosciences 169 (2001) 89±107

From Proposition 2 we derive the following result about the functions Qm ; m P 0: Proposition 3. fQm g1 mˆ0 is a sequence of polynomials satisfying the following recurrence relation: Qm‡1 …x† ˆ

1 ‰…m ‡ ix†Qm …x† ‡ x…1 ‡ x†Q0m …x†Š; m‡1

m P 0:

…39†

Proof. Using (36) we express the recurrence relation (38) in terms of the functions Qm as follows: i 1 h …1 ‡ x†iÿmÿ1 Qm‡1 …x† ˆ m…1 ‡ x†iÿm Qm …x† ‡ …i ÿ m†x…1 ‡ x†iÿmÿ1 Qm …x† ‡ x…1 ‡ x†iÿm Q0m …x† ; m‡1 which implies (39) for m P 1: In view of (37) relation (39) also holds for m ˆ 0. Together with (37) this implies that every function Qm ; m P 0; is a polynomial of degree m. Proposition 3 is proved. Our next statement provides the required formula for the polynomials Qm . For i > m; it o€ers a signi®cant computational advantage over formula (32) for the polynomials Pm : Proposition 4. Qm …x† ˆ

  m  X mÿ1 i‡kÿ1 kˆ1

mÿk

k

xk ;

m P 1:

…40†

Proof. Denote by Rm ; m P 1; the polynomial in the right-hand side of (40). For m ˆ 1; using (37) we have R1 …x† ˆ ix ˆ Q1 …x†: Therefore, it suces to show that polynomials Rm ; m P 1; satisfy the recurrence relation (39). For a ®xed m P 1; denote 1 ‰…m ‡ ix†Rm …x† ‡ x…1 ‡ x†R0m …x†Š m ‡ 1" #     m  m  X X 1 mÿ1 i‡kÿ1 k mÿ1 i‡kÿ1 kÿ1 …m ‡ ix† : ˆ x ‡ x…1 ‡ x† kx mÿk k mÿk k m‡1 kˆ1 kˆ1

U …x† :ˆ

We have to check that the polynomial U …x† ˆ Pm‡1 V …x† :ˆ Rm‡1 …x† ˆ kˆ1 vk xk ; where   m‡1  X m i‡kÿ1 k x: Rm‡1 …x† ˆ mÿk‡1 k kˆ1

Pm‡1 kˆ1

…41† uk xk coincides with the polynomial

From (41) and(42) we ®nd that u1 ˆ …mi ‡ i†=…m ‡ 1† ˆ i ˆ v1 : Also,     i‡m i‡mÿ1 i‡m um‡1 ˆ ˆ ˆ vm‡1 : m m‡1 m‡1 Next, for 2 6 k 6 m; we have by (41)

…42†

L.G. Hanin / Mathematical Biosciences 169 (2001) 89±107

uk ˆ

m‡k m‡1



mÿ1 mÿk



i‡kÿ1 k

 ‡

i‡kÿ1 m‡1



mÿ1 mÿk‡1



99



i‡kÿ2 : kÿ1

…43†

Observe that     mÿk‡1 mÿ1 m ˆ : mÿk mÿk‡1 m Similarly,         kÿ1 k m i‡kÿ2 i‡kÿ1 mÿ1 and ˆ : ˆ mÿk‡1 kÿ1 k mÿk‡1 m i‡kÿ1 Therefore, we continue (43) to ®nd through a simple calculation that    m i‡kÿ1 vk ˆ ˆ vk ; 2 6 k 6 m: mÿk‡1 k This completes the proof of Proposition 4. The following result solves Problem 1 formulated in Section 1. Theorem 1. Let V be a non-negative integer valued r.v. the p.g.f. of which has form (23) with some natural number i (in particular, this is the case for V ˆ M and V ˆ L). Then for the distribution of V we have  a i  b m  D  ; m P 0; …44† Qm pm ˆ c a bc where D ˆ ad ÿ bc; Q0 …x† ˆ 1, and polynomials Qm , m P 1, are given in (40). Proof. Since uV …1† ˆ 1; we set z ˆ 1 in (23) to ®nd that a ÿ b ˆ c ÿ d; hence a ˆ b ‡ c ÿ d: Also, setting z ˆ 1 in (30) we obtain that D ˆ …b ÿ d†…d ÿ c†; then bc ‡ D ˆ d…b ‡ c ÿ d† ˆ ad: Now invoking formulas (33), (36) and (45) we ®nd that  i  m  iÿm      m   b d D D a i b D ˆ ; 1‡ Qm Qm pm ˆ d c bc bc c a bc

…45†

m P 0:

Theorem 1 is proved. Remark 2. Every r.v. V with p.g.f. (23), where c; d 6ˆ 0 and c 6ˆ d; can be represented as the sum of i independent copies of a r.v. having generalized geometric distribution. In the particular case of geometric distribution …a ˆ 0†; the distribution of r.v. V is negative binomial. Therefore, distribution (44) can (and will) be called generalized negative binomial.

100

L.G. Hanin / Mathematical Biosciences 169 (2001) 89±107

4. Limit theorem for the iterated birth and death process The following theorem establishes the form of the limiting distribution of r.v. M as i ! 1 and s ! 0. Parameters n P 1, s > 0; k > 0 and m P 0 are assumed to be ®xed. For m ˆ 0; k ! 0, this theorem gives the classical Poisson limit of a sequence of binomial distributions. Theorem 2. Suppose that i ! 1 and s ! 0 so that there exists a limit h :ˆ lim isn ;

0 < h < 1:

…46†

Then the distribution of r.v. M converges to the probability distribution p ˆ fpm g1 mˆ0 given by pm ˆ eÿA qÿm rm ;

m P 0;

…47†

where qˆ

k ÿ am ; k…1 ÿ a†

…48†

Aˆ

h…q ÿ 1† h…k ÿ m† ˆ nÿ1 ; an q a …k ÿ am†

…49†

r0 ˆ 1, rm ˆ

 k m  X mÿ1 d kˆ1

mÿk

k!

;

m P 1;

…50†

and 2

dˆ

2

h…q ÿ 1† h…k ÿ m† ˆ nÿ2 : an q ka …1 ÿ a†…k ÿ am†

…51†

Proof. It follows from (12) and (48) that under conditions of the theorem x ! q. Note that qÿ1ˆ

a…k ÿ m† > 0; k…1 ÿ a†

whence q > 1. Recall formulas (10) and (18) to ®nd that a ! ÿq;

b ! ÿ1;

c ! ÿq:

…52†

Also, by (24), (46), (52) and (51) iDM h…q ÿ 1†2 ! ˆ d: bc an q Observe that according to (44)   m  b DM pm ˆ p0 Qm ; m P 0: bc a

…53†

…54†

L.G. Hanin / Mathematical Biosciences 169 (2001) 89±107

101

First, we compute the limit for p0 ˆ …a=c†i . Note that by (18) a xÿ1 ˆ 1 ‡ ln n : c l ÿx Therefore, p0 ˆ …a=c†i ! eÿA ;

where A ˆ lim iln

xÿ1 h…q ÿ 1† ˆ ; x ÿ ln an q

…55†

compare with (49). Next, invoking Proposition 4, we write  k ÿ1     X k m  DM j iDM mÿ1 1 Y 1‡ ; m P 1: ˆ Qm m ÿ k k! jˆ0 i bc bc kˆ1 Passage to limit in this equation with (53) taken into account yields    k m  X DM mÿ1 d ! ˆ rm ; m P 1: Qm m ÿ k k! bc kˆ1

…56†

Finally, in view of (52), b=a ! 1=q. Combining this with (55) and (56) we obtain from (54) the required form (47) of pm for m P 1: It remains to show that p is a proper probability distribution, that is, 1 X pm ˆ 1: …57† mˆ0

To this end, we will compute p.g.f. up for the sequence p ˆ fpm g1 mˆ0 . By (47) and (50) we have for jzj < q "  #  m 1 1 1  m X m  X X X z z dk m ÿ 1 up …z† ˆ pm zm ˆ eÿA rm ˆ eÿA 1 ‡ m ÿ k q q k! mˆ0 " mˆ0 kˆ1 # mˆ1     m 1 1 k X d X mÿ1 z ˆ eÿA 1 ‡ q k! mˆk m ÿ k kˆ1 "    j #   k 1 1 X 1 dz X z k‡jÿ1 ÿA ; ˆe 1‡ j k! q q kˆ1 jˆ0

…58†

where j :ˆ m ÿ k. Observe that the last sum represents the binomial series ÿk  j  1  X z z k‡jÿ1 ˆ 1ÿ ; jzj < q; k P 1: j q q jˆ0 Therefore, using (51) we continue (58):  k     1 X 1 dz h…q ÿ 1† dz ÿA up …z† ˆ e ˆ exp exp k! q ÿ z an q qÿz kˆ0     h…q ÿ 1† …q ÿ 1†z h…q ÿ 1† 1 ÿ z ˆ exp ÿ ˆ exp ÿ :
1ÿ an q qÿz an qÿz Setting here z ˆ 1 we establish (57). This completes the proof of Theorem 2.

…59†

102

L.G. Hanin / Mathematical Biosciences 169 (2001) 89±107

Remark 3. For the iterated pure birth process (m ˆ 0) , the limiting distribution of r.v. M under conditions of Theorem 2 is given by (56) with qˆ

1 ; 1ÿa

Aˆ

h anÿ1

;

dˆ

h anÿ2 …1

ÿ a†

; where a ˆ eÿks :

In the critical case k ˆ m > 0 we have qˆ1‡

1 ; ks

Aˆ

h 1 ‡ ks

and d ˆ

h : ks…1 ‡ ks†

Remark 4. In view of formulas (4) Eq. (59) leads to the following expressions for the expectation and variance of the limiting distribution p of r.v. M: Ep ˆ

h an

and Var p ˆ

h…q ‡ 1† h‰k…2 ÿ a† ÿ amŠ ˆ : an …q ÿ 1† an‡1 …k ÿ m†

The same results arise when we pass to limit in formulas (25) and (26) for the expectation and variance of r.v. M when i ! 1 and s ! 0 under condition (46). By contrast to r.v. M, the limiting distribution of r.v. L under the same conditions is a Poisson distribution. Theorem 3. Let p~m , m P 0, be the distribution of r.v. L. Then under the conditions of Theorem 2, ~m

~A

p~m ! eÿA

m!

;

m P 0;

h where A~ ˆ nÿ1 : a

…60†

Proof. Proceeding from (21) we easily ®nd that  i aln …x ÿ 1† nÿ1 ! eÿh=a : p~0 ˆ 1 ‡ n 1 ÿ x ÿ s ‡ al

…61†

Next, it follows from (20), (27) and (46) that sÿ1

1 b~ ! qÿ1 a~

and s

iDL h…q ÿ 1† ! : ~ anÿ1 b~ c

Passage to limit as s ! 0 and i ! 1 in the equality !m   !m k   m  kÿ1  X ~ DL j iDL b~ b m ÿ 1 smÿk Y ÿ1 Qm 1‡ ; ˆ s s ~c ~c mÿk k! jˆ0 i a~ a~ b~ b~ kˆ1

…62†

m P 1;

with (62) taken into account leads together with (61) to the desired conclusion (60). Remark 5. It follows from (28) and (29) that under conditions of Theorem 2, EL and Var L ~ that is, to the mean and variance of the Poisson distribution P …A†. ~ converge to A,

L.G. Hanin / Mathematical Biosciences 169 (2001) 89±107

103

5. Numerical example Theorems 2 and 3 lead to the following natural question: What are the rates of convergence of the distributions of r.v.s M and L to the limiting distributions (47) and (60), respectively? Sample calculations described below were designed to shed some light on this problem for r.v. L. A more thorough theoretical and numerical investigation will be a matter of another publication. It follows from Theorem 3 that the distribution PL of r.v. L can be approximated by the Poisson ~ with distribution P …h† h~ :ˆ isn =anÿ1 :

…63†

These two distributions were compared graphically and by means of computing the probability metric m ~ h ~ ~ :ˆ sup Pr…L ˆ m† ÿ eÿh …64† d…PL ; P …h†† m! mP0 in the case when no cell death occurs between exposures (m ˆ 0) and for the number of fractions n ˆ 30. The common expected value h~ of the two distributions was taken to be equal to 2; in other words, it was assumed that on the average two clonogens are alive after delivery of the last fraction of radiation. ~ were graphed for a ˆ 0:9 and i ˆ 105 ; 107 . The corresponding The distributions PL and P …h† value s of the survival probability per fraction was computed from (63). The exact distributions of ~ are represented on Fig. 1 by lines 1 r.v. L for i ˆ 105 ; 107 and the Poisson approximation P …h† (solid), 2 (dotted) and 3 (dashed), respectively. For easier visualization, all three discrete distri~ provides a butions were smoothed. Fig. 1 shows that, for large i, the Poisson distribution P …h† good ®t to the distribution PL . The distance (64) was computed for a ˆ 0:7; 0:8; 0:9; 0:95 and i ˆ 10N with N ˆ 5; 6; 7; 8. Re~ is sults of this calculation given in Table 1 con®rm that, for large i, the Poisson distribution P …h† close to PL . Furthermore, we conclude that the distance (64) monotonically tends to 0 when i ! 1 for ®xed a and monotonically decreases as a ! 1 for ®xed i. Observe that in the limiting case a ˆ 1 …k ˆ 0†, that is, in the absence of cell proliferation and spontaneous cell death, (64) represents the distance between the binomial distributions B…i; sn † and its Poisson approximation P …h†; h ˆ isn . According to the Law of Rare Events (see e.g. [14, pp. 279±281]), in this case d…PL ; P …h†† 6 is2n ˆ h~2 =i. Therefore, the distance (64) between the generalized negative binomial ~ for a < 1 is larger than between the corredistribution PL and the Poisson distribution P …h† n sponding binomial distribution B…i; s † and P …h† in the limiting case a ˆ 1. Table 1 ~ ˆ 30; m ˆ 0; h~ ˆ 2† The distance between the exact distribution PL and the approximating Poisson distribution P …h†…n i

a ˆ 0:7

a ˆ 0:8

a ˆ 0:9

a ˆ 0:95

105 106 107 108

0.177 0.143 0.119 0.097

0.122 0.098 0.080 0.067

0.062 0.049 0.040 0.033

0.031 0.025 0.020 0.019

104

L.G. Hanin / Mathematical Biosciences 169 (2001) 89±107

Fig. 1. Exact distribution PL for i ˆ 105 (solid line 1), i ˆ 107 (dotted line 2), and the approximating Poisson distri~ (dashed line 3) …n ˆ 30; m ˆ 0; a ˆ 0:9; h~ ˆ 2† bution P …h†

6. Discussion Our results solve the main problem raised in [1] and followed up in [6,11]. According to Theorem 1, the exact distribution of the number L of surviving clonogens belongs to the family (44) of generalized negative binomial distributions. This family does not contain binomial distributions. However, as discovered in Theorem 3, the distribution of r.v. L still converges under the natural condition (46) to the Poisson distribution (60). Therefore, the distribution of r.v. L can ~ with h~ ˆ isn =anÿ1 . Note that h~ is not equal to be approximated by the Poisson distribution P …h† n h ˆ is unless a ˆ 1, that is, k ˆ m. This explains why the Poisson distribution P …h† was rejected in [1,6,11] on the basis of numerical computations which made the authors think that any Poisson approximation is impossible. We observe also that the Poisson approximation P …h† is valid not only in the absence of cell proliferation and death (k ˆ m ˆ 0) but also in the more general case when their rates are equal (k ˆ m). It is interesting to mention that, although the distributions of r.v.s M and L are quite close when the interdose interval s is small, the limiting distribution (47) of r.v. M under the condition (46) is not Poisson and never reduces to it. The model of fractionated radiation cell survival discussed in this work depends on six parameters n; s; i; s; k; m, the last four of which are unobservable. More importantly, the number of surviving clonogens is unobservable as well. To make use of the computed distribution of the size of the surviving fraction of clonogens, one has to relate it to an observable endpoint, like the time to tumor recurrence or death caused by a speci®c recurrent cancer. In what follows, we will

L.G. Hanin / Mathematical Biosciences 169 (2001) 89±107

105

focus on the ®rst possibility. Speci®cally, let T be the time to tumor recurrence counted from the moment of delivery of the last fraction. It is assumed that r.v. T is absolutely continuous. There are two competing models of post-treatment tumor development. According to the clonogenic model introduced in [21], a recurrent tumor arises from a single clonogenic cell. Every surviving clonogen can be characterized by a latent time (called the progression time) during which it could potentially propagate into a detectable tumor. It is assumed additionally that progression times of surviving clonogens are i.i.d. r.v.s. Suppose that the number L of surviving clonogens is equal to m. Then for the observed time T to tumor recurrence we have T ˆ min1 6 j 6 m Tj ; where Tj P m ~ p is the progression time of the jth clonogen. Let uL …z† ˆ 1 mˆ0 m z be the p.g.f. of r.v. L: It follows from the assumptions of the clonogenic model that the conditional survival function FT jLˆm …t† :ˆ Pr…T > t j L ˆ m†; t P 0; of r.v. T is given by FT jLˆm ˆ F1m ; where F1 is the common survival function of the progression times of surviving clonogens. Then FT …t† ˆ

1 X

FT jLˆm …t† Pr…L ˆ m† ˆ

mˆ0

1 X mˆ0

p~m F1m …t† ˆ uL …F1 …t††;

t P 0:

…65†

Alternatively, the cumulative model proceeds from the assumption that all surviving clonogenic cells contribute to the resulting recurrent tumor, see for example [5]. Assume additionally that their contributions are identically distributed (but not necessarily independent). For m P 0; denote by rm the conditional hazard function of r.v. T given that L ˆ m: Then  Z t   FT jLˆm …t† ˆ exp ÿ rm …u† du : t P 0: 0

By our assumptions, rm ˆ mr; where r is the common hazard function associated with the contributions of the surviving clonogens to the resulting detectable tumor. Then FT jLˆm ˆ F2m ; m P 0; where F2 is the survival function corresponding to the hazard rate r, and hence FT …t† ˆ

1 X mˆ0

p~m F2m …t† ˆ uL …F2 …t††;

t P 0:

…66†

Comparing (65) and (66) we conclude that, surprisingly, the distribution of the time T to tumor recurrence within both models of post-treatment tumor development has the same form FT …t† ˆ uL …F…t††;

t P 0;

…67†

where F is the survival function of an absolutely continuous non-negative r.v. Suppose, in particular, that r.v. L follows the Poisson distribution P …A†. Then we derive easily from (67) that FT …t† ˆ eÿAF …t† ;

t P 0;

…68†

where F :ˆ 1 ÿ F is the corresponding cumulative distribution function. Although this form of the recurrence time distribution was originally suggested in [21] in conjunction with the clonogenic model of post-treatment tumor development, it proved to be instrumental for statistical estimation of the probability p of tumor cure (also termed the cure rate) in many other settings [7,22±26]. Observe that according to (68), p ˆ eÿA ; and the positivity of p is due to the fact that p~0 ˆ Pr…L ˆ 0† > 0:

106

L.G. Hanin / Mathematical Biosciences 169 (2001) 89±107

Signi®cance of formula (67) is two-fold. First, it suggests that knowledge of the entire distribution of the number of surviving clonogens (not only of the tumor control probability p~0 ) is critical for developing biologically motivated post-treatment survival models. Second, assuming some parametric form for the function F and using for the distribution of r.v. L either its exact form (44) or the Poisson approximation (60) one can in principle estimate unknown quantities i; s and especially the all-important kinetic parameters k and m from the observed times to tumor recurrence. Making use of the exact distribution of r.v. L seems to be especially promising because it incorporates the known dose D per fraction through survival probability s, and this allows to develop a regression model accounting for the variation of the dose D and other clinically important covariates. Pursuing this line of reasoning presupposes, indeed, investigation into identi®ability properties of the resulting survival model.

Acknowledgements The author is thankful to Drs A.Yu. Yakovlev (University of Utah) and M. Zaider (Memorial Sloan-Kettering Cancer Center) for their insightful comments on various issues raised in this paper, and to Dr A. Zorin for his help in producing Fig. 1. The author is also grateful to the referees for their valuable suggestions.

References [1] S.L. Tucker, H.D. Thames, J.M.G. Taylor, How well is the probability of tumor cure after fractionated irradiation described by Poisson statistics? Radiat. Res. 124 (1990) 273. [2] W.S. Kendal, A closed-form description of tumour control with fractionated radiotherapy and repopulation, Int. J. Radiat. Biol. 73 (1998) 207. [3] T.E. Harris, The Theory of Branching Processes, Springer, Berlin, 1963. [4] P. Jagers, Branching Processes with Biological Applications, Wiley, London, 1975. [5] M. 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. [6] S.L. Tucker, Modelling the probability of tumor cure after fractionated radiotherapy, in: M.A. Horn, G. Simonett, G. Webb (Eds.), Mathematical Models in Medical and Health Sciences, Vanderbilt University, Nashville, 1999, p. 1. [7] A.Y. Yakovlev, A.D. Tsodikov, Stochastic Models of Tumor Latency and their Biostatistical Applications, World Scienti®c, Singapore, 1996. [8] J. Deasy, Poisson formulas for tumor control probability with clonogen proliferation, Radiat. Res. 145 (1996) 382. [9] A.Yu. Yakovlev, Comments on the distribution of clonogens in irradiated tumors (Letter to the Editor), Radiat. Res. 134 (1993) 117. [10] L.G. Hanin, A.Y. Yakovlev, L.V. Pavlova, Biomathematical Problems in Optimization of Cancer Radiotherapy, CRC, Boca Raton, 1994. [11] S.L. Tucker, J.M.G. Taylor, Improved models of tumour cure, Int. J. Radiat. Biol. 70 (1996) 539. [12] L.G. Hanin, M. Zaider, A.Yu. Yakovlev, Distribution of the number of clonogens surviving fractionated radiotherapy: a long-standing problem revisited, Int. J. Radiat. Biol., to appear. [13] M. Zaider, L.G. Hanin, A.D. Tsodikov, M.J. Zelefsky, S.A. Leibel, A.Yu. Yakovlev, A survival model for fractionated radiotherapy with an application to prostate cancer, preprint. [14] H.M. Taylor, S. Karlin, An Introduction to Stochastic Modeling, 3rd Ed., Academic Press, San Diego, 1998.

L.G. Hanin / Mathematical Biosciences 169 (2001) 89±107 [15] [16] [17] [18] [19] [20] [21]

[22] [23] [24] [25] [26]

107

D.G. Kendall, On the generalized birth and death process, Ann. Math. Stat. 19 (1948) 1. D.R. Cox, H.D. Miller, The Theory of Stochastic Processes, Chapman and Hall, New York, 1984. S. Karlin, A First Course In Stochastic Processes, Academic Press, New York, 1966. K.B. Athreya, P.E. Ney, Branching Processes, Springer, New York, 1972. S. Asmussen, H. Hering, Branching Processes, Birkh auser, Boston, 1983. M. Zaider, G.N. Minerbo, Tumor control probability: a formulation applicable to any protocol of dose delivery, Phys. Med. Biol. 45 (2000) 279. A.Yu. Yakovlev, B. Asselain, V.-J. Bardou, A. Fourquet, T. Hoang, A. Rochefordiere, A.D. Tsodikov, A simple stochastic model of tumor recurrence and its application to data on premenopausal breast cancer. in: B. Asselain, M. Boniface, C. Duby, C. Lopez, J.-P. Masson, J. Tranchefort (Eds.), Biometrie et Analyse de Donnees SpatioTemporelles, vol. 12, Rennes, France: Societe Francßaise de Biometrie, ENSA, 1993, p. 66. A.D. Tsodikov, B. Asselain, A. Fourquet, T. Hoang, A..Y. Yakovlev, Discrete strategies of cancer post-treatment surveillance. Estimation and optimization problems, Biometrics 51 (1995) 437. B. Asselain, A. Fourquet, T. Hoang, A.D. Tsodikov, A.Y. Yakovlev, A parametric regression model of tumor recurrence: an application to the analysis of clinical data on breast cancer, Statist. Prob. Lett. 29 (1996) 271. A.D. Tsodikov, A proportional hazards model taking account of long-term survivors, Biometrics 54 (1998) 138. A.D. Tsodikov, Asymptotic eciency of a proportional hazards model with cure, Stat. Prob. Lett. 39 (1998) 237. A.D. Tsodikov, M. Loe‚er, A.Yu. Yakovlev, A cure model with time-changing risk factor: an application to the analysis of secondary leukemia. A report from the International Database on Hodgkin's disease, Stat. Medicine 17 (1998) 27.