A stochastic model of tumor response to fractionated radiation: limit theorems and rate of convergence

Mathematical Biosciences 191 (2004) 1–17 www.elsevier.com/locate/mbs

A stochastic model of tumor response to fractionated radiation: limit theorems and rate of convergence Leonid G. Hanin

a,b,*

a

b

Department of Mathematics, Idaho State University, Pocatello, ID 83209-8085, USA Department of Biostatistics and Computational Biology, University of Rochester, Rochester, NY 14642, USA Received 12 February 2003; received in revised form 13 January 2004; accepted 29 April 2004 Available online 2 July 2004

Abstract The iterated birth and death Markov process is defined as an n-fold iteration of a birth and death Markov process describing kinetics of certain population combined with random killing of individuals in the population at moments s1 ; . . . ; sn with given survival probabilities s1 ; . . . ; sn . A long-standing problem of computing the distribution of the number of clonogenic tumor cells surviving an arbitrary fractionated radiation schedule is solved within the framework of iterated birth and death Markov process. It is shown that, for any initial population size i, the distribution of the size N of the population at moment t P sn is generalized negative binomial, and an explicit computationally feasible formula for the latter is found. It is shown that if i ! 1 and sn ! 0 so that the product is1    sn tends to a finite positive limit, the distribution of random variable N converges to a probability distribution, which for t ¼ sn turns out to be Poisson. In the latter case, an estimate of the rate of convergence in the total variation metric similar to the classical Law of Rare Events is obtained.
2004 Elsevier Inc. All rights reserved. Keywords: Clonogenic tumor cell; Fractionated irradiation; Generalized geometric distribution; Generalized negative binomial distribution; Iterated birth and death Markov process; Law of Rare Events; Limiting distribution; Poisson distribution; Probability generating function; Rate of convergence

*

Address: Department of Mathematics, Idaho State University, Pocatello, ID 83209-8085, USA. Tel.: +1-208 282 3293; fax: +1-208 282 2636. E-mail address: [email protected] (L.G. Hanin). 0025-5564/$ - see front matter
2004 Elsevier Inc. All rights reserved. doi:10.1016/j.mbs.2004.04.003

2

L.G. Hanin / Mathematical Biosciences 191 (2004) 1–17

1. Introduction The present work solves, under realistic biological assumptions, the following long-standing problem in radiation oncology, see e.g. [1]: To find the distribution of the number of clonogenic tumor cells surviving a given schedule of fractionated irradiation. Solving this problem is critical for developing quantitative methods of assessment of the efficacy of radiation cancer treatment and designing optimal schedules of cancer radiotherapy. The clonogenic model of post-treatment tumor development introduced in [2] provides a link between the distribution of the number of surviving clonogenic cells and the distribution of the observed time to tumor recurrence. This makes solving the problem at hand a crucial step in developing biologically motivated models of post-treatment survival. For an at length discussion of biomedical significance of the above problem and previously developed approaches to its solution, see [3–5] and references therein. We will proceed from the following widely accepted model of tumor population kinetics, see e.g. [6,7]. A tumor initially comprising a non-random number i of clonogenic cells is exposed to a fractionated radiation schedule consisting of n instantaneously delivered doses D1 ; D2 ; . . . ; Dn administered at times s1 ; s2 ; . . . ; sn , where 0 6 s1 < s2 <    < sn . It is assumed that every clonogen survives exposure to the dose Dk with the same probability sk , 1 6 k 6 n, given that it survived the previous exposures, and independently of other clonogens. Observe that it is not necessarily supposed that all tumor cells are clonogenic. We also assume that irradiated tumor cells are killed instantaneously at the moment of exposure to radiation. According to a convention commonly accepted in radiation biology, a tumor cell is killed if it is incapable of producing a viable clone. Between the exposures, surviving clonogens proliferate or die spontaneously independently of each other. The corresponding time-dependent rates of their birth bðtÞ and spontaneous death dðtÞ (defined in the sense of Markov processes), may depend on the sequence of doses D1 ; . . . ; Dn and are assumed to be non-negative and locally integrable. When necessary, these functions will be thought of as defined beyond the interval ½s1 ; sn ; their values for t < s1 and t > sn will then represent the respective birth and death rates before and after irradiation. Spontaneous death of clonogens may be due to any mechanism of their removal from the population of clonogenic cells including necrosis, apoptosis, and entering differentiation pathway. Finally, it is assumed that spontaneous cell death between radiation exposures is an instantaneous event, and that both descendants of a clonogenic cell are clonogenic. According to the model, the size of the population of clonogenic cells between exposures follows birth and death Markov process, and at the moment of exposure to dose Dk is subject to random killing with the survival probability sk , 1 6 k 6 n. The combined stochastic process will be called in what follows iterated birth and death Markov process, compare with [3]. It is well-known that tumor response to radiation depends critically on intracellular damage repair processes. They operate on the time scale of minutes to hours while, in the case of fractionated radiation, the inter-dose intervals typically range from one to several days. Such time intervals are long enough for the accomplishment of transient processes of inactivation and recovery of damaged cells, so that survival probability of an irradiated cell depends, conditional on its survival under previous exposures, only on the radiation dose D thus making a single survival probability s ¼ sðDÞ accountable for the resultant effect of damage repair processes.

L.G. Hanin / Mathematical Biosciences 191 (2004) 1–17

3

