Might stochasticity and sampling variation be a possible explanation for variation in clonal population survival curves

mechanisimsof ageing and development ELSEVIER SCIENCE IRELAND

Mechanisms of Ageing and Development 73 (1994) 223-248

Might stochasticity and sampling variation be a possible explanation for variation in clonal population survival curves Matthew Witten University of Texas System, Centerfor High Performance Computing, Baleones Research Center, 1.154 CMS, 10100Burnet Road, Austin, TX 78753-4497, USA (Received 16 July 1990; revision received 11 July 1992; accepted 13 August 1993)

Abstract Biogerontological survival analysis attempts to understand, through the use of mathematical and computer models, how biological and environmental processes affect the dynamics of survival. The survival model parameters are assumed to reflect an average or mean response to some intervention. Further, these parameters are usually assumed to be constant over the time course of the experiment and across the elements of the experimental cohort. In this paper, we introduce stochastic (random) features to the survival curve parameters and we observe how this might affect our interpretation of the biology; as reflected in the estimates of the model parameters. In particular, we provide a possible explanation for variation in parameter estimates within sample populations drawn from a population of genetic clones.

Key words." Survival; Stochastic; Random; Diffusion process; First passage time; Mortality rate.

I. Introduction T h e analysis o f survival, in biological p o p u l a t i o n s [1-9,15,19], r o u t i n e l y involves fitting survival d a t a to various families o f curves w h o s e f u n c t i o n a l f o r m S(a) is 0047-6374/94/$00.00 © 1994 Elsevier Science Ireland Ltd. All rights reserved. SSDI 0047-6374(94) 01393-M

Matthew Witten / Mech. Ageing Der 73 (1994) 223-24~;

224

obtained by solving the general survival differential equation dS(a) da

-

A(a)S(a)

(1.1)

S(O) = I

where S ( a ) is the probability of surviving until at least age a (the fraction surviving until at least age a), and A(a) is the mortality rate (the instantaneous per capita rate of dying in the interval (a, a + Aa) given that you have survived until age a). For a more exhaustive discussion and details see Finch et al. [7], Lawless [10], and Witten [11-15] and the references contained therein. Different mortality rate functions A(a) give rise to different traditional families of survival curves. The mortality rate A(a) given by A(a) = constant = h(~

(1.2)

gives rise to the classical exponential survival curve. For a mortality rate A(a), given by A(a) = h o e ~'",

(1.3)

we obtain the two parameter Gompertzian survival curve [13]. A detailed discussion of the theory may be found in Witten [11-15], the references therein contained, and in Gompertz [16-18]. A detailed discussion of the biological application, particularly as applied to the Gompertzian survival model, may be found in Finch et al [7] (and the numerous references therein contained). Experimental survival data are routinely fit to such analytical survival models, in an effort to discern how survival/longevity might be related to the biodynamics of the organism [6,12-15,19-21]. Changes in the model parameters are often considered to be related to changes in the biology of the organism as reflected by changes in a particular experimental protocol [7,8] (Table 1). Fig. 1 of this paper illustrates sample survival curves from Johnson [8]. From the dynamics of popula-

Table 1 Table of error estimates for

lh)

Sample size (n)

Mean estimate for h. (day)- i

Standard deviation of estimates for h 0 (day)

8 16 32 64 128 256 512 1024

0.00944383 (I.I)0473795 0.00334812 0.00264936 0.00229755 0.00216367 0.0019381)7 0.00202570

0.0153899 0.0058977 0.00332969 0.00195955 0.00123923 0.000828517 (I.000534468 0.000384621

225

Matthew Witten / Mech. Ageing Deu. 73 (1994) 223-248

I00 • ee ~0

80

ee

o o

8a

• A 0

eeeO

4O o_ 20

wild type

fer-15 (b26) age-I (hx 542) 20

0

I0

20

,

30 fer-15(b26)

I

40

i) a ;

50

,I

60

age-I + age-l(hx542)

oS

fer-15(b26) i

J

~

.=

I0

l

.o0~ 0 & CL •

,, 40

00

0

• aa~.~

g 8O =~ 60 20

30

~ I•

o ~.

,

10 I00

eO

~° ,

t>

a

g

6O

o ti I

20

30 0

I0 20 Age (days)

30

40

50

60

Fig. 1. Sample survival curves from Johnson [8]. Survival curves of wild-type age-1 mutant and heterozygous strains of C. elegans. In each case, four identical populations of hermaphrodites were separated at three days of age and were monitored in mass cultures of 50 worms each until all worms had died. Survival analyses were performed except that survivors were counted every 24 hours (A) Survival curves for N2 (Bristol) hermaphrodites; m e a n life span of the four cultures were 20.5, 23.0, 20.9, and 21.0 days; culture 2 was found to be significantly longer lived than the other 3 cultures (p < 0.05) and was not used in the further analyses described in Table 1 of Johnson [8]. (B) Survival curves for DH26 [fer-15(b26)] hermaphrodites obtained from M. Klass; m e a n lifespans were 23.8, 20.7, 23.2, and 23.8 days and were not significantly different from each other. (C) Survival curves for TJ411 hermaphrodites; this strain carries fer-15(b26) age-l(hx542); m e a n lifespans were 34.3, 38.0, 37.6, and

36.6 days and were not significantlydifferent from each other. (D) Survivalcurves for F1 hermaphrodites and from the cross TJ411 and DH26 males; mean lifespans were 15.0, 18.7, 19.1, and 17.8 days and were not significantly different from each other.

tion survival, we wish to make inferences about the underlying biology of aging of the individual in the population [22]. In this paper, we present an alternative way of looking at variation in survival parameter estimates. We consider the possibility that one or more of the survival model parameters is, in reality, a stochastic or random variable and not constant over the life of the individual or the term of the experiment. The assumption of a non-constant mortality rate p a r a m e t e r has been recently shown to be biologically plausable [23,24]; though this is disputed in Riggs and Millecchia [40]. We will examine the theoretical consequences of such an assumption and we will show how to derive various important survival and mortality rate parameters and functions. We will then close with a discussion of the biological interpretation of the non-constant mortality rate formalism.

226

Matthew Witten /Mech. Ageing Det'. 73 (1994)223-248

2. Preliminary modeling For the purposes of our discussion, we will consider the Gompertzian survival model. We do this for two reasons. First, because the Gompertz model is a commonly used [7,8,11-15] survival distribution, and the second because it contains the exponential survival distribution as a limiting case. This analysis can be performed using any other distribution (Weibull, Gompertz-Makeham, etc.). We make no case for the validity of the Gompertz model as an adequate model for describing survival over the lifespan. This is discussed in numerous other papers [11-15,23,24,40]. From Eq. 1.3, we observe that the classical, deterministic, Gompertzian mortality rate A(a) may be re-expressed as the solution of the following differential equation dh(a)

