- Email: [email protected]
Control of a general lethal growth process
MATHEMATICAL
Control
BIOSCIENCES
13,235-252
of a General Lethal Growth
ELLEN A. CHERNIAVSKY
AND
HOWARD
235
(1972)
Process* M. TAYLOR
Center for Environmental Quality Management and Department of Operations Research, Cornell University, and
Ithaca, New
York
Department
Ithaca, New
York
of Operations
Communicated
Research,
Cornell
University,
by Richard Bellman
ABSTRACT It is assumed that a host carries a particle population, initially of size N,, that increases according to a general birth process. Treatments that remove a random number of particles, depending on the dosage level, may be administered. The treatment history and size of the particle population contribute to the death rate of the host. Both when the size of the population is observable and when it is not, the expected lifetime of the host is derived as a function of the frequency and dosage level of treatments. In some special cases the frequency and dosage level that maximize the host’s expected lifetime are given.
1. INTRODUCTION
AND SUMMARY
We consider a particle population whose size, in the absence of intervention, evolves stochastically as a terminating pure birth process. These particles are carried by, and in turn influence, a host, and termination of the birth process corresponds to death of the host. But we intervene in the process, modifying its evolution at intervals, by administering treatments that remove or destroy a random number of the particles present. Unfortunately, these treatments entail some danger to the host. The problem is to balance the benign and malignant effects so as to arrive at a frequency and intensity schedule for treatments that will maximize the expected lifetime of the host. Neuts [2] considered such a model under the assumptions that (i) the rate of growth of the particle population is proportional to its current size (a Yule process); (ii) the particles contribute to the death rate of the host in direct proportion to the particle population size; (iii) the random * This research was conducted under the auspices of the Center for Environmental Quality Management, Cornell University, 1969-1970, and was supported in part by PHS grant 1 TO1 ES 00130-03 and in part by NSF grant GK-21460. Copyright 0 1972 by American Elsevier Publishing Company, Inc.
236
ELLEN A. CHERNIAVSKY AND HOWARD M. TAYLOR
number of particles destroyed by treatment has a binomial distribution; and (iv) the host is given treatments of constant dosage at fixed and equal intervals of time. We extend Neuts’s work by allowing the evolution of the particle population to be a general pure birth process under general termination parameters. In addition, we allow a general distribution for the number of particles destroyed through treatment. More important, we consider not only periodic treatments, but introduce the concept of treating whenever the particle population reaches a pre-specified size. In a later section we consider a third type of policy in which treatments are given at equal intervals until the host either dies or is cured. In all these cases we derive general expressions for the expected host lifetime. Finally, in a variety of special cases, we further develop, analyze, and simplify these general formulas. 2.
THE MATHEMATICAL
MODEL
We begin by assuming that at time t = 0 there are N,, particles present. If at any instant t there are i particles present, then we suppose that the probability of one new particle being born in (t, t + dt) is Izi dt where {,$, i 2 O> are arbitrary nonnegative parameters, called the birth parameters of the process. We postulate that if the host is alive at time t, the probability that he dies during (t, t + dt) depends only on the number of particles present at time t and on the number of particles present at time t and on the number, dosages, and times of previous treatments, which we summarize as a treatment history h. Specifically, we assume that a host carrying i particles and having a previous treatment history h has a probability e,(h) dt of dying in (t, t + dt). For example, Neuts considers 0,(/z) = c + ip + vQ(a, p) where (T is the death rate of an uninfected individual, ip is the death rate due to the particle population, assumed proportional to the number of particles present, and the last term corresponds to the long-range aftereffects of v earlier treatments, each of dosage p and administered at intervals of length a. Intervention involves treating the host, and we allow several levels of treatment, indexed by p. A host with history h survives a treatment of level p with probability Il/(p, h), assumed known. If he survives, a random number of particles is destroyed. Let B(k, p, i) be the probability that out of a population of k particles, i particles survive a treatment at level p. The global sequence of events begins with N,, particles present in a host having a treatment history h,,, where ho signifies the absence of previous treatments. The population grows stochastically as a pure birth process with birth parameters {&} and termination parameters (ei(h,)>, until either the host dies or a treatment is initiated. In the latter case, suppose
231
LETHAL GROWTH PROCESS CONTROL
that the particle population has grown to size N(, at the time treatment is begun, and that the treatment is at level p. With probability 1 - $(p, ho) the host dies and the process terminates. Otherwise, the host survives and begins a new era, this time with N1 particles present and a history of hl, where N1 is drawn from the distribution B(N& p, .) and h, is the host’s updated history. The particle population now increases stochastically from N1 and the process repeats. The host’s life consists of a random number of such periods punctuated by treatments, during each of which the particle process grows undisturbed as a pure birth process. Throughout, we will refer to such periods as growth eras, or more simply, as eras. Each such era is a period of evolution in a pure birth process, whose starting state and parameters may differ from those of other eras. Each era other than the first begins just after a nonfatal treatment. Each era other than the last ends just before a treatment. Either the host dies in the last era, or survives the last era but dies in the ensuing treatment. It follows from this simple description that if we can compute two things about each era, we can compute the expected lifetime of the host. We need to know (i) P(i,j; h), the conditional probability that the host survives an era and ends it carryingj particles, given that he began the era with i particles and a treatment history h; and (ii) L(i, h), the expected length of an era, given that the host began it with i particles and a treatment history h. Letf(i, h) be the expected future lifetime, from the beginning of some era, of a host that carries i particles and has a treatment history of h. We want to computef(N,, h,). To simplify the final statement, let R(i, k; 4,) =
HP,
h,)
C j3i
W, j;
kJNj,
P,
4.
(1)
Then R(i, k; h,) is the conditional probability that the host survives an era and the ensuing treatment and that after the treatment he carries k particles, given that he began the era with i particles and a treatment history h,. Set R(‘)(i, k; h,) = R(i, k; h,)
and then for I > 1, define R”+“(i,
k; h,) = c R(‘)(i,j;
h,)R(j,
k; h,,,).
(2)
jB0
Then R”)(i, k; h,) is the conditional probability that the host survives 1 eras and the treatment following each era and that after the lth treatment he carries k particles, given that he began the first era with i particles and a history h,. Therefore ‘&oR(l)(No,j, h,)L(j, h,) is the expected length of the Zth era, given that the host began the first era with N,, particles and a
238
ELLEN
A. CHERNIAVSKY
AND HOWARD
M. TAYLOR
history ho. (The length of the Zth era is taken as 0 if the host fails to reach it.) The expected total life is the sum, over 1, of the expected life in the Ith era. Thus
SW0,h0)= JWO, h0) + I=1 f, f W~o,j; hoMj, h). j=O Our next task is to compute ted lifetime from (3).
P(i,j;
(3)
h) and L(i, h), and in turn, the expec-
3. THE GENERAL FORMULA Let X(t) be the size of the lethal particle population at time t. During an era this population evolves as a terminating pure birth process (X(t); 0 < t -c i) where [ is the lifetime of the process and corresponds to the lifetime of the host. Under this assumption, the stochastic process is governed by two sets of parameters: the birth parameters {Ai} and the termination parameters {e,}. (Since we are concerned here only with the process during a fixed era, we have suppressed, for the moment, the dependence on the history of previous treatments. Later, of course, when several eras forming a lifetime are simultaneously under discussion, we will again have to make explicit this dependence on the history.) Let P,(i, j) be the probability that such a process is in statej at time t, given that it began in state i at time t = 0. We mean, of course, P,(i,j)
= Pr{[ > t
and
X(t) =,ilX(O)
= i>,
(4)
so that Pr([ We postulate
> tlX(0)
= i} = x P,(i,j). jai
(5)
for i = 0, 1, . . . that (i) P,(i, i + 1) = hit + O(t), (ii) Pr([ E (0, t)lX(O) = i} = Oit + o(t),
(6)
(iii) P,(i, i) = 1 - (pi + Oi)t + O(t), (iv) P,(i, j) = a,,, where o(t) represents remainder terms of order less than t as t --f 0 and aij = 1 if i = j and 0 otherwise. Under these postulates, the birth parameters {AJ govern the growth of the particle population and the termination or mortality parameters {Oi} govern the probability of the host dying. In most cases ok will be nondecreasing in k; the more particles present, the more likely the death of the host. Note that the birth and mortality parameters do not depend on the age of the host, so that the process is stationary in this sense.
239
LETHAL GROWTH PROCESS CONTROL
Alternatively, a terminating pure birth process with parameters {jei} and {ei} may be characterized as a stochastic process on the states 0, I, . . . where the sojourn time in state i is an exponentially distributed random variable with parameter li + (Ii, all sojourn times being statistically independent, and in which, at the end of a visit to state i the process terminates with probability 0,/(2, + 0,) or else moves to state i + 1, which occurs with the remaining probability %J(ii + ei). Again these terminal decisions are statistically independent and independent of the sojourn times. From this description there follows a simple but useful construction of such a process. Let {Si} be independent exponentially distributed random variables with parameters E.i + ei. Let {Ii} be independent random variables, independent of the {Si> and having the distribution 1
with probability
(0
with probability
Ii =
Ai
;li +
ei’
ei ~~+ ei
integer i. Set
Fix a nonnegative
Vi(i) = 0, while forj
> i set vi(j)
Finally,
= Si + ZiSi+l + . . . + ZiZj+r x
,.
X
Zj_*Sj-1.
(7)
set i
=
v,(co) = 1
(8)
jai
and for t E [0, 0 define x(Z) =
.i
if and only if
V,(j) < t < Vi(.j + 1).
(9)
(We take n::! Z, = 1 in Eq. (8).) Then X(t) is a terminating pure birth process starting from i and governed by the parameters {ii} and {Oi}. We are still considering a fixed era, and the treatment history remains suppressed in our notation. The V,(j) may be interpreted as the random length of time until either the host dies or the particle population reaches size j, given that the era began with i particles present. A terminating pure birth process is a Markov process. In our case the transition function is P,(i,J)
= Pr{t < 5
and
X(t) =,ilX(O)
= i},
as given in Eq. (4). (When we again exhibit that the termination parameters {ei} depend on the previous treatment history 11, we will write
240
ELLEN A. CHERNIAVSKY
AND HOWARD
M. TAYLOR
Note that for t > 0
P,(i, j; h) for this quantity.)
C P,(i,j) < 1. j ai In fact, Q,(t), defined to be the cumulative distribution process lifetime [, is given by
Q,(t) = 1 - Jgif’t(I,_d
(10) function
of the
(11)
= Pr([ < tlX(0) = i}. So far we have described how the process evolves in the absence of intervention, or during the eras between treatments. An era ends with death of the host or when a treatment is begun, say at time T (which may be fixed or random, as yet unspecified). What we want to compute for each era, assuming the era begins with a particle population level i and a treatment history h, is (extending the notation in (4) to possibly random times) P,(i,j; h) = Pr{T < i and X(T) =jlX(O) = i} (12) and L(i, h) = E[min{ [, T} IX(O) = i] (13) where the history h appears on the right through the mortality parameters (0J = {Gi(h)} governing c. To continue the computations we need more information on T, the treatment time, and in particular. we need to know the rule determining T. We consider two types of treatment rules, corresponding to two alternative assumptions about the particle process. First we may assume that the number of the particles carried by the host is known to the observer at all times. In the case, called the observable case, it seems natural to initiate treatment whenever the number of particles reaches or exceeds a prespecified critical level M. That is, when the process is observable, we consider rules of the form T = min{t:
0 < t < 5
and
A’(t) 2 M}
(14)
where we agree to set T = co whenever the set in braces is empty. Alternatively, we may assume that nothing is known about the current level of the particle process, called the unobservable case. The only information available then is the passage of time, which forces the use of Jixed-time rules in which T is set at a fixed prespecified number. Of course, reality probably lies somewhere between these two extremes. The particle process may give some evidence of its presence but not reveal its exact size. The Observable
Case under Reach-State
Rules
Let A4 be the critical state determining the reach-state rule. We may suppose, without a significant loss in generality, that the process starts
241
LETHAL GROWTH PROCESS CONTROL
below M. Then, by the definition of T,whenever we treat before the host dies, that is, whenever T < [, we have X(T) = M. It follows from (9) that P,(i,
M; h) = Pr{T < ilX(O) = i}
(15)
and that P,(i,j; h) = 0 forj # M. Continuing, we reach state M, that is, T < 5, if and only if Zi = Zi+l = ...= ZM_l= 1 in the process representation described by Eq. (7)-(9). Since the Zi’S are independent, the probability is easily evaluated to be
Similarly,
from Eq. (7)-(9) min{T,
[> = V,(M),
(17)
so that Z(i, A) = which evaluates
E[J’i(Wl,
W3)
to M-l
L(i, h) =
c
(1, +
em)-’
(19)
m=i
using the independence of the Zi’s from the Sj’s and that E[S,] = (Aj + (3,)-l. (Note: Throughout, a product over an empty set is interpreted as unity.) Set A(M, P, h) =
f
B(M, P, j)P,(.L M; h)
j=O
=
P?~)~~l kcj (A, +fik(hJ
j;. wf,
(20)
Then A(M,p, 12) is the probability that a host with treatment survives any era other than the first. From (1) and (2) R”‘(i, k; ho) = $(P, h)P,(i,
history
h
M; ho)
l-l x
1
“gl
ACM>
P, k,)t4~,
h)
1
Wf,
P, k).
(21)
From Eq. (3) we may now write an explicit expression for the expected lifetime of the host. The subscript M indicates the observable case under a reach-state rule in which the host is treated when the particle population first reaches M. The expected lifetime is
242
ELLEN A. CHERNIAVSKY
j=O
AND HOWARD M. TAYLOR
m=j
(22) The Unobservable
Case under Fixed Time Rules
Suppose we treat at time T, a prespecified constant. We need to compute P(i, j; h) and L(i, h) for such rules, but since T is fixed, we have L(i,
h) =E[min(i,
T}IX(O)
= i]
T
Pr{i > t]X(O) = i} dt
=
s0
s1 =:
j& P,(4j; h) dt
=
(23)
(1 - Qi(t)} dt.
s
Thus, we need only compute P(i, j; h) = P,(i, j; h) for arbitrary fixed T = t. This is done in the following lemma, which generalizes the formula of Feller [I] to allow termination. LEMMA 1.
with parameters
Let (X(t); 0 < t < i) be a terminating pure birth process {Ai} and {oil subject to %i + Oi # lj + Oj for i # j. Set
Aij(rn) = 1
Then P,(i, j) = Pr{i PAi, j) = i]
> t and X(t) = jlX(0) Aij(m)
expC-(2,
+
for
i=j
for
i < j.
(24)
= i} is given by for
~&I
i < j.
(25)
m=i
X(t) = ilX(0) Proof. The case i = j is easy, since Pr{[ > t and i} = Pr{S, > t} = exp[-(I., + Oi)t]. Thus, suppose i < j. Let G,(i, j) = ,f Aij(m) exp[-(1, m=i
+ e,)t],
i
=
243
LETHAL GROWTH PROCESS CONTROL
The Laplace
transform
is m
Gr(i,j)
=
e-“G,(i,
j) dt
s0
where the second right-hand side is the partial fraction expansion of the third. If G;*(i, j) is the Laplace transform of the derivative; dG,(i,J)/dt, then G;*( i, j) = tG*( i, j), since for i < j. G,(i,j) = 0 Because of this correspondence
it is sufficient
to show
P;*(i, j) = Gi*(i, j)
But, using (9), we have = Pr{ vi(j) < t < Vi(j + I)}
P,(i,.j)
= Pr{ Vi(j) < t} - Pr(V,(j
+ I) < t)
and P;*(i,j)
=
J e-c’ dP,(i,j)
=
Jem@
J e-c’ dPr{I/,(j + 1) < t}
< t} -
dPr{K(j)
@evHV(j)lI - E{ev[-Wj = V*(i,j, 5) - V*(i, j + 1, 5)
+ I)]}
= where
I/*(4 j, 0 = E{exp[- OXj>~>
= =
+
E{eXp[-_(Si E(exp[ - @,I)
8.
=ii+Oi+~x
[
lii+
+
+
Ii
X ’ . .
X
Zj-_ZSj-_l)]}
+ * * * + Ii+l x ... x zj-2sj-l)l)
f&,“*(i +L.L 5)1 1 I
;li
4
ei + 5
’ . .
X E(exp[-~li(Si+1
ii + 8i
= Ai +
ZiSi+l +
Ai + Oi + 5
V*(i + 1, j, 5).
244
ELLEN A. CHERNIAVSKY
AND HOWARD
M. TAYLOR
Then
= I. + :, + 1
$V*(i + Lj, 5) - V*(i + Lj + 1, 01,
1
which iterates to P;*(i, j) = V*(i, j, 5) -
V*(i, j + 1, 5)
5 KZ;=f I, = nj,=& as was to be shown.
+ 8, + 5)’
From (23) we may compute
Again Eq. (25) and (26) could be substituted into Eq. (l)-(3) to give an explicit general expression for the expected lifetime of the host. In theory, we could then maximize with respect to the treatment level p and fixed treatment interval T. In practice, the resulting expression, as in the observable case under reach-state rules, is too complex in general for any but the most rudimentary numerical exploration. We therefore proceed to some special cases of interest in which further analysis can be made. 4.
SPECIAL
CASES UNDER
FIXED-TIME
RULES
The major difficulty in carrying out computations under fixed-time rules lies in the successive matrix-type multiplications required to compute R(‘)(i, k; A,). (The corresponding computations under reach-state rules are much simpler, as shown by Eq. (21).) In this section, we analyze some special cases in which it is possible to reduce the computational effort materially, provided we require only approximate answers. Throughout this section we will suppose that the mortality parameters e,(h) depend on the history h linearly in the number of previous treatments, where the size and timing of these treatments is reflected in the proportionality constant. Specifically, we will take, as does Neuts, Oi(hJ = d + il3 + ?Zy(U,p) (27) where o is the death rate of an uninfected individual, itl is the death rate due to the particle population, assumed proportional to the number of particles present, and the last term corresponds to the long-range aftereffects of the n earlier treatments, each of dosage p and administered at intervals of length a. We also assume, as does Neuts, that I&J, h,) = Ii/(p).
245
LETHAL GROWTH PROCESS CONTROL
It follows from Eq. (25) that P,(i, j; h,) = P,(i, j, h,) x exp(-
nyt)
(28)
and then from (2) that R(i, k; A,,) = R(i, k; h,) x exp(-nyT). Let R(O)(i,j) = a,,, the Kronecker
R”‘(i, j) = c R”-“(i,
k)R(k, j; ho).
k
LEMMA2.
(29)
delta, and for I > 1 let (30)
If e,(h,) = tli + ny, then
R”)(i, k; h,) = R(‘)(i, k) x exp{ -yT[Zn Proof. As mentioned (31) into (2)to get R(‘+l)(i, k; h,) = c
+ l(f -
1)/2]}.
(31)
in (29), the result is true for I = 1. Substitute
R”‘(i, j)R(j, k; h,+J
j>O
= R(‘+‘)(i, k) x exp
x exp
i [ -yT
-yT i
II (I+ ’ 11 In + ‘q
c
(1 + 1)n + ~
l>l
2
which completes an induction. Case 1 (Neuts’s Model). Here we suppose that 1, = m;l for m = 0, 1,. . .; that 0,(/z,) is as specified in (27), suppressing for the moment the dependence of y on (a, p); and that the particle treatment survival follows a binomial distribution B(n, p, i) = j,(nn; After substitution,
@l
-
O
PY,
Eq. (25) at ho reads
ew{-Cc+ (2 + Wit>
jei(l- exp[-(l
(32)
+ Qt]}j-i,
and from 1x1 < 1, we get C P,(i,j;
ho)
jai =
exp( - at)
From Eq. (24), defining Aij(m)
exp[-(A 1 -
[i/(n
+
+ O>t]
e)] x [i - exp[ -(A + e)t]
i.
(33)
Aij(m), we have
= (- l)m-i
‘-I
(j-
l)!i
i!(m - i)!(j - m!)’
(34)
246
ELLEN A. CHERNIAVSKY AND HOWARD M. TAYLOR
We may now use Eq. (33) in (23) or (34) in (26) to obtain tive expressions for L(i, h): T L(i,
It,)
=
s0
evC-(a
nr)tl
+
exp[-(2 ’
1 -
+ Qt]
[A/Q + 0)] x [l - exp[-(A j-i
1
= IJ
+
m
1
zz 0
+
ny
+
[c + ny
x(1 - exp{ -[c
+
+
e)m ( 1 + e > ny
+
(A
i)!(j -
+
(m
e)m]T})I
dt
-
i)!
1
JTI)!
+ (2 + O>m]T})
(m/2)-i il(_l)m-i
I (2
+ Qt]]
(j-l)!i(-I)“-’
i!(m -
ny + (1 + O)m
x (1 - exp{ m=i 4Y
the two alterna-
1
_[ 2 (&T”],
(35)
where I,,, is the modified Bessel function. At least in principle, we could substitute this expression for L(i, 71,) into the basic formula for the expected lifetime given by Eq. (3). In practice, the resulting expressions are too complex for analysis, the main difficulty being that the expression for R”‘(i, k; h,) does not simplify. It is possible to obtain relatively simple formulas approximating the expected lifetime starting from approximations to L(i, h,). First note from L(i, 12,) = E[min{i,
T}IX(O) = i]
that T x Pr{{ > TIX(0)
=
i;
h,}
6
L(i,
h,)
<
T.
(36) where the probability on the left is given in Eq. (33). The probability is of the form ctli and so is the coefficient of T on the right, for C = c( = 1. Thus we try simplification of Eq. (3) under the substitution ccdT for L(i, A,). Where W(i,j) is defined as in (30), we have
C
R(‘)(i,j)d
=
j =
f
C
k=O
jai
&
fYi,j;
P(i,j;
hoXv
h,)B(j,
+
p,
k)ctk
dj
/
= (ccp + q)i exp( -oT) ’
I
exp[--(A +
e>T]
1 - [up + q] x [A/i + e] x (1 - exp[-(R.
= exp(-07)
+ e>T]} Y (37)
241
LETHAL GROWTH PROCESS CONTROL
where 4=1--P, a = q(A + 0) exp[ -(A + WI, b = ~(1, + (3) exp[ - (A + WI, c = 2 + fil- ql{l - exp[-(3. + O)T]}, and d = --@[I
- exp{-(1.
+ QT}].
Thus, define M,(z)
= z
and
M,+,(z)
= M,
’
Then 1 R(‘)(i,j)aj j It follows by induction
= exp(-aT)M,(a)‘.
(38)
that T R”‘(i, j)a’ = exp( - loT)M,(cr)‘.
Finally,
from Eqs. (3) and (31)
fW0, ho) z ?’ f I=0
(‘3$(P, h.)} n=O
x exp
y(Z -
1); -
where c = 1 = c( for an upper bound
I>
la T
x cM,(c()~O,
(39)
and c = exp( - NT) and
exp[ -(A + e)z-1 cd = 1 - [n/(2 + 811 x { 1 - exp[ -(A + O)T]> = M1(1) for a lower bound. Case 2 (Growth /Yz= 0, 1,. . .; that and does not vary parameters is often rate. As in Case 1,
by Immigration). In this case we suppose &, = i for is, the rate of growth in the particle process is constant with the particle population size. Such a choice of used to model growth by immigration at a constant we keep f3,(h,) = d + if3 + ny
and assume the particle treatment survival bution. After substitution, Eq. (25) at k, reads j-i
1
P,(i, j; h,) =
$1 i
- exp(-et)]
x I
fo!lows
exp[ -(a
the binomial
distri-
+ L + ny + &)t] (j -
i)!
(40)
248
ELLEN A. CHERNIAVSKY
AND HOWARD M. TAYLOR
and C P,(i,j;
h,) = exp[-(o
+ ;1 + nr)t]
jai
x exp &l - exp( -et)] 1 From this, L(i, h,) may be computed from Eq. (26) by using =
~
(-ly-i
(42)
(j- m)!(m - i)!’
0
(41)
from Eq. (23), or it may be computed
j-i
Aij(m)
x [exp( -&Ii.
Again, in principle, we could substitute the resulting exact expression for L(i, h,) into the basic formula for the expected lifetime, but again as before, we turn to approximations resulting from
T x Pr{i > TlX(0) = i; h,} < L(i, h,) < T. Let us substitute CCQTfor L(i, A,). We have T R(‘)(i,j)crj
= C P(i, j; h,)(crp + 4)’ jai
=
exp[-(a
+ 1)Tl[B(~)]‘exp{B(~)
where B(U) = (CX~+ q) exp(B[B,(cr)]. Then by induction
0T).
Set
a [exp(OT) - l]} 0 B,(U) = B(U) and Bk+l(c() =
C R(‘)(i, j)czj = exp[ - Z(a + n)T] j x exp Finally,
a
10
e [exp(eT)
-
I]
i
B,(o))
x
k=l
[4(41’.
from Eq. (3) and (31) we have
x exp
a
i(>
e [exp(OT)
-
11 i Il.(a)) x Cw41N”; n=l
and we use c = M = 1 for an upper bound c = exp $1
- exp(-BT)]}
x exp[-(o
and + A)T]
and
CI = exp(-&)
for a lower bound. 5.
A SPECIAL
CASE UNDER
In this section meters
REACH-STATE
we will examine i_, = mA,
RULES
the observable
m = 0, 1, . . .,
case under
the para-
249
LETHAL GROWTH PROCESS CONTROL
and m = 0, %AkJ = and a binomial distribution $(p), and this probability expected lifetime is M-l
m = 1,2, . . .,
,
for treatment survival. We assume $(p, h,) = does not depend on the previous history. The
m(il+ e)(A&er-No+D
fM(NclT hl) = mzN,
xgvxPY(&J++No
x
P +
(1 -
PI,
=y
:
-
e
[(I- P&-$)-J~-)
where
The infinite
p>“-’
sum can be computed
explicitly awP)Cw
= A + 1 - $(p)(p
fMw)
+ (1 - p)[i/(A
to give an expected lifetime of + @I” -No
-t e)]”
- (1 - p)[A./(A + e)]“}’
where M-l
A=
1 A c m=Nom(n + 0) ( 2 + 8)
m-No
.
We expect that the assumptions on the termination parameters made in this section would be similar in most real cases to Neuts’s assumptions 8, = g + em, since CJis small relative to the typical value of f?m. 6.
OTHER
SPECIAL
MODELS
The fixed-time rules have the unfortunate characteristic of repeating treatment indefinitely even if the particle population has been eradicated. This is not a problem under reach-state rules, but the latter assume the particle process is completely observable, not always a realistic assumption. A model midway between these two assumptions can be created if we assume that whether or not the host is completely cured can be determined; but otherwise, the size of the particle population is unknown. We suppose that if a treatment at level p is initiated, the host is cured with probability g(p) and has an expected cured lifetime of l/o. If he is not cured, but survives the treatment, he is treated again, this process being continued 17
250
ELLEN A. CHERNIAVSKY
AND HOWARD
M. TAYLOR
until he either dies or is cured. Of course, realistically, we permit a certain interval to elapse between treatments. Here we assume this interval is small in comparison to the expected cured lifetime, and can thus be ignored. Let c(p) be the probability that the host survives the series of treatments at level p. He survives the first treatment and is cured with probability g@@(p). With probability Il/(p)(l - g(p)) he survives the first treatment, but is not cured and must be treated again. This leads to the recursion
C(P)= S(PMP)+ $(P)(l - dPMJ) or C(P) =
g(P)*(P) 1 -
ti(P)(l
-
S(P))’
(43)
Under our assumptions, his expected lifetime is c(p)/a. To find the optimal treatment level p, we need merely equate the derivative c’(p) to zero. The resulting equation for the optimal p = p* is dP*w(P:Y
+ d(P*MP*)
= dw)[$(P”r.
Not treating at all is equivalent to p = 0, and to prevent this as a solution we impose the condition c(O+) = 0. In words, we are assuming that when the level of treatment is small, the corresponding probability of ever being cured is correspondingly small. Apparently, little is known that points to appropriate choices for the g and $ functions. We will examine several purely illustrative examples. In all these cases t&p) decreases as p increases, which any reasonable such function should do. (i) The linear case. Suppose g(p) = p and $(p) = 1 - a - p for 0 < p < 1 - a, where a is a fixed constant subject to 0 < a < 1. The optimal p is p* = a* - a. (ii) The power case. Suppose g(p) = p but $(p) = 1 - pa when CIis a fixed constant subject to 0 < a < 1. The optimal p is p* = (1 - ~l)l/~. (iii) The exponential case. Suppose g(p) = P and $(P> = a exp(~PI, where a and b are fixed positive constants. The optimal p is approximately p* = (I - a2)‘12 - (1 - a) A more general model is as follows. We suppose with probability g(p). If complete cure is not to the original size No = the result of a treatment.)
(ba) ’ with many of the same desirable characteristics that treatment at level p cures the individual,
obtained, the particle population size is reduced N. (We may suppose this original size was itself As before, we assume a cured individual has an
LETHAL GROWTH
PROCESS CONTROL
251
expected lifetime of I/a. By analyzing the possibilities, as was done in arriving at Eq. (l), we see that the expected future lifetime satisfies f(N,
h,) = UN,
h,) + C& P(N,J’;
bJMP9 h”)
x (1 - dP))fW,
h+d + y
3
1
[
which may be iterated
to give I)
.
(44)
1 This equation is valid under both fixed-time and reach-state rules, provided the appropriate formulas for L(N, /I,) and P(N, j; hk) are used. The substitution of values as used in Sections 4 and 5 gives relatively simple results. In the fixed-time case we need only know
Q,,(f) = Pr{i > M(O) = N; 41, whence f(N,
ho) =
s
'QoW dt+ f 'ff[Il/(P,
h)Qk(T)l(l
-
g(P))‘-’
l=lk=O
0
x[(l-_g(p))j;Q,(O
dr + F].
(45)
In Neuts’s model, our Case 1 of Section 4, Q,(t) may be obtained from Eqs. (33) and (28), and in our Case 2, from Eqs. (41) and (28). Let us suppose that Q,(t) = Q(t) and I,&, h,) = I&); that is, that the history of previous treatments does not affect the current parameters. Then Eq. (45) becomes f(N) or
=
‘Q(t) dt + Q(W(P) il - dpN’-09 s0
s’
Q(t)
f(N) = O
+ ‘+
1
dt + Q(Mu>Cs(pY~l
1 - Q(W(P)C~ - g(p)] -
(46)
In view of the relative simplicity of Q(t), (46) represents a highly tractable expression. In many situations, it may be more natural to take Q(t) as the data for the problem, rather than compute Q(t) from more basic assumptions. Q(t) is the upper tail probability in the failure distribution, assuming no corrective action, and this may be directly estimable, as is usually the case in reliability work.
252
ELLEN A. CHERNIAVSKY
AND HOWARD
M. TAYLOR
REFERENCES 1 W. Feller, On the integro-differential equations of purely discontinuous Markoff processes, Trans. Amer. Math. Sot. 48, 488 (1940); 53, 477 (1945) [Errata]. 2 M. F. Neuts, Controlling a lethal growth process, M&z. Biosci. 2, 41 (1968). 3 M. Zelen, Application of exponential models to problems in cancer research, .I. Roy. Statist. Sot. A129, 368 (1966).
Control
BIOSCIENCES
13,235-252
of a General Lethal Growth
ELLEN A. CHERNIAVSKY
AND
HOWARD
235
(1972)
Process* M. TAYLOR
Center for Environmental Quality Management and Department of Operations Research, Cornell University, and
Ithaca, New
York
Department
Ithaca, New
York
of Operations
Communicated
Research,
Cornell
University,
by Richard Bellman
ABSTRACT It is assumed that a host carries a particle population, initially of size N,, that increases according to a general birth process. Treatments that remove a random number of particles, depending on the dosage level, may be administered. The treatment history and size of the particle population contribute to the death rate of the host. Both when the size of the population is observable and when it is not, the expected lifetime of the host is derived as a function of the frequency and dosage level of treatments. In some special cases the frequency and dosage level that maximize the host’s expected lifetime are given.
1. INTRODUCTION
AND SUMMARY
We consider a particle population whose size, in the absence of intervention, evolves stochastically as a terminating pure birth process. These particles are carried by, and in turn influence, a host, and termination of the birth process corresponds to death of the host. But we intervene in the process, modifying its evolution at intervals, by administering treatments that remove or destroy a random number of the particles present. Unfortunately, these treatments entail some danger to the host. The problem is to balance the benign and malignant effects so as to arrive at a frequency and intensity schedule for treatments that will maximize the expected lifetime of the host. Neuts [2] considered such a model under the assumptions that (i) the rate of growth of the particle population is proportional to its current size (a Yule process); (ii) the particles contribute to the death rate of the host in direct proportion to the particle population size; (iii) the random * This research was conducted under the auspices of the Center for Environmental Quality Management, Cornell University, 1969-1970, and was supported in part by PHS grant 1 TO1 ES 00130-03 and in part by NSF grant GK-21460. Copyright 0 1972 by American Elsevier Publishing Company, Inc.
236
ELLEN A. CHERNIAVSKY AND HOWARD M. TAYLOR
number of particles destroyed by treatment has a binomial distribution; and (iv) the host is given treatments of constant dosage at fixed and equal intervals of time. We extend Neuts’s work by allowing the evolution of the particle population to be a general pure birth process under general termination parameters. In addition, we allow a general distribution for the number of particles destroyed through treatment. More important, we consider not only periodic treatments, but introduce the concept of treating whenever the particle population reaches a pre-specified size. In a later section we consider a third type of policy in which treatments are given at equal intervals until the host either dies or is cured. In all these cases we derive general expressions for the expected host lifetime. Finally, in a variety of special cases, we further develop, analyze, and simplify these general formulas. 2.
THE MATHEMATICAL
MODEL
We begin by assuming that at time t = 0 there are N,, particles present. If at any instant t there are i particles present, then we suppose that the probability of one new particle being born in (t, t + dt) is Izi dt where {,$, i 2 O> are arbitrary nonnegative parameters, called the birth parameters of the process. We postulate that if the host is alive at time t, the probability that he dies during (t, t + dt) depends only on the number of particles present at time t and on the number of particles present at time t and on the number, dosages, and times of previous treatments, which we summarize as a treatment history h. Specifically, we assume that a host carrying i particles and having a previous treatment history h has a probability e,(h) dt of dying in (t, t + dt). For example, Neuts considers 0,(/z) = c + ip + vQ(a, p) where (T is the death rate of an uninfected individual, ip is the death rate due to the particle population, assumed proportional to the number of particles present, and the last term corresponds to the long-range aftereffects of v earlier treatments, each of dosage p and administered at intervals of length a. Intervention involves treating the host, and we allow several levels of treatment, indexed by p. A host with history h survives a treatment of level p with probability Il/(p, h), assumed known. If he survives, a random number of particles is destroyed. Let B(k, p, i) be the probability that out of a population of k particles, i particles survive a treatment at level p. The global sequence of events begins with N,, particles present in a host having a treatment history h,,, where ho signifies the absence of previous treatments. The population grows stochastically as a pure birth process with birth parameters {&} and termination parameters (ei(h,)>, until either the host dies or a treatment is initiated. In the latter case, suppose
231
LETHAL GROWTH PROCESS CONTROL
that the particle population has grown to size N(, at the time treatment is begun, and that the treatment is at level p. With probability 1 - $(p, ho) the host dies and the process terminates. Otherwise, the host survives and begins a new era, this time with N1 particles present and a history of hl, where N1 is drawn from the distribution B(N& p, .) and h, is the host’s updated history. The particle population now increases stochastically from N1 and the process repeats. The host’s life consists of a random number of such periods punctuated by treatments, during each of which the particle process grows undisturbed as a pure birth process. Throughout, we will refer to such periods as growth eras, or more simply, as eras. Each such era is a period of evolution in a pure birth process, whose starting state and parameters may differ from those of other eras. Each era other than the first begins just after a nonfatal treatment. Each era other than the last ends just before a treatment. Either the host dies in the last era, or survives the last era but dies in the ensuing treatment. It follows from this simple description that if we can compute two things about each era, we can compute the expected lifetime of the host. We need to know (i) P(i,j; h), the conditional probability that the host survives an era and ends it carryingj particles, given that he began the era with i particles and a treatment history h; and (ii) L(i, h), the expected length of an era, given that the host began it with i particles and a treatment history h. Letf(i, h) be the expected future lifetime, from the beginning of some era, of a host that carries i particles and has a treatment history of h. We want to computef(N,, h,). To simplify the final statement, let R(i, k; 4,) =
HP,
h,)
C j3i
W, j;
kJNj,
P,
4.
(1)
Then R(i, k; h,) is the conditional probability that the host survives an era and the ensuing treatment and that after the treatment he carries k particles, given that he began the era with i particles and a treatment history h,. Set R(‘)(i, k; h,) = R(i, k; h,)
and then for I > 1, define R”+“(i,
k; h,) = c R(‘)(i,j;
h,)R(j,
k; h,,,).
(2)
jB0
Then R”)(i, k; h,) is the conditional probability that the host survives 1 eras and the treatment following each era and that after the lth treatment he carries k particles, given that he began the first era with i particles and a history h,. Therefore ‘&oR(l)(No,j, h,)L(j, h,) is the expected length of the Zth era, given that the host began the first era with N,, particles and a
238
ELLEN
A. CHERNIAVSKY
AND HOWARD
M. TAYLOR
history ho. (The length of the Zth era is taken as 0 if the host fails to reach it.) The expected total life is the sum, over 1, of the expected life in the Ith era. Thus
SW0,h0)= JWO, h0) + I=1 f, f W~o,j; hoMj, h). j=O Our next task is to compute ted lifetime from (3).
P(i,j;
(3)
h) and L(i, h), and in turn, the expec-
3. THE GENERAL FORMULA Let X(t) be the size of the lethal particle population at time t. During an era this population evolves as a terminating pure birth process (X(t); 0 < t -c i) where [ is the lifetime of the process and corresponds to the lifetime of the host. Under this assumption, the stochastic process is governed by two sets of parameters: the birth parameters {Ai} and the termination parameters {e,}. (Since we are concerned here only with the process during a fixed era, we have suppressed, for the moment, the dependence on the history of previous treatments. Later, of course, when several eras forming a lifetime are simultaneously under discussion, we will again have to make explicit this dependence on the history.) Let P,(i, j) be the probability that such a process is in statej at time t, given that it began in state i at time t = 0. We mean, of course, P,(i,j)
= Pr{[ > t
and
X(t) =,ilX(O)
= i>,
(4)
so that Pr([ We postulate
> tlX(0)
= i} = x P,(i,j). jai
(5)
for i = 0, 1, . . . that (i) P,(i, i + 1) = hit + O(t), (ii) Pr([ E (0, t)lX(O) = i} = Oit + o(t),
(6)
(iii) P,(i, i) = 1 - (pi + Oi)t + O(t), (iv) P,(i, j) = a,,, where o(t) represents remainder terms of order less than t as t --f 0 and aij = 1 if i = j and 0 otherwise. Under these postulates, the birth parameters {AJ govern the growth of the particle population and the termination or mortality parameters {Oi} govern the probability of the host dying. In most cases ok will be nondecreasing in k; the more particles present, the more likely the death of the host. Note that the birth and mortality parameters do not depend on the age of the host, so that the process is stationary in this sense.
239
LETHAL GROWTH PROCESS CONTROL
Alternatively, a terminating pure birth process with parameters {jei} and {ei} may be characterized as a stochastic process on the states 0, I, . . . where the sojourn time in state i is an exponentially distributed random variable with parameter li + (Ii, all sojourn times being statistically independent, and in which, at the end of a visit to state i the process terminates with probability 0,/(2, + 0,) or else moves to state i + 1, which occurs with the remaining probability %J(ii + ei). Again these terminal decisions are statistically independent and independent of the sojourn times. From this description there follows a simple but useful construction of such a process. Let {Si} be independent exponentially distributed random variables with parameters E.i + ei. Let {Ii} be independent random variables, independent of the {Si> and having the distribution 1
with probability
(0
with probability
Ii =
Ai
;li +
ei’
ei ~~+ ei
integer i. Set
Fix a nonnegative
Vi(i) = 0, while forj
> i set vi(j)
Finally,
= Si + ZiSi+l + . . . + ZiZj+r x
,.
X
Zj_*Sj-1.
(7)
set i
=
v,(co) = 1
(8)
jai
and for t E [0, 0 define x(Z) =
.i
if and only if
V,(j) < t < Vi(.j + 1).
(9)
(We take n::! Z, = 1 in Eq. (8).) Then X(t) is a terminating pure birth process starting from i and governed by the parameters {ii} and {Oi}. We are still considering a fixed era, and the treatment history remains suppressed in our notation. The V,(j) may be interpreted as the random length of time until either the host dies or the particle population reaches size j, given that the era began with i particles present. A terminating pure birth process is a Markov process. In our case the transition function is P,(i,J)
= Pr{t < 5
and
X(t) =,ilX(O)
= i},
as given in Eq. (4). (When we again exhibit that the termination parameters {ei} depend on the previous treatment history 11, we will write
240
ELLEN A. CHERNIAVSKY
AND HOWARD
M. TAYLOR
Note that for t > 0
P,(i, j; h) for this quantity.)
C P,(i,j) < 1. j ai In fact, Q,(t), defined to be the cumulative distribution process lifetime [, is given by
Q,(t) = 1 - Jgif’t(I,_d
(10) function
of the
(11)
= Pr([ < tlX(0) = i}. So far we have described how the process evolves in the absence of intervention, or during the eras between treatments. An era ends with death of the host or when a treatment is begun, say at time T (which may be fixed or random, as yet unspecified). What we want to compute for each era, assuming the era begins with a particle population level i and a treatment history h, is (extending the notation in (4) to possibly random times) P,(i,j; h) = Pr{T < i and X(T) =jlX(O) = i} (12) and L(i, h) = E[min{ [, T} IX(O) = i] (13) where the history h appears on the right through the mortality parameters (0J = {Gi(h)} governing c. To continue the computations we need more information on T, the treatment time, and in particular. we need to know the rule determining T. We consider two types of treatment rules, corresponding to two alternative assumptions about the particle process. First we may assume that the number of the particles carried by the host is known to the observer at all times. In the case, called the observable case, it seems natural to initiate treatment whenever the number of particles reaches or exceeds a prespecified critical level M. That is, when the process is observable, we consider rules of the form T = min{t:
0 < t < 5
and
A’(t) 2 M}
(14)
where we agree to set T = co whenever the set in braces is empty. Alternatively, we may assume that nothing is known about the current level of the particle process, called the unobservable case. The only information available then is the passage of time, which forces the use of Jixed-time rules in which T is set at a fixed prespecified number. Of course, reality probably lies somewhere between these two extremes. The particle process may give some evidence of its presence but not reveal its exact size. The Observable
Case under Reach-State
Rules
Let A4 be the critical state determining the reach-state rule. We may suppose, without a significant loss in generality, that the process starts
241
LETHAL GROWTH PROCESS CONTROL
below M. Then, by the definition of T,whenever we treat before the host dies, that is, whenever T < [, we have X(T) = M. It follows from (9) that P,(i,
M; h) = Pr{T < ilX(O) = i}
(15)
and that P,(i,j; h) = 0 forj # M. Continuing, we reach state M, that is, T < 5, if and only if Zi = Zi+l = ...= ZM_l= 1 in the process representation described by Eq. (7)-(9). Since the Zi’S are independent, the probability is easily evaluated to be
Similarly,
from Eq. (7)-(9) min{T,
[> = V,(M),
(17)
so that Z(i, A) = which evaluates
E[J’i(Wl,
W3)
to M-l
L(i, h) =
c
(1, +
em)-’
(19)
m=i
using the independence of the Zi’s from the Sj’s and that E[S,] = (Aj + (3,)-l. (Note: Throughout, a product over an empty set is interpreted as unity.) Set A(M, P, h) =
f
B(M, P, j)P,(.L M; h)
j=O
=
P?~)~~l kcj (A, +fik(hJ
j;. wf,
(20)
Then A(M,p, 12) is the probability that a host with treatment survives any era other than the first. From (1) and (2) R”‘(i, k; ho) = $(P, h)P,(i,
history
h
M; ho)
l-l x
1
“gl
ACM>
P, k,)t4~,
h)
1
Wf,
P, k).
(21)
From Eq. (3) we may now write an explicit expression for the expected lifetime of the host. The subscript M indicates the observable case under a reach-state rule in which the host is treated when the particle population first reaches M. The expected lifetime is
242
ELLEN A. CHERNIAVSKY
j=O
AND HOWARD M. TAYLOR
m=j
(22) The Unobservable
Case under Fixed Time Rules
Suppose we treat at time T, a prespecified constant. We need to compute P(i, j; h) and L(i, h) for such rules, but since T is fixed, we have L(i,
h) =E[min(i,
T}IX(O)
= i]
T
Pr{i > t]X(O) = i} dt
=
s0
s1 =:
j& P,(4j; h) dt
=
(23)
(1 - Qi(t)} dt.
s
Thus, we need only compute P(i, j; h) = P,(i, j; h) for arbitrary fixed T = t. This is done in the following lemma, which generalizes the formula of Feller [I] to allow termination. LEMMA 1.
with parameters
Let (X(t); 0 < t < i) be a terminating pure birth process {Ai} and {oil subject to %i + Oi # lj + Oj for i # j. Set
Aij(rn) = 1
Then P,(i, j) = Pr{i PAi, j) = i]
> t and X(t) = jlX(0) Aij(m)
expC-(2,
+
for
i=j
for
i < j.
(24)
= i} is given by for
~&I
i < j.
(25)
m=i
X(t) = ilX(0) Proof. The case i = j is easy, since Pr{[ > t and i} = Pr{S, > t} = exp[-(I., + Oi)t]. Thus, suppose i < j. Let G,(i, j) = ,f Aij(m) exp[-(1, m=i
+ e,)t],
i
=
243
LETHAL GROWTH PROCESS CONTROL
The Laplace
transform
is m
Gr(i,j)
=
e-“G,(i,
j) dt
s0
where the second right-hand side is the partial fraction expansion of the third. If G;*(i, j) is the Laplace transform of the derivative; dG,(i,J)/dt, then G;*( i, j) = tG*( i, j), since for i < j. G,(i,j) = 0 Because of this correspondence
it is sufficient
to show
P;*(i, j) = Gi*(i, j)
But, using (9), we have = Pr{ vi(j) < t < Vi(j + I)}
P,(i,.j)
= Pr{ Vi(j) < t} - Pr(V,(j
+ I) < t)
and P;*(i,j)
=
J e-c’ dP,(i,j)
=
Jem@
J e-c’ dPr{I/,(j + 1) < t}
< t} -
dPr{K(j)
@evHV(j)lI - E{ev[-Wj = V*(i,j, 5) - V*(i, j + 1, 5)
+ I)]}
= where
I/*(4 j, 0 = E{exp[- OXj>~>
= =
+
E{eXp[-_(Si E(exp[ - @,I)
8.
=ii+Oi+~x
[
lii+
+
+
Ii
X ’ . .
X
Zj-_ZSj-_l)]}
+ * * * + Ii+l x ... x zj-2sj-l)l)
f&,“*(i +L.L 5)1 1 I
;li
4
ei + 5
’ . .
X E(exp[-~li(Si+1
ii + 8i
= Ai +
ZiSi+l +
Ai + Oi + 5
V*(i + 1, j, 5).
244
ELLEN A. CHERNIAVSKY
AND HOWARD
M. TAYLOR
Then
= I. + :, + 1
$V*(i + Lj, 5) - V*(i + Lj + 1, 01,
1
which iterates to P;*(i, j) = V*(i, j, 5) -
V*(i, j + 1, 5)
5 KZ;=f I, = nj,=& as was to be shown.
+ 8, + 5)’
From (23) we may compute
Again Eq. (25) and (26) could be substituted into Eq. (l)-(3) to give an explicit general expression for the expected lifetime of the host. In theory, we could then maximize with respect to the treatment level p and fixed treatment interval T. In practice, the resulting expression, as in the observable case under reach-state rules, is too complex in general for any but the most rudimentary numerical exploration. We therefore proceed to some special cases of interest in which further analysis can be made. 4.
SPECIAL
CASES UNDER
FIXED-TIME
RULES
The major difficulty in carrying out computations under fixed-time rules lies in the successive matrix-type multiplications required to compute R(‘)(i, k; A,). (The corresponding computations under reach-state rules are much simpler, as shown by Eq. (21).) In this section, we analyze some special cases in which it is possible to reduce the computational effort materially, provided we require only approximate answers. Throughout this section we will suppose that the mortality parameters e,(h) depend on the history h linearly in the number of previous treatments, where the size and timing of these treatments is reflected in the proportionality constant. Specifically, we will take, as does Neuts, Oi(hJ = d + il3 + ?Zy(U,p) (27) where o is the death rate of an uninfected individual, itl is the death rate due to the particle population, assumed proportional to the number of particles present, and the last term corresponds to the long-range aftereffects of the n earlier treatments, each of dosage p and administered at intervals of length a. We also assume, as does Neuts, that I&J, h,) = Ii/(p).
245
LETHAL GROWTH PROCESS CONTROL
It follows from Eq. (25) that P,(i, j; h,) = P,(i, j, h,) x exp(-
nyt)
(28)
and then from (2) that R(i, k; A,,) = R(i, k; h,) x exp(-nyT). Let R(O)(i,j) = a,,, the Kronecker
R”‘(i, j) = c R”-“(i,
k)R(k, j; ho).
k
LEMMA2.
(29)
delta, and for I > 1 let (30)
If e,(h,) = tli + ny, then
R”)(i, k; h,) = R(‘)(i, k) x exp{ -yT[Zn Proof. As mentioned (31) into (2)to get R(‘+l)(i, k; h,) = c
+ l(f -
1)/2]}.
(31)
in (29), the result is true for I = 1. Substitute
R”‘(i, j)R(j, k; h,+J
j>O
= R(‘+‘)(i, k) x exp
x exp
i [ -yT
-yT i
II (I+ ’ 11 In + ‘q
c
(1 + 1)n + ~
l>l
2
which completes an induction. Case 1 (Neuts’s Model). Here we suppose that 1, = m;l for m = 0, 1,. . .; that 0,(/z,) is as specified in (27), suppressing for the moment the dependence of y on (a, p); and that the particle treatment survival follows a binomial distribution B(n, p, i) = j,(nn; After substitution,
@l
-
O
PY,
Eq. (25) at ho reads
ew{-Cc+ (2 + Wit>
jei(l- exp[-(l
(32)
+ Qt]}j-i,
and from 1x1 < 1, we get C P,(i,j;
ho)
jai =
exp( - at)
From Eq. (24), defining Aij(m)
exp[-(A 1 -
[i/(n
+
+ O>t]
e)] x [i - exp[ -(A + e)t]
i.
(33)
Aij(m), we have
= (- l)m-i
‘-I
(j-
l)!i
i!(m - i)!(j - m!)’
(34)
246
ELLEN A. CHERNIAVSKY AND HOWARD M. TAYLOR
We may now use Eq. (33) in (23) or (34) in (26) to obtain tive expressions for L(i, h): T L(i,
It,)
=
s0
evC-(a
nr)tl
+
exp[-(2 ’
1 -
+ Qt]
[A/Q + 0)] x [l - exp[-(A j-i
1
= IJ
+
m
1
zz 0
+
ny
+
[c + ny
x(1 - exp{ -[c
+
+
e)m ( 1 + e > ny
+
(A
i)!(j -
+
(m
e)m]T})I
dt
-
i)!
1
JTI)!
+ (2 + O>m]T})
(m/2)-i il(_l)m-i
I (2
+ Qt]]
(j-l)!i(-I)“-’
i!(m -
ny + (1 + O)m
x (1 - exp{ m=i 4Y
the two alterna-
1
_[ 2 (&T”],
(35)
where I,,, is the modified Bessel function. At least in principle, we could substitute this expression for L(i, 71,) into the basic formula for the expected lifetime given by Eq. (3). In practice, the resulting expressions are too complex for analysis, the main difficulty being that the expression for R”‘(i, k; h,) does not simplify. It is possible to obtain relatively simple formulas approximating the expected lifetime starting from approximations to L(i, h,). First note from L(i, 12,) = E[min{i,
T}IX(O) = i]
that T x Pr{{ > TIX(0)
=
i;
h,}
6
L(i,
h,)
<
T.
(36) where the probability on the left is given in Eq. (33). The probability is of the form ctli and so is the coefficient of T on the right, for C = c( = 1. Thus we try simplification of Eq. (3) under the substitution ccdT for L(i, A,). Where W(i,j) is defined as in (30), we have
C
R(‘)(i,j)d
=
j =
f
C
k=O
jai
&
fYi,j;
P(i,j;
hoXv
h,)B(j,
+
p,
k)ctk
dj
/
= (ccp + q)i exp( -oT) ’
I
exp[--(A +
e>T]
1 - [up + q] x [A/i + e] x (1 - exp[-(R.
= exp(-07)
+ e>T]} Y (37)
241
LETHAL GROWTH PROCESS CONTROL
where 4=1--P, a = q(A + 0) exp[ -(A + WI, b = ~(1, + (3) exp[ - (A + WI, c = 2 + fil- ql{l - exp[-(3. + O)T]}, and d = --@[I
- exp{-(1.
+ QT}].
Thus, define M,(z)
= z
and
M,+,(z)
= M,
’
Then 1 R(‘)(i,j)aj j It follows by induction
= exp(-aT)M,(a)‘.
(38)
that T R”‘(i, j)a’ = exp( - loT)M,(cr)‘.
Finally,
from Eqs. (3) and (31)
fW0, ho) z ?’ f I=0
(‘3$(P, h.)} n=O
x exp
y(Z -
1); -
where c = 1 = c( for an upper bound
I>
la T
x cM,(c()~O,
(39)
and c = exp( - NT) and
exp[ -(A + e)z-1 cd = 1 - [n/(2 + 811 x { 1 - exp[ -(A + O)T]> = M1(1) for a lower bound. Case 2 (Growth /Yz= 0, 1,. . .; that and does not vary parameters is often rate. As in Case 1,
by Immigration). In this case we suppose &, = i for is, the rate of growth in the particle process is constant with the particle population size. Such a choice of used to model growth by immigration at a constant we keep f3,(h,) = d + if3 + ny
and assume the particle treatment survival bution. After substitution, Eq. (25) at k, reads j-i
1
P,(i, j; h,) =
$1 i
- exp(-et)]
x I
fo!lows
exp[ -(a
the binomial
distri-
+ L + ny + &)t] (j -
i)!
(40)
248
ELLEN A. CHERNIAVSKY
AND HOWARD M. TAYLOR
and C P,(i,j;
h,) = exp[-(o
+ ;1 + nr)t]
jai
x exp &l - exp( -et)] 1 From this, L(i, h,) may be computed from Eq. (26) by using =
~
(-ly-i
(42)
(j- m)!(m - i)!’
0
(41)
from Eq. (23), or it may be computed
j-i
Aij(m)
x [exp( -&Ii.
Again, in principle, we could substitute the resulting exact expression for L(i, h,) into the basic formula for the expected lifetime, but again as before, we turn to approximations resulting from
T x Pr{i > TlX(0) = i; h,} < L(i, h,) < T. Let us substitute CCQTfor L(i, A,). We have T R(‘)(i,j)crj
= C P(i, j; h,)(crp + 4)’ jai
=
exp[-(a
+ 1)Tl[B(~)]‘exp{B(~)
where B(U) = (CX~+ q) exp(B[B,(cr)]. Then by induction
0T).
Set
a [exp(OT) - l]} 0 B,(U) = B(U) and Bk+l(c() =
C R(‘)(i, j)czj = exp[ - Z(a + n)T] j x exp Finally,
a
10
e [exp(eT)
-
I]
i
B,(o))
x
k=l
[4(41’.
from Eq. (3) and (31) we have
x exp
a
i(>
e [exp(OT)
-
11 i Il.(a)) x Cw41N”; n=l
and we use c = M = 1 for an upper bound c = exp $1
- exp(-BT)]}
x exp[-(o
and + A)T]
and
CI = exp(-&)
for a lower bound. 5.
A SPECIAL
CASE UNDER
In this section meters
REACH-STATE
we will examine i_, = mA,
RULES
the observable
m = 0, 1, . . .,
case under
the para-
249
LETHAL GROWTH PROCESS CONTROL
and m = 0, %AkJ = and a binomial distribution $(p), and this probability expected lifetime is M-l
m = 1,2, . . .,
,
for treatment survival. We assume $(p, h,) = does not depend on the previous history. The
m(il+ e)(A&er-No+D
fM(NclT hl) = mzN,
xgvxPY(&J++No
x
P +
(1 -
PI,
=y
:
-
e
[(I- P&-$)-J~-)
where
The infinite
p>“-’
sum can be computed
explicitly awP)Cw
= A + 1 - $(p)(p
fMw)
+ (1 - p)[i/(A
to give an expected lifetime of + @I” -No
-t e)]”
- (1 - p)[A./(A + e)]“}’
where M-l
A=
1 A c m=Nom(n + 0) ( 2 + 8)
m-No
.
We expect that the assumptions on the termination parameters made in this section would be similar in most real cases to Neuts’s assumptions 8, = g + em, since CJis small relative to the typical value of f?m. 6.
OTHER
SPECIAL
MODELS
The fixed-time rules have the unfortunate characteristic of repeating treatment indefinitely even if the particle population has been eradicated. This is not a problem under reach-state rules, but the latter assume the particle process is completely observable, not always a realistic assumption. A model midway between these two assumptions can be created if we assume that whether or not the host is completely cured can be determined; but otherwise, the size of the particle population is unknown. We suppose that if a treatment at level p is initiated, the host is cured with probability g(p) and has an expected cured lifetime of l/o. If he is not cured, but survives the treatment, he is treated again, this process being continued 17
250
ELLEN A. CHERNIAVSKY
AND HOWARD
M. TAYLOR
until he either dies or is cured. Of course, realistically, we permit a certain interval to elapse between treatments. Here we assume this interval is small in comparison to the expected cured lifetime, and can thus be ignored. Let c(p) be the probability that the host survives the series of treatments at level p. He survives the first treatment and is cured with probability g@@(p). With probability Il/(p)(l - g(p)) he survives the first treatment, but is not cured and must be treated again. This leads to the recursion
C(P)= S(PMP)+ $(P)(l - dPMJ) or C(P) =
g(P)*(P) 1 -
ti(P)(l
-
S(P))’
(43)
Under our assumptions, his expected lifetime is c(p)/a. To find the optimal treatment level p, we need merely equate the derivative c’(p) to zero. The resulting equation for the optimal p = p* is dP*w(P:Y
+ d(P*MP*)
= dw)[$(P”r.
Not treating at all is equivalent to p = 0, and to prevent this as a solution we impose the condition c(O+) = 0. In words, we are assuming that when the level of treatment is small, the corresponding probability of ever being cured is correspondingly small. Apparently, little is known that points to appropriate choices for the g and $ functions. We will examine several purely illustrative examples. In all these cases t&p) decreases as p increases, which any reasonable such function should do. (i) The linear case. Suppose g(p) = p and $(p) = 1 - a - p for 0 < p < 1 - a, where a is a fixed constant subject to 0 < a < 1. The optimal p is p* = a* - a. (ii) The power case. Suppose g(p) = p but $(p) = 1 - pa when CIis a fixed constant subject to 0 < a < 1. The optimal p is p* = (1 - ~l)l/~. (iii) The exponential case. Suppose g(p) = P and $(P> = a exp(~PI, where a and b are fixed positive constants. The optimal p is approximately p* = (I - a2)‘12 - (1 - a) A more general model is as follows. We suppose with probability g(p). If complete cure is not to the original size No = the result of a treatment.)
(ba) ’ with many of the same desirable characteristics that treatment at level p cures the individual,
obtained, the particle population size is reduced N. (We may suppose this original size was itself As before, we assume a cured individual has an
LETHAL GROWTH
PROCESS CONTROL
251
expected lifetime of I/a. By analyzing the possibilities, as was done in arriving at Eq. (l), we see that the expected future lifetime satisfies f(N,
h,) = UN,
h,) + C& P(N,J’;
bJMP9 h”)
x (1 - dP))fW,
h+d + y
3
1
[
which may be iterated
to give I)
.
(44)
1 This equation is valid under both fixed-time and reach-state rules, provided the appropriate formulas for L(N, /I,) and P(N, j; hk) are used. The substitution of values as used in Sections 4 and 5 gives relatively simple results. In the fixed-time case we need only know
Q,,(f) = Pr{i > M(O) = N; 41, whence f(N,
ho) =
s
'QoW dt+ f 'ff[Il/(P,
h)Qk(T)l(l
-
g(P))‘-’
l=lk=O
0
x[(l-_g(p))j;Q,(O
dr + F].
(45)
In Neuts’s model, our Case 1 of Section 4, Q,(t) may be obtained from Eqs. (33) and (28), and in our Case 2, from Eqs. (41) and (28). Let us suppose that Q,(t) = Q(t) and I,&, h,) = I&); that is, that the history of previous treatments does not affect the current parameters. Then Eq. (45) becomes f(N) or
=
‘Q(t) dt + Q(W(P) il - dpN’-09 s0
s’
Q(t)
f(N) = O
+ ‘+
1
dt + Q(Mu>Cs(pY~l
1 - Q(W(P)C~ - g(p)] -
(46)
In view of the relative simplicity of Q(t), (46) represents a highly tractable expression. In many situations, it may be more natural to take Q(t) as the data for the problem, rather than compute Q(t) from more basic assumptions. Q(t) is the upper tail probability in the failure distribution, assuming no corrective action, and this may be directly estimable, as is usually the case in reliability work.
252
ELLEN A. CHERNIAVSKY
AND HOWARD
M. TAYLOR
REFERENCES 1 W. Feller, On the integro-differential equations of purely discontinuous Markoff processes, Trans. Amer. Math. Sot. 48, 488 (1940); 53, 477 (1945) [Errata]. 2 M. F. Neuts, Controlling a lethal growth process, M&z. Biosci. 2, 41 (1968). 3 M. Zelen, Application of exponential models to problems in cancer research, .I. Roy. Statist. Sot. A129, 368 (1966).