In the above model, the population of tumor cells is thought of as homogeneous with regard to both sensitivity of clonogenic cells to the killing action of radiation and to the pattern of their spontaneous death and proliferation. A more general case when the population of clonogenic tumor cells comprises several sub-populations with different radiosensitivities s1 ; . . . ; sn and/or different birth and spontaneous death rates bðtÞ and dðtÞ can also be easily accommodated: since these sub-populations evolve independently, the distribution of the total size of the surviving population of clonogens is the convolution of the corresponding distributions for its homogeneous sub-populations. A detailed discussion of the effects of tumor population heterogeneity on the probability of tumor extinction at the end of treatment (also termed tumor control probability, TCP) is given in [7]. In the particular case where doses of radiation and inter-dose intervals are equal (D1 ¼ D2 ¼    ¼ Dn ¼ D and sk ¼ ðk 1Þs, 1 6 k 6 n), and birth and death rates b and d are constant, the problem at hand was solved within an approach based on iterated birth and death Markov processes in [3]. It was shown that the distribution of the number N ðtÞ of surviving clonogenic cells at times t ¼ sn ¼ ðn 1Þs and t ¼ ns belongs to the family of generalized negative binomial distributions (see Section 3), and an explicit computationally feasible formula for the distributions from this family was obtained. A survival model based on the findings of [3] was successfully applied to statistical analysis of data on post-treatment recurrence of prostate cancer [5]. For a detailed discussion of a wide range of biomedical and statistical implications of the work [3], the reader is referred to [4]. The points made there promise even a larger impact for a far more general and realistic model studied in the present work. Application of this mechanistic model of tumor response to radiation coupled with a parametric model of post-treatment survival may allow for estimation of some of the unobservable kinetic parameters of the model from post-treatment survival data. The rationale for the extension of results obtained in [3] to arbitrary schedules of fractionated radiation is twofold. First, the most commonly used schedules of fractionated dose delivery consist of daily irradiation with equal or variable (usually escalating) doses on business days followed by weekend breaks. Second, the search for optimal schedules of fractionated irradiation presupposes variation of fractional doses and inter-dose intervals. Furthermore, consideration of variable birth and death rates is equally important for the following reasons: (1) A large fraction of tumors is detected and treated beyond the initial exponential phase of their progression. For such tumors, the birth and death rates depend on tumor size and thereby are functions of time. (2) It is well-known that exposure to ionizing radiation induces blocking of irradiated cells in various (most notably, G2 and M) phases of the mitotic cycle, see e.g. [8–10]. After remaining dormant for some time, such cells either continue to proliferate or disintegrate. This causes prolongation of the life cycle of surviving cells and leads to a complex dependence of birth and spontaneous death rates on time. (3) During the first weeks of radiation treatment the effective clonogen doubling time is relatively long after which it becomes much shorter, see [11] and references therein. Direct application of the mathematical techniques developed in [3] for obtaining an explicit formula for the probability generating function (p.g.f.) of the number of surviving clonogens to

4

L.G. Hanin / Mathematical Biosciences 191 (2004) 1–17

arbitrary schedules of fractionated irradiation and time-dependent birth and death rates leads to prohibitive mathematical difficulties. Because of this, in the present work we take up a different approach to studying compositions of fractional linear p.g.f.s. Using this approach, we will show that in the general case the distribution of the state N ðtÞ of iterated birth and death Markov process at time t P sn is still generalized negative binomial, and compute its parameters (Theorem 1). The results of [3] will then lead to a computationally feasible formula for the distribution PrðN ðtÞ ¼ kÞ, k P 0, of random variable (r.v.) N ðtÞ. Even in the simplest case of homogeneous iterated birth and death Markov process, our model depends on the unknown and unobservable parameters i, b, d as well as on parameters built into the dose-response function s ¼ sðDÞ. Estimation of these parameters from remote endpoints such as time to tumor recurrence or cancer specific death requires an adequate model that relates the time to endpoint to the distribution of the size of surviving population of clonogens [3,4]. Such estimation problem does not seem feasible unless the sample size is extremely large. This calls for asymptotic formulas for the distribution of the number of surviving clonogens that would make computation of the distribution easier and solving associated statistical problems of parameter estimation more tractable. The initial number i of clonogenic cells in a tumor is typically very large (1 cm3 of a solid tumor contains about 109 cells [12]; also, a clinically detectable tumor is estimated to contain at least 105 clonogenic cells probably ranging up to 109 cells or even more [11,13]). Additionally, the number n of exposures (typically ranging from 20 to 40) and the doses D1 ; . . . ; Dn (usually from 1 to 2 Gy) are selected in such a way as to ensure that the overall survival probability of a clonogenic cell, which equals sn :¼ s1 s2    sn , is very small. In the case of homogeneous iterated birth and death Markov process with equal doses and equal inter-dose intervals, such limiting distribution was found in [3]. Specifically, it was shown that as i ! 1 and sn ! 0, where s ¼ sðDÞ is the cell survival probability per fraction of radiation, in such a way that isn ! h for some h with 0 < h < 1, the distribution of r.v. N ðsn Þ converges to a Poisson distribution P ðAÞ with parameter A ¼ h eðb dÞsn . This invites obtaining, in the spirit of the classical Law of Rare Events (see e.g. [14, pp. 279–289]), a theoretical estimate of the total variation distance between the distribution of r.v. Nðsn Þ and the Poisson distribution P ðisn eðb dÞsn Þ. A few sample values of this distance were computed numerically in [3]. In the present work, we extend the classical Law of Rare Events, that provides an estimate for the total variation distance between a binomial distribution and the corresponding Poisson distribution, to the family of generalized negative binomial distributions (Theorem 2). This will allow us to identify conditions under which a sequence of generalized negative binomial distributions has a Poisson limit (Theorem 3) and apply this result to the distribution of r.v. N ðsn Þ. A distinct advantage of the Poisson approximation is that it reduces significantly the number of model parameters to be estimated from clinical or epidemiological data. Finally, in Theorem 4 we identify another family of probability distributions that are limiting for generalized negative binomial distributions. In particular, this applies to the distribution of r.v. N ðtÞ for t > sn . In what follows we will be constantly referring to a few basic facts about branching stochastic processes, in particular, birth and death Markov processes, and their p.g.f.s. Refs. [14, Chapter 3, Sections 8–9, and Chapter 6, Section 3], [15, Chapter 1, Sections 1–7] and [16, Chapter 7] serve as an ample source of information on these topics.

L.G. Hanin / Mathematical Biosciences 191 (2004) 1–17

5