= pA

(2.1) A(O) = h~ F o r t h e case where p = 0, we observe that equation (2.1a) reduces to the following equation

dA(a) d~ - 0

(2.2)

which may be solved to yield Eq. 1.2. The constant mortality rate equation (1.2) gives rise to the exponential survival distribution

S(a) = e -t....

(2.3)

For the case p ~ 0, we obtain A(a) as given in Eq. 1.3. Replacing Eq. 1.3 into Eq. 1.1 gives rise to the Gompertzian survival distribution

S(a)=exp{~(1-e"")} See Fig. 6 of Ref. 12 for an illustration of the dynamics of of the parameters h 0 and p (see also Ref. 25 by Finch).

(2.4)

S(a) for various values

3. Stochastic perturbation in the mortality rate model Models of the type specified by Eq. 11 (Gompertz, Weibull, Makeham, etc.) are

models in the mean. That is, they assume a homogeneous (identical) population of indiLqduals. This assumption is used to develop the following mathematical construction. Let S(a) and A(a) be as defined in the introduction. Let us assume that our population has an initial population of N ( 0 ) + N O individuals. From standard mass action arguments, it follows that the number of individuals at age a + Aa (assuming Aa --* 0 and N(a) --* ~c; thereby avoiding the problem of having fraction-

Matthew Witten/ Mech. Ageing Dec' 73 (1994) 223-248

227

al individuals in the population) is given by N ( a + Aa) = N ( a ) - h ( a ) A a N ( a )

(assuming no immigration, emigration, reproduction, or fusion/fission). Dividing both sides of the equation through by No, we obtain the associated survival mass action equation S(a + ha) = S(a) - h(a)haS(a)

Observe that this construction yields equation (1.1), as Aa ~ 0. However, this construction also assumes that all of the individuals in the population are equivalent to an average individual and that this population o f average individuals obeys a simple bulk mass transport (survival mass) equation where the rate behavior is replaced by a bulk average behavior over the aggregate clonal population. Clearly, in any population experiment, this type of behavior will not necessarily be the case. The probability of death, at a specific age a will vary, about some mean value, from individual to individual and, perhaps, over time [22-24]. We now consider one possible means to include such a variation. Within the context of the Gompertzian survival model, h0 is often considered to be the age-independent mortality rate coefficient and p the age-dependent mortality rate coefficient [7,25,26]. Suppose that p the age-dependent mortality rate coefficient, which we normally assume to be a constant, in fact is not a constant [23]. Rather, let us assume that we may express p as a stochastic (random) function of the form p ( a ) = Po + oo(a)

(3.1a)

where o~(a) is a mean zero, variance 2E (E is considered the noise parameter) Gaussian white noise process and P0 is a deterministic constant. That is, we are perturbing the constant value bulk mortality rate coefficient P0 with white noise (this could also be performed with colored noise) which is a function of the age of the organisms within a particular population. Other noise forms may be studied as well. In this formulation, we are assuming that p(a) is a function that varies randomly about the mean value P0. That is, the expectation value of the white noise process is E[ w(a)] = 0

(3.1b)

Additionally, we assume that the correlation of the stochastic process o~(a) is given by E[ w ( a 1) oo(a 2 )] = 2 E6(a 2 - a 1)

(3.1c)

leading to a variance o-2(a) = 2e. For E -- 0 there is no stochasticity and Eq. 3.1a reduces to the standard deterministic model (Eq. 2.1). Additional discussion may be found in Karmeshu and Gupta [27], and Witten [28]. An assumption of stochasticity in p is not unreasonable from a biological perspective. We have already observed that P0 is a bulk average value. Hence, it is

Matthew Witten / Mech. Ageing Det'. 73 (1994) 223-248

228

a mean about which individual values may vary. Alternatively, we might assume that p is constant between individual populations (clones) but that it varies stochastically rather than as a drift over time. It is this particular, perspective that we wish to investigate. In point of fact, it is not unreasonable that a mortality rate p a r a m e t e r might fluctuate over time; even if it does so imperceptibly [23]. Only in a truly homogeneous clonal population of individuals might a model in which there was no random variation apply [22]. And, as we shall see, the slightest variation or randomness introduces new concepts into the interpretation of mortality rate curves. In order to examine the effects of stochasticity on mortality rate, and later on survival, we will need to examine the solution of the tollowing general first order stochastic differential equation

dX(t) dt = f ( X ) + g ( X ) w ( t )

(3.2)

where w(t) is a mean 0, variance o-2(t), Gaussian white noise process, f and g are known functions, and X(t) is the variable of interest. Witten [28] has considered the specific case in which

f ( X ) = poX

(3.3a)

g(X) = x This particular case is equivalent to our current model discussion when p is given as in Eq. 3.1a and where we relabel X ~ A and t - a. That is, dA(a)

da

- p°A + Am(a)

(3.3b)

