Bulletin of Mathematical Biology, Vol. 41, pp. 305-324 Pergamon Press Ltd. 1979. Printed in Great Britain 0 Soaely for Mathematical Biology

MATHEMATICAL APPROACHES OF CANCER CHEMOTHERAPY

0007_4985/79/0501aO5

$02.00/O

TO OPTIMIZATION

n S. ZIETZ and C. NICOLINI Division of Biophysics, Department of Biophysics-Physiology, Temple University Health Sciences Center, Philadelphia, PA 19140, U.S.A.

This paper uses optimal control theory in conjunction with a Gompertzian type model for cellular growth to determine the optimal method of administering cycle non-specific chemotherapy or more generally the optimal durations of treatment and rest periods during chemotherapy. The performance criteria employed to determine the relative merits of the therapy include not only the destruction of malignant cells, but also the sparing of a critical normal tissue. Since these criteria are at odds with one another, the solutions are found which satisfy the Pareto optimality conditions.

1. Introduction. Cancer is. one of the most dreaded diseases of our time, killing approximately 350,000 individuals a year in the United States alone. One of the difficulties that we face in treating the wide variety of cancers is that almost all cancers metastasize. Metastasis is the process by which cells break away from the primary tumor site and spread through various pathways like the blood and lymph systems to distant sites of the body where they take up residence and form secondary growths. Once metastasis has occurred, treatment of cancer by localized methods like surgery and radiation therapy is no longer curative since these methods cannot destroy the multitude of microscopic metastatic growths, each of which have the capability of growing to life-threatening size. Therefore, in such instances, systemic therapy must be undertaken. The most successful of such treatments to date is chemotherapy in which patients are given drugs which are supposed to destroy the tumor cells. However, these drugs also destroy normal cells, especially the fast growing regenerative tissues like the bone marrow and gut epithelium. Thus, in administering chemotherapy the physician is dealing with a two-sided sword. In designing a chemotherapeutic protocol he must not only consider how well the drugs will destroy the malignant tissue, but also 305 11

306

S. ZIETZ

AND

C. NICOLINI

how best he can protect vital normal tissue and minimize other toxic reactions of the treatment. Often it is not obvious how “best” to administer chemotherapy in the face of conflicting criteria. Such treatment must always play off the benefits of treatment versus the risk of complications. Mathematical modeling of chemotherapy can help in exploring various strategies of treatment. In recent years various mathematical descriptions of cellular growth and the effect of chemotherapeutic agents on the growth have been developed (Aroesty et al., 1974; Nicolini et al., 1975, 1977; Zietz, 1977). Using the models it is possible to simulate therapy and to potentially make predictions of which protocol might be better than others. However, using these models, it is also possible to try to optimize the treatments. That is, out of all possible treatments (with a confined group of possible strategies for treatment), predict the “best” one. Of course, the “best” must be made into a precise mathematical statement encompassing all the relevant criteria. Also, it is conceivable that the answers may be model dependent, that is, the models that seem to be equally plausible given our present biological knowledge might yield completely different answers. In this case we must go back to basic biology and do the critical experiments to distinguish between the models. However, if all plausible models give approximately the same answers, then we can feel a little safer that we have captured the essential aspects of the optimization problem. Our approach to the optimization of cancer chemotherapy is to make judicious use of experiments, mathematical models and optimization techniques to explore optimal treatment strategies. Our modeling of chemotherapy takes place at many levels depending on the biological question under consideration. We feel that the mathematical complexity of the model to be considered must be directly related to the phenomena under consideration. Although a complex mathematical model might yield good quantitative predictions between experiment and theory often a less complex one will equally account for important experimental observations. Also, as the complexity of a model increases, the more intractable it becomes; and often the complexity causes great difficulty in trying to optimize for the best treatment. There have been several mathematical attempts to apply optimal control theory to design of treatment protocols. Noteworthy in this respect are Almquist and Banks (1976) and Swan and Vincent (1977), which besides containing excellent work cite much of the relevant published literature. In this work, we examine several simple cellular kinetic models of cell cycle non-specific chemotherapy * to determine the optimal method of *The models proposed let a patient rest between

more generally treatment.

apply

to the problem

of how often to treat and how long to

MATHEMATICAL

APPROACHES

TO OPTIMIZATION

OF CANCER

CHEMOTHERAPY

307