2. Birth and death Markov process and generalized geometric distribution Consider a cell population that starts at time t ¼ t0 from a single cell and which kinetics is governed by a birth and death Markov process with birth Rand death rates bðtÞ and dðtÞ, t respectively. Let d :¼ b d be the net birth rate and DðtÞ ¼ t0 dðuÞ du be the corresponding cumulative rate. It was shown by Kendall [17] (see also [7] and [15, Chapter 5, Section 7]) that p.g.f. of the size N ðtÞ of the population at time t P t0 equals uðzÞ ¼ 1

að1 zÞ ; 1 þ bð1 zÞ

ð1Þ

where a and b are functions of t given by Z t DðtÞ DðtÞ aðtÞ ¼ e and bðtÞ ¼ e bðuÞ e DðuÞ du;

t P t0 :

ð2Þ

t0

Observe that bP0

and 0 < a 6 1 þ b:

ð3Þ

In fact, the only inequality that requires a proof is a 6 1 þ b. Using (2) we rewrite it in the form Z t DðtÞ e þ bðuÞ e DðuÞ du P 1; t P t0 : ð4Þ t0

For t ¼ t0 , this is obvious. Next, the derivative of the left-hand side of (4) for almost all t P t0 equals e DðtÞ dðtÞ þ bðtÞ e DðtÞ ¼ dðtÞ e DðtÞ P 0, which completes the proof. Formula (1) can be rewritten in the form m 1 1  X a az a a b uðzÞ ¼ 1 zm : þ ¼1 þ 1 þ b ð1 þ bÞð1 þ b bzÞ 1 þ b ð1 þ bÞ2 m¼1 1 þ b Therefore, for any numbers a and b that satisfy conditions (3), formula (1) represents p.g.f. of a non-negative integer-valued r.v. X such that PrðX ¼ 0Þ ¼ 1 r

and PrðX ¼ mÞ ¼ rð1 qÞqm 1 ; m P 1;

ð5Þ

where r¼

a 1þb

and q ¼

b : 1þb

ð6Þ

Clearly, 0 < r 6 1 and 0 6 q < 1. Conversely, p.g.f. of every r.v. X , that has distribution (5) with parameters r 2 ð0; 1 and q 2 ½0; 1Þ, equals 1 X 1 r ðq rÞz ; ð7Þ qm 1 zm ¼ uX ðzÞ ¼ 1 r þ rð1 qÞ 1 qz m¼1 and hence can be represented in the form (1), where parameters a and b are given by r q and b ¼ ; a¼ 1 q 1 q and satisfy conditions (3).

6

L.G. Hanin / Mathematical Biosciences 191 (2004) 1–17

The distribution (5) will be referred to as generalized geometric distribution and denoted Gðr; qÞ. Thus, formula (1) with parameters a, b subject to (3) gives, along with (7), a general form of p.g.f. of generalized geometric distribution. In particular, the state at any given time of a birth and death Markov process starting from a single cell follows a generalized geometric distribution. Two particular cases of the generalized geometric distribution are worth mentioning: r ¼ 1 represents a plain geometric distribution GðqÞ and corresponds to a pure birth Markov process, while q ¼ 0 leads to a Bernoulli r.v. X that takes values 0 and 1 with probabilities 1 r and r, respectively and corresponds to a pure death Markov process. Note also that distribution Gðr; qÞ can be viewed as a mixture of d0 , the degenerate distribution at 0, and GðqÞ : Gðr; qÞ ¼ ð1 rÞd0 þ rGðqÞ. P.g.f. (7) is fractional linear, that is, has the form uðzÞ ¼

a bz : c dz

ð8Þ

We will assume that such p.g.f. u is non-degenerate in the sense that D :¼ ad bc 6¼ 0. It is easy to see that degenerate fractional linear p.g.f. corresponds to the distribution d0 . Conversely, it is straightforward to verify that every non-degenerate fractional linear p.g.f. is necessarily p.g.f. of a generalized geometric distribution, and hence is of the form (7) with r 2 ð0; 1 and q 2 ½0; 1Þ. Thus, formulas (1) and (7) provide equivalent two-parametric representations of non-degenerate fractional linear p.g.f.s (8). We will alternate freely between these three representations making use of the advantages that each of them offers. Observe that if p.g.f. of r.v. X is given by (1) then u0 ð1Þ ¼ a and u00 ð1Þ ¼ 2ab. Therefore, E X ¼ u0 ð1Þ ¼ a

and Var X ¼ u00 ð1Þ þ u0 ð1Þ ½u0 ð1Þ2 ¼ að2b a þ 1Þ;

ð9Þ

where E and Var stand for the expectation and variance, respectively.

3. Generalized negative binomial distribution If the birth and death Markov process described in Section 2 starts from a non-random number i of cells then, due to the fact that the sizes of populations emerging from initial cells are independent and identically distributed (i.i.d.), the total size N ðtÞ at time t is the sum of i i.i.d. r.v.s with generalized geometric distribution. For the plain geometric distribution GðqÞ, such sum follows negative binomial distribution NBði; qÞ. This motivates introducing the following class of probability distributions. Definition 1. Generalized negative binomial distribution NBði; r; qÞ is the distribution of the sum of i i.i.d. r.v.s having generalized geometric distribution Gðr; qÞ. If p.g.f. of the underlying generalized geometric distribution is represented in the form (1) then for p.g.f. of the sum N of its i i.i.d. copies we have 

i að1 zÞ uN ðzÞ ¼ 1 : 1 þ bð1 zÞ

ð10Þ

L.G. Hanin / Mathematical Biosciences 191 (2004) 1–17

7

In particular, from (9) we obtain E N ¼ ia

and Var N ¼ iað2b a þ 1Þ:

ð11Þ

