- Email: [email protected]
A stochastic model of tumor growth
A Stochastic Model of Tumor Growth*
FLOYD
B. HANSON
Depurrment
of Muthemutics,
Unioer.sit), of Illinois ut Chicago,
Chicago, I/knots 60680
CHARLES Deparment
TIER of Muthematics,
Unioerstty of Illtnois at Chicugo.
Chicago, Illinois 60680
Ud
Department of Engrneering Sciences und Applied Muthemutics Northwestern Unroersi!v, Eomston, Illinois 60201 Received 9 November
1981
ABSTRACT A stochastic
model
for tumor
growth
is derived
continuous-time, density dependent branching deterministic part. For the diffusion process,
as a diffusion
approximation
of a
process with a Gompertz growth law as the the conditional probabilities of extinction
reaching a size c, and doubling are computed along with the expected time of these events. The results are given in terms of integrals, which are evaluated by numerical methods that account
for the logarithmic
singularity
introduced
by the Gompertz
growth
law. When the
variance of the branching process is small compared to the deterministic term, simplified asymptotic expressions are given using methods that are modified for the Gompertzian logarithm. The results are used to find the probability of implant take in limiting-dilution assay experiments and the probability of a tumor becoming detectable in carcinogenesis.
1.
INTRODUCTION
Tumor growth models describe the evolution of the size of a tumor which is assumed to originate from an initial transformed, or progenitor, cell or from a multicell implant. We shall use the number of tumor cells, N(l), as a measure of tumor size. The number of tumor cells can usually be estimated indirectly using measurements of volume, weight, or chemical markers of the
*Research supported by the National Science Foundation and MCS 81-01698 and D.O.E. Grant ACOZ-78 ERO-4650. MATHEMATICAL
BIOSCIENCES
61:73-100
QElsevier Science Publishing Co., Inc., 1982 52 Vanderbilt Ave., New York. NY 10017
under Grants
(1982)
MCS 79-O I7 18
73
0025-5564/82/07073
+ 28$02.75
FLOYD B. HANSON AND CHARLES TIER
74
tumor. General deterministic
tumor growth models are of the form
(1.1) where g is the specific growth rate. If g is constant, the growth is density independent and exponential. When g = g(N), i.e. depends upon N, the growth is density dependent. Usually the specific growth rate g decreases as N increases. This can result from factors such as the immunological response of the host or the inability of solid tumors to efficiently supply nutrients and remove wastes as their size increases. A widely used deterministic tumor growth law is the Gompertz growth law
(1.2) Here K is the largest tumor size and 1/b is the length of time required for the specific growth rate to decrease by a factor of l/e, i.e. the e-folding time. The solution of the Gompertz equation (1.2) with initial value N(0) is
N(t) = N(0) exp [ In &$1-e?],
(1.3)
whose graph is a sigmoid curve that approaches K as t gets large. In the literature, the parameters ln[ K/N(O)] are sometimes replaced by the constants A/b, where A is the instantaneous growth rate at the first observation. In Figure 1, we display a graph of the solution (1.3) using the parameters for a Fortner plasmacytoma 1 in hamsters adapted from data given in Table 4 of [ 121. We have modified the data by converting the adjusted tumor weight to cell number using the estimate that 1 mg contains approximately IO6 cells. In addition, we have changed the initial value N(0) to place the experimental initial value at time zero instead of the value obtained using the method of least squares. This only results in a translation of the curve in time. From (1.3) we can compute the time required for the tumor to pass from the initial size N(0) to c> N(O), which we give as
(1.4)
The first-passage-time formula (1.4) is valid only if N(0) -C C-C K. We illustrate T using (1.4) as a function of N(0) in Figure 5, which appears in Section 3. We notice that T,, becomes very large as N(0) + 0.
STOCHASTIC
TUMOR
75
GROWTH
FIG. I, The tumor size N(t) given by the Gompertz
growth
b = 0.107, K = 2.93X lOlo, and N(0) = 1.84X IO’, adapted displays
the sigmoid
shape characteristic
of density
dependent
law
(1.3)with parameters
from [12. Table 41. The graph growth
and approaches
K as
1-w.
The tumor doubling time T,, the length of time required for a tumor of size N(0) to reach size 2N(O), can be found from (1.4) by setting c = 27?(O), which gives
T,=tln
ln[ VW)1 ln[ K/2 N(O)]
1’ O
Here T, depends upon the initial size N(0). The tumor doubling time goes to 0 as N(0) goes to 0, indicating a rapid initial growth, and increases with increasing N(O), indicating a slowing of the growth rate (see Figure 6 in Section 3). The doubling time becomes unbounded as N(0) + K/2. This should be compared with the doubling time for exponential growth, which is constant for all sizes. The works of Laird [6], Burton [l], Simpson-Herren and Lloyd [ 121, Sullivan and Salmon [ 151, and Steel [ 141 demonstrate the applicability of the Gompertz growth law to tumor growth. Their results are based upon curve fitting with actual data. When tumor data involves a wide range of sizes, the Gompertz law gives a better data fit than other growth laws such as the
76
FLOYD
B. HANSON
AND CHARLES
TIER
logistic model [14]. There seems to be no strong physical or biological argument to explain why the Gompertz model fits actual tumor data best. Thus the frequent use of the Gompertz model is empirical, based upon data fitting. The Gompertz model is deterministic in that the tumor size at any time t is given by the function (1.3). In real situations, tumor cells are subject to irregular growth due to random events (see Figure 2). The irregular growth can result in tumor sizes that are different than those predicted by the deterministic models. The irregular growth of tumors appears to be the rule rather than the exception. For example, deterministic tumor models alone cannot explain tumor implantation experiments in which only a certain proportion of the implants take or grow. To account for irregular growth, stochastic models of tumor growth have been introduced, particularly in models of carcinogenesis. Iverson and Arley [4] describe the growth of a transformed cell, a progenitor of a tumor, by a pure linear birth process. In their model, the probability of a birth is constant, which is analogous to a constant specific growth rate and hence a density independent model. Neyman and Scott [lo] used a linear birth-anddeath process to describe tumor growth. Their model used constant birth and death probabilities and hence is also density independent. A stochastic model for the growth of solid tumors was developed by Wette, Katz, and Rodin [18] based upon physical characteristics of the tumor. This leads to a density dependent stochastic process for the mean tumor size. Their process does not seem to be related to the Gompertz model, but has similar qualitative aspects.
G E 2
I -
.Olo Time
100. in days
200.
FIG. 2. Examples of irregular growth,of tumors from experimental data sketched from [14, Figure 1.51. Curve A is for an R39 sarcoma. The other curves are for the growth of primary BMBA-induced mammary carcinomas (for details see [ 141).
STOCHASTIC
TUMOR
77
GROWTH
Smith and Tuckwell [ 131 consider a random Gompertz law in which the parameter b is replaced by a Gaussian white-noise process. The resulting stochastic process is analyzed, and the mean and variance of the logarithm of the population size are computed. The model of Smith and Tuckwell does not allow for tumor extinction, and the exact meaning of the noise in the parameter h is uncertain. Recently, Prajneshu [ 1 I] analyzed a model similar to that of Smith and Tuckwell with heredity effects. Dubin [2] formulates a density dependent birth-and-death process to describe tumor growth subject to immunological response. The density dependence is due to a quadratic factor (nonlinearity) in the transition probability of the death of a tumor cell. The deterministic part of Dubin’s model is similar to the logistic growth law. Dubin uses various approximations to obtain information about the process. Swan [ 161 describes a method for obtaining the exact solution to Dubin’s model in his review of tumor models. We propose a stochastic model of tumor growth which is the diffusion limit of a continuous-time, density dependent branching process. In Section 2, we give a brief derivation of the branching process and its corresponding diffusion limit. We find that the infinitesimal mean and variance of the limiting diffusion process are M(x)
= bxln( K/x),
V(x) =ax.
(1.6)
Here x is a scaled population size, and the infinitesimal mean, which is the deterministic part of the process, is the Gompertz growth law (1.2). The variance arises from the stochasticity of the branching process. We have chosen a specific branching process that has a Gompertz growth law as its deterministic part in order to be consistent with empirical fits to actual data. We have no physical or biological motivation for this branching process. We feel that this model is useful, since it involves the Gompertz growth law, density dependence and allows for tumor extinction. In Section 3, we analyze this model by computing the conditional probabilities of the tumor becoming extinct and of the tumor reaching some upper boundary point, say c. The probability of reaching c will be interpreted as the probability of a tumor becoming detectable or the probability of take of a tumor implant. In addition, we shall compute the expected time for the above conditional events along with the expected doubling time. The results are computed exactly in terms of integrals by solving problems involving the backward Kolmogorov equation, described in Appendix 1. The integrals cannot be evaluated exactly, so asymptotic approximations are found when the deterministic part is large relative to the variance of the branching process. We measure this relation through the nondimensional parameter
78
FLOYD
B. HANSON
AND CHARLES
TIER
R = 2 bK/a B 1. In Appendix 2, we describe the nondimensional parameter R and a method for estimating the size of the parameter a in the variance of the branching process. The approximations are used to obtain useful information about the stochastic model. These approximations have been used in ecological models by Ludwig [8, 91 and by Hanson and Tier [3, 171. When R is not large, the integrals in the exact result are computed numerically. The numerical computations are compared with the asymptotic results when R Z+1 as a check of accuracy. Both the asymptotic and the numerical computations require special mathematical treatment due to the logarithm in the Gompertz growth law (1.6). These approximations are described in Appendix 3. Finally, in Section 4, we apply the results to. two applications in cancer research. First, we compute the probability of take in limiting-dilution assay experiments. The results are compared with a widely used theory. Next, we use the model for the growth phase in carcinogenesis. We compute the ultimate probability of finding a detectable tumor and the expected time from exposure to a carcinogen until the appearance of a detectable tumor.
2.
THE DIFFUSION
APPROXIMATION
We shall now heuristically derive a diffusion approximation to a continuous-time, density dependent branching process for the number of tumor cells, N(t). Here N(t) is a discrete-valued variable. The details of a rigorous proof of the convergence of density dependent branching processes to diffusion processes is given by Lipow [7]. Our interest is not in the rigorous proof of convergence, but in the analysis of the limiting diffusion process. Let the number of tumor cells at time t + h be given by N(t)
iv(t+h)=
2
B,,
(2.1)
t=l
where h is small. Here B, is a random variable equal to the number of cells resulting from tumor cell i, per unit time in (t, t + h). B, = 0 if the cell dies, = 1 if there is no change, and = 2 if division occurs. In density dependent models, B, is a function of the current tumor size N(t) as well as h. We shall assume that the random variables B, are independent and identically distributed with mean and variance
E[B,IN(t)=n]
=m,(n;h),
(2.2) Var[BfIN(t)=n]
=o,(n;h).
STOCHASTIC
TUMOR
GROWTH
79
Using (2.1) with (2.2), we can compute the expected value and variance of thechangeinthetumorsize,AN(t)=N(t+h)-N(t),as
E[ANIN(t)=n]=E =n[m,(n;h)-1]
(2.3)
and
Var[AN(N(t)=n]
=E[[AN-~(~,-~)]‘IN(~)=~]
=‘[/j,iB,-m~)]~] =i i &B,-m,)@,-d] r=i
]=I
= nv,( n; h).
(2.4)
Here we have used the independence of the random variables B,, i.e. Covar[ B,, B,] = 0, to simplify (2.4). To formally obtain a diffusion approximation to the random process defined by (2.1), we must transform the discrete-variable process into a continuous one. A continuous variable is defined by X( t’) = N( Lt’)/L, where L is some large reference tumor size, such as the detectable size of the tumor, and t’ = t/L is a new time variable. The diffusion approximation also requires that small changes in X(t’) occur in small time intervals, i.e., At’=h/L is small and m,(n;h)=l+[g(h/L)+o(l)]h/L is near replacement value. Also the “offspring” variance is assumed to be of the form v,( n; h) = [a + o(l)]h. The diffusion approximation is then obtained as At’- 0. The limiting diffusion process X(t) is described by its infinitesimal mean and variance,
w”)=A~y
0
E[AXjX(t’)=x] At’
’
Var[AXIX(t’)=x] V(x) = ,‘,im At’ ‘-0 Introducing
the coupled
scalings x = n/L
and t’= t/L
(2.5) . in (2.3) and (2.4),
80
FLOYD
then letting At’-
B. HANSON
AND CHARLES
TIER
0, we obtain
M(x) = xdx), v(x) = ux.
(2.6)
We note that the derivation of the infinitesimal mean M(x) is analogous to the derivation of the deterministic model (1.1) except that the scaling from discrete to continuous variables is usually ignored. Finally, if we choose M(x)=bxln(K/x), V(x) = ax,
(2.7)
then we have a stochastic Gompertz model with the stochasticity due to the randomness of the branching process. We note that a in (2.7) could also depend on x. In Section 3, we shall analyze the diffusion approximation and evaluate some important probabilistic quantities for it. 3.
MATHEMATICAL
ANALYSIS
We are interested in studying the diffusion process for the scaled tumor size X(t), which is defined in the interval [0, c]. Here 0 represents tumor extinction and c has different interpretations depending upon the application. For example, c could represent the size at which a tumor becomes detectable or the size at which an implant is said to have taken. Of particular interest are the conditional probabilities of extinction or reaching c, the conditional probability of doubling, and the expected times of these events. To determine these quantities, we need to study problems involving the backward Kolmogorov operator for the process X(t), which is Lu=;V(x)u,,+M(x)u,. Here M(x) and V(x) are the infinitesimal mean and variance, which are derived in Section 2 and given by M(x)
(3.1) respectively,
= bxln( K/x),
V(x) = CIX.
(3.2)
In the above model, tumor extinction is possible, since x = 0 is a singular exit boundary point. The boundary point x = c is a regular boundary point (see [5]). In addition, we assume that c K K, i.e., c is much smaller than K, which can be interpreted for example as the size of a detectable tumor being much smaller than the largest size to which the tumor can grow, K. This appears to be a reasonable assumption for most tumors.
STOCHASTIC
CONDITIONAL
TUMOR
GROWTH
81
PROBABILITIES
Let P,(n) be the conditional probability of a tumor reaching x = c before extinction, given an initial size x. Then, as is shown in Appendix 1, P,(x) satisfies ~V(x)P:‘(x)+M(x)PC~(x)=O, PC(O) = 0,
(3.3)
P,(c)=l.
The solution of (3.3) with M(x) and V(x) given by (3.2) is x
_2hMK/r)+
t
I 1exp1 P<(X)= O( _2hlln(K/t)+r dr’ I Jexp1 a
dt
(3.4)
a
0
If P,(x) is the conditional probability of reaching extinction, given x and that the tumor has not reached x = c, then it satisfies (3.3) but with boundary conditions
PO(O) = 1,
P,(c)=O.
(3.5)
The solution of (3.3) with (3.5) is
(3 6)
In general the integrals in (3.4) and (3.6) cannot be evaluated explicitly. In such cases, we must either evaluate the integrals numerically or use an asymptotic approximation technique to simplify them. We shall describe asymptotic approximations of (3.4) and (3.6) valid when the deterministic part is large relative to the variance, i.e. R = 2 hK/a B 1. These results will then be compared with numerical computations. The details of the asymptotic approximations and the numerical computations are given in Appendix 3. In addition we assume that c< K. We find that for R B 1, (3.4) becomes
P<(X)-
I-exp(-x/x,,) I- exp( - c/x,,)
’
(3.7)
FLOYD
82
B. HANSON
AND CHARLES
TIER
where x,, is a critical tumor size that is the smaller solution of the transcendental equation (3.8) For R za 1, the smaller root of (3.8) is approximated
by
x,, N -In R < 1, K R
(3.9)
The critical tumor size x,, plays an important qualitative role in the approximation of the conditional probabilities. The conditional probability P,.(x) given by (3.7) can be simplified further to give P,(x)
-1 -exp(-
x/x,,),
x,,(x(c.
(3.10)
We feel that (3.10) is useful in applications such as detection of tumors and take probabilities of implants, since for many tumors x,, a CK K. The probability of tumor extinction, P,(x), can be approximated using (3.7), since Pa(x) = 1 - PC(x), which simplifies as in (3.10) to pa(x)-cxp(-x/x,,),
FIG. 3.
The conditional
P~.num, and asymptotically
with c = 1 X 10”.
probability f,(x) evaluated using (3.10) with label P,,,,,,
x,,
(3.11)
numerically from (3.4) with label for the case described in Figure I
STOCHASTIC
TUMOR
83
GROWTH
In Figure 3, we illustrate the conditional probability of exit at x = c using the parameters for the Fortner plasmacytoma described above and in Appendix 2. The probability PC(x) is computed numerically from the exact solution (3.4) and by using the asymptotic approximation (3.7). For this value of R the numerical and asymptotic results are in close agreement. For a smaller value of R the agreement between the numerical and asymptotic results should not be as good. The conditional probability of doubling, P,(x), is the probability that a tumor with initial size x doubles before extinction. This can be found by evaluating PC(c/2) using (3.7) for R z~ 1 and then replacing c by 2x to get
P,(x) -
1
l+exp(-x/x,,).
(3.12)
Here x represents the initial tumor size, and x,, is the critical tumor size. Thus, for large R, the conditional probability of doubling is between f and 1 (see Figure 4). We note that P,(O) = 4, because for small x both extinction and doubling have nearly equal probabilities, as in a simple random walk.
FIG. 4. The probability P D.““rnTand asymptotically ure I.
of doubling P,(x), evaluated numerically using (3.7) with label using (3.12) with label PD.,,,, for the case described in Fig-
84
FLOYD
B. HANSON
AND CHARLES
TIER
The larger the value of R, the more rapidly this approaches 1 as the initial size increases. For sufficiently large x, the tumor is almost certain to double before extinction. This is very much like the deterministic models, since the effect of the randomness is reduced for larger tumor sizes. EXPECTED
TIMES
Another quantity useful in describing tumor growth is the expected, or mean, time for the conditional events described above. These are the expected time for the tumor to reach c before extinction and the expected doubling time given that the tumor does not go extinct. The expected time of the conditional events is computed from a conditional process which includes only realizations of the original process satisfying the conditioning. Since the conditional process is defined by excluding the realizations not satisfying the conditioning, the total probability of the conditional event is less than one. This is assuming that the discarded realizations have a nonzero probability. Thus, the conditional process must be normalized to make the conditional event a proper random event. Once a proper random event is defined, the expected time of the conditional process is computed as the mean time of the realizations. As described in Appendix 1, these expected times can be computed using the backward Kolmogorov operator (3.1). Let E/f,.] be the expected time for the tumor to reach c, excluding tumors that become extinct. It is given by
E[t,.l = ~,(x)/~,.(x), where M,.(x) satisfies (see Appendix
1)
f~(x)M:‘(x)+M(x)M~(x)= M.(O) = 0,
(3.13)
M,.(c)=O.
-P,(x),
(3.14) (3.15)
Here P,(x) is the conditional probability of reaching c before extinction and is given by (3.4). Using the solution of (3.14) and (3.15) with M(x) and V(x) given by (3.2) then the expected time for conditional exit at x = c can be found from (3.13) and is given by
(3.16)
STOCHASTIC
TUMOR
GROWTH
85
where for convenience
=
R[xln(K/x)+xl K
(3.17)
.
The expected doubling time E[t,] can be found using (3.16) with the initial value x = c/2. Again we shall consider the doubling from x to 2x, so we replace c by 2x in (3.16) to get
(3.18) This result should be compared with the doubling time for the deterministic Gompertz growth law given by (1.5) with N(0) the initial size instead of x. The expressions for the expected times are complicated and cannot be evaluated explicitly. Asymptotic expansions of the multiple integrals (3.16) and (3.18) are also difficult, due to the presence of the growing and decaying exponential terms. Instead we shall find deterministic approximations for E[t,.] and E[t,] using (3.13) and (3.14). When R B 1 and x is not near 0, the second-derivative term in (3.14) is small and PC(x ) - 1 from (3.10). Thus, (3.14) can be written approximately as M(x)kqx)-
-1,
(3.19)
with
M,.(c)=O.
(3.20)
Here we have retained the boundary condition at x = c from (3.15), since our approximation is only valid away from 0. The solution of (3.19) with (3.20) is the deterministic passage time T, given by (1.4), so using (3.13), E[t,.]-M,.(x)-T,.. Similarly,
the expected doubling
time can be approximated
E[t,l -To. Here T, is the deterministic
doubling
time given by (1.5).
(3.21) by (3.22)
86
FLOYD
R. HANSON
AND CHARLES
TIER
Since tumor data are quite variable, the conditions for the above approximations are not always satisfied. For other cases, the multiple integrals (3.16) and (3.18) must be evaluated numerically. The details of the numerical computations are given in Appendix 3. In Figure 5, we graph E[ t,.] and T,. for the Fortner plasmacytoma described above. We notice that E[ t,.] and T,. are close for x away from 0 as expected from (3.21). It is important to note that for small x, E[t,.] approaches a finite value while T is unbounded. The graphs of E[t,] and TD are displayed in Figure 6 for the Fortner plasma-
0
012
014
0:6
018
X, CEL LS/lO” FIG. 5. The deterministic first-passage time q given by (I .4), and the expected time for conditional exit at X= c, E[t,], given by (3.16). for the case described in Figure I with (‘= I x IO’“.
STOCHASTIC
TUMOR
GROWTH
FIG. 6. The deterministic doubling time T, given by (1.5). and the expected time E[ tn] given by (3. IS), for the case described in Figure I,
87
doubling
cytoma. Again, they are in close agreement for x away from 0. Near x = 0, E[t,] has a finite slope, while the slope of To is unbounded. If instead c> K/2, then E[ tn] would be bounded for all values of x, while To is unbounded as x --t K/2. The numerical calculations of Figures 5 and 6 show that the deterministic approximations (3.21) and (3.22) are suitable for R large and x away from 0. This is also true for the parameter values in Appendix 2. 4.
APPLICATIONS
We shall now apply the results in Section 3 to two specific applications in cancer research. First we shall consider limiting-dilution assay experiments and the probability of a successful implant take. Next we shall consider the problem in carcinogenesis of finding the probability and expected time of tumor occurrence after exposure to a carcinogen. LIMITING-
DILUTION
ASSA Y
A limiting-dilution assay is an assay method used to determine the growth and colony growing characteristics of tumors (see [ 141). It can also be used to determine the surviving fraction of treated tumors. The method consists of
88
FLOYD
B. HANSON
AND
CHARLES
TIER
either plating or transplanting into an experimental environment a multicell suspension (implant) derived from a particular tumor. The experiment is repeated using implants of the same size, and the total number of colonies is observed and used as a measure of the proportion or probability of implant take. An implant is said to have taken if it reaches some prescribed size c within a certain time interval. The probability of take should depend upon the growth characteristics of the original tumor. To determine the surviving fraction of a treated tumor, the experiment is performed on both the treated and the untreated tumors. For each tumor, the number of cells in an implant required for a 50% take probability is determined. The ratio of the number of cells required for a 50% take probability in the untreated to that in the treated tumor is the surviving fraction. The 50% take probability is chosen arbitrarily as a standard. The results are then extrapolated to the treatment of tumors in patients. If one uses only the deterministic Gompertz growth law (1.2) or (1.3) as a model for the limiting-dilution assay experiments, it predicts that all implants, no matter how small, will be successful, though the time of take will vary. Thus, the probability of take equals 1. To account for the experimental results, it has been proposed (see [ 141 for a review) that each tumor contains only a small fraction of cells that can be the progenitor of a colony (i.e. clonogenic cells). Thus for an implant to be successful, it must contain at least one clonogenic cell that proliferates. To account for the stochastic effect in the experiments, it is assumed that the probability of obtaining an implant with no clonogenic cells is exponentially distributed, as e-“‘, where m is the average number of clonogenic cells per implant. Thus, the probability of obtaining an implant with one or more clonogenic cells is I- epm, and hence the probability of take is Pr {implant take} = I- epm =
1 -
e-/n’.
(4.1)
In the second equality, the parameter N is the total number of cells in the implant, and f is the fraction of clonogenic cells. The stochasticity only enters in the probability of selecting clonogenic cells in the implant. We should note that it is in general not possible to distinguish between clonogenic and nonclonogenic cells. In the stochastic Gompertz growth law described in Section 2, the stochasticity is due to the branching process, with each cell division giving rise to a random number of proliferating descendents. The stochastic effects are most evident when the implant sizes are small, as in these experiments. The probability of take of an implant is the conditional probability of reaching some prescribed upper limit c, i.e. P,(x). The conditional probabil-
STOCHASTIC
TUMOR
GROWTH
89
ity PC(x) is given by (3.4), or asymptotically and c >>x,, , the probability of take is
for R B 1 by (3.7). When R B 1
P,.(x)-l-exp(-x/x,,). Here x is the initial size of the implant
(4.2)
and x,, is the smaller root of (3.8): (4.3)
The value of x,, depends upon the growth characteristics of the tumor, including b, K, and the stochasticity a. The values of b and K are estimated from empirical data fits using the Gompertz growth law and represent average values over a range of tumor growth. The additional parameter a measures the variance of the birth schedule in the branching process (see Appendix 2). We note that P<(x) represents the ultimate probability of the tumor size reaching x = c before extinction. Finally, we compute the surviving fraction for a treated tumor using the limiting-dilution assay. We shall assume that the growth characteristics of the treated and untreated tumors are different. Thus, the parameters in the Gompertz growth laws of these tumors will be different. For the untreated tumor we shall assume the probability of take is given by (4.2) while for treated tumor it is
P,(x) = 1-exp(
- x/x’,,)
(4.4)
where x& satisfies (4.3) with the treated-tumor growth parameters. We compute the surviving fraction by taking the ratio of the numbers of cells from the untreated and the treated tumor required for a 50% take probability. This gives surviving fraction = x,, /x&.
(4.5)
Thus, the surviving fraction reduces to the ratio of the two critical tumor sizes when R za 1. These critical sizes are functions of the growth parameters for each of the tumors and also of the stochasticity. The formula (4.5) represents the ultimate surviving fraction. If the tumor can eventually repair the damage caused by treatment and regain the growth characteristics of the untreated tumor, then the surviving fraction will ultimately be 1. CA RCINOGENESIS
The development of a tumor is believed to be a two-phase process. The first is the transformation phase, in which normal cells are transformed into progenitors of a tumor. The transformation may result from the contact of
90
FLOYD
B. HANSON
AND CHARLES
TIER
normal cells with carcinogens and is generally believed to occur in a series of stages. In the second or growth phase, the transformed cell grows into a detectable tumor. We shall compute the ultimate probability of finding a detectable tumor given the exposure of normal cells to a carcinogen, and the expected time from exposure to a carcinogen until the appearance of a detectable tumor. Stochastic models have been proposed for both the phases of carcinogenesis (see [ 191). We shall use the stochastic Gompertz model outlined in Section 2 for the growth phase. Let q = Pr {normal cell is ultimately
transformed}
,
(4.6)
which may result from a multistage transformation process. If p is the ultimate probability of a normal cell becoming a detectable tumor, then P = #C(l)
-4[1-w(-l/xcr)l,
R>l.
(4.7)
Here PC(1) is the conditional probability the tumor reaches x = c given that the initial number of cells was 1 (the single transformed cell). It is given by (3.4) or approximately by (3.10). The value of c is the smallest detectable tumor size. Thus, in a collection of J identical, normal cells, the ultimate probability of finding at least one cell that will give rise to detectable tumor is found by assuming that the probability of finding k tumors, Pk, is binomial (see [ 191): J! pk= k!(J-
k)!
pk(l-py.
For J large and p small, this can be approximated
by the Poisson distribution
Pk - e -Jp( Jp)k,‘k!. Thus, the probability
of at least one tumor is 1 - PO, or approximately Pr{ at least 1 tumor} - 1 - eCJP.
(4.8)
The tumor time t = t, + t, is the total time from exposure to a carcinogen until the appearance of a detectable tumor. Here t, is the transformation time and r2 is the growth time. Since the density function for t is a convolution of the densities of t, and t2, then
E[t]=E[t,]+E[t,].
(4.9)
91
STOCHASTIC TUMOR GROWTH
Here the random time variables t, t,, and t, all refer to times of conditional events. Thus E[t] represents the average time from exposure to a carcinogen until the appearance of a detectable tumor given that the normal cell eventually gives rise to a tumor. The expected time E[ t2] is the expected time for the growth phase given that the tumor becomes detectable, and is the same as E[t,.] given by (3.13) and (3.16).
APPENDIX
1.
BACKGROUND
PROBABILITY
In this appendix, we give intuitive derivations of the equations for the conditional probabilities and the expected times described in Section 3. Let q(x, y, t) be the transition probability density function for the diffusion process X(t) = y given that X(0) =x. Then q satisfies the forward equation
4r+J,.M=o
OC.Y
(Al.l)
where
J[ql(xTY,t) = -Iv(Y)ql v+WY)q is the probability equation
flux. The density
function
q,=L[ql=tV(x)q,,+M(x)q,, The initial condition
for the backward
(Al .2)
q also satisfies the backward
ocxcc.
and forward equations
(Al .3) is
limq(x,y,t)=6(y-xx).
(Al .4)
r-0
The coefficients M(x) and V(x) are the infinitesimal mean and variance, respectively; see (2.11). The boundary conditions for (Al.l) and (A1.3) depend upon the particular process (see [5]). Let r,(x) be the exit time of X(t) from (0, c), defined as te( x) = inf {t 1X(t) Then the distribution given by
not in (0, c)] X(0) = x in (0, c)}.
of t, is the probability
(Al .5)
that X(t) is not in (0, c) and is
(Al .6)
92
FLOYD
with a density function,
B. HANSON
the time derivative
AND
of the distribution
CHARLES
TIER
(A1.6), given
by (Al .7) Using (Al. 1) and (Al .2), the density function ~(x,t)=
limJ_Y - c
can be written
,ti_moJ=u~(x,t)-po(x,f)
(Al .8)
and hence has two components: pa(x, t) for exits at x = 0, and p&x, t) for exits at x = c. Since q satisfies the backward equation (Al .3), then so do p( x, t), pO( x, t), and p,,( x, t), which are derived as
p,..,=J[q,l(x,c,t)=J[~[qll Here we have interchanged
differentiation
=L[J[qll =L[JJcl. and integration.
Similarly,
P,=L[Pl, PO,* =L[Pol.
The ultimate probability
of exit from either 0 or c is
p(x) = jg"P(X, t) dt,
(Al .9)
which also satisfies the backward equation (A 1.3), i.e. LP( x) = 0. The proper boundary conditions are P(0) = 1 and P(c) = 1, since a process with initial points 0 or c is outside the interval. These absorbing boundary conditions imply that the ultimate probability of exit is certain, i.e. P(x) = 1. Similarly, the ultimate probability of conditional exit at x = 0 is (A1.lO) and at x = c is (Al.ll) Each satisfies the backward
equation
L[P,l(x) = L[P,lb) =o,
(A1.12)
STOCHASTIC
TUMOR
GROWTH
93
which is Equation (3.3). The boundary conditions for P,(x) are [see (3.5)] P,(O) = 1 and P,(c) = 0, since a process that starts at x = 0 is outside (0, c), while if it starts at x = c it can never exit at x = 0. Using a similar argument, we obtain the boundary conditions for PC(x) which are given in (3.3). The expected time of exit at either x = 0 or c is the first moment of p( x, t) and is given by (A1.13) which after applying the backward P(x) = 1, we find to satisfy
operator,
integrating
LIT]=~rtLIP1dr=jo”~~,dr=tp]~-IXPdt=
by parts, and using
-P(x)=
-1.
0
(A1.14) The boundary conditions are T(0) = T(c) = 0, since it takes no time for the process to exit if its initial point is either x = 0 or c. The expected time of the conditional exits are more complicated, since the probability of conditional exit is in general less than one, and hence conditional exits are not proper random events. To define the conditional exit at x = c to be a proper random event, we normalize its density to get p,.( x, t)/P<( x). This now integrates to 1. The expected time for conditional exit at x = c is then defined by
fk(*. E[ t,.] =
t) dt
p,.(x)
(Al .15)
Unfortunately, this does not satisfy a simple equation. Instead, it is easy to show using the same steps as on T(x) that the conditional moment M,.(x) defined by
w(x)
=pPc(x,t)dt=p,.(x)E[t,.]
(~1.16)
satisfies
Hence, the expected found using (3.14).
time of conditional
exit at x = c satisfies (3.13) and is
94 APPENDIX
FLOYD
2.
ESTIMATE
B. HANSON
OF THE DIFFUSION
AND CHARLES
PARAMETER
TIER
u
We could estimate the parameter a from tumor measurements, using the variance Var[ X] of the process. Var[ X] can be exactly calculated in principle from the density for X(t), but the density must be obtained from the forward Kolmogorov equation with a great deal of effort. Hence we employ an Omstein-Uhlenbeck approximation of the density about the maximum, deterministic tumor size K, by letting y = K( I- z) with z small and using it in the forward equation (Al. 1):
(A2.1)
We note that the dimensionless “signal-to-noise ratio,” R = 2 hK/a, naturally arises when y or x is scaled by K, and t is scaled by 1/b. Setting qt = 0 in (A2. l), we can find a quasi-stationary density (see [5]), which represents the density of those tumors which are not yet extinct. The normalized, quasi -stationary density is approximately q^.(R/27r)“2exp[-Rz2/2],
(A2.2)
with mean K and variance K2/R ‘1’ for the variable y. Hence the coefficient of variation, or relative standard deviation, is CV-
l/R’/‘.
(A2.3)
We have also computed formulas for the CV from Ludwig’s ray method [8, 91 by expanding the exponent of the density about the deterministic trajectory. These computations give the CV as a function of y and agree with (A2.3) when y = K. However, we will report these details in another investigation. Simpson-Herren and Lloyd [ 121 in their Table 4 give average, regression adjusted weights and their standard deviations for Fortner plasmacytoma 1 harvested from hamster hosts. Coefficients of variation computed from their standard deviations and adjusted weights range from 0.224 to 0.930 with a mean of 0.523. The measurements closest to K have the smaller CVs and have an extrapolated value close to 0.300. Given b = O.l07/day and K = 29,300 mg from their data and with CV= 0.300, we find that a rough estimate of a is 600 mg/day = 6 X lo9 cells per day. We note that although the “noise” a is large in an absolute sense, its relative magnitude is small when properly measured in terms of the dimensionless group that actually appears in the theory, i.e. the reciprocal of the “signal-to-noise ratio,” 1/R = a /2 bK = 0.0957. Therefore a small-noise approximation is more appropriate, and we consider this in the next appendix.
STOCHASTIC
TUMOR
APPENDIX
3.
95
GROWTH
ASYMPTOTIC AND NUMERICAL APPROXIMATIONS
In this appendix we derive special asymptotic and numerical approximations for the conditional probabilities and the expected times described in Section 3. The stochastic model with the Gompertz growth law involves a logarithm, which leads to nonanalytic behavior. The standard methods used in asymptotic and numerical analyses are not applicable directly, and hence we use a modified analysis based upon a transformation to remove the logarithm. APPROXIMATION
OF THE CONDITIONAL
PROBABILITIES
We seek an asymptotic approximation of P<(x) given by (3.4) when R = 2 bK/u B 1. Rewriting (3.4) using (3.17), we get
(A3.1)
where G(x)
=
Rx[ln(K/x)+ll K
To eliminate
the logarithm,
we let y = ln( K/x),
P,(x) =
J-“,
s
(A3.2)
.
so that (A3.1) becomes
(A3.3)
)
eHcs)ds
Y
H(y)=-y-G(KeeY)=
(A3.4)
-y-Rep’(y+l),
and Y = ln( K/c). In the spirit of Laplace’s method for the asymptotic expansion of exponential-type integrals, we seek a maximum of the integrands in (A3.3). The critical point y,, of H(y) satisfies H’(.K,) = - 1+ Ry,,evA maximum
(A3.5)
IL> = 0.
occurs at the largest root in (A3.5), which is y,,=lnR+lny,,
-lnR+ln(lnR),
p-
+co.
(~3.6)
96
FLOYD
B. HANSON
AND CHARLES
TIER
Converting (A3.5) and (A3.6) back to the variable x, we obtain (3.8) and (3.9) respectively, for x,,. Testing the second derivative of H(y) at y,,, we find that H”(y,.,) = - 1+ l/y,, - - 1 as R - + 00. Hence, although y,, is the location of a wide maximum, it is not sharp enough to allow the use of Laplace’s method. To obtain the final result, we use the transformation z = ln(x,,/x) on (A3.3) to make it easier to examine its behavior near the maximum of the exponent. The transformation leads to
(A3.7)
where
h(z)=
-z-ePi[l+F],
(A3.8)
and Z = ln( x0/c). Here L = ln( K/x,,), large, we may approximate h(z) by h(z)This approximation
and L ZJ1 when R B 1. With
-z-eeZ,
LBl.
L
(A3.9)
will be valid when e-‘(z+l)/Lal.
(A3.10)
However, when the approximation (A3.9) is applied to (A3.7) the result is uniformly valid for 0
P,.(x)-
1-exp(-e-‘)
(A3.11)
l-exp(-e-Z)
Transforming back to x, we obtain (3.7) or when c/x,,< 1 the further simplified form (3. IO). The results for the probability of doubling, PD( x), are analogous and given by (3.12). QUASI
-DETERMINISTICAPPROXIMATION
OF THE TIME TO EXIT
FROM
(0, c)
A deterministic approximation of E[ t,] given in (Al. 13) by T,, is less suitable because E[ t,] ---r+ 0 as x + +O, while Tc does not behave this way. However, we have found that a “quasi-deterministic” approximation,
-qt,l -
p,.T,,
(A3.12)
STOCHASTIC
FIG. 7. deterministic
TUMOR
97
GROWTH
The expected
exit time E[t,]
given by
(Al. 13) compared
with the quasi-
7; P, given by (A3.12).
works reasonably well, although it consistently underestimates the exit time. E[ r,] is compared with P,.T,. in Figure 7. In summary, we have found good asymptotic approximations for the probabilities of exit (detection) and doubling; and we have found deterministic approximations (3.21) and (3.22) and a quasi-deterministic approximation (A3.12) for the exit or doubling times where more accurate asymptotic approximations are very difficult to compute. NUMERICAl.
TREATMENT
OF LOGARITHMIC
INTEGRALS
A general rule of numerical integration schemes is that the order of the local error depends upon some power of the mesh size h, and upon a corresponding derivative of the functionf( x), such as O( h”+ ‘f(“)(x)). Here n - 1 is the highest polynomial that can be integrated exactly by the numerical scheme. However, when the function has a nearby discontinuity, the actual local error is not of the general form. Instead, n = 0 must be used in the error estimate in almost all cases. Since our problem involves In x, which is not well behaved near x = 0, a special treatment is necessary to improve the local error to better than O(h). We shall transform variables to remove
98
FLOYD
B. HANSON
AND CHARLES
TIER
the logarithm and obtain better-behaved integrals. This will allow us to use the Gaussian integration rules with their high degree of polynomial precision. We shall first analyze P<(x). We let y = ln(x/t) in both integrals in (A3.1), which we rewrite as
xF,(x) 4(c) ’
P,.(x) = ~
(A3.13)
where
F,(x) =/,, e-‘fo(Y;x>dy, fa(y;x)=exp[-G(xe-‘)I. The integrals in (A3.13) now have smooth integrands but are singular in their domain. They can be integrated without problems using a Gauss-Laguerre integration rule [i.e. Gaussian on (0, co) with weight factor e -“I. Some results are shown in Figure 3 for our sample parameters. The probability of doubling, PD( x), can also be evaluated using the same formula. Some results for PO are shown in Figure 4. A similar technique can be used on the integral formula for E[t,.] in (3.16) except that a “catastrophic cancellation” can occur between the two integrals. Even for moderate values of the “signal-to-noise ratio” R, the two integral terms in (3.16) can be almost the same size, and their difference leads to a great loss of significant digits. Fortunately, if we integrate by parts to change the order of integration in (3.16) and then rearrange the result so that positive integrals are added rather than subtracted, we obtain the useful form EL1
,] L
=
2[Q(c)-
Q(x)]
AQ(c)Q(x) 2
+AQ(c)
x Q2b)eGrV'dy
Jo
c.Q(~)~~(~)
Y
[QWQb)l
sx
dy,
tA3.14j
Y
where Q(x)=~xepG(‘)dt
and
Q(c)-Q(x)=/Ce-G(‘Ldt. X
The integral over (0, x) is treated like those in (A3.13), using the logarithm transformation and Gauss-Laguerre rules. The integrals over (x, c) or (y, c) as in Q(c) - Q(y) are still affected by the singularity of In x for small x, even though x is excluded from the domain. Hence on (x, c) [and similarly on (y, c)], we use the transformation
,=1_21n(c/~) ln(c/x)
’
99
FLOYD B. HANSON AND CHARLES TIER
and then use a Gauss-Legendre integration rule [i.e. a Gaussian (- 1,l) with weight 11. Upon transforming the integrals, we obtain
xF,(x; E[t,.] =In(
c)&(x)/F,(x)+
cln( t)F,(x;
rule on
c)/4
f) M,(c) (A3.15)
where
F,(x; c) = j-;,f,(s; x; c) h, f,(s; x; c) = exp{g(s;x; g(s;x;c)= F,(x) f2(t.
x)
-In =im =
c>-G( cexp[g(s;x; c)])},
I-s
(f 12,
e-‘f2( t; x) dt,
FoZ(xe-“*)
x>’ F,(x; c) =j--lh(s; x; c> ds, 2f,(t/2
f3(s;x; c) = Cl-s)exp{g(s; x; c)+G(cedg(s; x; c)l)}
XF,(cexp[g(s;x; c)l)F,(cexp[g(s;x; c)l; c). The integrandsf, are fairly well behaved, so that the numerical integration is straightforward. For the multiple integrals, the inner integrals of functions must be evaluated at each of the nodes of the outer integrals. Some results for E[ t,.] are shown in Figure 5. The numerical computation of the expected doubling time E[ tD] proceeds directly from (A3.15) with c set equal to 2x, which gives
E[ r,] = (xln2)
F,(x;2x)F2(x)/F,(x)+(ln2)F,(x;2x)/2 a&,(2x) (~3.16)
Certain minor simplifications can be made [e.g., g(s; x; 2x) = (ln2)(1- s)/2 is independent of x], but the primary reduction was made in deriving (A3.15). In Figure 6 are shown some results for E[t,]. The numerical decomposition of E[t,] is similar to that of E[ t,.] in (A3.15), and we omit the formula here. Sample results are presented in
FLOYD
100
B. HANSON
AND CHARLES
TIER
Figure 7. The steeper slope near x = 0 than near x = c indicates that the expected time for a tumor to go extinct (remission) is much longer than the expected time to reach x = c (detection of tumor or death of host). REFERENCES I
A. L. Burton, Rate 30:157-176 (1966).
2
N. Dubin,
growth
of solid
tumors
as a problem
A Stochnstic Model for Immunological
and Approximations.
Lect. Notes in Biomath..
of diffusion,
Growth
Feedback in Carcmogenesis:
Ana!vsrs
Vol. 9, Springer,
3
F. B. Hanson and C. Tier, An asymptotic singular diffusion arising in population (1981).
4
S. Iversen and N. Arley, On the mechanism Microbial. Stand. 27:773-803 (1950).
5
S. Karlin and H. M. Taylor, York, 1981.
6
A. K. Laird, The dynamics
7
C. Lipow, Limiting diffusions for population Appl. Probab. 14: 14-24 (1977).
8 9
D. Ludwig,
Stochastic Population Theories, Springer,
D. Ludwig, 17:605-640
Persistence (1975).
IO
J. Neyman
and E. Scott, Statistical
Prajneshu,
12
(1979). L. Simpson-Herren experimental
13
15 16 17 I8 I9
of experimental
of tumor growth,
of dynamical
Gompertz
and
solution to the first passage time problem for biology, SIAM J. Appl. Math. 53:89-l 17 Acra. Path.
Br. J. Cancer 28:490-502 size dependent
New York,
processes,
aspects
of carcinogenesis.
of the problem
model with heredity Kinetics
J.
1914.
perturbations,
Sratisiics und Probability,
New
(1964).
branching
systems under random
H. H. Lloyd,
tumor systems,
carcinogenesis,
SIAM
Rev.
in Fi/rh
Univ. of Calif. Press,
effect, J. Math. Biol. 8: l89- I96
parameters
Cancer Chemo. Reps. 54: 143-174
and
growth
curves
for
(1970).
C. E. Smith and H. C. Tuckwell, Some stochastic growth processes. Problems in Bio/ogv - Victoria Conference, Lecture Notes in Biomath., New York,
14
A stochastic
1976.
A Second Course in Srochastrc Processes, Academic,
Berkeley Symposium on Mathematical Berkeley, 1967, pp. 745-776. I1
New York,
in Muthemarical Vol. 2. Springer,
1974.
G. G. Steel, Growth Kinetics of Tumors, Clarendon. Oxford, 1977. P. W. Sullivan and S. E. Salmon, Kinetics of tumor growth and regression multiple myeloma, J. C/in. Invest. 51: 1697-1708 (1972).
in IgG
G. W. Swan, Some current mathematical topics in cancer research, Univ. Microfilms Inter.. Ann Arbor, 1977. C. Tier and F. B. Hanson, Persistence in density dependent stochastic populations, Mafh. Biosci. 53:89-II7 (1981). R. Wette, I. N. Katz, and E. Y. Rodin, Stochastic processes for solid tumor kinetics, Math. Biosci. 19:231-255 (1974). A. Whittemore and J. B. Keller, Quantitative theories of carcinogenesis, SIAM Rev. 2O:l-30 (1978).
FLOYD
B. HANSON
Depurrment
of Muthemutics,
Unioer.sit), of Illinois ut Chicago,
Chicago, I/knots 60680
CHARLES Deparment
TIER of Muthematics,
Unioerstty of Illtnois at Chicugo.
Chicago, Illinois 60680
Ud
Department of Engrneering Sciences und Applied Muthemutics Northwestern Unroersi!v, Eomston, Illinois 60201 Received 9 November
1981
ABSTRACT A stochastic
model
for tumor
growth
is derived
continuous-time, density dependent branching deterministic part. For the diffusion process,
as a diffusion
approximation
of a
process with a Gompertz growth law as the the conditional probabilities of extinction
reaching a size c, and doubling are computed along with the expected time of these events. The results are given in terms of integrals, which are evaluated by numerical methods that account
for the logarithmic
singularity
introduced
by the Gompertz
growth
law. When the
variance of the branching process is small compared to the deterministic term, simplified asymptotic expressions are given using methods that are modified for the Gompertzian logarithm. The results are used to find the probability of implant take in limiting-dilution assay experiments and the probability of a tumor becoming detectable in carcinogenesis.
1.
INTRODUCTION
Tumor growth models describe the evolution of the size of a tumor which is assumed to originate from an initial transformed, or progenitor, cell or from a multicell implant. We shall use the number of tumor cells, N(l), as a measure of tumor size. The number of tumor cells can usually be estimated indirectly using measurements of volume, weight, or chemical markers of the
*Research supported by the National Science Foundation and MCS 81-01698 and D.O.E. Grant ACOZ-78 ERO-4650. MATHEMATICAL
BIOSCIENCES
61:73-100
QElsevier Science Publishing Co., Inc., 1982 52 Vanderbilt Ave., New York. NY 10017
under Grants
(1982)
MCS 79-O I7 18
73
0025-5564/82/07073
+ 28$02.75
FLOYD B. HANSON AND CHARLES TIER
74
tumor. General deterministic
tumor growth models are of the form
(1.1) where g is the specific growth rate. If g is constant, the growth is density independent and exponential. When g = g(N), i.e. depends upon N, the growth is density dependent. Usually the specific growth rate g decreases as N increases. This can result from factors such as the immunological response of the host or the inability of solid tumors to efficiently supply nutrients and remove wastes as their size increases. A widely used deterministic tumor growth law is the Gompertz growth law
(1.2) Here K is the largest tumor size and 1/b is the length of time required for the specific growth rate to decrease by a factor of l/e, i.e. the e-folding time. The solution of the Gompertz equation (1.2) with initial value N(0) is
N(t) = N(0) exp [ In &$1-e?],
(1.3)
whose graph is a sigmoid curve that approaches K as t gets large. In the literature, the parameters ln[ K/N(O)] are sometimes replaced by the constants A/b, where A is the instantaneous growth rate at the first observation. In Figure 1, we display a graph of the solution (1.3) using the parameters for a Fortner plasmacytoma 1 in hamsters adapted from data given in Table 4 of [ 121. We have modified the data by converting the adjusted tumor weight to cell number using the estimate that 1 mg contains approximately IO6 cells. In addition, we have changed the initial value N(0) to place the experimental initial value at time zero instead of the value obtained using the method of least squares. This only results in a translation of the curve in time. From (1.3) we can compute the time required for the tumor to pass from the initial size N(0) to c> N(O), which we give as
(1.4)
The first-passage-time formula (1.4) is valid only if N(0) -C C-C K. We illustrate T using (1.4) as a function of N(0) in Figure 5, which appears in Section 3. We notice that T,, becomes very large as N(0) + 0.
STOCHASTIC
TUMOR
75
GROWTH
FIG. I, The tumor size N(t) given by the Gompertz
growth
b = 0.107, K = 2.93X lOlo, and N(0) = 1.84X IO’, adapted displays
the sigmoid
shape characteristic
of density
dependent
law
(1.3)with parameters
from [12. Table 41. The graph growth
and approaches
K as
1-w.
The tumor doubling time T,, the length of time required for a tumor of size N(0) to reach size 2N(O), can be found from (1.4) by setting c = 27?(O), which gives
T,=tln
ln[ VW)1 ln[ K/2 N(O)]
1’ O
Here T, depends upon the initial size N(0). The tumor doubling time goes to 0 as N(0) goes to 0, indicating a rapid initial growth, and increases with increasing N(O), indicating a slowing of the growth rate (see Figure 6 in Section 3). The doubling time becomes unbounded as N(0) + K/2. This should be compared with the doubling time for exponential growth, which is constant for all sizes. The works of Laird [6], Burton [l], Simpson-Herren and Lloyd [ 121, Sullivan and Salmon [ 151, and Steel [ 141 demonstrate the applicability of the Gompertz growth law to tumor growth. Their results are based upon curve fitting with actual data. When tumor data involves a wide range of sizes, the Gompertz law gives a better data fit than other growth laws such as the
76
FLOYD
B. HANSON
AND CHARLES
TIER
logistic model [14]. There seems to be no strong physical or biological argument to explain why the Gompertz model fits actual tumor data best. Thus the frequent use of the Gompertz model is empirical, based upon data fitting. The Gompertz model is deterministic in that the tumor size at any time t is given by the function (1.3). In real situations, tumor cells are subject to irregular growth due to random events (see Figure 2). The irregular growth can result in tumor sizes that are different than those predicted by the deterministic models. The irregular growth of tumors appears to be the rule rather than the exception. For example, deterministic tumor models alone cannot explain tumor implantation experiments in which only a certain proportion of the implants take or grow. To account for irregular growth, stochastic models of tumor growth have been introduced, particularly in models of carcinogenesis. Iverson and Arley [4] describe the growth of a transformed cell, a progenitor of a tumor, by a pure linear birth process. In their model, the probability of a birth is constant, which is analogous to a constant specific growth rate and hence a density independent model. Neyman and Scott [lo] used a linear birth-anddeath process to describe tumor growth. Their model used constant birth and death probabilities and hence is also density independent. A stochastic model for the growth of solid tumors was developed by Wette, Katz, and Rodin [18] based upon physical characteristics of the tumor. This leads to a density dependent stochastic process for the mean tumor size. Their process does not seem to be related to the Gompertz model, but has similar qualitative aspects.
G E 2
I -
.Olo Time
100. in days
200.
FIG. 2. Examples of irregular growth,of tumors from experimental data sketched from [14, Figure 1.51. Curve A is for an R39 sarcoma. The other curves are for the growth of primary BMBA-induced mammary carcinomas (for details see [ 141).
STOCHASTIC
TUMOR
77
GROWTH
Smith and Tuckwell [ 131 consider a random Gompertz law in which the parameter b is replaced by a Gaussian white-noise process. The resulting stochastic process is analyzed, and the mean and variance of the logarithm of the population size are computed. The model of Smith and Tuckwell does not allow for tumor extinction, and the exact meaning of the noise in the parameter h is uncertain. Recently, Prajneshu [ 1 I] analyzed a model similar to that of Smith and Tuckwell with heredity effects. Dubin [2] formulates a density dependent birth-and-death process to describe tumor growth subject to immunological response. The density dependence is due to a quadratic factor (nonlinearity) in the transition probability of the death of a tumor cell. The deterministic part of Dubin’s model is similar to the logistic growth law. Dubin uses various approximations to obtain information about the process. Swan [ 161 describes a method for obtaining the exact solution to Dubin’s model in his review of tumor models. We propose a stochastic model of tumor growth which is the diffusion limit of a continuous-time, density dependent branching process. In Section 2, we give a brief derivation of the branching process and its corresponding diffusion limit. We find that the infinitesimal mean and variance of the limiting diffusion process are M(x)
= bxln( K/x),
V(x) =ax.
(1.6)
Here x is a scaled population size, and the infinitesimal mean, which is the deterministic part of the process, is the Gompertz growth law (1.2). The variance arises from the stochasticity of the branching process. We have chosen a specific branching process that has a Gompertz growth law as its deterministic part in order to be consistent with empirical fits to actual data. We have no physical or biological motivation for this branching process. We feel that this model is useful, since it involves the Gompertz growth law, density dependence and allows for tumor extinction. In Section 3, we analyze this model by computing the conditional probabilities of the tumor becoming extinct and of the tumor reaching some upper boundary point, say c. The probability of reaching c will be interpreted as the probability of a tumor becoming detectable or the probability of take of a tumor implant. In addition, we shall compute the expected time for the above conditional events along with the expected doubling time. The results are computed exactly in terms of integrals by solving problems involving the backward Kolmogorov equation, described in Appendix 1. The integrals cannot be evaluated exactly, so asymptotic approximations are found when the deterministic part is large relative to the variance of the branching process. We measure this relation through the nondimensional parameter
78
FLOYD
B. HANSON
AND CHARLES
TIER
R = 2 bK/a B 1. In Appendix 2, we describe the nondimensional parameter R and a method for estimating the size of the parameter a in the variance of the branching process. The approximations are used to obtain useful information about the stochastic model. These approximations have been used in ecological models by Ludwig [8, 91 and by Hanson and Tier [3, 171. When R is not large, the integrals in the exact result are computed numerically. The numerical computations are compared with the asymptotic results when R Z+1 as a check of accuracy. Both the asymptotic and the numerical computations require special mathematical treatment due to the logarithm in the Gompertz growth law (1.6). These approximations are described in Appendix 3. Finally, in Section 4, we apply the results to. two applications in cancer research. First, we compute the probability of take in limiting-dilution assay experiments. The results are compared with a widely used theory. Next, we use the model for the growth phase in carcinogenesis. We compute the ultimate probability of finding a detectable tumor and the expected time from exposure to a carcinogen until the appearance of a detectable tumor.
2.
THE DIFFUSION
APPROXIMATION
We shall now heuristically derive a diffusion approximation to a continuous-time, density dependent branching process for the number of tumor cells, N(t). Here N(t) is a discrete-valued variable. The details of a rigorous proof of the convergence of density dependent branching processes to diffusion processes is given by Lipow [7]. Our interest is not in the rigorous proof of convergence, but in the analysis of the limiting diffusion process. Let the number of tumor cells at time t + h be given by N(t)
iv(t+h)=
2
B,,
(2.1)
t=l
where h is small. Here B, is a random variable equal to the number of cells resulting from tumor cell i, per unit time in (t, t + h). B, = 0 if the cell dies, = 1 if there is no change, and = 2 if division occurs. In density dependent models, B, is a function of the current tumor size N(t) as well as h. We shall assume that the random variables B, are independent and identically distributed with mean and variance
E[B,IN(t)=n]
=m,(n;h),
(2.2) Var[BfIN(t)=n]
=o,(n;h).
STOCHASTIC
TUMOR
GROWTH
79
Using (2.1) with (2.2), we can compute the expected value and variance of thechangeinthetumorsize,AN(t)=N(t+h)-N(t),as
E[ANIN(t)=n]=E =n[m,(n;h)-1]
(2.3)
and
Var[AN(N(t)=n]
=E[[AN-~(~,-~)]‘IN(~)=~]
=‘[/j,iB,-m~)]~] =i i &B,-m,)@,-d] r=i
]=I
= nv,( n; h).
(2.4)
Here we have used the independence of the random variables B,, i.e. Covar[ B,, B,] = 0, to simplify (2.4). To formally obtain a diffusion approximation to the random process defined by (2.1), we must transform the discrete-variable process into a continuous one. A continuous variable is defined by X( t’) = N( Lt’)/L, where L is some large reference tumor size, such as the detectable size of the tumor, and t’ = t/L is a new time variable. The diffusion approximation also requires that small changes in X(t’) occur in small time intervals, i.e., At’=h/L is small and m,(n;h)=l+[g(h/L)+o(l)]h/L is near replacement value. Also the “offspring” variance is assumed to be of the form v,( n; h) = [a + o(l)]h. The diffusion approximation is then obtained as At’- 0. The limiting diffusion process X(t) is described by its infinitesimal mean and variance,
w”)=A~y
0
E[AXjX(t’)=x] At’
’
Var[AXIX(t’)=x] V(x) = ,‘,im At’ ‘-0 Introducing
the coupled
scalings x = n/L
and t’= t/L
(2.5) . in (2.3) and (2.4),
80
FLOYD
then letting At’-
B. HANSON
AND CHARLES
TIER
0, we obtain
M(x) = xdx), v(x) = ux.
(2.6)
We note that the derivation of the infinitesimal mean M(x) is analogous to the derivation of the deterministic model (1.1) except that the scaling from discrete to continuous variables is usually ignored. Finally, if we choose M(x)=bxln(K/x), V(x) = ax,
(2.7)
then we have a stochastic Gompertz model with the stochasticity due to the randomness of the branching process. We note that a in (2.7) could also depend on x. In Section 3, we shall analyze the diffusion approximation and evaluate some important probabilistic quantities for it. 3.
MATHEMATICAL
ANALYSIS
We are interested in studying the diffusion process for the scaled tumor size X(t), which is defined in the interval [0, c]. Here 0 represents tumor extinction and c has different interpretations depending upon the application. For example, c could represent the size at which a tumor becomes detectable or the size at which an implant is said to have taken. Of particular interest are the conditional probabilities of extinction or reaching c, the conditional probability of doubling, and the expected times of these events. To determine these quantities, we need to study problems involving the backward Kolmogorov operator for the process X(t), which is Lu=;V(x)u,,+M(x)u,. Here M(x) and V(x) are the infinitesimal mean and variance, which are derived in Section 2 and given by M(x)
(3.1) respectively,
= bxln( K/x),
V(x) = CIX.
(3.2)
In the above model, tumor extinction is possible, since x = 0 is a singular exit boundary point. The boundary point x = c is a regular boundary point (see [5]). In addition, we assume that c K K, i.e., c is much smaller than K, which can be interpreted for example as the size of a detectable tumor being much smaller than the largest size to which the tumor can grow, K. This appears to be a reasonable assumption for most tumors.
STOCHASTIC
CONDITIONAL
TUMOR
GROWTH
81
PROBABILITIES
Let P,(n) be the conditional probability of a tumor reaching x = c before extinction, given an initial size x. Then, as is shown in Appendix 1, P,(x) satisfies ~V(x)P:‘(x)+M(x)PC~(x)=O, PC(O) = 0,
(3.3)
P,(c)=l.
The solution of (3.3) with M(x) and V(x) given by (3.2) is x
_2hMK/r)+
t
I 1exp1 P<(X)= O( _2hlln(K/t)+r dr’ I Jexp1 a
dt
(3.4)
a
0
If P,(x) is the conditional probability of reaching extinction, given x and that the tumor has not reached x = c, then it satisfies (3.3) but with boundary conditions
PO(O) = 1,
P,(c)=O.
(3.5)
The solution of (3.3) with (3.5) is
(3 6)
In general the integrals in (3.4) and (3.6) cannot be evaluated explicitly. In such cases, we must either evaluate the integrals numerically or use an asymptotic approximation technique to simplify them. We shall describe asymptotic approximations of (3.4) and (3.6) valid when the deterministic part is large relative to the variance, i.e. R = 2 hK/a B 1. These results will then be compared with numerical computations. The details of the asymptotic approximations and the numerical computations are given in Appendix 3. In addition we assume that c< K. We find that for R B 1, (3.4) becomes
P<(X)-
I-exp(-x/x,,) I- exp( - c/x,,)
’
(3.7)
FLOYD
82
B. HANSON
AND CHARLES
TIER
where x,, is a critical tumor size that is the smaller solution of the transcendental equation (3.8) For R za 1, the smaller root of (3.8) is approximated
by
x,, N -In R < 1, K R
(3.9)
The critical tumor size x,, plays an important qualitative role in the approximation of the conditional probabilities. The conditional probability P,.(x) given by (3.7) can be simplified further to give P,(x)
-1 -exp(-
x/x,,),
x,,(x(c.
(3.10)
We feel that (3.10) is useful in applications such as detection of tumors and take probabilities of implants, since for many tumors x,, a CK K. The probability of tumor extinction, P,(x), can be approximated using (3.7), since Pa(x) = 1 - PC(x), which simplifies as in (3.10) to pa(x)-cxp(-x/x,,),
FIG. 3.
The conditional
P~.num, and asymptotically
with c = 1 X 10”.
probability f,(x) evaluated using (3.10) with label P,,,,,,
x,,
(3.11)
numerically from (3.4) with label for the case described in Figure I
STOCHASTIC
TUMOR
83
GROWTH
In Figure 3, we illustrate the conditional probability of exit at x = c using the parameters for the Fortner plasmacytoma described above and in Appendix 2. The probability PC(x) is computed numerically from the exact solution (3.4) and by using the asymptotic approximation (3.7). For this value of R the numerical and asymptotic results are in close agreement. For a smaller value of R the agreement between the numerical and asymptotic results should not be as good. The conditional probability of doubling, P,(x), is the probability that a tumor with initial size x doubles before extinction. This can be found by evaluating PC(c/2) using (3.7) for R z~ 1 and then replacing c by 2x to get
P,(x) -
1
l+exp(-x/x,,).
(3.12)
Here x represents the initial tumor size, and x,, is the critical tumor size. Thus, for large R, the conditional probability of doubling is between f and 1 (see Figure 4). We note that P,(O) = 4, because for small x both extinction and doubling have nearly equal probabilities, as in a simple random walk.
FIG. 4. The probability P D.““rnTand asymptotically ure I.
of doubling P,(x), evaluated numerically using (3.7) with label using (3.12) with label PD.,,,, for the case described in Fig-
84
FLOYD
B. HANSON
AND CHARLES
TIER
The larger the value of R, the more rapidly this approaches 1 as the initial size increases. For sufficiently large x, the tumor is almost certain to double before extinction. This is very much like the deterministic models, since the effect of the randomness is reduced for larger tumor sizes. EXPECTED
TIMES
Another quantity useful in describing tumor growth is the expected, or mean, time for the conditional events described above. These are the expected time for the tumor to reach c before extinction and the expected doubling time given that the tumor does not go extinct. The expected time of the conditional events is computed from a conditional process which includes only realizations of the original process satisfying the conditioning. Since the conditional process is defined by excluding the realizations not satisfying the conditioning, the total probability of the conditional event is less than one. This is assuming that the discarded realizations have a nonzero probability. Thus, the conditional process must be normalized to make the conditional event a proper random event. Once a proper random event is defined, the expected time of the conditional process is computed as the mean time of the realizations. As described in Appendix 1, these expected times can be computed using the backward Kolmogorov operator (3.1). Let E/f,.] be the expected time for the tumor to reach c, excluding tumors that become extinct. It is given by
E[t,.l = ~,(x)/~,.(x), where M,.(x) satisfies (see Appendix
1)
f~(x)M:‘(x)+M(x)M~(x)= M.(O) = 0,
(3.13)
M,.(c)=O.
-P,(x),
(3.14) (3.15)
Here P,(x) is the conditional probability of reaching c before extinction and is given by (3.4). Using the solution of (3.14) and (3.15) with M(x) and V(x) given by (3.2) then the expected time for conditional exit at x = c can be found from (3.13) and is given by
(3.16)
STOCHASTIC
TUMOR
GROWTH
85
where for convenience
=
R[xln(K/x)+xl K
(3.17)
.
The expected doubling time E[t,] can be found using (3.16) with the initial value x = c/2. Again we shall consider the doubling from x to 2x, so we replace c by 2x in (3.16) to get
(3.18) This result should be compared with the doubling time for the deterministic Gompertz growth law given by (1.5) with N(0) the initial size instead of x. The expressions for the expected times are complicated and cannot be evaluated explicitly. Asymptotic expansions of the multiple integrals (3.16) and (3.18) are also difficult, due to the presence of the growing and decaying exponential terms. Instead we shall find deterministic approximations for E[t,.] and E[t,] using (3.13) and (3.14). When R B 1 and x is not near 0, the second-derivative term in (3.14) is small and PC(x ) - 1 from (3.10). Thus, (3.14) can be written approximately as M(x)kqx)-
-1,
(3.19)
with
M,.(c)=O.
(3.20)
Here we have retained the boundary condition at x = c from (3.15), since our approximation is only valid away from 0. The solution of (3.19) with (3.20) is the deterministic passage time T, given by (1.4), so using (3.13), E[t,.]-M,.(x)-T,.. Similarly,
the expected doubling
time can be approximated
E[t,l -To. Here T, is the deterministic
doubling
time given by (1.5).
(3.21) by (3.22)
86
FLOYD
R. HANSON
AND CHARLES
TIER
Since tumor data are quite variable, the conditions for the above approximations are not always satisfied. For other cases, the multiple integrals (3.16) and (3.18) must be evaluated numerically. The details of the numerical computations are given in Appendix 3. In Figure 5, we graph E[ t,.] and T,. for the Fortner plasmacytoma described above. We notice that E[ t,.] and T,. are close for x away from 0 as expected from (3.21). It is important to note that for small x, E[t,.] approaches a finite value while T is unbounded. The graphs of E[t,] and TD are displayed in Figure 6 for the Fortner plasma-
0
012
014
0:6
018
X, CEL LS/lO” FIG. 5. The deterministic first-passage time q given by (I .4), and the expected time for conditional exit at X= c, E[t,], given by (3.16). for the case described in Figure I with (‘= I x IO’“.
STOCHASTIC
TUMOR
GROWTH
FIG. 6. The deterministic doubling time T, given by (1.5). and the expected time E[ tn] given by (3. IS), for the case described in Figure I,
87
doubling
cytoma. Again, they are in close agreement for x away from 0. Near x = 0, E[t,] has a finite slope, while the slope of To is unbounded. If instead c> K/2, then E[ tn] would be bounded for all values of x, while To is unbounded as x --t K/2. The numerical calculations of Figures 5 and 6 show that the deterministic approximations (3.21) and (3.22) are suitable for R large and x away from 0. This is also true for the parameter values in Appendix 2. 4.
APPLICATIONS
We shall now apply the results in Section 3 to two specific applications in cancer research. First we shall consider limiting-dilution assay experiments and the probability of a successful implant take. Next we shall consider the problem in carcinogenesis of finding the probability and expected time of tumor occurrence after exposure to a carcinogen. LIMITING-
DILUTION
ASSA Y
A limiting-dilution assay is an assay method used to determine the growth and colony growing characteristics of tumors (see [ 141). It can also be used to determine the surviving fraction of treated tumors. The method consists of
88
FLOYD
B. HANSON
AND
CHARLES
TIER
either plating or transplanting into an experimental environment a multicell suspension (implant) derived from a particular tumor. The experiment is repeated using implants of the same size, and the total number of colonies is observed and used as a measure of the proportion or probability of implant take. An implant is said to have taken if it reaches some prescribed size c within a certain time interval. The probability of take should depend upon the growth characteristics of the original tumor. To determine the surviving fraction of a treated tumor, the experiment is performed on both the treated and the untreated tumors. For each tumor, the number of cells in an implant required for a 50% take probability is determined. The ratio of the number of cells required for a 50% take probability in the untreated to that in the treated tumor is the surviving fraction. The 50% take probability is chosen arbitrarily as a standard. The results are then extrapolated to the treatment of tumors in patients. If one uses only the deterministic Gompertz growth law (1.2) or (1.3) as a model for the limiting-dilution assay experiments, it predicts that all implants, no matter how small, will be successful, though the time of take will vary. Thus, the probability of take equals 1. To account for the experimental results, it has been proposed (see [ 141 for a review) that each tumor contains only a small fraction of cells that can be the progenitor of a colony (i.e. clonogenic cells). Thus for an implant to be successful, it must contain at least one clonogenic cell that proliferates. To account for the stochastic effect in the experiments, it is assumed that the probability of obtaining an implant with no clonogenic cells is exponentially distributed, as e-“‘, where m is the average number of clonogenic cells per implant. Thus, the probability of obtaining an implant with one or more clonogenic cells is I- epm, and hence the probability of take is Pr {implant take} = I- epm =
1 -
e-/n’.
(4.1)
In the second equality, the parameter N is the total number of cells in the implant, and f is the fraction of clonogenic cells. The stochasticity only enters in the probability of selecting clonogenic cells in the implant. We should note that it is in general not possible to distinguish between clonogenic and nonclonogenic cells. In the stochastic Gompertz growth law described in Section 2, the stochasticity is due to the branching process, with each cell division giving rise to a random number of proliferating descendents. The stochastic effects are most evident when the implant sizes are small, as in these experiments. The probability of take of an implant is the conditional probability of reaching some prescribed upper limit c, i.e. P,(x). The conditional probabil-
STOCHASTIC
TUMOR
GROWTH
89
ity PC(x) is given by (3.4), or asymptotically and c >>x,, , the probability of take is
for R B 1 by (3.7). When R B 1
P,.(x)-l-exp(-x/x,,). Here x is the initial size of the implant
(4.2)
and x,, is the smaller root of (3.8): (4.3)
The value of x,, depends upon the growth characteristics of the tumor, including b, K, and the stochasticity a. The values of b and K are estimated from empirical data fits using the Gompertz growth law and represent average values over a range of tumor growth. The additional parameter a measures the variance of the birth schedule in the branching process (see Appendix 2). We note that P<(x) represents the ultimate probability of the tumor size reaching x = c before extinction. Finally, we compute the surviving fraction for a treated tumor using the limiting-dilution assay. We shall assume that the growth characteristics of the treated and untreated tumors are different. Thus, the parameters in the Gompertz growth laws of these tumors will be different. For the untreated tumor we shall assume the probability of take is given by (4.2) while for treated tumor it is
P,(x) = 1-exp(
- x/x’,,)
(4.4)
where x& satisfies (4.3) with the treated-tumor growth parameters. We compute the surviving fraction by taking the ratio of the numbers of cells from the untreated and the treated tumor required for a 50% take probability. This gives surviving fraction = x,, /x&.
(4.5)
Thus, the surviving fraction reduces to the ratio of the two critical tumor sizes when R za 1. These critical sizes are functions of the growth parameters for each of the tumors and also of the stochasticity. The formula (4.5) represents the ultimate surviving fraction. If the tumor can eventually repair the damage caused by treatment and regain the growth characteristics of the untreated tumor, then the surviving fraction will ultimately be 1. CA RCINOGENESIS
The development of a tumor is believed to be a two-phase process. The first is the transformation phase, in which normal cells are transformed into progenitors of a tumor. The transformation may result from the contact of
90
FLOYD
B. HANSON
AND CHARLES
TIER
normal cells with carcinogens and is generally believed to occur in a series of stages. In the second or growth phase, the transformed cell grows into a detectable tumor. We shall compute the ultimate probability of finding a detectable tumor given the exposure of normal cells to a carcinogen, and the expected time from exposure to a carcinogen until the appearance of a detectable tumor. Stochastic models have been proposed for both the phases of carcinogenesis (see [ 191). We shall use the stochastic Gompertz model outlined in Section 2 for the growth phase. Let q = Pr {normal cell is ultimately
transformed}
,
(4.6)
which may result from a multistage transformation process. If p is the ultimate probability of a normal cell becoming a detectable tumor, then P = #C(l)
-4[1-w(-l/xcr)l,
R>l.
(4.7)
Here PC(1) is the conditional probability the tumor reaches x = c given that the initial number of cells was 1 (the single transformed cell). It is given by (3.4) or approximately by (3.10). The value of c is the smallest detectable tumor size. Thus, in a collection of J identical, normal cells, the ultimate probability of finding at least one cell that will give rise to detectable tumor is found by assuming that the probability of finding k tumors, Pk, is binomial (see [ 191): J! pk= k!(J-
k)!
pk(l-py.
For J large and p small, this can be approximated
by the Poisson distribution
Pk - e -Jp( Jp)k,‘k!. Thus, the probability
of at least one tumor is 1 - PO, or approximately Pr{ at least 1 tumor} - 1 - eCJP.
(4.8)
The tumor time t = t, + t, is the total time from exposure to a carcinogen until the appearance of a detectable tumor. Here t, is the transformation time and r2 is the growth time. Since the density function for t is a convolution of the densities of t, and t2, then
E[t]=E[t,]+E[t,].
(4.9)
91
STOCHASTIC TUMOR GROWTH
Here the random time variables t, t,, and t, all refer to times of conditional events. Thus E[t] represents the average time from exposure to a carcinogen until the appearance of a detectable tumor given that the normal cell eventually gives rise to a tumor. The expected time E[ t2] is the expected time for the growth phase given that the tumor becomes detectable, and is the same as E[t,.] given by (3.13) and (3.16).
APPENDIX
1.
BACKGROUND
PROBABILITY
In this appendix, we give intuitive derivations of the equations for the conditional probabilities and the expected times described in Section 3. Let q(x, y, t) be the transition probability density function for the diffusion process X(t) = y given that X(0) =x. Then q satisfies the forward equation
4r+J,.M=o
OC.Y
(Al.l)
where
J[ql(xTY,t) = -Iv(Y)ql v+WY)q is the probability equation
flux. The density
function
q,=L[ql=tV(x)q,,+M(x)q,, The initial condition
for the backward
(Al .2)
q also satisfies the backward
ocxcc.
and forward equations
(Al .3) is
limq(x,y,t)=6(y-xx).
(Al .4)
r-0
The coefficients M(x) and V(x) are the infinitesimal mean and variance, respectively; see (2.11). The boundary conditions for (Al.l) and (A1.3) depend upon the particular process (see [5]). Let r,(x) be the exit time of X(t) from (0, c), defined as te( x) = inf {t 1X(t) Then the distribution given by
not in (0, c)] X(0) = x in (0, c)}.
of t, is the probability
(Al .5)
that X(t) is not in (0, c) and is
(Al .6)
92
FLOYD
with a density function,
B. HANSON
the time derivative
AND
of the distribution
CHARLES
TIER
(A1.6), given
by (Al .7) Using (Al. 1) and (Al .2), the density function ~(x,t)=
limJ_Y - c
can be written
,ti_moJ=u~(x,t)-po(x,f)
(Al .8)
and hence has two components: pa(x, t) for exits at x = 0, and p&x, t) for exits at x = c. Since q satisfies the backward equation (Al .3), then so do p( x, t), pO( x, t), and p,,( x, t), which are derived as
p,..,=J[q,l(x,c,t)=J[~[qll Here we have interchanged
differentiation
=L[J[qll =L[JJcl. and integration.
Similarly,
P,=L[Pl, PO,* =L[Pol.
The ultimate probability
of exit from either 0 or c is
p(x) = jg"P(X, t) dt,
(Al .9)
which also satisfies the backward equation (A 1.3), i.e. LP( x) = 0. The proper boundary conditions are P(0) = 1 and P(c) = 1, since a process with initial points 0 or c is outside the interval. These absorbing boundary conditions imply that the ultimate probability of exit is certain, i.e. P(x) = 1. Similarly, the ultimate probability of conditional exit at x = 0 is (A1.lO) and at x = c is (Al.ll) Each satisfies the backward
equation
L[P,l(x) = L[P,lb) =o,
(A1.12)
STOCHASTIC
TUMOR
GROWTH
93
which is Equation (3.3). The boundary conditions for P,(x) are [see (3.5)] P,(O) = 1 and P,(c) = 0, since a process that starts at x = 0 is outside (0, c), while if it starts at x = c it can never exit at x = 0. Using a similar argument, we obtain the boundary conditions for PC(x) which are given in (3.3). The expected time of exit at either x = 0 or c is the first moment of p( x, t) and is given by (A1.13) which after applying the backward P(x) = 1, we find to satisfy
operator,
integrating
LIT]=~rtLIP1dr=jo”~~,dr=tp]~-IXPdt=
by parts, and using
-P(x)=
-1.
0
(A1.14) The boundary conditions are T(0) = T(c) = 0, since it takes no time for the process to exit if its initial point is either x = 0 or c. The expected time of the conditional exits are more complicated, since the probability of conditional exit is in general less than one, and hence conditional exits are not proper random events. To define the conditional exit at x = c to be a proper random event, we normalize its density to get p,.( x, t)/P<( x). This now integrates to 1. The expected time for conditional exit at x = c is then defined by
fk(*. E[ t,.] =
t) dt
p,.(x)
(Al .15)
Unfortunately, this does not satisfy a simple equation. Instead, it is easy to show using the same steps as on T(x) that the conditional moment M,.(x) defined by
w(x)
=pPc(x,t)dt=p,.(x)E[t,.]
(~1.16)
satisfies
Hence, the expected found using (3.14).
time of conditional
exit at x = c satisfies (3.13) and is
94 APPENDIX
FLOYD
2.
ESTIMATE
B. HANSON
OF THE DIFFUSION
AND CHARLES
PARAMETER
TIER
u
We could estimate the parameter a from tumor measurements, using the variance Var[ X] of the process. Var[ X] can be exactly calculated in principle from the density for X(t), but the density must be obtained from the forward Kolmogorov equation with a great deal of effort. Hence we employ an Omstein-Uhlenbeck approximation of the density about the maximum, deterministic tumor size K, by letting y = K( I- z) with z small and using it in the forward equation (Al. 1):
(A2.1)
We note that the dimensionless “signal-to-noise ratio,” R = 2 hK/a, naturally arises when y or x is scaled by K, and t is scaled by 1/b. Setting qt = 0 in (A2. l), we can find a quasi-stationary density (see [5]), which represents the density of those tumors which are not yet extinct. The normalized, quasi -stationary density is approximately q^.(R/27r)“2exp[-Rz2/2],
(A2.2)
with mean K and variance K2/R ‘1’ for the variable y. Hence the coefficient of variation, or relative standard deviation, is CV-
l/R’/‘.
(A2.3)
We have also computed formulas for the CV from Ludwig’s ray method [8, 91 by expanding the exponent of the density about the deterministic trajectory. These computations give the CV as a function of y and agree with (A2.3) when y = K. However, we will report these details in another investigation. Simpson-Herren and Lloyd [ 121 in their Table 4 give average, regression adjusted weights and their standard deviations for Fortner plasmacytoma 1 harvested from hamster hosts. Coefficients of variation computed from their standard deviations and adjusted weights range from 0.224 to 0.930 with a mean of 0.523. The measurements closest to K have the smaller CVs and have an extrapolated value close to 0.300. Given b = O.l07/day and K = 29,300 mg from their data and with CV= 0.300, we find that a rough estimate of a is 600 mg/day = 6 X lo9 cells per day. We note that although the “noise” a is large in an absolute sense, its relative magnitude is small when properly measured in terms of the dimensionless group that actually appears in the theory, i.e. the reciprocal of the “signal-to-noise ratio,” 1/R = a /2 bK = 0.0957. Therefore a small-noise approximation is more appropriate, and we consider this in the next appendix.
STOCHASTIC
TUMOR
APPENDIX
3.
95
GROWTH
ASYMPTOTIC AND NUMERICAL APPROXIMATIONS
In this appendix we derive special asymptotic and numerical approximations for the conditional probabilities and the expected times described in Section 3. The stochastic model with the Gompertz growth law involves a logarithm, which leads to nonanalytic behavior. The standard methods used in asymptotic and numerical analyses are not applicable directly, and hence we use a modified analysis based upon a transformation to remove the logarithm. APPROXIMATION
OF THE CONDITIONAL
PROBABILITIES
We seek an asymptotic approximation of P<(x) given by (3.4) when R = 2 bK/u B 1. Rewriting (3.4) using (3.17), we get
(A3.1)
where G(x)
=
Rx[ln(K/x)+ll K
To eliminate
the logarithm,
we let y = ln( K/x),
P,(x) =
J-“,
s
(A3.2)
.
so that (A3.1) becomes
(A3.3)
)
eHcs)ds
Y
H(y)=-y-G(KeeY)=
(A3.4)
-y-Rep’(y+l),
and Y = ln( K/c). In the spirit of Laplace’s method for the asymptotic expansion of exponential-type integrals, we seek a maximum of the integrands in (A3.3). The critical point y,, of H(y) satisfies H’(.K,) = - 1+ Ry,,evA maximum
(A3.5)
IL> = 0.
occurs at the largest root in (A3.5), which is y,,=lnR+lny,,
-lnR+ln(lnR),
p-
+co.
(~3.6)
96
FLOYD
B. HANSON
AND CHARLES
TIER
Converting (A3.5) and (A3.6) back to the variable x, we obtain (3.8) and (3.9) respectively, for x,,. Testing the second derivative of H(y) at y,,, we find that H”(y,.,) = - 1+ l/y,, - - 1 as R - + 00. Hence, although y,, is the location of a wide maximum, it is not sharp enough to allow the use of Laplace’s method. To obtain the final result, we use the transformation z = ln(x,,/x) on (A3.3) to make it easier to examine its behavior near the maximum of the exponent. The transformation leads to
(A3.7)
where
h(z)=
-z-ePi[l+F],
(A3.8)
and Z = ln( x0/c). Here L = ln( K/x,,), large, we may approximate h(z) by h(z)This approximation
and L ZJ1 when R B 1. With
-z-eeZ,
LBl.
L
(A3.9)
will be valid when e-‘(z+l)/Lal.
(A3.10)
However, when the approximation (A3.9) is applied to (A3.7) the result is uniformly valid for 0
P,.(x)-
1-exp(-e-‘)
(A3.11)
l-exp(-e-Z)
Transforming back to x, we obtain (3.7) or when c/x,,< 1 the further simplified form (3. IO). The results for the probability of doubling, PD( x), are analogous and given by (3.12). QUASI
-DETERMINISTICAPPROXIMATION
OF THE TIME TO EXIT
FROM
(0, c)
A deterministic approximation of E[ t,] given in (Al. 13) by T,, is less suitable because E[ t,] ---r+ 0 as x + +O, while Tc does not behave this way. However, we have found that a “quasi-deterministic” approximation,
-qt,l -
p,.T,,
(A3.12)
STOCHASTIC
FIG. 7. deterministic
TUMOR
97
GROWTH
The expected
exit time E[t,]
given by
(Al. 13) compared
with the quasi-
7; P, given by (A3.12).
works reasonably well, although it consistently underestimates the exit time. E[ r,] is compared with P,.T,. in Figure 7. In summary, we have found good asymptotic approximations for the probabilities of exit (detection) and doubling; and we have found deterministic approximations (3.21) and (3.22) and a quasi-deterministic approximation (A3.12) for the exit or doubling times where more accurate asymptotic approximations are very difficult to compute. NUMERICAl.
TREATMENT
OF LOGARITHMIC
INTEGRALS
A general rule of numerical integration schemes is that the order of the local error depends upon some power of the mesh size h, and upon a corresponding derivative of the functionf( x), such as O( h”+ ‘f(“)(x)). Here n - 1 is the highest polynomial that can be integrated exactly by the numerical scheme. However, when the function has a nearby discontinuity, the actual local error is not of the general form. Instead, n = 0 must be used in the error estimate in almost all cases. Since our problem involves In x, which is not well behaved near x = 0, a special treatment is necessary to improve the local error to better than O(h). We shall transform variables to remove
98
FLOYD
B. HANSON
AND CHARLES
TIER
the logarithm and obtain better-behaved integrals. This will allow us to use the Gaussian integration rules with their high degree of polynomial precision. We shall first analyze P<(x). We let y = ln(x/t) in both integrals in (A3.1), which we rewrite as
xF,(x) 4(c) ’
P,.(x) = ~
(A3.13)
where
F,(x) =/,, e-‘fo(Y;x>dy, fa(y;x)=exp[-G(xe-‘)I. The integrals in (A3.13) now have smooth integrands but are singular in their domain. They can be integrated without problems using a Gauss-Laguerre integration rule [i.e. Gaussian on (0, co) with weight factor e -“I. Some results are shown in Figure 3 for our sample parameters. The probability of doubling, PD( x), can also be evaluated using the same formula. Some results for PO are shown in Figure 4. A similar technique can be used on the integral formula for E[t,.] in (3.16) except that a “catastrophic cancellation” can occur between the two integrals. Even for moderate values of the “signal-to-noise ratio” R, the two integral terms in (3.16) can be almost the same size, and their difference leads to a great loss of significant digits. Fortunately, if we integrate by parts to change the order of integration in (3.16) and then rearrange the result so that positive integrals are added rather than subtracted, we obtain the useful form EL1
,] L
=
2[Q(c)-
Q(x)]
AQ(c)Q(x) 2
+AQ(c)
x Q2b)eGrV'dy
Jo
c.Q(~)~~(~)
Y
[QWQb)l
sx
dy,
tA3.14j
Y
where Q(x)=~xepG(‘)dt
and
Q(c)-Q(x)=/Ce-G(‘Ldt. X
The integral over (0, x) is treated like those in (A3.13), using the logarithm transformation and Gauss-Laguerre rules. The integrals over (x, c) or (y, c) as in Q(c) - Q(y) are still affected by the singularity of In x for small x, even though x is excluded from the domain. Hence on (x, c) [and similarly on (y, c)], we use the transformation
,=1_21n(c/~) ln(c/x)
’
99
FLOYD B. HANSON AND CHARLES TIER
and then use a Gauss-Legendre integration rule [i.e. a Gaussian (- 1,l) with weight 11. Upon transforming the integrals, we obtain
xF,(x; E[t,.] =In(
c)&(x)/F,(x)+
cln( t)F,(x;
rule on
c)/4
f) M,(c) (A3.15)
where
F,(x; c) = j-;,f,(s; x; c) h, f,(s; x; c) = exp{g(s;x; g(s;x;c)= F,(x) f2(t.
x)
-In =im =
c>-G( cexp[g(s;x; c)])},
I-s
(f 12,
e-‘f2( t; x) dt,
FoZ(xe-“*)
x>’ F,(x; c) =j--lh(s; x; c> ds, 2f,(t/2
f3(s;x; c) = Cl-s)exp{g(s; x; c)+G(cedg(s; x; c)l)}
XF,(cexp[g(s;x; c)l)F,(cexp[g(s;x; c)l; c). The integrandsf, are fairly well behaved, so that the numerical integration is straightforward. For the multiple integrals, the inner integrals of functions must be evaluated at each of the nodes of the outer integrals. Some results for E[ t,.] are shown in Figure 5. The numerical computation of the expected doubling time E[ tD] proceeds directly from (A3.15) with c set equal to 2x, which gives
E[ r,] = (xln2)
F,(x;2x)F2(x)/F,(x)+(ln2)F,(x;2x)/2 a&,(2x) (~3.16)
Certain minor simplifications can be made [e.g., g(s; x; 2x) = (ln2)(1- s)/2 is independent of x], but the primary reduction was made in deriving (A3.15). In Figure 6 are shown some results for E[t,]. The numerical decomposition of E[t,] is similar to that of E[ t,.] in (A3.15), and we omit the formula here. Sample results are presented in
FLOYD
100
B. HANSON
AND CHARLES
TIER
Figure 7. The steeper slope near x = 0 than near x = c indicates that the expected time for a tumor to go extinct (remission) is much longer than the expected time to reach x = c (detection of tumor or death of host). REFERENCES I
A. L. Burton, Rate 30:157-176 (1966).
2
N. Dubin,
growth
of solid
tumors
as a problem
A Stochnstic Model for Immunological
and Approximations.
Lect. Notes in Biomath..
of diffusion,
Growth
Feedback in Carcmogenesis:
Ana!vsrs
Vol. 9, Springer,
3
F. B. Hanson and C. Tier, An asymptotic singular diffusion arising in population (1981).
4
S. Iversen and N. Arley, On the mechanism Microbial. Stand. 27:773-803 (1950).
5
S. Karlin and H. M. Taylor, York, 1981.
6
A. K. Laird, The dynamics
7
C. Lipow, Limiting diffusions for population Appl. Probab. 14: 14-24 (1977).
8 9
D. Ludwig,
Stochastic Population Theories, Springer,
D. Ludwig, 17:605-640
Persistence (1975).
IO
J. Neyman
and E. Scott, Statistical
Prajneshu,
12
(1979). L. Simpson-Herren experimental
13
15 16 17 I8 I9
of experimental
of tumor growth,
of dynamical
Gompertz
and
solution to the first passage time problem for biology, SIAM J. Appl. Math. 53:89-l 17 Acra. Path.
Br. J. Cancer 28:490-502 size dependent
New York,
processes,
aspects
of carcinogenesis.
of the problem
model with heredity Kinetics
J.
1914.
perturbations,
Sratisiics und Probability,
New
(1964).
branching
systems under random
H. H. Lloyd,
tumor systems,
carcinogenesis,
SIAM
Rev.
in Fi/rh
Univ. of Calif. Press,
effect, J. Math. Biol. 8: l89- I96
parameters
Cancer Chemo. Reps. 54: 143-174
and
growth
curves
for
(1970).
C. E. Smith and H. C. Tuckwell, Some stochastic growth processes. Problems in Bio/ogv - Victoria Conference, Lecture Notes in Biomath., New York,
14
A stochastic
1976.
A Second Course in Srochastrc Processes, Academic,
Berkeley Symposium on Mathematical Berkeley, 1967, pp. 745-776. I1
New York,
in Muthemarical Vol. 2. Springer,
1974.
G. G. Steel, Growth Kinetics of Tumors, Clarendon. Oxford, 1977. P. W. Sullivan and S. E. Salmon, Kinetics of tumor growth and regression multiple myeloma, J. C/in. Invest. 51: 1697-1708 (1972).
in IgG
G. W. Swan, Some current mathematical topics in cancer research, Univ. Microfilms Inter.. Ann Arbor, 1977. C. Tier and F. B. Hanson, Persistence in density dependent stochastic populations, Mafh. Biosci. 53:89-II7 (1981). R. Wette, I. N. Katz, and E. Y. Rodin, Stochastic processes for solid tumor kinetics, Math. Biosci. 19:231-255 (1974). A. Whittemore and J. B. Keller, Quantitative theories of carcinogenesis, SIAM Rev. 2O:l-30 (1978).