(Various forms of related discussion may be found in Keiding [29]; Kiester and Barakat [30]; Levins [31]; Lewonton and Cohen [32]; Capocelli and Ricciardi [33]; May [34]; and Turelli [35]). This extension of the deterministic mortality rate model to a stochastic mortality rate model gives rise to an interesting concept in survival analysis. No longer may we ask for the mortality rate A(a) at a specified age a. As A(a) is now a stochastic variable, we must ask for the probability density distribution function PA of having a certain mortality rate A at age a. We denote this as follows PA( A, a; h 0) - (the probability density distribution of having mortality rate A, at age a, given that h(j was the initial mortality rate}

(3.4)

The random nature of p(a) has forced us to consider the fact that, at a given age a and for a given initial mortality rate h 0, it is possible to have any one of a number of mortality rates. The likelihood of occurrence of any given value of the mortality rate is based upon a probability distribution pA(A, a; h0), which we shall momentarily determine. Observe that this formulation also implies that the M R D T (mortality rate doubling time) [7,20] will also be a stochastic variable with mean value ( l n 2 ) / p 0. Such a formulation is, as we shall see, quite testable biologically.

Matthew Witten/ Mech. Ageing Deb'. 73 (1994)223-248

229

2.0

1.5 o - 3.0

o"- 0.25

1.0

/, O0

03

06

09

1.2

1.5

1.8

2.!

2.4

2.7

Fig. 2. In this figure we illustrate log-normal probability density functions with p. = 0 and E= 0.25, 1.5, and 3.0. Figure is reproduced from Lawless [35] where

exp[

4. Solving for Px(h, a; h o) and related variables It is possible to solve for the probability density function p , (as based upon Eq. 2.1 of the mortality rates [27,28]. It is given by

p,(A, a; h 0)

2a A ~

{

exp -

~ea- ~

'°a]

A > O.

(4.0)

For fixed age a, Eq. 4.0 is a log-normal distribution in the mortality rate A. That is, the mortality rate values A, at a fixed age a will be log-normally distributed. Fig. 2 illustrates sample log-normal distributions (see Lawless [10] for details on log-normal distribution functions). O u r conjecture concerning stochasticity in p is biologically testable. Consider an experiment in which a large n u m b e r (say n = 1000) of survival colonies of at least 10 000 organisms are watched. T h e n u m b e r 10 000 is arbitrarily chosen so that small sample issues at the tail of the distribution will not be crucial to the analysis. F r o m each of the colonies, it is possible to c o m p u t e a mortality rate curve Aj(a). The n mortality rate curves Al(a) , A2(a) . . . . . Aj(a) . . . . . A,(a) will be distributed, in some distribution, at each fixed age a. That is, if we fix the value of age at a*, and we then examine the mortality rate values Aj(a*) for j = 1. . . . . n, these values, we

230

Matthew Witten / Mech. Ageing Dev. 73 (1994) 223 248

conjecture, will be distributed in a log-normal distribution of the form given by Eq. 4.0 (Fig. 2). Since Eq. 4.0 is a parameterized probability density distribution and since we have a set of values Ai(a) i = 1, 2 . . . . . n (for each fixed value of a), we can fit Pa(A, a; h 0) by standard statistical techniques to the data; thereby estimating Pc~, h0, e. In addition, should such an hypothesis be true, one can then examine the effects of diet, temperature, or other factors upon the parameter estimates. And, as such, assess whether or not contributions to changes in the original deterministic parameters may actually be attributed to changes in E (the level of internal randomness) a n d / o r P0. For example, if we consider the work of Johnson [2,8], we might inquire as to whether the age-1 gene contributes to changing the values of P0 the age-related mortality rate coefficient a n d / o r E, the internal randomness of the biological system, for a Gompertz fit. It is important and necessary that we observe that E is not related to the error of the estimates for each of the model parameters. These error estimates are related to yet another set of variances crp~ and O'e2 associated with the particular choice of fitting and estimation procedure used by the investigator. The value e measures the level of randomness associated with p(a). And, as such, is different from the errors of the fixed constant parameter estimates. We can easily derive additional information concerning the distribution properties of Pa. For example, we can show that the modal or maximum mortality rate (at age a) is given by m l ( a ) = h~e (~'''- zE')"

(4.1)

Further, the median mortality rate is m2(a) = ho et'~'"

(4.2)

Hence, the modal value of A, at a given age a, is constant if ( P0 = 2ae), decreasing if P0 < 2ae, and increasing if P0 > 2ae. The rth moment of the mortality rate A is given by E[A', a; hi,] =h'~;e '"~''~*rE')

r = 1,2 .....

(4.3)

From this it follows that the expected or average value (r = 1) of the mortality rate, at age a is just E[ A] = hoe "~ t,,,4 2,,)

(4.4a)

Observe that, as e ~ 0 (no noise), E [ A ] ~ h o e ~.... which is the deterministic exponential mortality rate distribution with which we began. The variance of A is given by Var [A] = h2e 2(2a~+~''')" (e 4~a: - 1)

(4.4b)

Fig. 3 a - c illustrates the relationship between the noise • level perturbation of the age-dependent Gompertz mortality rate coefficient p and that effect upon the median, expected, and modal mortality rate distributions. In Fig. 3a we illustrate a three dimensional diagram of the mortality rate probability as given in Eq. 4.0. As

2/h0

MRD

2 MRD

M°dal/~J

0

0.5/hO

l.O/hO

1.5/h0

epsilon - 0.05 rho0 squared

3 MRD h0

age

Media

epsilon

mortality rate 2 MKD 3 ~ age ~

~

h0

odal / ~ S h

0.05 rho0 squazed

4 MRD

=

0

0.5/h0

1.0/h0

1.5/h0

0

10 hO 0 mortality fete

5 ~4

h0

eps~ior. = 6.61 rno0 squared

h0

mortality

rate

15 h0

Fig. 3. This figure illustrates the relationship between the noise o-2 perturbation of the age-dependent Gompertz mortality rate coefficient p, and its subsequent effect on the median, expected, and modal mortality rate distributions. In part A we illustrate a three dimensional diagram of the mortality rate probability as given in Eq. 4.0. As can be readily seen, for any fixed age a, there is a log-normal distribution of mortality rates possible. This three dimensional surface is then projected into a two dimensional mortality rate vs. age plane. Observe that for very small noise (a) the expected, median, and modal mortality rate distributions are all nearly equivalent. This is what we would expect if there was little variation in the Gompertz age-dependent coefficient. Parts B and C illustrate the effect of increasing the level of noise in the age-dependent mortality rate p.

0

0.5/h0

1.0/h0

1.5/h0

ep$i~on = 6.01 rno0 square~

4~ Oo

~O

.<

o~

232

Matthew Witten/ Mech. Ageing Del: 7,+(1994) 223-248

can be readily seen from this figure, for any fixed age a, there is a log-normal distribution of mortality rates possible. This three dimensional surface is then projected onto a two dimensional mortality rate vs. age plane. Observe that for very small noise (Fig. 3a) the expected, median, and modal mortality rate distributions are all nearly equivalent. This is what we would expect if there was little variation in the age-dependent mortality rate coefficient. Fig. 3b,c illustrates the effect of increasing the level of noise in the age-dependent mortality rate coefficient p. Given that the mortality rate A is now a random function of age, we may also ask the following question: What is the probability that the mortality rate a lies in some given range or interval at a given age a? That is, we wish to compute the probability that A lies in the interval [A~, A2] at age a. We denote this cumulative distribution function by Fa(Al:, a; h0), where A~

(4.5)

FA(AI2, a; h o) = f -pa(X, a; htl)dx Ai

If we s e t A 1 = 0 and A: = A*, Eq. 4.5 may be integrated to obtain the following closed form for the cumulative distribution function Fa, the cumulative distribution of mortality rate values A in the interval [0, A* ] given an initial mortality rate h.. ? J~(a*, a; h~,) = # e r f

(4.6)

2azd7

where err(x) is the error function with argument x. 5. First passage and the stochastic mortality rate function k(a)

We may also ask the following question: What is the probability that A(a) will exceed the critical value A~ for the first time, in the age interval (a, a + Aa), given A(0) = h0? This is equivalent to asking at what age a will the mortality rate A exceed a particular critical value A~ for the first time. Mathematically, this is called a first passage problem. We are interested in determining the functional form of FA(Ac, a , h0), the required probability density distribution for first passage.; FA is the probability density function that the mortality rate a exceeds the critical value a,,, for the first time, at age a, given an initial mortality rate h.. Karmeshu and Gupta [27], and Witten [28] have demonstrated that this density function is given by

ra(A,, a; hi,)

[ h° ~ ( ln[T") expl2 a ~

In

poa

--7

1

and that it is also a log normally distributed variable for fixed values of a.

(5.0)

Matthew Witten/ Mech. Ageing Dev. 73 (1994) 223-248

233

The most likely age of first passage (the age at which the critical value is most likely to be exceeded for the first time) is given by the mode of Eq. 5.0.

12aE

P0 = 0

--3+ V/((3)2--(~)21n(-~c°))

(5.1)

portO

We can now consider two problems. The first, what is the first passage distribution for the mortality rate Ac = h 0 + 7 In (N) which is the mortality rate at the time t N for which the survival distribution (Gompertz) takes on the value SOn) = 1 / N . The second problem is to consider the question of the first passage distribution for the mortality rate Ac = 2h 0, which is the mortality rate at the time t = M R D T the mortality rate doubling time [7]. The M R D T is given by M R D T = (In 2)/3,. Substitution of this into the Gompertzian mortality rate equation (1.3) yields the critical mortality rate Ac = 2h 0. Substitution of this critical mortality rate into Eq. 5.0 and after some algebraic manipulations, we arrive at the following equation for the M R D T first passage time which is given by

Fa(A c = 2 h 0 , a ; h o ) = a/yRD~/( T_..7_8rrea2)l exp -

[MRDT-a]21 'I[ 8Y2 " a 2 t]

~

(5.2a)

Introducing the change of variables t/= a / M R D T allows us to transform to a dimensionless age first passage distribution given by

1

M R D T X Fx(Ac = 2h0, a; h 0) =

, exp

(

[1 - a ] 2

7) 8ca 2

(5.2b)

This equation is plotted, for different values of the noise (epsilon ~) in Fig. 4. The formula, for t N is identical, save that we replace M R D T with t N. In fact, we can make the general statement that, if t c, is the time at which the Gompertz mortality rate takes on the value Ac, then if we use this particular value Ac in Eq. 5.0, we will end up with the following equation for the first passage distribution

1

GXFa(a~,a;ho)=g

t

/

[1 -a] 2

/[8,1.1.~e2]exp~ (8~(~2 ] Vl-

-I

3`2J

(5.2b)

Matthew Witten / Mech. Ageing Det~ Z~ (1994)223-248

234

epsiton 2.5

i

= 0.1 " g a m m a " 2 i

i

i

epsilon = 0 . 0 1 " g a m m a " 2 epsilon = 0 . 1 " g a m m a " 2 epsilon = 1 . 0 " g a m m a * ' 2

=o-

- ..... ....

1.5

c~ E (.9

.,

,. ,,, .,

£3

,, i' ", /

rr

i 0.5

"'..

/ / J

0 0

"

i

0.5

) 1

1,5 agetMRDT

2

2,5

3

Fig. 4. In this figure we illustrate the graph of Eq. 5.2b (the dimensionless age first passage distribution for the mortality rate doubling time MRDT as based upon a Gompertz model). Each of the curves illustrates a different amount of noise.

w h e r e w e n o w define fi = alto. This allows us to consider any first passage time (median, mean, etc.) and the associated mortality as our first passage variables of interest; thereby generalizing the result. Let us s u m m a r i z e our discussion up to this point. W e have e x a m i n e d the effect of introducing stochastic variation into the G o m p e r t z mortality rate parameter p - - t h e age-dependent mortality rate coefficient. W e have s e e n that this introduces a n e w level of complexity into the classical mortality rate analysis of survival curves in that w e n o w have probability distributions of mortality rates at each age. Thus, w e can no longer specify the mortality rate at a given age a. Rather, w e must ask for the probability of having a certain mortality rate at age a. It follows, as a natural c o n s e q u e n c e of this construction, that w e would expect to see a variation in survival curves even though the populations are s u p p o s e d to be clonal. W e now consider what happens if w e introduce stochasticity into a survival m o d e l . 6. Stochasticity in the survival m o d e l In the previous sections, w e considered the effect of r a n d o m perturbation of the G o m p e r t z slope p p a r a m e t e r and h o w that would affect our vision of the mortality rate. W e n o w consider h o w stochastic perturbation can affect a deterministic

Matthew Witten /Mech. Ageing Deu 73 (1994) 223-248

235

survival model. In particular, let us now consider the simple constant rate mortality model giving rise to exponential or wild-type survival, as given by Eq. 2.3. This survival model is just

= -h°S

(6.0)

S(0) = 1 If we assume that h0, the age-independent mortality rate coefficient [7,25,26] is a random variable of the form

ho(a) = (ho) o + oJ(a)

(6.1)

satisfying E[ho(a)] = (ho) o then it follows, directly from previous discussion that the probability density distribution for S is given by

1 { (lnS+(ho)oa) 2} ps(S, a; 1) = 2aS2v~-~-~ exp 8~a2

(6.2)

We can no longer specify the probability of surviving until a given age a. Rather, we must ask for the probability of having a certain probability of survival. That is ps(S, a; 1) is the probability density function for the probability distribution describing the probability S of surviving to age a given S(0) = 1. It follows, from our previous discussion, that the mode of (6.2) is given by ml(a)

= e ( - h ° - 2")a = e -(h°+ 2a')a

(6.3)

and the rth moment of the survival distribution

E[S r]

=

e ra(-h°-rae) r = 1, 2 ....

(6.4)

From Eq. 6.4 it follows that the expected (average) survival value, at any age is just

E[S] = e a(-h°-a')

(6.5)

Observe that as ¢ ~ 0 (the case of no noise), we arrive at the standard deterministic exponential survival model. In Fig. 5a-c we illustrate the effect of increasing the noise pertu[batlon of the age-independent Gompertz mortality rate coefficient h 0 in the simple exponential survival model. Fig. 5a illustrates the result of small noise perturbation or variation. Again, we illustrate the three dimensional probability of survival surface and its projection into the survival vs age plane. As we would expect, the mean (expected) and the modal survival distribution values are nearly identical. As we increase the level of noise (Fig. 5b,c) we see these two distributions become increasingly separated. We may also compute the probability that 0 < S < S* at age a. This is given by the cumulative probability distribution Fs(S*, a; 1)

Fs(S*, a; 1) = _

ps(x, a; 1) dx 2 erf[lnS* + h 0 a ]] ~aa~--~

(6.6)

age

0

5

15

20

age

2 hO

probability density

C

~d!

~

epsilon = 0.01 h0 squared

hO

Expected

epsilon = 0.05 h0 squared

8 1

15 I0

probabilitydensity 20

B

0.2

0.6 0.4 survival

O.8 i

epsilon = 0.02 h0 squared

1

Fig. 5. Illustrates the effect of increasing the noise perturbation of the age-independent Gompertz mortality rate coefficient h o in the simple exponential survival model given in Eq. 6.0. Part A illustrates the result of small noise perturbation or variation. Again, we illustrate the three dimensional probability of survival surface and its projection into the survival vs. age plane. As we would expect, the mean lexpected) and the modal survival distribution values are nearly identical. As we increase the level of noise (B,Ct we see these two distributions become increasingly separated.

0

5

15

20

probability density

A

4~ ~e

~D

g

2

Matthew Witten/ Mech, Ageing Dec. 73 (1994) 223-248

237

h0"Gamrna(a,S;1) gam(a,s)

h0*age

3.5

4

I o.o.

45

0.1 '

5

Fig. 6. We illustrate the first passage probability density distribution for the case in which the Gompertz parameter h 0 (the age-independent mortality rate) is a stochastic variable. That is, we illustrate the first passage probability density distribution for the probability that the survival distribution will be less than a critical amount S by a given age a. The deterministic case would be an exponential survival model. This is clearly indicated by the exponential behavior of the location of the peak ih the illustrated surface. Observer, however, that for any fixed value of age, there is a distribution in possible survival values.

As before, we may ask for the probability density furfedon describing the probability that S(a) will exceed the value So, for the first time, in the interval (a, a + An) given S(0) = 1. Because survival functions are decreasing functions over time, this question is not correctly posed. What we really wish to ask is the following question. What is the probability density function describing the probability that the survival is less than or equal to the value S~ for the first time in a given infinitesmal time interval. The previously described first passage problem may be correctly posed in terms of S~ and that yields following equation as the required answer s

Iln(Sc)l

[

Fs(Sc, a; 1 ) = 2a2 2v~--~e e x p ] -

(In(So) +h°a)2 / 8ca2 -)

(6.7)

We illustrate this probability density distribution in Fig. 6. Eq. 6.7 is also of use to us in that we can determine the probability density distribution for population extinction. That is, if a population starts off with N individuals, then the value S = 1 / N is the survival value at which only one individual is alive. That is, for S < 1 / N the population has become extinct. It follows that when we substitute

238

M a t t h e w Witten / M e c h . A g e i n g D e v . 73 (1994) 2 2 3 - 2 4 8 1.2 S S S S

il

= 0.01 = 0.05 =0.10 =0.15

-..... ...... .......

i 7, 0.8

ilia: 03

0.6

iil/::

E

il

0.4

i 0.2

\

:i i :

I "

\ ".

/// ]

jl 2

4

8

8

10

age h 0

Fig. 7. In this figure we illustrate the distribution of first passage times as a function of age, for different present m a x i m u m lifetimes. Observe that, as we get closer to the upper percentiles of the survival distribution, the first passage time for the m a x i m u m lifespan becomes more and more G a u s s i a n looking.

S s = 1/N into Eq. 6.7 we obtain the probability density function describing the probability that the survival function S(a) will be less than 1/N in the time interval (a, a + Aa) given S(0)= 1. Thus, the density function is given by the following equation

Fs(1,a;l)=

ln(~}lexp 2a2 21/2-~

-

(6.8) 8ca 2

Observe that, for fixed sample size N, this is a Gaussian distribution in age a. Fig. 7 illustrates simulated survival distributions and the distribution of maximum lifespans as a function of sample population size. From these figures it is clear that the simulation results are Gausslan in a for fixed population sample size N. And, it is also clear that this model predicts the simulated behavior. This is of interest in that it shows that stochastic variation in the h 0 Gompertz parameter and the random effects of pure sampling may be indistinguishable. This would make it very difficult to assign variation in parameter values to either sampling variation effects and/or internal genetic/environmental effects. Details on the simulation algorithm appear in section 7 below.

Matthew Witten/Mech. Ageing Dev. 73 (1994) 223-248

239

7. Biological implications and conclusions Variation in parameter estimates for survival models, such as the Gompertz or Weibull survival model, may come from a variety of sources. Variation may be due to heterogeneity in the experimental population. Hence samples drawn on such a population are likely to have some variation in parameter estimates due to component sample variation. The very act of sampling of a population can carry with it a component of variation, even if the population is known to be genetically pure [36]. Additionally, sample size can play into the amount of variation seen [37]. Additional contributors to variation include environmental factors. Diet, existence of environmental pathogens, exercise, and a host of other factors can and do affect lifespan and alter survival; subsequently affecting parameter estimates from sample to sample. Coupled to these contributing factors is the effect of individual frailty and/or vitality in this face of environmental challenge. In this paper we propose an alternative, though not mutually exclusive explanation for variation in survival model parameter estimates. In particular, we assume that the model parameters, such as h 0 and p in the Gompertz model, are not constant over the period of the experiment (lifetime); an already illustrated fact [23]. Rather, we assume that they are randomly perturbed about some deterministic or mean value. The result of such an assumption is that the mortality rate (hazard function) is no longer deterministic. Instead, for any fixed time (or age value), there is a probability density distribution of possible mortality rates that might occur. This distribution is given by Eq. 4.0 (for a particular example). This result has particular bearing on the issue of variation of parameter estimates in a genetically pure colony which is sampled into small subsets. The results of this theoretical analysis demonstrate that, all things being equal, even a genetically pure colony would exhibit variation in the survival model parameter estimates. The amount of this variation, in the case of white noise variation would be + 24 where 24 is the variance of the noise function ~o(a). Such a result implies that we will, at best, know a distribution of model parameter values; even in a genetically pure stock. This hypothesis is clearly experimentally testable. A large population of genetic clones is sampled into a number of smaller samples. For each sample, the survival model parameters are estimated. Given that all influences are controlled for, any parameter variation is intrinsic systemic noise. Probability distribution models can be fitted to the mortality rate values and a variance estimate can be obtained. One way to examine the effects of sample size on survival model parameter estimates is to simulate samples drawn from a known population. The algorithm for this discussion is illustrated in Fig. 9. First, a survival model is picked. In the case of the upcoming discussion, we will make use of the previously discussed two parameter Gompertz model. This approach is not limited to the Gompertz survival model alone. It may be readily applied to any parametric survival model. Having picked the survival model, we next pick the actual values of the survival parameters h 0 and P0- For the purposes of our analysis, we chose values close to realistic values of known biological populations. In particular, we chose h 0 = 0.002 per day

i

5O

I

800

J 1000 11~

1200

13Q0

1400

1500

150

fSO

0600

50

2OO

2OO

100

250

25O

D

5o

lOO

150

20O

250

3OO

350

400

300

UZO 10

E;

3OO

age at Oe~h. ~

~ze TO

35~

sL~n y~" ~

agl (days}

s x ~ w a rnaxrnum age IJ oea~, ~

35O

4OO

0600

50

100

150

2GO

250

300

350

400

I

700

I

800

s~!~

8o0

1000

~ m

900

' 12~

1400

A

death, sample ~ e IO0

1200 age (Oays)

age at

11too

death, sarape s~ze lOO

1000 age (trays)

san~'u a maxtrnum age at

1600

13Ioo

]~ig. 8. In this series of four figures we illustrate the results of simulating random sampling effects, based upon a Gompertzian distribution and using the ad libitum and diet restricted data of Yu et al. [38]. The actual lifespan data is available from the GAIA MultiSpecies Survival Database (see Ref. 43 for details). Two sample sizes, for each diet protocol, were simulated (N=I 0 and N=I00). For each sample size N, 1000 samples of that size were generated and the maximum lifespan (maximum age at death) distribution was calculated. This distribution is illustrated in part A (N=10, ad libitum rats, part B (N=100, ad libitum rats, part C/N=10. diet restricted rats), and part D (N=I00, diet restricted ratst. Observe the similarity in shape of these distributions to the predicted distributions based upon stochastic perturbation of h 0 initial mortality rate parameter.

A

4~ ~e

241

Matthew Witten /Mech. Ageing Dev. 73 (1994)223-248

Weibull Define survival model to be studied

~--~

Gompertz ] .-- Other

ho = 0.002

Fix parameter values for given survival model

,l Fix number of samples and size of samples

7=2"31

I

.•.------• 1,000 samples [ ~__~

size of sample 8,16,32,..1024 I

Generate lifespan samples and write to disk for storage

Compute statistics

Fig. 9. In this figure we illustrate the simulation diagram used to generate the various samples for the Gompertz parameter estimation simulation. The investigator initially picks his model and parameter values (Gompertz Weibull, whatever). The number of samples and the different sample sizes to be generated are then input. The computer generates the requisite samples and writes them to files. These files contain lifespan samples of the requested sample size and it contains an input number of repeats (or drawings) at that sample size. The statistics are then calculated and are summarized in Table 1 of this paper.

Matthew Witten/Mech. Ageing Deu. 73 (1994) 223-248

242

C) r-

n=8

E to

(~ 13.. 0.0075 N t:: CL

E

o £9

0.0050 n-16

>

n-32

%°, E 0.0025

n=512 DO-O--

-0-

-

-

0

.

.

.

.

~

.

.

.

.

n=1024 ' •

.

c

'G 0.0000 <:

I

0

250

,

,

,

,

I

500

,

,

,

,

I

.

.

.

.

750

I

1000

.

.

.

.

I

1250

Size Of Sample (n) Fig. 10. An illustration of the variation of the estimates of the Gompertzian survival parameter h o (day)- l (also called the age-independent mortality rate coefficient) as a function of sample size when 1000 samples of each size are generated. Observe that the true coefficient estimate should be h 0 = 0.002. The open circles represent the exact answer and the filled black diamonds represent the simulation mean value.

and P0 = 2.31 per day as based upon the original ad l i b i t u m / d i e t restricted rats groups of Yu and Masoro [38]. Several tens of thousands of simulated samples of various sizes drawn from a Gompertz distribution typical of laboratory rat lifespan data were generated and subjected to both nonlinear regression and maximum likehood analysis techniques; providing samples of the parametric estimate distributions. Nonlinear regression is a more sophisticated methodology for fitting survival curves than is frequently used in the gerontological literature. It is also routinely assumed to be much more accurate than log-transformations of the survival curve and subsequent use of linear regression techniques. Details of the problems associated with the various methodologies of parameter estimation in survival curves may be found in Eakin et al. [39]. Simulations results indicate that the Gompertz parameter h0 is surprisingly hard to estimate from small sample datasets (n < 256). In the following discussion, it should be noted that we chose values of the Gompertz model parameters close to that of known published literature for rat (Fischer 344) populations. Fig. 10 illustrates the actual and estimated values of h 0 for a simulation in which 1000 samples of size n = 8, 16, 32, 64, 128, 256, 512, and 1024 are drawn from a known distribution h 0 = 0.002 per day. The open circles represent the actual value of the parameter and the black diamonds represent the mean over 1000 samples of the

Matthew Witten/ Mech. Ageing Dev. 73 (1994)223-248

243

n=8

0 r-

== 0.0150 > ¢.. ¢0 0.)

0.0100 LL

n=16 0,)

1:3

0.0050 =32

t~ 10 t,--

2 5n ~ 6 _ . _ _

0.0000

|

i

i

i

I

25O

,

,

l

=

~I I

n=512 •

500

n=1024 o

o

,

I

.

.

.

.

750 Size Of Sample (n)

--I

-'~" .

1000

.

.

.

I

1250

Fig. 11. A n illustration of the variation in the standard deviation of the estimate for the m e a n h 0. Observe that there is order of magnitude error as large as sample sizes of n = 512 (significantly larger than most diet-restriction experiments).

estimated h 0. Fig. 11 illustrates the standard deviation of the estimate for h 0 in the same population. We summarize the data in Table 1. Summary of simulation results for 1000 populations of different sample size n = 8, 16 ..... 1024. For this simulation, h 0 was fixed at a known value of h 0 = 0.002 per day. Observe that for sample sizes as large as n = 128 (greater than most of the published rat literature) the standard deviation in the estimates for h 0 is of the same order as the value of the estimate itself. In Fig. 12 we illustrate the effect of changing the sample parameter h 0. The same experiment was run for a h 0 set at double and at half the original value of h 0 = 0.002 per day. The same dynamics can be seen for these values of h 0 as well. Fig. 13 illustrates the sample variances for all three cases. We also examined the estimates for P0 the age-dependent per capita mortality rate coefficient (the slope of the log-linear mortality rate line). Fig. 14 illustrates the simulation results for the estimate of P0 over the same sample sizes as used for estimating h 0. Observe the large amount of fluctuation in the estimates (over 1000 samples) for the P0 parameter. The true population value is P0 = 2.31 per day. Clearly, as the size of the sample population becomes large, we obtain a more accurate estimate for P0. However, observe the wide fluctuation in estimates for samples even as large as n = 512. Interestingly, the sample variances of the

Matthew Witten / Mech. Ageing Det~ 73 (1994) 223-248

244

sample 0.016

I

i

means I

of

h

0

!

h 0 = 0.002 h

0

=

0.004

~---

h-O = 0 . 0 0 1 -o

0.014

0.012

0.01

0.008

0.006

0.004

.................

0

t

"

.............

° ...................

l 200

0

+ ........................................

2

; ..........................................

l 400

I 600

1 800

I 1000

1200

Fig. 12. I n this f i g u r e w e r e p e a t t h e e x p e r i m e n t i l l u s t r a t e d in F i g u r e [10] w i t h t h e a d d i t i o n a l p r o v i s o t h a t w e a r e r e p o r t i n g t h e r e s u l t s u s i n g h a l f a n d t w i c e t h e v a l u e o f F i g u r e [11] h 0 = 0 . 0 0 2 p e r d a y as t h e e x a c t v a l u e . A s c a n b e s e e n , t h e r e s u l t s a r e n e a r l y i d e n t i c a l in t h e i r b e h a v i o r w i t h r e s p e c t to t h e e s t i m a t e o f t h e m e a n h 0.

sample variances J i

0.025

of

h

0

,

I

l* 0 = 0 . 0 0 2 h-O = 0 . 0 0 4 -+-h-O = 0 . 0 0 1 O 0.02

0.015

0005t 0.01

"

0

I 0

"-----i- . . . . . . . .

I 200

+

....... i . . . .. ... ... . . . .e ~. .... ... . . . . . . . . .I. . . . . . . . .- .. .. ... ... .. . . . . t-. . .. ... ... .. ... ... .. . . . . . . . . 400

600

800

T ~*~-

1000

1200

Fig. 13. T h i s f i g u r e i l l u s t r a t e s t h e s a m p l e v a r i a n c e s f o r t h e e x p e r i m e n t r u n in Fig. 12. A g a i n , w e s e e n e a r l y i d e n t i c a l b e h a v i o r in t h e v a r i a n c e o f t h e e s t i m a t e s .

Matthew Witten/Mech. Ageing Deu. 73 (1994) 223-248

sample ,

1.2

245

variances ,

of

gamma ,

h 0 = 0.002 h-0 = 0.004

-e--~--

h-O = 0.001 -o--

0.8

0.6

0.4

0.2

I 200

A 400

I 600

I 800

I000

1200

Fig. 14. We used the same datasets to estimate the value of P0 the age-dependentmortality rate (slope of the log-mortality rate curve). In this figure we see that the estimates vary quite wildly with sample size, particularly with small sample sizes (as large as n = 128).

sample

means

of

gamma

2.35 ~; 2.34

h 0 = 0.002 -¢"-h-O 0.004 -+-'-" h-O = o . o o i o =

~"

ii 2.33

2.32

2.31

l!"i.... Q.. I o-÷. "'-.. 4-"" "-... - ' . .

2.3

2.29

2.28

2.27

2.26

2.25

I

I

I

I

I

200

400

600

800

1000

1200

Fig. 15. Observe, however, that the variation in the sample variances, as a function of sample size, is very uniform and predictable. Further, it is extremely tight for all cases of h 0.

246

Matthew Witten / Mech. Ageing Der 73 (1994) 223-248

estimates for P0 are very tight, even for small population samples. This is illustrated in Fig. 15. These simulation results indicate that estimates of Gompertz model parameters, particularly for small samples sizes may be very spurious and unreliable. As we have seen, even for large numbers of sample populations (n = 1000), the estimates have order of magnitude standard deviation. This calls to question the value of using any of the currently used methods for estimating Gompertz parameters for experimental populations. 8. The GAlA Multispecies Survival Database M. Witten, R. Shouman, T. Eakin, Y. Qi, and G. Liu, The GAlA Multispecies Surt~ual Database, is an effort to organize multispecies genetic and survival data into a database structure that is publicly available. The database is currently accepting data submissions. Details may be requested by writing to the following address: GAIA Multispecies Survival Project, UT System Center for High Performance Computing, BRC, 1.154 CMS, 10100 Burnet Road, Austin, TX 78758-4497 USA. Or, you may call 1-800-262-2472 and ask for a Survival Project package to be sent to you. You may fax us a request at 1-512-471-2445. Lastly, you may mail an electronic request to the database server at [email protected] with the subject line containing the word ADMIN or HUMAN and the body of your message containing your postal mail address. The GAIA Project is supported, in part by a grant from SUN MicroSystems Corporation, a grant from the University of Texas System Center For High Performance Computing, a grant from IMSL Corporation, Inc. (now Visual Numerics Corporation), and grant NIH 1 RO1 AGl1079-01. The author's electronic mail address is: mwitten @chpc.utexas.edu 9. Acknowledgements This research was partially supported by a grant from the Nathan and Margaret Shock Aging Research Foundation. I would like to thank Frank Guess for his many discussions on survival theory, Richard Greenberg for his reading of earlier drafts of this manuscript and for his statistical mentorship and his ongoing interest in this research project, and Kenneth Cooke for interesting me in the stochastic aspects of population models [41, 42]. I would like to thank Tom Johnson for giving me permission to reproduce his figures and for his continued support of these efforts. I would also like to thank David Harrison, Richard Miller, and Devendra Dubey for their time spent in discussion on experimental verification of the theoretical ideas presented in this work. I would like to thank Radey Shouman for generating the three dimensional figures and for performing the supercomputer simulations discussed in this manuscript (supported by a grant from the Nathan and Margaret Shock Aging Research Foundation). I would also like to thank him for proofreading the manuscript. I would like to thank Tim Eakin for his meticulous rereading of

Matthew Witten/ Mech. Ageing Dev. 73 (1994) 223-248

247

the final version of this manuscript and for his many excellent and clarifying comments. Lastly, I would like to thank the referee for his comments and suggestions. 10. References 1 D.B. Friedman and T.E. Johnson, A mutation in the age-1 gene in- C. elegans lengthens life and reduces hermaphroditic fertility. Genetics, 118 (1988) 75-86. 2 T.E. Johnson, Aging can be genetically dissected into component processes using long-lived lines of C. elegans. Proc. Natl. Acad. Sci. USA, 84 (1987) 3777-3781. 3 J. Heicklen and E. Brown, Increase in life expectancy for mice fed with diethylhydroxylamine (DEHA). J. Gerontol., 42 (1987) 674-680. 4 B.P. Yu, E.J. Masoro, and C.A. McMahan, Nutritional influences on aging of Fischer 344 rats: I. Physical, metabolic, and longevity characteristics. J. Gerontol., 40 (1985) 657-670. 5 K.E. Cheney, R.K. Liu, R.E. Leung, M.R. Mickey, and R.L. Walford, Survival and disease patterns in C57BI/6J mice subjected to undernutrition. Exp. Gerontol., 15 (1980) 237-258. 6 R. Gelman, A. Watson, E. Yunis, and R.M. Williams, Genetics of survival in mice: Sub-regions of the major histocompatibility complex, (1989), preprint. 7 C.E. Finch, M.C. Pike, and M. Witten, Slow increases of the Gompertz mortality rate during aging in humans also occur in other mammals and birds and suggest that short lived rodents and birds diverged from ancestors with long potential lifespans and slow senescence. Science 249 (1990) 902-905. 8 T.E. Johnson, Age-1 mutants of Caenohabditis elegans prolong life by modifying the Gompertz rate of aging. Science, 249 (1990) 908-912. 9 D.E.L. Promislaw, Senescence in natural populations of mammals: A comparative study. Evolution, 45(8) (1991) 1869-1887. 10 J.F. Lawless, Statistical Models and Methods For Lifetime Data. John Wiley and Sons, N.Y. (1982). 11 M. Witten, Reliability theoretic methods and aging: Critical elements, hierarchies, and longevity interpreting biological survival curves, (in) Molecular Biology of Aging (eds.). A. Woodhead, A. Blackett, and A. Hollaender (Plenum Press, N.Y., 1985). 12 M. Witten, A return to time, cells, systems and aging: Relational and reliability theoretic approaches to the study of senescence in living systems. Mech. Aging and Dev., 2 (1984) 323-340. 13 M. Witten, A return to time, cells, systems, and aging: III. Gompertzian models of biological aging and some possible roles for critical elements. Mech. Aging and Dev., 32 (1985) 141-177. 14 M. Witten, A return to time, cells, systems, and aging: IV. Further thoughts on Gompertzian survival dynamics-The neonatal years. Mech. Aging and Dev., 33 (1986) 177-190. 15 M. Witten, a return to time, cells, systems, and aging: V. Further thoughts on Gompertzian survival dynamics-The geriatric years. Mech. Aging and Dev., 46 (1988) 175-200. 16 B. Gompertz, A sketch on the analysis and the notation applicable to the value of life contingencies. Philosophical Transactions of the Royal Society, 110 (1820) 214-294. 17 B. Gompertz, On the nature of the function expressive of the law of human mortality and on a new mode of determining the value of life contingencies. Philosophical Transactions of the Royal Society, 115 (1825) 513-585. 18 B. Gompertz, A supplement to two papers published in the transactions of the Royal Society, On the science connected with human mortality; the one published in 1820, and the other in 1825. Philosophical Transactions of the Royal Society, 52 (1862) 511-559. 19 J.E. Riggs, Longitudinal Gompertzian analysis of adult mortality in the U.S. 1900-1986. Mech. Aging and Dev., 54 (1990) 235-247. 20 M. Witten and W. Satzer, Gompertz survival model parameters: Estimation and Sensitivity. Appl. Math. Letters 5(1) (1992) 7-12. 21 R. Lestienne, On the thermodynamical and biological interpretation of the Gompertzian mortality rate distribution. Mech. Aging and Dev., 42 (1988) 197-214

248 22

23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38

39 40 41 42

Matthew Witten / Mech. Ageing Dev. 73 (1994) 223-248 F. Guess and M. Witten, A population of exponentially distributed lifespans cannot lead to Gompertzian or Weibull (with increasing mortality rate) survival dynamics. Bull, Math. Biol., 50 (1988) 187 192. A.L. Ekonomov, C.L. Rudd, and A.J. Lomakin, Actuarial aging rate is not constant with the human lifespan. Gerontol., 35 (1989) 113-120. F.A. Lints, Rate of living theory revisited. Gerontol., 35 (1989) 36-57. C.E. Finch, Longetity, Senescence, and the Genome, (Chicago University Press, Chicago, 1990). M. Witten, Re-examining the Gompertzian model of aging, Institute for Mathematics and its Applications, University of Minnesota. IMA Preprint Series 483 (1989), Minneapolis, MN. Karmeshu and C.K. Gupta, A one-compartmental model with stochastic parameter. Bull. Math. Biol., 43 (1981) 503-512. M. Witten, The diffusion process approach to compartmental stochastic models: A mathematical note. Bull. Math. Biol., 45(3) (1983) 425-430. N. Keiding, Extinction and exponential growth in random environments. Theor. Pop. Biol., 8 (1975) 49-63. A. Ross Kiester and R. Barakat, Exact solutions to certain stochastic differential equation models of population growth. Theor. Pop. Biol., 6 (1974) 199-216. R. Levins, The effect of random variations of different types in population growth. Proc. Nat. Acad. Sci., 62 (1969) 1061-1065. R.C. Lewontin and D. Cohen, On population growth in a randomly varying environment. Proc. Nat. Acad. Sei., 62 (1969) 1056-1060. R.M. Capocelli and L. Ricciardi, A diffusion model for population growth in a random environment. Theor. Pop. Biol., 5 (1974) 28-41. R.M. May, Stability in randomly fluctuation versus deterministic environments. Amer. Natur., 107 (1973) 621-650. M. Turelli, Random environments and stochastic calculus. Theor. Pop. Biol., 12 (1977) 140-178. R. Shouman and M. Witten, Survival estimates and sample size: What can we conclude. J. Gerontol., (1994) in press. R. Shouman and M. Witten, Studying the biology and the evolution of aging: Sample size issues in survival model parameter estimation. Bull. Math. Biol., (1993) submitted. B.P. Yu, E.J. Masoro, I. Murata, H.A. Bertrand, and F.T. Lynd, Lifespan study of SPF Fisher 344 male rats fed Ad Libitum or restricted diets: Longevity, growth, lean body mass, and disease. J. Gerontol., 37 (1982) 130-141. T. Eakin, R. Shouman, Y. Qi, G. Liu, and M. Witten, Estimating survival model parameters: problems and insights. J. Gerontol. (1994) in press. J.E. Riggs and R.J. Millecchia, Mortality among the elderly in the U.S. 1956-1987; Demonstration of the upperbound to Gompertzian mortality. Mech. Aging and Dev., 62 (1992), 191-199. K.L. Cooke, R. Elderkin, and M. Witten, Harvesting procedures with management policy in iterative density-dependent population models. Natural Resource Modeling, 2 (1988) 383 420. K.L. Cooke and M. Witten, One-dimensional linear and logistic harvesting models. Math. Modeling, 7 (1986) 301-340.