Recall that in the case of birth and death Markov process fNðtÞgt P t0 with Nðt0 Þ ¼ i, parameters a and b are specified in (2). Explicit formulas for the generalized negative binomial distribution were obtained in [3]. Let X be a r.v. with such distribution. Denote pk :¼ PrðX ¼ kÞ, k P 0, and let  i a bz uX ðzÞ ¼ c dz be p.g.f. of r.v. X . We have p0 ¼ uX ð0Þ ¼ ða=cÞi . Raising the formula   a bz b D 1 ¼ 1  c dz d bd z c=d to power i and differentiating k times at z ¼ 0 we find that  j  i  k X i   b d D i kþj 1 ; k P 1: pk ¼ j k d c bc j¼1 For large i, this formula is barely computationally feasible and also untoward for the purpose of asymptotic analysis as i ! 1. A much more convenient alternative formula was derived in [3, Theorem 1]. Specifically, it was shown that  a i  b k  D  Qk ; k P 0; ð12Þ pk ¼ c a bc where Q0 ðxÞ ¼ 1, and polynomials Qk for k P 1 are given by   k  X k 1 iþj 1 j x ; k P 1: Qk ðxÞ ¼ k j j j¼1 Representing the underlying p.g.f. (8) in the form (7) we rewrite (12) as follows:  k  rð1 qÞ  i q r pk ¼ ð1 rÞ Qk ; k P 0: 1 r q r

ð13Þ

ð14Þ

Alternatively, the distribution of r.v. X can be expressed through (6) in terms of parameters a and b:  i  k   a b a a Qk ; k P 0: ð15Þ pk ¼ 1 1þb 1þb a ð1 þ bÞðb aÞ

4. Composition of fractional linear p.g.f.s Unlike formulas (7) and (8), parameterization (1) of fractional linear p.g.f.s lends itself easily to computing compositions. If p.g.f.s u1 and u2 are given by

8

L.G. Hanin / Mathematical Biosciences 191 (2004) 1–17

uj ðzÞ ¼ 1

aj ð1 zÞ ; 1 þ bj ð1 zÞ

j ¼ 1; 2;

then a straightforward calculation yields ðu1  u2 ÞðzÞ ¼ 1

a1 a2 ð1 zÞ : 1 þ ða2 b1 þ b2 Þð1 zÞ

ð16Þ

In particular, if wðzÞ ¼ 1 sð1 zÞ is p.g.f. of a Bernoulli r.v. with parameter s, 0 6 s 6 1, and u is given by (1) then ðu  wÞðzÞ ¼ 1

asð1 zÞ : 1 þ bsð1 zÞ

ð17Þ

For parameters of p.g.f. (16) we have in view of (3): a1 a2 > 0, a2 b1 þ b2 P 0, and a1 a2 6 ð1 þ b1 Þa2 6 1 þ a2 b1 þ b2 : Thus the set of p.g.f.s defined in (1), (3) forms a semigroup under composition which is isomorphic to the semigroup of all pairs ða; bÞ satisfying conditions (3) with the semigroup operation  given by ða1 ; b1 Þ  ða2 ; b2 Þ ¼ ða1 a2 ; a2 b1 þ b2 Þ. In other words, the set of non-degenerate fractional linear p.g.f.s is closed under composition. Therefore, the composition of any number of p.g.f.s of generalized geometric distributions belongs to the same class. An immediate extension of (16) by induction leads to the following explicit formula for the composition of any number of fractional-linear p.g.f.s (or equivalently, for the product of any number of 2 · 2 matrices with equal sums of entries in the two rows). Proposition 1. Let uj ðzÞ ¼ 1

aj ð1 zÞ ; 1 þ bj ð1 zÞ

1 6 j 6 n:

Then ðu1  u2      un ÞðzÞ ¼ 1

a1  a2    an ð1 zÞ ; 1 þ ða1 b1 þ a2 b2 þ    þ an bn Þð1 zÞ

where ak :¼ akþ1    an ; 1 6 k 6 n 1;

and

an :¼ 1:

Remark. For the n-fold composition uðnÞ of p.g.f. (1) with itself we have uðnÞ ðzÞ ¼ 1

an ð1 zÞ n

1 þ bðaa 1 1Þ ð1 zÞ

;

ð18Þ

with an obvious modification for the case a ¼ 1. A formula equivalent to (18) was obtained in [3] in the context of iterated homogeneous birth and death Markov processes through a different parameterization of fractional linear p.g.f. in which one parameter was defined as the non-trivial

L.G. Hanin / Mathematical Biosciences 191 (2004) 1–17

9

fixed point of the p.g.f. and the other as the expected value of the underlying r.v., see also [7], [15, p. 9], [16, pp. 296–298], [18, p. 7] and [19, p. 57]. 5. Distribution of the state of iterated birth and death Markov process In this section we identify the distribution of the state at any time t P sn of the iterated birth and death Markov process described in Section 1. We begin with a general argument that will show that this distribution is generalized negative binomial. Consider iterated birth and death Markov process that starts at t ¼ 0 from i clonogenic cells, as described in Section 1. Disregarding exposures to radiation at times s1 ; . . . ; sn denote by uk p.g.f. of the number of clonogens at time sk given that only one clonogen is present at time sk 1 , 1 6 k 6 n, where s0 :¼ 0. Similarly, let ft be p.g.f. of the number of clonogens at time t given one clonogen at time sn . Finally, denote by wk ðzÞ ¼ 1 sk ð1 zÞ p.g.f. of the Bernoulli distribution that describes random killing of clonogenic cells with the survival probability sk , 1 6 k 6 n. A standard branching process argument leads to the following formula for p.g.f. Ut of the state N ðtÞ at time t P sn of the resulting iterated birth and death Markov process: i

Ut ¼ ðu1  w1  u2  w2      un  wn  ft Þ : P.g.f.s uk , wk , 1 6 k 6 n, and ft are fractional linear and non-degenerate. Therefore, as shown in Section 4, r.v. NðtÞ follows a generalized negative binomial distribution. An explicit formula for its p.g.f. Ut is given below. Theorem 1 2 Ut ðzÞ ¼ 41 where DðtÞ ¼

Rt 0



1 þ eDðtÞ r1

R s1 0

f ðuÞdu þ r2

R s2 s1

s1 s2  sn e

3i

DðtÞ

ð1 zÞ 5; R sn Rt f ðuÞdu þ   þ rn sn 1 f ðuÞdu þ sn f ðuÞdu ð1 zÞ