administering treatment. In these models it is assumed that the pharmacodynamics of the agents does not play an important role. In future works, we will present the results of using other more realistic models and the comparison between the predicted response to treatment and the actual in viuo response. But to begin our discussion of cell cycle non-specific chemotherapy it is reasonable to observe some facts about the administration of chemotherapy. Chemotherapy is usually given on a cyclical basis. The rationale behind such protocols are simple: during the treatment intervals, both normal and malignant cells are killed. During the treatment free intervals, the normal tissue (and unfortunately the tumors) have a chance to recover and replenish their depleted numbers. The goal, of course, is to eradicate the tumor while keeping the normal cells at a reasonable level. Such a protocol has a hope of achieving its objectives if (1) tumor cell kill is greater than normal cell kill or (2) normal cell growth is faster than tumor cell growth so that during the treatment free intervals, the normal cells recover faster, thus allowing us to treat before the tumor has grown to its pre-treatment size. Clearly, the success or failure of such cyclic chemotherapy lies in the knowledge of how long to treat and how long to let the patient recover. If we treat for too long or if we do not let the patient rest long enough between treatments, extreme toxicity will occur due to depletion of the body’s normal cells, especially the bone marrow. However, if we don’t treat long enough, or wait too long between treatments, we allow the tumor to recover and may even permit the tumor to increase in size despite the treatment. Thus, the timings of when to give treatment must take into account the relative rates of growth of the tumor and critical normal tissues and also the killing efliciency of the chemotherapeutic protocol used. 2. A Mathematical Model. To try to study the above questions, it is first important to get a model of tumor and normal cell growth which approximates the actual growth of these cells. The Gompertzian function (Laird, 1969) has been successfully used to describe the growth of both tumorous and normal populations of cells both in vitro and in vitio. Therefore, it will be used as the description of cellular growth. Second, it is important to understand how drugs interact and kill cells. It is observed that a given concentration of drug kills a given fraction of the cells, not a given number. Thus, to a first approximation, drug killing obeys first order kinetics. By this, we do not necessarily imply that the logarithm of cell killing vs. drug concentration curves are linear. All we assume is that there is a function k(c) of the concentration of the drug which gives the fraction of the population killed.

308

S. ZIETZ

AND C. NICOLINI

Of course, when dealing with a chemotherapeutic protocol, we are not dealing with the administration of just one agent, but with an entire range of treatment modalities. We must not think of administering a concentration of a drug, but of treating a patient with a certain level of “aggressiveness”. By aggressiveness, we are thinking of a global version of drug concentration. Thus, a more aggressive protocol would kill more cells than a less aggressive one. To get a handle on the lethal effects of drugs, let us examine a typical survival curve. Most dose-response curves are S shaped (on a logarithmic dosage scale) (Goodman and Gilman, 1974) since low concentrations of drugs do not affect the cells while at high concentrations (extreme aggressive treatment). a point of diminishing returns is reached probably due to resistant cells. In the middle region, the response is usually monotonic and, in our approximations, we will assume it to be linear. With the above assumptions, our model of cell growth becomes: N=N[k-plnN]-u(t)yN(t)=N[k--fiInN--yu(t)],

(I)

