Generalization of the Critical Volume NTCP Model in the Radiobiology

Alexander Bulinski a , * * T h i s w o r k i s p a r t i a l l y s u p p o r t e d b y t h e R F B R g r a n t 03 − 01 − 00724 , b y t h e g r a n t 1758.2003.1 o f S c i e n t i f i c S c h o o l s a n d b y t h e I N T A S g r a n t 03 − 51 − 5018 . , Andrei Khrennikov b

November 27, 2006

a   Dept. of Mathematics and Mechanics of the Moscow State University, Moscow 119992, Russia E-mail: bulinski@mech.math.msu.su b   Int. Center for Mathematical Modelling in Physics and Cognitive Sciences, MSI, University of Växjö, S-35195, Sweden E-mail: Andrei.Khrennikov@msi.vxu.se
Abstract
A generalization of the well known critical volume NTCP model is proposed to take into account dependence of the functional subunits of irradiated organ (or tissue). A new statistical version of the CLT is established to analyze the corresponding random fields.
AMS classification: 60F05, 62E20, 62G15, 62P10.
Key words: NTCP, dependence conditions, random fields, CLT, statistical version of the CLT. 1. Introduction The problem of finding optimal radiation doses for organs or tissues in therapy of cancer belongs to the principal ones in the modern Radiobiology (see, e.g., [26] , [16] , [18] , [19] , [28] , [30] , [23] and references therein).
The complexity of this problem is related to nondeterministic character of the oncological therapy results.
The aim of this paper is to study stochastic models for collective effects in the behaviour of the irradiated cells. We provide a generalization of the well-known critical volume (CV) normal tissue complication probability (NTCP) model to comprise a concept of functional subunits (FSUs) dependence.
Such a model is beyond the scope of [9] and for its investigation new limit theorems are required. Note that here the links between Probability and Geometry are stipulated by the dependence structure of a random field under consideration which is governed by the configuration of a graph used as a parameter set. Moreover, it seems natural from the biological view-point to assume some dependence in the collective performance of cells (and FSUs). The description of the FSUs behaviour by means of non-binary random variables is considered as well. Other important biological response models and further research directions are tackled in the last Section. 2. Accuracy of the CV NTCP model We recall the basic critical volume model (see [15] , [20] ) and after that consider more carefully its framework. The organ (or tissue) modelled is assumed to be composed of independent FSUs and it is supposed that complications in its functioning arise only if sufficiently many FSUs (”the functional reserve”) are destroyed. More precisely (see, e.g., [27] ), the assessment of the impact of irradiation is divided into two stages. The first one is the reaction of the cells forming an FSU which gives the probability of death or survival for the FSU. The second one is a probability to have sufficient number of FSUs survived the irradiation in order to maintain the organ's functionality (it means the complications absence probability). Note that the serial and tumor control models are special cases when the functional reserve consists of one or all the FSUs respectively. One speaks also of the critical element models. The last special models are natural for organs having the one-dimensional shape (e.g., spinal chord). For other type of organs it is desirable to consider the ”parallel architecture”.
For the sake of simplicity we start with a model of organ consisting of n   similar FSUs. Introduce independent identically distributed (i.i.d.) random variables (r.v's) X 1 , … , X n   defined on some probability space ( Ω , ℱ , P )   such that
P ( X i = 1 ) = p , P ( X i = 0 ) = q w h e r e q = 1 − p , 0 < p < 1 . (1)
More exactly consider an array of r.v's X i = X i ( D )   where D   is a positive parameter (the irradiation dose). Assume that X i ( D )   shows the state of the i   -th FSU after irradiation of dose D   ( i = 1 , … , n )   .
Namely, the random event { X i = 1 }   means that the i   -th FSU is killed and { X i = 0 }   corresponds to the case that the i   -th FSU survives. As usual the argument ω ∈ Ω   is omitted and we write X i ( D )   (or simply X i   ) instead of X i ( ω , D )   .
Set S n ( D ) = ∑ i = 1 n X i ( D )   . In other words consider a random variable equal to the number of FSUs killed due to irradiation dose D   . Thus NTCP is P ( S n ( D ) ≥ L )   with the threshold L   being some positive integer.
Using the convergence rate estimate in the central limit theorem (see, e.g., [10] , p. 323) one has
sup − ∞ < x < ∞ | P ( S n ( D ) ≥ x ) − 1 + Φ ( z ) | ≤ c / n p ( D ) q ( D ) (2)
where z = ( x − n p ( D ) ) / n p ( D ) q ( D )   , p ( D )   and q ( D )   appear in  1 for X i = X i ( D )   , positive constant c ≤ 0.7975   and the c.d.f. of a standard normal law Φ ( z ) = 1 2 π ∫ − ∞ z e − u 2 2 d u , z ∈ R .   Thus if n p ( D ) q ( D )   is large enough, 1   for any threshold x   (possibly depending on n   and p ( D )   ) the following relation holds
P ( S n ( D ) ≥ x ) ≈ 1 − Φ ( z ) (3)
where z   is introduced above and the proximity is evaluated by the right hand side of  2 . It is well known that the power − 1 2   of n   in  2 is the best possible (see, e.g., [22] for that and also about non uniform estimates of the convergence rates in CLT for independent summands).
Consequently, for a given value γ ∈ ( 0 , 1 )   (close enough to 1), using the table of function Φ   one can find z = z γ   as the unique root of the equation Φ ( z ) = γ   . Then the probability of complications is approximately 1 − γ   (with exactly specified boundaries) and the threshold
x γ = n p ( D ) + ( n p ( D ) q ( D ) ) 1 / 2 z γ . (4)
Evidently, x γ = x γ ( n , p ( D ) )   .
The distribution of r.v. S n   has atoms at the points 0 , … , n   . Therefore one could choose an integer threshold L γ = [ x γ ]   where [ ⋅ ]   stands for the integer part of a real number. In this case a trivial estimate for the continuity module of the function Φ   shows that | P ( S n ( D ) ≥ L γ ) − 1 + Φ ( z ) | ≤ 1 n p ( D ) q ( D ) ( c + 1 2 π ) .   We observe that if n p ( D ) q ( D )   is large enough then a search for an integer threshold L γ   is not important. Note in passing that for the number of ”successes” S n   in the Bernoulli scheme (of n   independent trials with probability p   for ”success”) one can apply (see, e.g., [31] ) a little bit different approximation using for 0 ≤ k ≤ m ≤ n   the relation P ( k ≤ S n ≤ m ) = Φ ( t 2 ) − Φ ( t 1 ) + q − p 6 2 π σ { ( 1 − t 2 ) e − t 2 / 2 } | t 1 t 2 + Δ   where σ = n p q   , t 1 = ( k − 1 2 − n p ) / σ   , t 2 = ( m + 1 2 − n p ) / σ   and the error term Δ   satisfies for σ ≥ 5   the inequality | Δ | ≤ ( 0.12 + 0.18 | p − q | ) σ − 2 + e − 3 σ / 2 .   Actually we deal with an equivalent description of the well-known critical volume model. Namely, suppose that the volume 2   of the irradiated organ is V   and let V i   represent the volume of the i   -th FSU ( i = 1 , … , n )   .
Clearly instead of S n ( D )   (the number of killed FSUs) we could consider the random damage volume V ~ = ∑ i : X i = 1 V i   . In the case when V i = V / n   , i = 1 , … , n   , one has
V ~ = V S n ( D ) n . (5)
Thus we come to description of the irradiation result in terms of the damage volume and one can specify the threshold v c   for P ( V ~ ≥ v c )   .
Formula  5 suggests that it is natural to introduce a threshold of the type 3   x = ϰ n   where ϰ ∈ ( 0 , 1 )   is the fraction of killed FSUs. Thus
ϰ = p + c ( p ( 1 − p ) ) 1 / 2 w h e r e p = p ( D ) , c = z γ n − 1 / 2 . (6)
Note that z γ ≥ 0   for γ ≥ 1 / 2   and consequently c ≥ 0   . Evidently for c = 0   (i.e. γ = 1 / 2   ) one has ϰ = p   . For each c > 0   the graph of a function ϰ = ϰ ( p )   has the following features. One can easily verify that ϰ ( p 1 ) = 1   for p 1 = 1 / ( 1 + c 2 )   , ϰ ′ ( 0 + ) = + ∞   , ϰ ′ ( 1 − ) = − ∞   and the concave function ϰ ( p )   attains its maximum ϰ * = 1 2 ( 1 + 1 + c 2 )   at the point p * = 1 2 ( 1 + 1 1 + c 2 )   . Moreover, for γ ≥ 1 / 2   one has 0 ≤ inf 0 ≤ p ≤ 1 ( ϰ ( p ) − p ) ≤ sup 0 ≤ p ≤ 1 ( ϰ ( p ) − p ) ≤ z γ / ( 2 n ) .   Note also that if p ( D ) ≥ p 1   then relation  6 can not be satisfied for any ϰ ∈ ( 0 , 1 )   .
On the other hand, given n ∈ N   , γ ≥ 1 / 2   (i.e. c = z γ / n   ) and ϰ ∈ ( 0 , 1 )   , there is a unique root p = p ¯   of equation  6 
p ¯ = ( ϰ + c 2 2 + c ( ϰ − ϰ 2 + c 2 4 ) 1 / 2 ) / ( 1 + c 2 ) . (7)
We remark also that p ( D )   should be nondecreasing function on ( 0 , ∞ )   and if p ( D )   is continuous then for any ϰ ∈ ( ϰ 1 , ϰ 2 )   where ϰ 1 = ϰ ( inf D > 0 p ( D ) )   and ϰ 2 = min { 1 , ϰ ( sup D > 0 p ( D ) ) }   there exists (unique if p   is strictly increasing) D ¯   such that p ( D ¯ ) = p ¯   .
Now we discuss the models providing p ( D )   . Assume that every FSU consists of n 0   cells. The surviving fraction of these cells after irradiation of dose D   is determined (see, e.g., [30] ) by S F ( D ) = exp { − α D }   where α > 0   is the radiosensivity 4   of the cells. This is a so-called single-hit model. Suppose that each cell of FSU behaves in the same manner as other ones. Usually one admits that an FSU can regenerate from a single surviving cell, which means it is disabled only when no cell survives. Thus the probability of killing an FSU due to irradiation of dose D   is
p ( D ) = ( 1 − e − α D ) n 0 . (8)
To obtain  8 one supposes that all n 0   cells in a FSU evolve independently of each other.
Now assume that every cell contains m   targets, each of them must be hit at least once to inactivate the cell. Then the probability that all targets of a cell will be hit at least once is ( 1 − e − α D ) m   . Thus for this multi-target model S F ( D ) = 1 − ( 1 − e − α D ) m , p ( D ) = ( 1 − e − α D ) m n 0 .   Most experimental survival curves have an initial slope whereas the multi-target/single-hit model predicts no initial slope. To have a more adequate description one uses the family of functions S F ( D ) = e − α D ( 1 − ( 1 − e − β D ) m ) , α > 0 , β > 0 .   Note (see, e.g., [16] ) that for S F ( D )   a linear quadratic (LQ) model is also widely used with S F ( D ) = e − ( α D + β D 2 ) , α > 0 , β > 0 .   Now we concentrate on the generalization of the CV model considered above.
3. Variant of the central limit theorem for dependent random fields Let X ( D ) = { X j ( D ) , j ∈ Z d }   ( d ≥ 1 )   be a family of random fields defined on a probability space ( Ω , ℱ , P )   for D > 0   . Employing instead of an integer lattice Z d   a parameter set T = δ Z d   , with δ > 0   , one can easily reformulate all the results for a family of FSUs assuming, e.g., that the k   -th FSU is a cube with a center at a point t k ∈ T   and with an edge length equal to δ   . In other words one can use a scale appropriate to the problem under consideration. Thus without loss of generality we restrict ourselves to the study of a random field X ( D )   on a lattice Z d   . Moreover, we can assume that the random variable X j ( D )   describes the state of the corresponding FSU after its irradiation of dose D   . This gives us a possibility to consider not only the death and survival of an FSU but also to consider the ”intermediate” states. Then the general (collective) effect of irradiation is represented by the following sum S ( U , D ) = ∑ j ∈ U X j ( D )   where U is a finite subset of Z d   . For a fixed D   we also write simply X j   and S ( U )   . For a finite set I ⊂ Z d   with cardinality | I |   introduce the σ   -algebra A ( I ) = σ { X j , j ∈ I }   , that is consider the σ   -field generated by a field X   over a set I   .
There are different methods (see, e.g., [3] , [13] ) to describe the dependence structure of a field X   . Here we use the maximum correlation coefficient for A ( I )   and A ( J )   over finite disjoint sets I , J ⊂ Z d   which is defined as follows
ρ ( I , J ) = sup { | c o r r ( ξ , η ) | : ξ ∈ L 2 ( Ω , A ( I ) , P ) , η ∈ L 2 ( Ω , A ( J ) , P ) } (9)
where c o r r ( ξ , η )   is the correlation coefficient for (nondegenerate, square integrable) real-valued random variables ξ   and η   measurable with respect to σ   -algebras A ( I )   and A ( J )   .
Assume that for all finite disjoint sets I , J ⊂ Z d   , some positive c 0   and λ   one has
ρ ( I , J ) ≤ c 0 | I | | J | ( d i s t ( I , J ) ) − λ (10)
where d i s t ( I , J ) = min { ∥ q − j ∥ : q ∈ I , j ∈ J } , ∥ z ∥ = max 1 ≤ k ≤ d | z k | , z ∈ Z d .   Remark 1. Employing condition  10 has the following motivation. In many stochastic models it is reasonable to assume that dependence between the random variables { X j , j ∈ I }   and { X j , j ∈ J }   is rather small if the distance between I   and J   is large enough. However, due to the paper [12] it was realized that for random fields (in contrast to stochastic processes corresponding to the case d = 1   ) one cannot, in general, measure the dependence between A ( I )   and A ( J )   only in terms of the distance between I   and J   . Namely in many situations the dependence between A ( I )   and A ( J )   could increase for sets I   and J   growing, e.g., in such a way that the distance between them is fixed. Dependence notions based on correlations are quite familiar in various domains of applied probability. Appearance of the factors | U |   , | V |   and parameters c , λ   in  10 is intended to account, in a qualitative sense, for the effect of possible increase of dependence between σ   -algebras A ( I )   and A ( J )   when I   and J   are growing so that the d i s t ( I , J )   is preserved. Besides,  10 implies the same (i.e. power-type) decrease of correlations when I   and J   are moved apart so that the distance between them tends to infinity. Of course, a simple and natural hypothesis of m   -dependence is a particular case of our condition. Recall that a random field X = { X j , j ∈ Z d }   is m   -dependent (with some m > 0 )   whenever A ( I )   and A ( J )   are independent if d i s t ( I , J ) ≥ m   . Thus, we include, in particular, a useful model of dependent nearest neighbouring FSUs. See also Remark 2 and Section 5.
First of all we establish the central limit theorem (CLT) with convergence rate for partial sums
S ( U n ) = ∑ j ∈ U n X j , n ∈ N , (11)
of multi-indexed dependent r.v's where summation is carried over the integer cubes U n = [ − n , n ] d ∩ Z d   , n ∈ N   .
Theorem 1 Let X ( D ) = { X j ( D ) , j ∈ Z d }   , D > 0   , be a family of strictly stationary random fields such that for some δ ∈ ( 0 , 1 ]   , c 2 + δ ( D ) > 0   and any D > 0  
E | X 0 ( D ) | 2 + δ ≤ c 2 + δ ( D ) . (12)
Assume that condition  10 holds for all fields X ( D )   with the same λ > 4 d ( 1 + δ ) / δ   and c 0   . Then there exists ν = ν ( d , λ , δ ) > 0   such that for each D > 0   and any n ∈ N  
sup x ∈ R | P ( ( S ( U n , D ) − | U n | E X 0 ( D ) ) / ( σ ( D ) | U n | 1 / 2 ) ≤ x ) − Φ ( x ) | ≤ A | U n | − ν (13)
where A = A 0 ( d , λ ) max { 1 , c 0 E X 0 2 ( D ) } max { 1 , c 2 + δ ( D ) / σ 2 + δ ( D ) }   and
σ 2 ( D ) = ∑ j ∈ Z d c o v ( X 0 ( D ) , X j ( D ) ) ≠ 0 . (14)
Proof is based on the classical blocks technique initiated by Bernstein, so we only indicate the main steps and concentrate in the next Section on a statistical version of this result.
For every n ∈ N   introduce p = p ( n ) = [ n α ]   and q = q ( n ) = [ n β ]   where 0 < β < α < 1   and [ ⋅ ]   stands for an integer part of a number. Consider k = k ( n ) = [ ( 2 n + 1 ) / ( 2 p + q ) ]   . Then one can write [ − n , n ] =   I 1 ∪ I 1 ′ ∪ … I k ∪ I k ′ ∪ I k ′ ′   where I m , I m ′ , I k ′ ′   are disjoint intervals of the form I m = [ a m , a m + 2 p ]   , I m ′ = ( a m + 2 p , a m + 2 p + q )   and I k ′ ′ = ( a k + 2 p + q , n ]   ( I k ′ ′   can be empty, a 1 = − n , a m ∈ [ − n , n ]   , m = 1 , … , k )   .
Set B i = I i 1 × … × I i d ∩ Z d   where i = ( i 1 , … , i d ) ∈ M n = { 1 , … , k } d   and let V n = ∪ i ∈ M n B i   .
It is easy to verify that for all n large enough
1 | U n | E ( S ( U n ) − E S ( U n ) − ( S ( V n ) − E S ( V n ) ) ) 2 ≤ ( 1 − | V n | | U n | ) ∑ j ∈ Z d | c o v ( X 0 , X j ) | ≤ 4 d n − γ v ( D ) (15)
where γ = min { 1 − α , α − β }   and the series
∑ j ∈ Z d | c o v ( X 0 ( D ) , X j ( D ) ) | = v ( D ) (16)
converges in view of  10 for λ > d   . Consequently,
E | S ( U n ) − | U n | E X 0 σ ( D ) | U n | 1 / 2 − S ( V n ) − | V n | E X 0 σ ( D ) | U n | 1 / 2 | ≤ 2 ( d v ( D ) ) 1 / 2 n − γ / 2 σ − 1 ( D ) . (17)
Lemma 1 Let X ( D ) = { X j ( D ) , j ∈ Z d }   be a wide-sense stationary random field such that  10 holds with some λ > d   . Then for all n ∈ N  
| v a r S ( U n ) | U n | − σ 2 ( D ) | ≤ a f ( n , d , λ ) (18)
where a = a 0 ( d , λ ) c 0 E X 0 2 ( D )   and f ( n , d , λ ) = { n d − λ , d < λ < d + 1 , ( 1 + ln n ) / n , λ = d + 1 , n − 1 , λ > d + 1 .  
Proof. One has | v a r S ( U n ) − σ 2 ( D ) | U n | | ≤ ∑ i ∈ U n ∑ j / ∈ U n | c o v ( X i , X j ) | = T 1 + T 2   where the sums T 1   and T 2   are taken respectively over the sets { i ∈ U n , j / ∈ U n , ∥ i − j ∥ } = r , r > n }   and { i ∈ U n , j / ∈ U n , ∥ i − j ∥ } = r , r ∈ { 1 , … , n } }   . Evidently, T 1 ≤ 2 d 3 d − 1 | U n | c 0 E X 0 2 ∑ r > n r d − 1 − λ ≤ 2 d 3 d − 1 c 0 E X 0 2 n d − λ | U n | / ( λ − d ) ,   T 2 ≤ ∑ r = 1 n ∑ n − r < ∥ i ∥ ≤ n , ∥ j − i ∥ = r | c o v ( X i , X j ) |   ≤ 2 d c 0 E X 0 2 ∑ r = 1 n r − λ ( ( 2 n + 1 ) d − ( 2 ( n − r ) + 1 ) d ) ( 2 r + 1 ) d − 1 ≤ 4 d 2 3 d − 1 c 0 ( 2 n + 1 ) d − 1 E X 0 2 ∑ r = 1 n r d − λ .   Using a trivial estimate ∑ r = 1 n r − μ ≤ 1 + ∫ 1 n x − μ d x , μ > 0 ,   we come to relation  18 . The Lemma is proved.
Set Y j = S ( B j ) − E S ( B j )   where B j   belongs to a collection of ”large” blocks, j ∈ M n = { 1 , … , k }   , k = k ( n )   . Clearly Y j = Y j ( p n , D )   . Introduce independent copies Z j   , j ∈ M n   , of random variables Y j   , j ∈ M n   .
Then it is easily seen that for any t ∈ R   , i 2 = − 1   and all n   large enough | E exp { i t ∑ j ∈ M n Y j } − E exp { i t ∑ j ∈ M n Z j } | ≤ 4 | M n | c 0 q − λ ( 2 p + 1 ) d | U n | ≤ 4 c 0 3 d n 2 d q − λ ,   E ( ∑ j ∈ M n Y j σ ( D ) | U n | 1 / 2 ) 2 ≤ v ( D ) / σ 2 ( D ) , E ( ∑ j ∈ M n Z j σ ( D ) | U n | 1 / 2 ) 2 ≤ v ( D ) / σ 2 ( D ) .   Thus
| E exp { i t ∑ j ∈ M n Y j σ ( D ) | U n | 1 / 2 } − E exp { ∑ j ∈ M n Z j σ ( D ) | U n | 1 / 2 } | ≤ min { 4 c 0 3 d n 2 d q − λ , 2 | t | v ( D ) / σ ( D ) } . (19)
Using Lemma 1 one can verify that for all n   large enough E ( ∑ j ∈ M n Z j σ ( D ) | U n | 1 / 2 − ∑ j ∈ M n Z j ( v a r ∑ j ∈ M n Z j ) 1 / 2 ) 2 ≤ ( 4 d n − γ + a σ − 2 ( D ) f ( p , d , λ ) ) 2   where γ   , a   and f   are the same as in  15 and  18 . Therefore for λ > d + 1   we can write
E | ∑ j ∈ M n Z j σ ( D ) | U n | 1 / 2 − ∑ j ∈ M n Z j ( v a r ∑ j ∈ M n Z j ) 1 / 2 | ≤ C 1 n − τ (20)
where C 1 = C 1 ( d , a ) max { 1 , σ − 2 ( D ) }   , τ = min { γ , α }   .
Now the Esseen inequality implies that for every T > 0   one has
P ( S ( U n , D ) σ ( D ) | U n | 1 / 2 ≤ x ) − Φ ( x ) | ≤ a 1 ∫ | t | ≤ T | E exp { i t S ( U n ) σ ( D ) | U n | 1 / 2 } − exp { − t 2 2 } t | d t + a 2 T − 1 (21)
where a 1   and a 2   are absolute positive constants. Applying the Berry–Esseen estimate of the convergence rate in the CLT for independent summands Z i , i ∈ M n   , with finite absolute moments of order 2 + δ   (see, e.g., [10] , p. 322), using  17 –  20 and estimating the integral in the right hand side of  21 as a sum of integrals ∫ | t | ≤ 1 / T   and ∫ 1 / T < | t | ≤ T   and finally taking T = b n ζ   with appropriately small ζ   and specified b > 0   we arrive at  13 . This completes the proof of Theorem 1.
Remark 2. There are many versions of the CLT for random fields under various dependence conditions (see, e.g., [2] , [1] , [3] , [11] , [7] ) and references therein). In the same manner we could use instead of the maximal correlation coefficient ρ   , e.g., the Rosenblatt-type mixing coefficient. We proved here the CLT with rate because it permits to establish the law of the iterated logarithm (announced in [8] ) under the dependence conditions of the type  10 . It is worth mentioning that to this end we need only arbitrary slow power-type estimate of the convergence rate in the CLT without specifying an exponent ν   in  13 . We do not provide here an explicit cumbersome expression for ν   . More restrictive mixing conditions than  10 , i.e. ρ ( I , J ) ≤ c 0 | I | | J | exp { − a d i s t ( I , J ) }   where a   and b   are some positive parameters, were recently used in [24] , [25] (see also the references therein) for CLT and LIL. We do not consider here growing subsets U n ⊂ Z d   more general than ”integer” cubes. For generalizations of this kind we refer to [1] , [7] .
4. Statistical version of the CLT There are two ways for applications of Theorem 1. Namely, if we believe in the model describing the stochastic behaviour of each FSU (see Section 1) then we can calculate E X 0   . However, the problem for dependent FSUs is the following one. Now we cannot claim (in general) that the variance of the sum S ( U n , D )   is equal to the sum of variances of summands. Thus for every D > 0   , in contrast to the CLT for the Bernoulli scheme, i.e.
S n ( D ) − n p ( D ) n p ( D ) ( 1 − p ( D ) ) → D Z ∼ N ( 0 , 1 ) a s n → ∞   discussed in Section 1 (here p ( D ) = E X 0   , Z   is a standard normal r.v.), the relation
S ( U n , D ) − | U n | E X 0 ( D ) σ ( D ) | U n | 1 / 2 → D Z ∼ N ( 0 , 1 ) a s n → ∞ (22)
contains an unknown function σ ( D )   . As usual → D   stands for weak convergence of random variables distributions.
Fortunately it is possible to construct a sequence of nonnegative statistical estimates C ^ ( U n , D )   for σ 2 ( D )   such that for any D > 0  
C ^ ( U n , D ) → P σ 2 ( D ) a s n → ∞ (23)
where → P   means the convergence in probability as usual. We employ here a family of consistent statistical estimates introduced in [5] for random fields 5   , for stochastic processes we refer to the paper [21] .
Then by virtue of  22 and  23 we come, for every D > 0   (if σ ( D ) 2 ≠ 0   ), to the formula
( C ^ ( U n , D ) | U n | ) − 1 / 2 ( S ( U n , D ) − | U n | E X 0 ( D ) ) → D Z ∼ N ( 0 , 1 ) a s n → ∞ . (24)
In other words a random normalization is used in the CLT. Consequently to determine (approximately) for a given value γ ∈ ( 0 , 1 )   the threshold x γ   we can apply the following analogue of formula  3 
P ( S n ≥ x ) ≈ 1 − Φ ( ξ ) (25)
where ξ = ( x − n E X 0 ( D ) ) / C ^ ( U n , D ) | U n |   , x ∈ R   , n ∈ N   . However, now in the right hand side of  25 there is a r.v. Φ ( ξ )   , i.e. we use Φ ( ξ )   as statisitical estimate for P ( S n < x )   . Note that we have used only the value E X 0 ( D )   provided by the model of stochastic behaviour for FSUs and we did not suppose here that the collective effect of the evolution of cells under irradiation is described by independent binary random variables.
Another way of using Theorem 1 is to construct approximate confidence intervals for the unknown mean value E X 0 ( D )   without hypotheses concerning the explicit formulas (discussed in Section 1) for distribution of random variables X j   , j ∈ U n   .
Thus in both cases it is desirable to establish the CLT for dependent random fields using random normalization.
Remark 3. As far as we know, in previous applications of the CLT to NTCP models for independent FSUs the question of convergence rate was not raised, so Section 1 covers this gap. However, the same question in case of dependent FSUs is more involved. We intend to investigate the accuracy of the proposed model in a special publication. One can consider Theorem 1 as the first step in this direction.
Moreover, we can obtain the power-type estimate in the CLT with random normalization. However, the rate of convergence will be slower than that for independent random summands. The effect of convergence rate in the CLT sensitivity to the dependence conditions was demonstrated for positively or negatively associated random fields in [4] .
For j ∈ U ⊂ Z d   ( 1 ≤ | U | < ∞ )   and b = b ( U ) > 0   set
K j ( b ) = { t ∈ Z d : ∥ j − t ∥ ≤ b } , Q j = Q j ( U , b ) = U ∩ K j ( b ) , (26)
C ^ ( U , D ) = 1 | U | ∑ j ∈ U | Q j | ( S ( Q j , D ) | Q j | − S ( U , D ) | U | ) 2 . (27)
Note that the averaged variables S ( Q j , D ) / | Q j |   arise for dependent summands (in contrast to the traditional estimates of variance used for independent observations).
Theorem 2 Let the conditions of Theorem  1 be satisfied. Let { U n } n ∈ N   be a sequence of ”integer” cubes, i.e. U n = [ − n , n ] d ∩ Z d   , n ∈ N   ( d ≥ 1 )   . Assume b ( U n ) = b n   where { b n } n ∈ N   is a sequence of positive integers such that
b n → ∞ , b n = o ( n ) a s n → ∞ . (28)
Then for every D > 0   relation  24 holds if σ ( D ) ≠ 0   .
Proof. The estimate C ^ ( U , D )   introduced by means of  27 and σ 2 ( D )   are invariant under the transformation X j ↦ X j − E X 0   , j ∈ Z d   . So, without loss of generality we can further on assume that E X 0 = 0   . Let ∥ ξ ∥ L   stand for the norm of a real-valued random variable ξ   in a space L = L 1 ( Ω , ℱ , P )   . For any U n ⊂ Z d   and n ∈ N   one has ∥ C ^ ( U n , D ) − σ 2 ( D ) ∥ L ≤ I 1 ( U n , D ) + I 2 ( U n , D ) + I 3 ( U n , D )   where I 1 ( U n , D ) = 1 | U n | ‖ ∑ j ∈ U n | Q j | { ( S ( Q j ) | Q j | − S ( U n ) | U n | ) 2 − ( S ( Q j ) | Q j | ) 2 } ‖ L ,   I 2 ( U n , D ) = 1 | U n | ‖ ∑ j ∈ U n 1 | Q j | ( S 2 ( Q j ) − E S 2 ( Q j ) ) ‖ L ,   I 3 ( U n , D ) = | 1 | U n | ∑ j ∈ U n 1 | Q j | E S 2 ( Q j ) − σ 2 ( D ) | .   Here and below S ( U ) = S ( U , D )   for U ⊂ Z d   and D > 0   . We have
| Q j | − 1 E S 2 ( Q j ) ≤ v ( D ) , j ∈ Z d , (29)
where v ( D )   is the same as in  16 . By virtue of condition  28 it is clear that I 1 ( U n , D ) ≤ | U n | − 3 E S 2 ( U n ) ∑ j ∈ U n | Q j | + 2 | U n | − 2 ∑ j ∈ U n E | S ( Q j ) S ( U n ) |  
≤ v ( D ) { | K 0 ( b n ) | | U n | − 1 + 2 | K 0 ( b n ) | 1 / 2 | U n | − 1 / 2 } → 0 as n → ∞ . (30)
For a fixed c > 0   introduce the functions
h 1 ( x ) = sign ( x ) min { | x | , c } , h 2 ( x ) = x − h 1 ( x ) , x ∈ R . (31)
Given a nonempty finite set Q ⊂ Z d   let S ¯ ( Q ) = S ( Q ) / | Q | .   Note that
I 2 ( U n , D ) ≤ ∑ p , m = 1 2 I 2 ( p , m ) ( U n , D ) (32)
where I 2 ( p , m ) ( U n , D ) = 1 | U n | ∥ ∑ j ∈ U n h p ( S ¯ ( Q j ) ) h m ( S ¯ ( Q j ) ) − E h p ( S ¯ ( Q j ) ) h m ( S ¯ ( Q j ) ) ∥ L .   For b , n ∈ N   introduce the sets T n ( b ) = { s ∈ U n : inf t ∈ ∂ U n ∥ s − t ∥ ≤ b }   where ∂ U n = { j ∈ U n : ∃ q / ∈ U n such that ∥ j − q ∥ = 1 }   . Put T n = T n ( b n )   , n ∈ N   , where b n   meet condition  28 . Due to  29 one has I 2 ( 1 , 2 ) ( U n , D ) + I 2 ( 2 , 1 ) ( U n , D ) + I 2 ( 2 , 2 ) ( U n , D )   ≤ 2 | U n | − 1 ∑ j ∈ U n ( 2 E | h 1 ( S ¯ ( Q j ) ) h 2 ( S ¯ ( Q j ) ) | + E h 2 2 ( S ¯ ( Q j ) ) )   ≤ 2 ( 2 E | h 1 ( S ¯ ( K 0 ( b n ) ) ) h 2 ( S ¯ ( K 0 ( b n ) ) ) | + E h 2 2 ( S ¯ ( K 0 ( b n ) ) ) + 3 | T n | | U n | − 1 v ( D ) )   ≤ 4 ( v ( D ) E ( S ¯ 2 ( K 0 ( b n ) ) 1 { | S ¯ ( K 0 ( b n ) ) | ≥ c } ) ) 1 / 2   + 2 E ( S ¯ 2 ( K 0 ( b n ) 1 { | S ¯ ( K 0 ( b n ) ) | ≥ c } )  
+ 6 | T n | | U n | − 1 v ( D ) (33)
where 1   is an indicator function.
Condition  29 implies that
| T n | | U n | − 1 → 0 as n → ∞ . (34)
It is easy to see that a family { S ¯ 2 ( K 0 ( b n ) ) } n = 1 ∞   is uniformly integrable. Consequently, taking into account  33 , for any ɛ > 0   we can find c = c ( ɛ )   such that for all n   large enough
I 2 ( 1 , 2 ) ( U n , D ) + I 2 ( 2 , 1 ) ( U n , D ) + I 2 ( 2 , 2 ) ( U n , D ) < ɛ , (35)
furthermore,
( I 2 ( 1 , 1 ) ( U n , D ) ) 2 ≤ | U n | − 2 ∑ j , t ∈ U n | cov ( h 1 2 ( S ¯ ( Q j ) ) , h 1 2 ( S ¯ ( Q t ) ) ) | . (36)
In view of  31 we obtain the inequalities | U n | − 2 ∑ j , t ∈ U n , ∥ j − t ∥ ≤ 4 b n | cov ( h 1 2 ( S ¯ ( Q j ) ) , h 1 2 ( S ¯ ( Q t ) ) ) |   ≤ 2 c 2 | U n | − 2 ∑ j , t ∈ U n , ∥ j − t ∥ ≤ 4 b n E h 1 2 ( S ¯ ( Q j ) )  
≤ 2 2 d + 1 c 2 | U n | − 1 | K 0 ( b n ) | v ( D ) . (37)
Now condition  10 with λ > d   entails the estimate | U n | − 2 ∑ j , t ∈ U n , ∥ j − t ∥ > 4 b n | cov ( h 1 2 ( S ¯ ( Q j ) ) , h 1 2 ( S ¯ ( Q t ) ) ) | ≤ c 0 c 4 | U n | − 2 ∑ j , t ∈ U n , ∥ j − t ∥ > 4 b n | Q j | | Q t | ∥ j − t − 2 b n ∥ − λ  
≤ c 0 c 4 d 4 d | U n | − 1 | K 0 ( b n ) | 2 ∑ r > 2 b n r d − 1 − λ . (38)
Taking into account  35 –  38 ,  28 and  10 with λ ≥ 2 d   , we verify that
I 2 ( U n , D ) → 0 as n → ∞ . (39)
Now observe that | U n | − 1 ∑ j ∈ U n | Q j | − 1 E S 2 ( Q j ) = | U n | − 1 | U n \ T n | | K 0 ( b n ) | − 1 E S 2 ( K 0 ( b n ) )   + | U n | − 1 ∑ j ∈ T n | Q j | − 1 E S 2 ( Q j ) .   According to  18 the following relation is valid | K 0 ( b n ) | − 1 E S 2 ( K 0 ( b n ) ) → σ 2 ( D ) as n → ∞ .   Due to  28 and  29 we conclude that
I 3 ( U n , D ) → 0 as n → ∞ . (40)
Relations  30 ,  39 and  40 yield  21 . The proof of Theorem 2 is complete.
5. Concluding Remarks Besides the concept of a functional reserve of an organ undergoing irradiation it seems desirable to use the models taking into account the geometrical configuration of survived FSUs (or/and cells). Here the concepts of random clusters appear naturally. In this regard we refer to a quite recent paper [29] (cf. [17] ).
The stochastic models of disordered structures (involving point random fields) could be applied also to describe the damage volumes of irradiated organ. Note that it is possible to describe the relations between various FSUs (cells) in terms of random graphs and study the models where some vertices (or edges) are destroyed at random. An interesting problem is to find the optimal dose of irradiation taking into account not only the complication probabilities but the balance of conditions for irradiated organ (tissue) and its normal environment. Moreover, it is important to consider non uniform irradiation, another problem is to study a population of non-identical patients (see, e.g., [32] ,[14] ). To conclude we mention a deep problem of constructing dynamical models describing the evolution of an irradiated organ in space and time.
Acknowledgements The authors thank Professor L.Klepper for drawing their attention to stochastic models in Radiobiology. A.Bulinski is grateful for invitation to the Växjö University where in August 2004 the results of this joint paper were obtained and he is grateful also to the Dept. of Probability Theory and Stochastic Models of the University Paris-VI where the final version of the paper was written.
References

  1. Bolthausen E., On the central limit theorem for stationary mixing random fields, Ann. Probab., 10(1982), 1047–1050.
  2. Bulinski A.V., Zhurbenko I.G., A central limit theorem for additive random functions, Theory Probab. Appl., 21(1976), 687–697.
  3. Bulinski A.V., Limit Theorems under Weak Dependence Conditiions, MSU, 1989 (in Russian).
  4. Bulinski A.V., On the convergence rates in the CLT for positevely or negatively dependent random fileds, In.: I.A.Ibragimov and A.Yu.Zaitsev (Eds.) Proc. of the Kolmogorov Semester Int. Euler Math. Inst., St.-Petersburg (02.03.1993 – 23.04.1993). Gordon and Breach, 1996, 3 – 14.
  5. Bulinski A.V., Vronski M.A., Statistical variant of the central limit theorem for associated random fields, Fund. Appl. Math., 2(1996), 891–908 (in Russian).
  6. Bulinski A.V., Statistical version of the central limit theorem for vector-valued random fields, Math. Notes, 76(2004), 455–464.
  7. Bulinski A.V., Shashkin A.P., Rates in the central limit hteorem for weakly dependent random variables, J. Math. Sci., 122(2004), 3343–3358.
  8. Bulinski A.V., A law of the iterated logarithm for random fields with power decay of correlations, Rev. Appl. Indust. Math., 11(2004), 503–504 (in Russian).
  9. Bulinski A.V., Stochastic models in microand macroworld, Proc. Semin. ”Time, Chaos and Mathematical Problems”, MSU, 3(2004), 231–254 (in Russian).
  10. Chow S.C., Teicher H., Probability Theory. Independence, Interchangebility, Martingales, Springer, New York etc., 1997 (3-d ed).
  11. Dedecker J., A central limit theorem for stationary random fields, Probab. Th. Rel. Fields, 110(1998), 397–426.
  12. Dobrushin R.L., Description of a random field by means of conditional probabilities and conditions of its regularity, Theory Probab. Appl., 13(1968), 197–225.
  13. Doukhan P., Mixing: Properties and Examples, Lecture Notes in Statistics, 85(1994), Springer–Verlag.
  14. Fenwick J.D., Nahum A.E., Series model volume effects in a population of non-identical patients: how low is low?, Phys. Med. Biol., 46(2001), 1815–1834.
  15. Jackson A., Kutcher G.J., Yorke E.D., Probability of radiation-induced complications for normal tissues with parallel architecture subject to non-uniform irradiation, Med. Phys., 20(1993), 613–625.
  16. Jones B., Dale R.G., Mathematical models of tumours and normal tissue response, Acta Oncol., 38(1999), 883–893.
  17. Klepper L.I., The occurrence of radiation complications in body organs and tissues, Med Tekh., 5(2000), 36–40 (in Russian).
  18. Lind B.K., Mavroidis P., Hyodynmaa S., Kappas C., Optimization of the dose level for a given treatment plan to maximize the complication-free tumor cure, Acta Oncol., 38(1999), 787–798.
  19. Moiseenko V., Battista J., Van Dyk J., Normal tissue complication probabilities: dependence on choice of biological model and dose-volume histogram reduction scheme, Int. J. Radiat. Oncol. Biol. Phys., 46(2000), 983–993.
  20. Niemierko A., Goitein M., Modeling of normal tissue response to radiation: the critical volume model, Int. J. Radiat. Oncol. Biol. Phys., 25(1993), 135–145.
  21. Peligrad M., Shao Q-M., Estimation of variance of partial sums of an associated sequence of random variables, Stoch. Proc. Appl., 56(1995), 307–319.
  22. Petrov V.V., Limit Theorems of Probability Theory: Sequences of Independent Random Variables, Oxford Univ. Press, 1995.
  23. Philippsens M.E., Pop L.A., Visser A.G., Schelleksns S.A., Van Der Kogel A.J., Dose-volume effects in rat thoracolumbar spinal chord: an evaluation of NTCP models, Int. J. Radiat. Oncol. Biol. Phys., 60(2004), 578–590.
  24. Schmuland B., Sun W., A central limit theorem and law of the iterated logarithm for random fields with exponential decay of correlations, Canad. J. Math. 56(2004), 209–222.
  25. Schmuland B., Sun W.A., The law of large numbers and law or the iterated logarithm for infinite dimensional intaracting diffusion processes, Inf. Dim. Anal., Quantum Probab. Rel. Topics (IDAQP), 6, 489–503.
  26. Smith A.R. (Ed.), Medical Radiology: Radiation Therapy Physics, Springer – Verlag, Berlin, 1995.
  27. Stavrev P., Stavreva N., Niemierko A., Goitein M., Generalization of a model of tissue response to radiation based on the idea of functional subunits and binomial statistics, Phys. Med. Biol., 46(2001), 1501–1518.
  28. Stavreva N., Stavrev P., Warkentin B., Fallone B.G., Derivation of the expressions for gamma 50 and D50 for different individual TCP and NTCP models, Phys. Med. Biol., 47(2002), 3591–3604.
  29. Thames H.D., Zhang M., Tucker S.L., Liu H.H., Dong L., Mohan R., Cluster models of dose-volume effects, Int. J. Radiat. Oncol. Biol. Phys., 59(2004), 1491–1504.
  30. van Luijk P., Delvigne T.C., Schilstra C., Schippers J.M., Estimation of parameters of dose-volume models and their confidence limits, Phys. Med. Biol., 48(2003), 1863–1884.
  31. Weisstein E.W., ”de Moivre–Laplace Theorem”. From MathWorld – A Wolfram Web Resource. http://mathworld.wolfram.com/deMoivre-LaplaceTheorem.html
  32. Yorke E.D., Modeling the effect of inhomogeneous dose distributions in normal tissues, Semin. Radiat. Oncol., 11(2001), 197–209.

1   Otherwise one has to use different approximations.

2   Possibly a length or an area, the interpretation depends on the model of an organ

3   We do not use in this paper the theory of large deviations for sums of r.v's.

4   One writes also S F ( D ) = exp { − D / D 0 }   where D 0   is called the mean lethal dose.

5   Vector-valued random fields satisfying other dependence conditions are studied in [6] .