dðuÞ du, f ðuÞ ¼ bðuÞ e DðuÞ , and rj :¼ sj  sjþ1    sn , 1 6 j 6 n.

Proof. Indeed, it is enough to consider the case i ¼ 1. Arguing by induction in n we first prove the required formula for t ¼ sn , that is, show that p.g.f. of the number of surviving clonogens immediately after the nth exposure is given by Usn ðzÞ ¼ 1

h

1 þ eDðsn Þ r1

R s1 0

s1 s2    sn eDðsn Þ ð1 zÞ i : Rs R sn f ðuÞ du þ r2 s12 f ðuÞ du þ    þ rn sn 1 f ðuÞ du ð1 zÞ

ð19Þ

For n ¼ 1, this is a combination of (17) and (2). Suppose that for some k, 1 6 k 6 n 1, we have Usk ðzÞ ¼ 1 where Ak ¼ rk;1

Z

s1 s2    sk eDðsk Þ ð1 zÞ ; 1 þ Ak eDðsk Þ ð1 zÞ

s1

f ðuÞ du þ rk;2 0

Z

s2

f ðuÞ du þ    þ rk;k s1

ð20Þ Z

sk

f ðuÞ du; sk 1

ð21Þ

10

L.G. Hanin / Mathematical Biosciences 191 (2004) 1–17

and rk;j :¼ sj  sjþ1    sk , 1 6 j 6 k. Then according to formulas (20), (16) and (2) R skþ1 dðuÞ du Dðsk Þ s1 s2    sk e e sk ð1 zÞ  R skþ1 R skþ1  Ru ðUsk  ukþ1 ÞðzÞ ¼ 1 R dðuÞ du dðuÞ du skþ1 dðvÞ dv Dðs Þ sk sk sk k e þe bðuÞ e du ð1 zÞ 1 þ Ak e sk ¼1

s1 s2    sk eDðskþ1 Þ ð1 zÞ h i : R skþ1 Dðs Þ kþ1 1þe Ak þ sk f ðuÞ du ð1 zÞ

Hence, from (17) and (21) we conclude that p.g.f. Uskþ1 ¼ Usk  ukþ1  wkþ1 has the desired representation (19) with n ¼ k þ 1. This establishes (19) for arbitrary n. Finally, the argument used above for computing the composition Usk  ukþ1 applies equally to the composition Ut ¼ Usn  ft , and yields the required formula for Ut . Theorem 1 is proved. h As shown in Theorem 1, p.g.f. Ut has the form (10) with a ¼ s1 s2    sn eDðtÞ

ð22Þ

and b¼e

DðtÞ



Z

s1

r1 0

f ðuÞ du þ r2

Z

s2

f ðuÞ du þ    þ rn s1

Z

sn

f ðuÞ du þ sn 1

Z

t

 f ðuÞ du ;

ð23Þ

sn

while for the corresponding extinction probability (or TCP) we have  i a p0 ¼ PrðN ðtÞ ¼ 0Þ ¼ Ut ð0Þ ¼ 1 : 1þb Furthermore, formulas (13) and (14) (or (15)) allow us to find explicitly the probabilities pk ¼ PrðN ðtÞ ¼ kÞ, k P 1, and proceeding from (11) one can compute the expectation and variance of r.v. N ðtÞ. Suppose that birth and death rates do not vary between exposures: bðtÞ ¼ bk , dðtÞ ¼ dk for sk < t < skþ1 , 0 6 k 6 n 1. Setting dk :¼ bk dk and Dsk :¼ skþ1 sk , 0 6 k 6 n 1, we find that, for t ¼ sn , Pn 1 a ¼ s1 s2    sn e k¼0 dk Dsk and b¼

Pn 1 n 1 X rkþ1 bk d Ds ð1 e dk Dsk Þ e j¼k j j : d k k¼0

Let, in particular, birth and death rates be constant: bðtÞ ¼ b and dðtÞ ¼ d, 0 6 t 6 sn . Then a ¼ s1 s2    sn edsn

n 1 X b and b ¼ edsn rkþ1 ðe dsk e dskþ1 Þ; d k¼0

where d ¼ b d. Finally, in the case where s1 ¼ s2 ¼    ¼ sn ¼ s we obtain

L.G. Hanin / Mathematical Biosciences 191 (2004) 1–17

n dsn

a¼s e

n 1 b dsn X and b ¼ e sn kþ1 ðe dsk e dskþ1 Þ: d k¼0

11

ð24Þ

6. Limit theorems and rate of convergence for generalized negative binomial distribution For probability measures P , Q on the Borel r-algebra BðRÞ, define the metric dðP ; QÞ ¼ supfjP ðAÞ QðAÞj : A 2 BðRÞg:

ð25Þ

It is easy to see that 2dðP ; QÞ equals the total variation of the signed measure P Q. The classical Law of Rare Events consists of the following estimate for the distance d between a binomial distribution Bði; pÞ and the Poisson distribution P ðipÞ: dðBði; pÞ; P ðipÞÞ 6 ip2 :

ð26Þ

In particular, this implies the following classical theorem about Poisson limit of a sequence of binomial distributions: if Xi 2 Bði; pi Þ, where ipi ! h as i ! 1, 0 < h < 1, then for each k 2 Zþ , PrðXi ¼ kÞ ! e h

hk k!

as i ! 1:

ð27Þ

The goal of this section is to extend relations (26) and (27) to generalized negative binomial distributions. Applications of these results to the iterated birth and death Markov process are discussed in the next section. We start, for the reader’s convenience, with the following well-known elementary statement. Proposition 2. Let n and g be two r.v.s, and A be a Borel subset of R. Then j Prðn 2 AÞ Prðg 2 AÞj 6 Prðn 6¼ gÞ: Proof. Consider events B ¼ ðn 2 AÞ and C ¼ ðg 2 AÞ. Clearly, for their symmetric difference we have BMC  ðn 6¼ gÞ. Therefore, j Prðn 2 AÞ Prðg 2 AÞj ¼ j PrðB n CÞ PrðC n BÞj 6 PrðB n CÞ þ PrðC n BÞ ¼ PrðBMCÞ 6 Prðn 6¼ gÞ; which completes the proof. h Theorem 2. dðNBði; r; qÞ; P ðirÞÞ 6 irðr þ qÞ. i

