B U L L E T I N OF ~r

BIOLOGY

VOLU~ 37, 1975

:MATHEMATICAL ANALYSIS OF CANCER C H E M O T H E R A P Y

[ ] SHUI-NAN CHUANGand HARRIS H. LLOYD Kettering.Meyer Laboratory, Southern Research Institute, Birmingham, Alabama 35205

This paper presents a mathematical analysis of a tumor model first proposed b y Skipper and Zubrod. The tumor model is comprised of three compartments, a proliferative compartment, a nonproliferative b u t viable compartment, and a dead compartment. By the suitable selection of functions describing loss of cells from the proliferative and nonproliferative compartments, the model is capable of describing tumor behavior during periods of growth and drug treatment. The loss functions during t r e a t m e n t are related to pharmacokinetic functions and m a y be chosen according to known drug properties. Tumor properties m a y be simulated b y the appropriate choice of cell cycle parameters. I t therefore seems feasible to simulate tumor behavior for scheduled treatment with chemotherapeutic agents. Another important result of this analysis is the derivation of a fraction labelled mitosis function which incorporates the nonproliferative compartments.

1. Introduction. I n recent years there has been considerable effort in applying the known growth characteristics of tumor cells toward improving the effectiveness of cancer treatment. It became clear that an extensive study of pharmacokinetics of antitumor agents, cell cycle kinetic analysis of tumor and sensitive normal tissues, drug-cell interactions and host toxicities is necessary in designing an optimal dose and schedule for cancer chemotherapy. Although some of the involved areas have been analyzed theoretically as mentioned in the following paragraphs, a comprehensive mathematical study for designing an optimal dosage regimen has not been developed. Among a variety of techniques for cell cycle analysis, the most successful and prevalent one is the fraction labelled mitosis (FLM) method of Quastler and Sherman (1959). The experimental FLM curve expresses the fraction-147

148

S H U I - N A N CHUA:NG A N D H A R R I S H. L L O Y D

number of labelled mitoses/number of mitoses--as a function of the time elapsed since the administration of a single short pulse of tritiated thymidine to the population. The effect of the pulse is to label the cells which are in S phase at the time of administration. Mathematical expressions for the theoretical FLM function have earlier been derived by Barrett (1966), Takahashi (1966, 1968), Truceo and Brockwell (1968), and MacDonald (1970). In these different approaches normal, log-normal or gamma distributions were chosen for the probability density functions of the phase durations, and optimization criteria, such as maximum likelihood, least squares, or Monte Carlo techniques were used to fit the theoretical FLM functions to experimental curves. All of these recent models for cell cycle analysis were concerned only with the proliferating cells. However, there is evidence (Lajtha et al., 1962) which suggests that tumors are composed of at least three types of cell populations-proliferating (consisting of the four cycle phases G1, S, G2, M), temporarily nonproliferatiug (Go, resting) and dead (D, permanently nonproliferating or destined to die) cells (Skipper, 1969). Proliferating cells are senstive to all anticancer agents, provided that drug concentration and duration of exposure are sufficient (Perry, 1971). Nonproliferating or resting cells in vitro exhibit varying degrees of sensitivity, depending on the nature or classification of the drug. I t will be presumed that the temporarily nonproliferating cells are capable of returning to the proliferative state at G1. The permanently nonproliferating cells, having lost the ability to divide, are of no concern to the chemotherapist. Nevertheless, their presence is reflected in physical determinations of tumor size (solid tumors) or cell count (leukemias). The above considerations initiated our investigation of the population growth kinetics of these three different types of cells in untreated tumor systems. Also, since the cells in the resting state and the four consecutive phases of the proliferating cycle respond differently to drug treatment, the representations of tumor behaviors during the treatment for each proliferating phase and resting state by separate equations will express explicitly the cell-cycle or cellcycle-phase specificity of antitumor drug in the model. Some theoretical papers in mathematical studies of predicting cell kill by drug treatments have been published (Himmelstein and Bisehoff, 1973; Jusko, 1973; Bisehoff et al., 1973), but the aforementioned concerns have not been done. Based on the model of Figure 1, we will consider in the following sections a theoretical FLM function, which is in terms of parameters of both proliferative and nonproliferative compartments, and the functions for the cell kinetics of the population distributed in: (1) the untreated tumor system where the loss functions are constants; (2) the tumor system under drug treatment where the loss functions are related to the drug concentration in the tumor site; and (3) the

MATHEMATICAL ANALYSIS OF CANCER CHEMOTHERAPY

149

tumor system during multiple-dose, scheduled treatments with chemotherapeutic agents. A description of computerized procedures for predicting the population fluctuations during the schedules is also given. In these analyses we assume t h a t the transit times of cells in the proliferating phases are random variables with independently distributed probability density functions.