where N(t)= the number of cells at time t given N(0) cells at time t=O, u(t) is the level of “aggressiveness” of our treatment at time t, and k, p, y are various parameters of the model. When no treatment is administered, the model assumes that the cell growth is Gompertzian. When u(t)= 1 [on some appropriate scale], we have N =N[k -/? In N - y]. Notice this equation is also Gompertzian. If we assume k - p In N - y < 0, then the population will decrease monotomically until N =O, this is until k -/?ln N -y =0 or N =exp(k -y//I). Observe that if y > k -/II In N and y 9 k, this equilibrium value of N is almost 0, which can be interpreted as representing the complete eradication of the tumor. Also in this case, the decrease in size is approximately exponential. Let us now consider the constraints on u(t). Since the treatments destroy cells, u(t)20 with u(t)=0 when no treatment is given. However, toxicity usually limits the aggressiveness of the treatment (respectively the concentration of a specific drug) that can be administered at any given time. For instance methotrexate, a highly successful agent used to treat many malignancies causes neurotoxicity when used in high doses. To incorporate these toxicities into the model we assume that at any given time there is a maximally aggressive protocol (respectively maximum concentration of drug) that can be given. The parameters in the model are set so that this maximum corresponds to u = 1. Therefore 0 5 u(t) 5 1. With the above assumptions, we can now ask, what is the best manner to administer the treatments. However, we have not yet defined what we mean by “best”. Suppose we focus our attention on the tumor cell

MATHEMATICAL

APPROACHES

TO OPTIMIZATION

OF CANCER

CHEMOTHERAPY

309

population, and one normal cell population, say bone marrow. We can then ask what is the best treatment protocol to minimize the tumor population, while maximizing the normal population at the end of a certain amount of time, say T days. The best will be defined in the Pareto sense, that is, we seek solutions whereby if we decrease the final tumor population the final normal population must also decrease. (See Appendix A). With the above assumptions, the control problem becomes: Problem 1. Let: xi(t) represent the number of normal cells present at time r, x2(r) represent the number of tumor cells present at time r, where x1 and x2 obey the following dynamics

~i(t)=X,[ki-BilnXi-riU(t)], where x,(O), ki,Bi,yi

i= 1,2

i= 1,2 are given constants and

u(r) is subject to the control constraint Oiu(r)s 1. Then given the treatment interval [0, T] we wish to maximize x1 (T) Notice the performance .form :

and

minimize x2 ( T )

criteria can be written in the following equivalent

s s T

maximize J1 =

2i (t)dt

0

T

minimize J, =

I,(r)dr.

0

Before proceeding to solve the above problem, some assumptions about the relative normal and tumor growth rates must be made. Since tumors are usually less inhibited in their growth by density effects than are normal cells, we assume that the /? for the tumor population is less than the /I for the normal population. As will be shown, this assumption is sufficient to completely characterize the form of the optimal treatment schedule. 3. Mathematical Solution to the Optimal Control Problem 1. In order to solve the Pareto Optimal control problem, we utilize the Pareto version of the Pontryagin Maximum Principle established by Yu and Leitmann (see

310

S. ZIETZ

AND C. NICOLINI

Appendix B). We therefore form the Hamiltonian: H= -cx,~~ +a&

+&x& = -~+,Ch

+/llx,[kt

-Prlnx,

-y,u(t)]

-82Inx2-y2Wl -PI

(2)

lnx, -wWl

which must be minimized subject to the differential constraints ~i=Xi[ki-pilnXi-Yiu(t)],

Ai=

_~+1y+1 i

cliB-;li[ki-Bi

Xi(O) = Xi0

In Xi-_iU(t)

1

+(-l)‘cci[ki-pil

i=1,2.

IlXi-YiU(t)]+il$i,

The requirement that the adjoint variables, be perpendicular terminal manifold tf = T implies that ,I.,(T) = 1, (T) = 0. Rearranging the Hamiltonian, we can write H=u(t)Cy,x,h

to the

-3L1)-~2~2(~2+3L2)1

+[-~a,xk,+a,/3,x,lnxI+a2x2k2-a2~2x21nx2+;llklxl -1I/3Ix,

lnx, +A2x2k2 -n2/32x2 lnx,].

Since u(t) enters linearly in the Hamiltonian, in order to minimize the Hamiltonian at each instant of time we need to choose 0 if s(t)>0 u(t)=

1 if s(t)p2, the optimul trajectory cannot switch from 0 to 1.

Recall the switching function is

and since eDiT> 0 for all t, the sign of S is determined by

Since /?r >b2, S(t) is an increasing function of t. Therefore if s^(t,)>O at any time t,, S(t) 20 for all t 2 to. However S(t)>0 implies that u =O, so if u =0 at any time, u =0 for all subsequent times, thus eliminating case 3. Q.E.D. The other three cases are not only mathematically, but also clinically feasible. Notice that by choosing a2 properly (putting emphasis on tumor

MATHEMATICAL

APPROACHES

TO OPTIMIZATION

OF CANCER

CHEMOTHERAPY

313

cell kill) we can arrange the S(0) ~0, so the initial stage of the protocol is treatment. To determine the actual optimal control sequence for a specific set of parameters we will now integrate the state equations (1). Let t, be the time of “switch” of protocol from u = 1 to u = 0. Note t, = 0 implies Case 1 while t,= T, implies Case 2. Thus the following analysis includes all possible dynamic behavior that the optimal treatment can exhibit. To find the optimal control u(t) notice that on the interval [O,t,] ii = Xi[ (ki - Ui) - Bi ln xi],

1,2,

i=

or x.(t)=x.(O)e-@exp I I

ki - Yi

(

,tE(O,t,).

gl-[l-ePD~r] I

(8)

>

Therefore k,-y,

7[1

x.(t,)=x.(O)e-Vsexp I

I

-e-4’s]

.

I

>

Now on [ts, T] li = Xi[ ki - fli In xi]

,

i=

1,.2,

or Xi(t)=Xi(ts)e-Bi”-zdeXp

(

$[l--e

-Bi(t-rs)],

i=

1,2.

I

Thus Xi(T)=Xi(ts)e-PicT-ts)eXp

(;[I

_,-sic-q

, i=1,2. >

1

Substituting the derived expressions for x1 (ts), x,(t,) and simplifying, we obtain

k

(equation 8) into (9)

Ee-Pi’[l

-ePirs]

which gives the values of the state variables at the terminal function of the parameters and the switching times.

(11)

time as a

314

S. ZIETZ

AND

C. NICOLINI

Recall that the switching time t, is characterized as the [unique] zero of the switching function. The switching function is zero when

Substituting obtain

the expressions

for x1(T)

and x,(T)

and rearranging

we

(P2-Pl)fs+;e-

(P2+al)T&tP2+B,)=

1

I

y,a,x, (0)e-PlT

In

w2x2(0,

epPJ

exp $1 ( !!?[I

exp (

_e-Plr] >

e(P*-P1)T

(11)

_e-PZT] >

!

which is a transcendental equation for the switching time t,. Given the parameters of the system, the length T of the total treatment interval and the initial numbers of cells, we can solve for the switching time. Varying c(~ and a2 can achieve various final values for the state variables. However, all solutions satisfy the necessary conditions for Pareto Optimality. However, a severe limitation to the usefulness of the model is that the only factors that entered into the performance criteria were the population levels at the end of the treatment. In the next section we will alter the performance criteria to take into account population levels during the treatment. 4. A Model With More Realistic Performance Criteria. In this section we will consider a slight generalization of the previous model which, although less tractable mathematically than the previous case, is more realistic in assessing the “cost” of treatment. Recall that the preceding model presumed there were two populations of cells that were important to monitor, that they obeyed Gompertzian dynamics, and the treatment objective was to maximize the number of normal cells while minimizing the number of malignant cells, at the end of a fixed treatment interval. It was shown that the general form of the optimal trajectory consisted of alternating periods of maximum treatment with periods of no treatment, although in that specific case there was at most one switch between these modalities. An explicit equation for the time of this switch was derived. However the preceding problem did not take into consideration that

MATHEMATICAL

APPROACHES

TO OPTIMIZATION

OF CANCER

CHEMOTHERAPY

315

there are costs involved in having too few normal or too many malignant cells at any time during treatment. For instance, if the normal leukocyte population falls below a critical level the patient requires special sterile room techniques since his resistance is low and thus he is extremely susceptible to infection. In fact, a major limitation to therapy seems to be toxicities caused by killing too many normal cells. Similarly, if the malignant population is allowed to get too large, the patient shows clinical signs of the disease. Thus a more realistic performance criteria would ask to keep the tumor size low and the normal population high, for as long as possible during treatment, while again achieving tumor cell kill. To incorporate these factors into the model, let us define two population levels M,, M2 and ask that during the entire treatment interval, the normal population level in excess of Ml is maximized and the tumor population level below M, is minimized. Mathematically the cost criteria become

s T

maximize

0

(x,(t)-Mr)dt

T

minimize s0 With the above modifications,

(x2(t)--Mz)dt.

our model now becomes

s s T

Maximize J 1 =

(xl(t)--MlW

0

and

T

Minimize J, =

(x2(f) - M,)dr,

0

subject to the dynamics ai=Xi[ki-plnXi_yiu(t)],

Xi(O)=XiO,

i-172,

and control constraints 0 5 u(t) 5 1 over the time interval [0, T]. As before we form the Hamiltonian

H=alWfl -~~(t))+a,(x,(t)-M~)+~,(t)x,(t)(k, . +~,(t)x,(t)Ck,-Bzlnx,-y,u(t)l

-PI lnx, -ylu(t)

316

S. ZIETZ

AND C. NlCOLlNl

and seek the control constraints

u(t) which minimizes H subject to the differential

liCx,Ck, -&lnxi-Yiu(t)],

Xi(O)=Xio, i=l,2,

where

and 1(T)=O, i= 1,2, since (A,, ,I,) is perpendic,ul;lr to the terminal manifold tJ= T. Rearranging the Hamiltonian we get

H= -u(t)Cy,l,(t)x,(t)+y,~,(t)x*(t)l f%(M, -xl(t))+Q(x,(t)--2) W)xi(Nh-Pl

lnxll+&C~2-B21nx,l.

As before, since the dynamics are linear in the control variable to minimize H, we need to choose 0

u(t)=

if

S(t)>0

1 if

S(t)