Proof. Consider the probability space X ¼ ½0; 1Þ with the Lebesgue measure. Points of X will be represented in the form x ¼ ðx1 ; . . . ; xi Þ. For each j ¼ 1; 2; . . . ; i, define r.v. nj : X ! R as follows: nj ðxÞ ¼ 0 if 0 6 xj < 1 r and nj ðxÞ ¼ 1 if 1 r 6 xj < 1. Clearly, Pir.v.s n1 ; n2 ; . . . ; ni are independent and have the distribution Bð1; rÞ, hence their sum n :¼ j¼1 nj follows the binomial distribution Bði; rÞ.

12

L.G. Hanin / Mathematical Biosciences 191 (2004) 1–17

Similarly, we will construct a sequence g1 ; g2 ; . . . ; gi of independent r.v.s on X that follow 1 generalized geometric distribution Gðr; qÞ. To this end, define a sequence fqm gm¼0 by setting qm :¼ 1 r þ rð1 qÞ

m X

qm 1 ¼ 1 rqm :

k¼1

For j ¼ 1; 2; . . . ; i, let gj ðxÞ ¼ 0, if 0 6 xj < q0 ¼ 1 r and gj ðxÞ ¼ m, if qm 1 6 xj < qm , m P 1. m 1 , m P 1. Note that for each j, Prðgj ¼ 0Þ ¼ 1 r and Prðgj ¼ mÞ ¼ qm qP m 1 ¼ rð1 qÞq i Therefore, in view of (5), gj 2 Gðr; qÞ, 1 6 j 6 i. Setting g :¼ j¼1 gj we conclude that g 2 NBði; r; qÞ. Observe that nj ðxÞ 6¼ gj ðxÞ if and only if q1 6 xj < 1. Hence, Prðnj 6¼ gj Þ ¼ 1 q1 ¼ rq: Therefore, applying Proposition 2 we have for all A  Zþ : j Prðn 2 AÞ Prðg 2 AÞj 6 Prðn 6¼ gÞ 6 Pr