'L5 ~ I

~6 ..

GI S G2I M A

II

Figure 1. Model for tumor growth and drug treatment. The proliferative compartment consists of the four phases G1, S, G2 and 1YL After each binary fission, (2-A) cells enter the nonproliferative compartment Go, and A cells continue in the proliferative cycle where 1 _< A _ 2. Cells may die either naturally or by drug treatment as determined by the functions Lt(t), i = 1-5. When leaving Gi, S, G2, M or Go, they enter compartment D, the dead cell compartment. Loss from the tumor site is determined by 16. Also, a certain proportion, p, of G0-cells may re-enter the proliferative cycle at G1 The algorithms and computer programs for simulating the cancer chemotherapeutic effects and fitting the derived FLY[ function (which incorporates the nonproliferative compartments) to experimental data to predict the cell kinetic parameters are being prepared.

2. Mathematical Model.

The mathematical model of a tumor shown in Figure 1 is based on the three-compartment system suggested b y Skipper and Zubrod (Skipper, 1969). The proliferating compartment consists of four phases (G1, S, G2, M) through which the cells pass in t h a t order. The cells in each phase m a y lose their ability to divide, due to natural factors or drug treatment, and convert to the permanently nonproliferating compartment D, in which the cells are dead or will die. A binary fission following mitotic phase is assumed with a 'distribution parameter' A/2, the expected fraction of mitoses returning to phase G1, and a 'decycling probability' 1 - (A/2), the expected fraction entering the resting state G0. A can have values between one and two during growth. The resting cells may revert to the proliferating cycle or die as a result of drug treatment or natural factors.

150

S H U I - N A N C H U A N G A N D H A R R I S H. L L O Y D

In Figure 1 and the following derivations, the subscripts 1, 2, 3, 4, 5 and 6 are related to G1, S, G2, M, Go and D respectively. The population density functions n~(t, a~) of the proliferating phases satisfy Von Foerster's population balance equation (Von Foerster, 1959) ~n~(t, a~) ~n~(t,ai) ~ § ~a-----~.= -[2~(a,) + L~(t)]n~(t, a~), i = 1, 2, 3, 4

(1)

where the a~, i = 1-4, represent cell age within each phase and the L~(t) are loss functions. The boundary conditions n~(t, O) = a~(t) (the birth rates) are

nl(t , O) = A

2~(a4)n4(t, a4) da 4 + pNs(t) ,=o

ni+ l(t, O)

(2)

~(a~)n~(t, a~) dai,

i = 1, 2, 3

|~0

where pNs(t ) is the rate at which G o cells revert to the proliferating state. The initial conditions (initial age distributions) axe n~(0, ai) = fl~(a,),

i = 1, 2, 3, 4.

(3)

In (1) and (2), the functions 2~(ai), i -- 1-3, are the transfer coefficients, and 2~(a4) is the generation coefficient. They are related to the transit time's probability density functions f~(a~) b y the expressions 2i(a,) = fi(a')

i ---- 1, 2, 3, 4

r

(4)

with r

fi(x) dx = exp

= 1 =0

If; -

]

2~(x) dx , =0

i--- 1 , 2 , 3 , 4 .

(5)

As mentioned in the introduction, a variety of functional forms can be fitted to the probability density functions f,(a,), i = 1-4, of the transit times which are assumed to be independently distributed. The loss functions L,(t), i = 1-4, in (1) and Ls(t ) in (6) are the rate constants for cells losing their ability to divide due to drug treatment or some other factor. In an untreated tumor system Lt(t ) -- 1,, i = 1-5, are small constant coefficients, while for tumors under chemotherapeutic treatment, they are functions of drug concentration C(t) at the tumor site. The drug-cell interaction relationship will be discussed later in Section 5. The term pNs(t ) in (2) is the rate at which cells revert to the proliferating cycle from the resting state. That is, p represents the rate constant for resting cells

M A T H E M A T I C A L A N A L Y S I S O F CAL~'CER C H E M O T H E R A P Y

151

entering phase GI. The expected number of cells Ns(t) in the resting state at time t is governed by the equation

d'-----i-- = - [ p + Ls(t)]Ns(t) + (2 - A)

;~4(a4)n~(t,a4) da4.

(6)

4=0

The initial value, Ns(O), must be chosen according to known tumor characteristics. Finally, the governing equation for the expected number of cells in compartment D is dNr(t) dt

5 1,N,(t) + ~ L,(t)N,(t)

(7)

$=1

where 16 is the rate constant for D cells leaving the tumor site. The initial condition/Vr(0 ) depends on tumor characteristics and will be discussed later. The expected number of cells in each phase of the proliferating cycle is obtained by integration of the population density functions n~(t, a~), i.e. N~(t) =

f

n~(t, a~) da~, i = 1, 2, 3, 4.

(8)

t=0

The totM population in the system at time t, then, would be 6

v(t) =

(9) i=I

3. Cell Kinetics of Tumor System. A single dose of an antitumor chemotherapeutic agent will provide cytotoxic concentrations for some interval, depending upon the characteristics and pharmacokinetics of the drug. However, for m a n y drugs fractionated doses provide greater therapeutic effectiveness with relatively less host toxicity. Although both normal and tumor celIs are affected during the interval(s) of cytotoxic concentrations, we will be concerned in this section only with the effect upon tumor cells. The following analysis is general in the sense of accommodating (or simulating) any treatment schedule with a drug of known (or assumed) properties against tumors with known (or assumed) kinetic parameters. Only gross drug characteristics are required, i.e. eel/cycle-specific or -nonspecifie, degree of effect upon GO cells and serum concentration with time. The analysis or simulation of treatment is accomplished over a series of time intervals. There is an interval of growth before treatment, an interval during which the drug level is above the threshold cytotoxio concentration, and an interval of regrowth during which the drug concentration is less than the minimum effective concentration.

152

S H U I - N A N C H U A N G A N D H A R R I S H. L L O Y D

During treatment with cytotoxic drugs, the loss functions L~(t), i = 1-5, for each proliferating phase and the resting state are functions of drug concentration at the tumor site (see Section 5). For the untreated tumor system or for tumor regrowth between treatment intervals, the loss functions L,(t) = I, i --- 1-5, and 16 in (1), (6) and (7) will be presumed independent of time and cell ages. There are six simultaneous differential equations for the cell population in six states. To separate the system (1)-(3) of the proliferating compartment from those of Go and D compartments, N~(t) must be expessed in terms of n4(t, a~) and substituted into the boundary condition nl(t, 0) of (2). Let U~(t) = exp

[f -

] foo f

L~(~-) d , , T=0

Then from (6) we have

[

No(t) = e-~tUs(t) N5(0 ) + (2 - A)

i = 1, 2, 3, 4, 5.

Ua(~.)

~0

(10)

]

,~4(a4)n~(~-,a4) da4 d~4=0

(11)

in terms ofn4(t, a4). With this substitution, (1)-(3) give the population density functions for the proliferating phases (Trucco, 1965).

Ia~(t - ai)~(ai) ui(Ui(t_)ai ), t >_ a~ hi(t, a,)

=