i [

! ðnj 6¼ gj Þ 6

j¼1

i X

Prðnj 6¼ gj Þ ¼ irq:

j¼1

In view of (25) this implies that dðNBði; r; qÞ; Bði; rÞÞ 6 irq: Finally, applying the Law of Rare Events (26) we find that dðNBði; r; qÞ; P ðirÞÞ 6 dðNBði; r; qÞ; Bði; rÞÞ þ dðBði; rÞ; P ðirÞÞ 6 irq þ ir2 ¼ irðr þ qÞ: Theorem 2 is proved. h Remark 1. In the case q ¼ 0, Theorem 2 coincides with the classical Law of Rare Events (26). Remark 2. Suppose that generalized geometric distribution Gðr; qÞ is parameterized by a, b related to r, q in accordance with (6). Then the result of Theorem 2 becomes      a b ia iaða þ bÞ d NB i; : ð28Þ ; ;P 6 1þb 1þb 1þb ð1 þ bÞ2 Remark 3. In the classical Law of Rare Events (26), the binomial and the approximating Poisson distributions have the same expectation ip. By contrast, the expectation of the generalized negative binomial distribution Bði; r; qÞ, which equals irð1 qÞ, is different from the expectation of its Poisson approximation P ðirÞ. Remark 4. Theorem 2 provides the following estimate for the error in the Poisson approximation of the tumor control probability p0 : jp0 e ir j 6 irðr þ qÞ:

L.G. Hanin / Mathematical Biosciences 191 (2004) 1–17

13

An analogue of (27) follows immediately from Theorem 2. Theorem 3. Let Xi 2 NBði; ri ; qi Þ, where iri ! h, 0 < h < 1, and qi ! 0 as i ! 1. Then for each k 2 Zþ , hk as i ! 1: k! It is curious to mention that if iri ! h, 0 < h < 1, and qi ! q, 0 < q < 1, then the sequence of distributions NBði; ri ; qi Þ converges to a (proper) probability distribution different from Poisson. PrðXi ¼ kÞ ! e h

Theorem 4. Let Xi 2 NBði; ri ; qi Þ, where iri ! h, 0 < h < 1, and qi ! q, 0 < q < 1, as i ! 1. Then for each k 2 Zþ , PrðXi ¼ kÞ ! e h ck qk where c0 ¼ 1, and  k  X k 1 lj ; ck ¼ k j j! j¼1

as i ! 1;

k P 1;

ð29Þ

ð30Þ

with l ¼ hðq 1 1Þ:

ð31Þ

Furthermore, the limiting distribution in (29) is proper. Proof. We refer to formula (14). Clearly, under conditions of Theorem 4, ð1 ri Þi ! e h as i ! 1, which establishes (29) for k ¼ 0. Observe that iri ð1 qi Þ ! hðq 1 1Þ ¼ l qi ri

as i ! 1:

ð32Þ

Rewriting polynomial Qk given by (13) in the form  j 1  k  X m k 1 1 Y 1þ Qk ðxÞ ¼ ðixÞj ; k P 1; k j j! m¼0 i j¼1 setting x ¼ ri ð1 qi Þ=ðqi ri Þ, and taking into account (32) and (30) we find that   ri ð1 qi Þ Qk ! ck as i ! 1; k P 1: qi ri Finally, ðqi ri Þ=ð1 ri Þ ! q as i ! 1, which completes the proof of (29). It remains to show that the limiting distribution in (29) is proper, that is, 1 X k¼0

e h ck qk ¼ 1:

ð33Þ

14

L.G. Hanin / Mathematical Biosciences 191 (2004) 1–17

To this end, we will compute its g.f. u. We have for jzj < q 1 : "  # 1 1 k  X X X k 1 lj k k h h uðzÞ ¼ e ck ðqzÞ ¼ e 1 þ ðqzÞ k j j! k¼0 k¼1 #j¼1 "   1 1 X lj X k 1 ¼ e h 1 þ ðqzÞk k j j! k¼j j¼1 # "  1 1  X ðlqzÞj X jþm 1 m h ¼e 1þ ðqzÞ ; m j! m¼0 j¼1

ð34Þ

where m :¼ k j. Observe that the last sum in (34) represents a binomial series:  1  X jþm 1 m j ðqzÞ ¼ ð1 qzÞ ; jzj < q 1 ; j P 1: m m¼0 Therefore, using (31) we continue (34) as follows:   j  1 X 1 lqz 1 z h uðzÞ ¼ e ¼ exp h : j! 1 qz 1 qz j¼0

ð35Þ

Setting here z ¼ 1 we find that uð1Þ ¼ 1 which establishes (33). This completes the proof of Theorem 4. h Remark 1. In the limiting case q ¼ 0, Theorem 4 reduces to Theorem 3. Remark 2. Formula (35) coincides with the one obtained by passage to limit in (10) as i ! 1 with a ¼ ri =ð1 qi Þ and b ¼ qi =ð1 qi Þ.

7. Application of limit theorems to iterated birth and death Markov process It was shown in Section 5 that the distribution of the state NðtÞ at time t of iterated birth and death Markov process is generalized negative binomial. According to Theorem 1 p.g.f. Ut of r.v. NðtÞ for t P sn has the form (10) with parameters a and b specified in (22) and (23). In particular, this applies to the number N ðtÞ of clonogenic tumor cells surviving fractionated radiation at the end of treatment (t ¼ sn ) or at any time t after treatment (t > sn ), see Section 1. Applicability of the approximation and limit theorems from Section 6 to the distribution of r.v. N ðtÞ depends on whether t ¼ sn or t > sn . Below we discuss these cases separately. A general premise underlying our analysis is that i ! 1 and sn ! 0 in such a way that the product is1 s2    sn is bounded (note that the survival probabilities s1 ; . . . ; sn 1 may or may not depend on i, and recall that s1 s2    sn is the overall survival probability of a cell in the absence of cell proliferation and spontaneous death so that is1 s2    sn is the expected number of cells surviving the entire irradiation schedule). It follows from these assumptions (see (22)) that a ! 0 and ia is bounded.

L.G. Hanin / Mathematical Biosciences 191 (2004) 1–17

(1) Suppose that t ¼ sn . In this case (23) becomes  Z s1 Z s2 Z Dðsn Þ b¼e s1    sn f ðuÞ du þ s2    sn f ðuÞ du þ    þ sn 0

s1

sn

15

 f ðuÞ du ;

sn 1

where f ðuÞ ¼ bðuÞ e DðuÞ . Obviously, b ! 0. Therefore, estimate (28) of the distance d between the distribution of r.v. N ðsn Þ and the Poisson distribution suggested by Theorem 2 tends to 0. Practical efficiency of this approximation is confirmed by our numerical example in the next section. If i ! 1 and sn ! 0 in such a way that is1    sn ! c, 0 < c < 1, then conditions of Theorem 3 are met which implies that the limiting distribution for r.v. N ðsn Þ is a Poisson distribution P ðhÞ with h ¼ c eDðsn Þ . For iterated homogeneous birth and death Markov process under conditions s1 ¼ s2 ¼    ¼ sn and sk ¼ ðk 1Þs, 1 6 k 6 n, this result was obtained in [3] directly on the basis of formulas (12) and (13). Rt (2) Suppose that t > sn . Then it follows from (23) that b ! b0 :¼ eDðtÞ sn f ðuÞ du, so that in view of (6), q ! b0 =ð1 þ b0 Þ. Observe that b0 > 0 if the birth rate b does not vanish almost everywhere on the interval ½sn ; t. In this case the Poisson approximation (28) is impractical unless ia is very small. Assuming b0 > 0 we find that if i ! 1 and sn ! 0 in such a way that is1    sn ! c, 0 < c < 1, then conditions of Theorem 4 are satisfied with h ¼ c eDðtÞ =ð1 þ b0 Þ and q ¼ b0 =ð1 þ b0 Þ. Therefore, the limiting distribution for r.v. N ðtÞ, t > sn , is given by (29). A very special case of this result can be found in [3, Theorem 2], see also [4]. Specifically, it was shown there that, for iterated homogeneous birth and death Markov processes with s1 ¼ s2 ¼    ¼ sn and sk ¼ ðk 1Þs, 1 6 k 6 n, the limiting distribution for r.v. NðnsÞ belongs to the two-parametric family identified in Theorem 4.

8. Numerical example To assess how well Poisson distribution approximates generalized negative binomial distribution and see how tight the upper bound for the distance between them provided by Theorem 2 is, the following numerical experiment designed in the framework of iterated birth and death Markov process was conducted. The initial size of clonogenic tumor cell population was taken in the form i ¼ 10m , where 5 6 m 6 9. Fractionated irradiation schedule consisted of n ¼ 30 equal fractions of D ¼ 2 Gy delivered at the same time on every business day starting Monday with weekend breaks so that the total duration of radiation treatment was 39 days. The survival 2 probability s ¼ sðDÞ was assumed to follow the linear-quadratic model, that is, sðDÞ ¼ e lD mD , see e.g. [20], where l and m are positive constants that can be interpreted in terms of cell radiosensitivity and damage repair capability, respectively. The values l ¼ 0:4 Gy 1 and m ¼ 0:003 Gy 2 adopted for the computation fall into the ranges reported in radiobiological literature [7]. For such parameters l and m, the survival probability s of clonogenic tumor cells per fraction of radiation equals 0.444. As suggested in [7], spontaneous birth and death rates b and d were set as follows: b ¼ ln 2=10 days 1 , d ¼ ln 2=20 days 1 . These values correspond to the tumor doubling time of 20 days. For each value of m, the distance d between the exact distribution PN of the number N of surviving clonogens at the end of treatment given by (15) and the approximating Poisson distribution P ðhÞ with h ¼ ia=ð1 þ bÞ was computed. Note that in the case under study

16

L.G. Hanin / Mathematical Biosciences 191 (2004) 1–17

Table 1 Sample values of the distance d ¼ dðPN ; P ðhÞÞ and its upper bound d^ i

d

d^

105 106 107 108 109

5.5179 · 10 7 5.5174 · 10 6 5.5127 · 10 5 5.4653 · 10 4 5.0133 · 10 3

5.5179 · 10 7 5.5179 · 10 6 5.5179 · 10 5 5.5179 · 10 4 5.5179 · 10 3

parameters a and b of the generalized negative binomial distribution for r.v. N are given by formulas (24). The values of the distance d and its upper bound d^ :¼ iaða þ bÞ=ð1 þ bÞ2 suggested in Theorem 2 for all considered initial population sizes i are given in Table 1. Results of the numerical calculation suggest that Poisson approximation to the distribution of the number of clonogens surviving a given schedule of fractionated radiation is quite efficient in the realistic range of model parameters and that the estimate of the distance between the two distributions obtained in Theorem 2 is very close to the true value of the distance.

9. Further discussion From a formal mathematical standpoint, the approach of the present work makes it possible to obtain an explicit formula for the distribution of the number of clonogenic tumor cells surviving exposure to continuous irradiation with arbitrary dose rate. Solving this problem presupposes the knowledge of the rate hðtÞ of radiation-induced cell death which is assumed to be either a given function of time or parametrically specified as a function of the accumulated total dose. This avenue was explored in [7], where p.g.f. Ut for the number of surviving clonogens at time t P sn and the corresponding TCP were computed using the approach of Kendall [17] in the case of continuous exposure. Specifically, it was shown in [7] that " #i SðtÞ eDðtÞ ð1 zÞ ; ð36Þ Ut ðzÞ ¼ 1 Rt 1 1 þ SðtÞ eDðtÞ 0 bðuÞ½SðuÞ e DðuÞ du ð1 zÞ Rt hðuÞ du where SðtÞ ¼ e 0 is the cell survival probability at time t under continuous irradiation in the absence of cell proliferation and spontaneous cell death, compare with Theorem 1. Mathematically, the cases of continuous and fractionated radiation can be viewed as limiting cases of each other under appropriate limit procedures. However, a rigorous justification of the passage to the limit in the distribution of the number of surviving clonogens in both cases is quite challenging. Arguing heuristically, the authors of [7] obtained a formula for the TCP in the case of fractionated irradiation by equal doses separated by equal time intervals and constant birth and death rates on the basis of their formula for TCP under continuous irradiation; their result agrees with the corresponding formula in [3]. Similarly, one can formally derive from (36) an expression for the p.g.f. of the number of surviving clonogens under arbitrary schedule of fractionated irradiation. This formula coincides with the one derived in the present work (see Theorem 1) so that the latter serves as a rigorous justification of the pertinent limiting procedure. More

L.G. Hanin / Mathematical Biosciences 191 (2004) 1–17

17

importantly, though, from a biological point of view, mechanistic incorporation of the allimportant processes of damage repair in the case of continuous irradiation requires general nonMarkovian age-dependent branching stochastic processes, for which a theoretical form of the distribution of the state of the process at a given time is not available. Whether the cumulative effect of damage repair kinetics on the distribution of the number of surviving clonogens can be adequately described or approximated within the framework of non-homogeneous birth and death Markov processes still remains to be seen.

Acknowledgement Research of the author was supported by NSF grant DMS-0109895.

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] A.Y. 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, Societe Francßaise de Biometrie, ENSA, Rennes, France, 1993, p. 66. [3] L.G. Hanin, Iterated birth and death process as a model of radiation cell survival, Math. Biosci. 169 (2001) 89. [4] L.G. Hanin, M. Zaider, A.Y. Yakovlev, Distribution of the number of clonogens surviving fractionated radiotherapy: a long-standing problem revisited, Int. J. Radiat. Biol. 77 (2) (2001) 205. [5] M. Zaider, M.J. Zelefsky, L.G. Hanin, A.D. Tsodikov, A.Y. Yakovlev, S.A. Leibel, A survival model for fractionated radiotherapy with an application to prostate cancer, Phys. Med. Biol. 46 (2001) 2745. [6] W.S. Kendal, A closed-form description of tumour control with fractionated radiotherapy and repopulation, Int. J. Radiat. Biol. 73 (1998) 207. [7] M. Zaider, G.N. Minerbo, Tumor control probability: a formulation applicable to any protocol of dose delivery, Phys. Med. Biol. 45 (2000) 279. [8] S. Okada, Radiation Biochemistry, in: Cells, vol. 1, Academic Press, New York, 1970. [9] B.G. Weiss, Perturbations in precursor incorporation into DNA of X-irradiated HeLA S3 cells, Radiat. Res. 48 (1971) 128. [10] J.B. Mitchell, J.S. Bedford, S.M. Bailey, Observations of the first postirradiation division of HeLa cells following continuous or fractionated exposure to c rays, Radiat. Res. 80 (1979) 186. [11] S.L. Tucker, J.M.G. Taylor, Improved models of tumour cure, Int. J. Radiat. Biol. 70 (1996) 539. [12] 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. [13] S.L. Tucker, Modeling 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. [14] H.M. Taylor, S. Karlin, An Introduction to Stochastic Modeling, 3rd Ed., Academic Press, San Diego, CA, 1998. [15] T.E. Harris, The Theory of Branching Processes, Springer, Berlin, 1963. [16] S. Karlin, A First Course In Stochastic Processes, Academic Press, New York, London, 1966. [17] D.G. Kendall, On the generalized ‘‘birth and death’’ process, Ann. Math. Stat. 19 (1948) 1. [18] K.B. Athreya, P.E. Ney, Branching Processes, Springer, New York, 1972. [19] S. Asmussen, H. Hering, Branching Processes, Birkh€ auser, Boston, 1983. [20] E.J. Hall, Radiobiology for the Radiologist, 4th Ed., Lippincott-Raven, 1994.