|R.ta. -- t~ L r''' '

' r

r

i = 1, 2, 3, 4.

U.ID, - t)

"'

(12)

t < a,

These r e l a t i o n s state t h a t a cell, w h i c h is i n phase i w i t h age a s < t a t t i m e t,

must have entered phase i at time t - a i and have survived to age a~, and a cell with age a, > t must have been in phase i with age a, - t at time zero and survived to age a t. I t is reasonable to assume exponential growth when the untreated tumor is small. Thus, the initial age distributions take the form

fl,(a,) = Kir

e-(C+',)a,,

i = 1, 2, 3, 4

(13)

where c is the specific growth rate of the tumor when in the early stage of exponential growth. For a later simulation interval, the initial age distributions fii(ai), i = 1-4, and initial populations N~(0) and N6(0 ) are equal to the population density functions of the proliferative phases and cell populations in compartments Go and D at the end of the previous interval, i.e.

fii(ai) = nt(t~, at), i = 1, 2, 3, 4t N~(O) = Nt(t~) , i 5, 6 where t~ is the final time of the previous interval.

(14)

MATHEMATICAL ANALYSIS OF CANCER C H E M O T H E R A P Y

153

The recursive equations for the birth rate functions at(t ) for the treatment interval are obtained as

al(t ) = A Ur

"~(t -- ~) f~(a~) dar 4=0

~=t r

t) f4(a4) da~ + pe-PtUs(O Ns(O )

a~(r - a~) A(a ) da4 A) I t eO~ U4(r) f ~ U----~T)LJa,=o ~ a4) J ~ =0 fl4(a~- r) da4~dr ]

+ (2 ~|

/

(15)

J

a~+i(t)= Ut(t)(~21=o~t(t--~f~(a,)d a, +

,=t r

t)

f~(a~)da,~ ,

i=

1,2,3

by substituting the population density functions n~(t, a~) of (12) into the boundary conditions of (2). Equation (15) can only be solved iteratively by numerical methods with a reasonably chosen starting function for any one of the a~(t). For the growth interval, a constant function with value at(0 ) of the initial birth rate will be chosen as the starting function while for the first treatment interval, the starting function will be a constant function with value a~(tf), the birth rate at the final time of the untreated interval. For subsequent treatment (growth) intervals, the starting function can be any one of the a~(t) of the preceding interval scaled by the ratio of the population sizes at the beginning of the current interval and the preceding interval. With the birth rate functions thus obtained, the expected cell populations in each state are N (t) = u,(t) { f l

a , ( t -- a,) ,=o Ui(t at) r

r ~ fl~(ait) r + ja,=t r t)

da~ da~}

[

Na(t) = e-~

Ns(0 ) + (2 - A)

x

{

fi4(a4- r)

N6(t) = e -z6t N6(O ) +

;

e ~ U4(T) --o

Udr)

f4(a4) da~

,=o U4(r

+ rt

i = 1, 2, 3, 4

fl

=0

da~}d1-]

e~ ~

]}

L,(~-)N,(r) d r Li=I

(16)

154

SHUI-NANCHUANG .4~7D HARRIS H. LLOYD

The total expected number of cells in the tumor is the sum of the cells in each compartment expressed b y (16). In (13), K~ are constants defined as the birth rate at time zero and depend on the size of the initial population. B y setting t = 0 in (15) they are related b y

g l = AK4L(c + 14) + pNs(O) ] K~+I Kif,(c + l~), i = 1,2,3 f

(17)

with one of them as a scaling factor obtained from experimental data as discussed later in the simulation procedures. Also, from (17) we see that the relative growth rate is determined uniquely as the nonnegative solution of the equation

1 = Afl(c + ll)L(c + 12)f3(c + la)L(c + 14) + pN~(o) K-q-T-"

(18)

However, the estimation of c must be performed iteratively with the estimation of the phase parameters offi(a~) b y a curve-fitting process of a theoretical F L M function to the observed data.

~. Theoretical -FLM .Function.

With the superscript asterisk (*) referring to quantities associated with labelled cells, the theoretical F L M function for any period of time can be derived from (1) and (2) with the initial age distributions

n*2(O, a2) = f12(a2) n~*(0, ai) = 0,

}

(19)

i = 1,3,4

and Ns*(0 ) = 0. That is, an expression for the expected labelled mitotic cell population N4*(t) can be obtained exactly as was done in the preceding section for the expected total mitotic cell population N4(t ). The F L M function is F L M (t) = N*(t) N4(t)

(20)

However, the observed F L M curve is taken from the untreated tumor system. In this case, fl~,(a2) is assumed to have the form of (13) and U~(t) = exp (-l~t), i=1-5. B y the Laplace transformation method, the reeursive integral (15) gives

a~(s)

_ b(~)

[L(c + l~) - L ( s

+ l~)]L(s + l~)

(21)

with

~9(s) = 1 -

[

l A + 8 + p + 15J/i(s + ll)L(8 + 12)L(8 + t3)L(s + l~) (22)

MATHEMATICAL ANALYSIS OF CANCER CHEMOTHERAPY

155

and the Laplace transform of N~*(t) can be obtained from (16), i.e.

~*(s) =- d*(s)~4(s + 14)

(23)

(~(s) = [1 - L(s)]/s.

(24)

Similarly, the analytic solutions for the Laplace transforms of a4(t) and N~(t) of the total mitotic cells are obtained as 1

a~(s) - b(s)

[{

pNa(0 )

S + p +15

+ [A+

• ] - I ~ ( s + zj) + i=1

+ s - ck~

p(2-A)] 1s 9 s + p -~-I5 ~ [ f 4 ( c

+ 1 4 ) - f~(s +14)]

}

[Z(~ + l,) - Z(s + l,)] 1-[ L(s + lj) ~,=i

1=~+1

[L(c + z~) -L(s + z~)]]

(25)

and

_~(s) = a~(s)q~(s + q) + s K- ~ [q~(c + q) - ~ ( s + Z~)].

(26)

The functions N~(t) and N~(t) in the time domain can be calculated by a computer program for numerical inversion of Laplace transforms (Simon et al., 1972). From the definition of FLM function and (23) and (26), the curve-fitting process, with a chosen optimization criterion will give the estimates of cycle phase parameters.

5. Drug Dependency of the Loss Functions. In the preceding section, the loss functions were considered only generally. A more detailed description of these functions is given here, so that the analysis presented above can be used to simulate the kinetics of populations treated by cytotoxic drugs. As mentioned before, for tumors under chemotherapeutic treatment the loss functions depend on the drug concentration C(t) in the tumor site. This is assumed under the cell-kill hypothesis of Skipper et al. (1964) that for most drugs a clear relationship exists between the drug dose and the extent of eradication of tumor cells within the limits of toxicity to the host, and that a given dose of a drug kills a constant fraction of cells, not a constant number, regardless of the cell numbers present at the time of therapy. With this hypothesis, the loss functions are expressed by the general form Flmax, i

C(t) > Cmax, t

L~(t) = ~l_ma~.,[C(t) - Cm~.,] + Z,[C~=.~ - C(t)], !

Ll~

C(t)