- Email: [email protected]
Computer models of anticancer drug interaction
Pharmac. Ther. Vol. 4, pp. 245-280, 1979. Pergamon Press Ltd.
Printed in Great Britain
Specialist Subject Editors: ALANC. SARTORELLI,WILLIAMA. CREASEYand JOSEPHR. BERTINO
COMPUTER
MODELS OF ANTICANCER INTERACTION
DRUG
R. C. JACKSON Laboratory for Experimental Oncology, Indiana University School of Medicine, Indianapolis, Indiana 46202, USA
and K. R. HARRAP Department of Biochemical Pharmacology, Institute of Cancer Research, Sutton, Surrey, UK
1. TYPES OF SYSTEM WHICH HAVE BEEN STUDIED USING COMPUTER MODELS Anticancer drugs, administered to a tumor-bearing host, may interact at a number of levels of biological organization. This report will consider four kinds of interaction. Firstly, p h a r m a c o d y n a m i c interactions: this category of interaction refers to mutual effects on the spatial disposition of two or more agents. Thus, one drug may alter the intestinal absorption, tissue distribution, transport across cellular membranes, or renal excretion of other agents. Secondly, inter-tissue metabolic interactions may arise, e.g. some drugs are only active after being chemically modified by hepatic microsomal enzymes, and such drugs would be antagonized by agents which inhibited those microsomal enzymes; conversely if one drug inhibits the conversion of another to an inactive product, then potentiation may occur. Thirdly, intracellular biochemical interactions may be expected. This category differs from the second kind of interaction, in that the latter concerns interactions involving drug activation or inactivation within organs of the host, and these effects modify the response to a drug or drug metabolite at a distant site; by intracellular biochemical interactions is meant the effect that two or more drugs may have on each other within the same cell. Such interactions fall broadly into two kinds: (a) an influence of one agent on activation or inactivation of another, e.g. an adenosine analog may be active in the nucleotide form, but because the nucleotide cannot cross the plasma membrane, it must be administered as the nucleoside; in this case another agent which is an adenosine kinase inhibitor may well be antagonistic for the first agent. (b) The second kind of biochemical interaction occurs when two or more compounds lead to opposing changes in biochemical pathways; for example, arabinosylcytosine (AraC) acts in many cells as the triphosphate form, which inhibits DNA polymerase, and this inhibition is competitively opposed by deoxycytidine triphosphate (dCTP); thus if another drug leads to elevated nuclear dCTP levels, it may serve to antagonize the cytotoxic action of AraC. Fourthly, drugs may interact in relation to their cyte~'inetic effects. If, for example, a particular compound is active specifically against S-phase cells and another agent blocks progression from G~ to S, then the second agent will antagonize the action of the first compound. The progression through the various stages of the cell cycle is, of course, determined by intracellular biochemical events; however, it is convenient to distinguish this category of drug interaction from the third category in that the dimension of time is here a critical feature. Mathematical models have been written to describe drug behavior at all of these levels, and many of these models are considered below. The models, in general, were written with the intention of understanding better the pharmacology of a specific drug in a particular experimental system. Models of the effects, on tumor and/or host, of JPT Vol. 4, No. 2--A
245
246
R. C. JACKSONand K. R. HARRAP
two or more drugs, which may include, or predict, interactions of each drug on the effects caused by the other drug(s), have not been widely used. Clearly other kinds of interactions between antitumor drugs are possible. For example, if one agent is immunogenic and another is immunosuppressive, they will interact at the level of the immune system. Such interactions are excluded from this discussion, not because they are unimportant, but because they have not been described in sufficient quantitative detail to form the basis of mathematical models. 2. THE NEED FOR MATHEMATICAL MODELS OF DRUG INTERACTION The need for these models reflects the fact that there are certain questions that cannot readily be answered without them. The reasons for this follow from the nature of antitumor drug-target interactions. (i) Antitumor drug-target interactions are complex. The biochemical make up of a cell involves many thousands of components. An anticancer drug will typically react at several points, and each such reaction may lead to a cascade of consequential effects. Consider, for example, the antimetabolites 5-fluorouracil (FU) and 5fluorodeoxyuridine (FUdR) (Fig. 1). These are generally considered to act as antineoplastic agents primarily after conversion to the level of the deoxynucleotide, FdUMP, which is a potent inhibitor of thymidylate synthetase. Figure 1 shows that there are at least three possible routes from FU (or FUdR) to FdUMP. Each of the fluorinated intermediates produced will probably influence to some extent the corresponding reaction with the normal uracil or thymine intermediate; these perturbations, in turn, will alter the rates of other reactions. (ii) Biological systems are highly interactive. Consider, for example, the highly simplified diagram of methotrexate (MTX)-FdUMP interactions (Fig. 2). Thymidylate synthetase (TS) is inhibited by FdUMP, but the inhibition is dependent upon methy-
TPi
I FiTP I
I F
P I
[ FUMP
l_
~1
I
FUR
I
I
Pip~~ R-I-
FUdR I
~dR_l)i
I-P
catabolism FIG. 1. Metabolism of 5-fluorouracil and its nucleoside and nucleotide derivatives. FU, 5-fluorouracil; dR-I-P, deoxyribose-l-phosphate; Pi, inorganic phosphate; FUdR. 5fluorouracildeoxyriboside; FUR, 5-fluorouridine; FUMP, 5-fluorouridine monophosphate; FUDP, 5-fluorouridin¢ diphosphate; FdUDP, 5-fluorodeoxyuridinediphosphate; FdUMP, 5-fluorodeoxyuridine monophosphate; FUTP, 5-fluorouridinetriphosphate; FCTP, 5-fluorocytidine triphosphate.
247
Computer models of anticancer drug interaction MTX I I I
,
MTX"/
~
I~lrines
CH2FHI'~ 4
dUMP r
FIt~. 2. The interaction of MTX and FdUMP. FH2, dihydrofolate: FH4, tetrahydrofolate: CH2FH4, N-S,10-methylenetetrahydrofolate; dUMP, deoxyuridine monophosphate; dTMP, thymidine monophosphate; FdUMP, 5-fluorodeoxyuridine monophosphate; MTX, metho-
trexate. lenetetrahydrofolate (CH2FH4), since FdUMP binds to the TS--CH2FH4 complex. The level of CH2FH4 depends upon the activity of dihydrofolate reductase, which is inhibited by MTX. Thus, MTX can modify the efficacy of FdUMP. However, the effectiveness of MTX, in turn, will depend upon the concentration of the competing substrate, FH2, which is produced by thymidylate synthetase. Thus, the presence of MTX influences the effectiveness of FdUMP, and vice versa. Such interactions are particularly complex where feedback regulation is involved; consider the combination of AraCTP with FdUMP. The latter drug leads to lower levels of cellular dTI'P, which would be expected to decrease the activation of ribonucleotide reductase for purine ribonucleoside diphosphate reduction, thus establishing a deficiency of substrates for DNA polymerase (Reichard, 1972). Now dCTP competes With AraCTP for binding to DNA polymerase, and as a competing substrate would be expected to accumulate behind the AraCTP block; this increased dCTP may then activate deoxycytidylate deaminase, leading to an increased concentration of dUMP. This action would reduce the inhibition by FdUMP, and tend to restore the dTTP level. What actually happens to the nucleotide concentrations in situations like this, of course, must be determined for the particular system in question by experiment. When the result is known, the question arises, can this particular result be explained by what is known about the biochemistry of these pathways? Do the known control mechanisms account for the observed behavior? Where the system under observation is highly interactive, human intuition is unable to keep track of the tangle of interrelationships, and a mathematical model becomes most desirable. A further reason why human intuition or simple non-mathematical models become inadequate is that beyond a certain level of complexity biochemical systems may exhibit the phenomenon of 'multistability', i.e. multiple alternative stable steady states may exist, and small changes in the system parameters may lead to abrupt changes in the levels of ~i particular enzyme or substrate; thus, a small quantitative change in concentration of a compound may trigger a qualitative change in behavior of the system. Prigogine and Babloyantz (1972), for example, discuss the induction of/3-galactosidase in E. coll. This 'is an all-or-none phenomenon'. Bacterial populations grown at low concentrations of inducer in a fixed medium are composed of individuals which are in either of two possible steady states; maximally induced, or noninduced. The bacteria can be shifted from one stable state to the other by transitory changes in the environment... the observed experimental c u r v e s . . , show a sudden jump in the quantity of the enzyme, while the inducer concentration changes relatively slowly'. These authors explain the properties of this system in terms of the interactions between the galactoside substrate, the/3-galactosidase enzyme, the galactoside permease, and the lac operon repressor protein. In such a system, where not only the quantitative aspects of control, but also the qualitative behavior, is affected by small differences in
248
R. C. JACKSONand K. R. HARRAP
concentration of nutrient or metabolite, a mathematical model is necessary to analyze and predict the system behavior. (iii) A further reason why mathematical models of drug interaction are desirable arises from the combinatorial nature of drug treatments. Assume that for a particular tumor twenty agents are known which have some cytotoxic effect. Then the number of combinations of four or fewer drugs which could be selected is (20!/16! + 20!/17! + 20!/18!+20!/19!)= 123,520. This makes the assumption that the order of administration is important, which in many instances is known to be the case. If we allow for different dose-levels and varying time intervals between administration of the successive drugs, it is possible to devise tens of millions of protocols of four or fewer drugs from a choice of twenty useful agents. Since, in practice, only a minute proportion of these protocols can be tested, even in experimental systems, it is interesting to consider the logical processes by which the present combination chemotherapy protocols have arisen. For the most part, these have been statistical; if drug A gave a 40 per cent remission rate for a particular tumor, and drug B gave a 30 per cent remission rate, and if the modes of action were completely independent, then A plus B would be expected to give a remission rate of [40+ (100- 40) x (30/100)] per cent = 58 per cent. This reasoning only applies if the modes of action of the two compounds are completely independent, and in view of the highly interactive nature of cellular biochemistry, this is unlikely to be the case. The likelihood of negative interactions can to some extent be minimized by deliberately choosing agents whose mechanisms of action are apparently unrelated, and perhaps for this reason some successful combinations have been designed consisting of, for example, an alkylating agent, an antimetabolite, and a mitotic spindle inhibitor. Other rules of thumb have been suggested, such as the principle of avoiding cumulative toxicity; thus, if drug A is particularly myelotoxic, drugs used in combination with it should preferably exert their limiting host toxicity in a tissue other than the bone marrow. In this way it may be hoped that the drug effects will be additive toward the tumor, but not toward any one host organ or tissue. Thus, it may be claimed that for the most part presently used clinical combination chemotherapy regimens have been devised primarily on statistical grounds using general principles based on pharmacological and clinical commonsense and intuition. Nevertheless, for the reasons stated above, when a system passes a certain level of complexity, human i0tuition, or even explicit non-mathematical models, become tools of uncertain reliability; it is at this point that quantitative models become desirable. Such models, of course, require detailed accurate data concerning the system to be modeled if their predictions are to be of any value. (iv) The remaining reasons for using mathematical models of drug interaction are that sometimes there are, for technical reasons, no experimental means of answering a question, and a computer simulation may in such cases be the next best thing to a real experiment. In this way we can attempt to answer questions about imaginary situations. Consider an example: the antimetabolite MTX, although a potent inhibitor of dihydrofolate reductase, also inhibits thymidylate synthetase, albeit more weakly, in every system so far examined. Some authorities have considered that this secondary effect is an important component of the drug's cytotoxicity, and others have disagreed. If two tumor cell lines existed which were identical in every respect, except that one had a thymidylate synthetase which was susceptible to MTX inhibition, and the other had a MTX-resistant thymidylate synthetase, then the question could be answered directly by considering the comparative effect of the drug on the two tumor lines. Since such an experimental system does not exist this approach is not possible. However, if a mathematical model of the system is written, containing detailed information concerning the kinetic properties of the relevant enzymes, including their inhibition by MTX, the following question can be asked, 'how does removal of the inhibitory effect of MTX on thymidylate synthetase affect the cytotoxicity predicted by the model?'. This kind of question, although imaginary, is not trivial, since answers to questions of this kind provide information which may be used to help in the design of new inhibitors.
Computer models of anticancer drug interaction
249
Biochemistry, as usually practiced, is largely an analytical subject. The various components of cells are purified, and their properties in isolation are studied in great detail. This analytical approach is necessary if we are to understand the nature of drug-receptor interactions, However, the tools of biochemistry have allowed less progress with synthetic aspects. Sometimes it has been possible to combine purified cellular components into functional reconstituted systems, but this 'synthetic biochemistry' is extremely difficult and is still far behind analytical biochemistry. This is an area where computer techniques can make a contribution by examining the possible ways in which the separate components of a system may fit together and interact. 3. PHARMACOKINETIC (PHARMACODYNAMIC) MODELING OF ANTITUMOR DRUGS The invention of the modern compartmental approach to the science of pharmacokinetics is often credited to Teorell (1937). This subject, dealing with the spatial distribution of drugs in different compartments of the body, and the temporal changes in these concentrations resulting from drug uptake, movement between compartments, and elimination from the body, necessarily relies to a great extent on mathematical models, and makes extensive use of computers. For drugs which are water soluble, and slow acting, and where complex nonlinear transport processes can be ignored, it is often found that a plot of log concentration vs time approximates a straight line. The mathematics of systems of this kind were developed by KriigerThiemer (1966a, b). Such models have been used to predict optimal dosage regimens for a number of agents (Kriiger-Thiemer, 1969). Important aspects of pharmacokinetic processes were developed by Bellman et ai. (1960) and Jacquez et al. (1960) who modeled concentration gradients in capillaries. The mathematical techniques used in compartmental analysis are the subject of a detailed monograph by Jacquez (1972). An early model of MTX action, incorporating some pharmacokinetic assumptions, was written by Werkheiser (1971). This model will be discussed in Section 7. The modern approach to detailed pharmacokinetic modeling of anticancer drugs was pioneered chiefly by Bischoff, Dedrick, Zaharko and their co-workers (Bischoff, 1975; Bischoff et al., 1970, 1971; Dedrick et al., 1970, 1973b; Zaharko et al., 1971). These studies simulated the distribution and excretion of MTX, an appropriate choice for such studies because it is not metabolized in man when administered at conventional doses. The compartmental model of Bischof et al. (1971) considered MTX distribution between blood plasma, liver, kidney, muscle and intestine, secretion into bile, absorption into gut, and losses into urine and feces (Fig. 3). Bischoff (1975) in a recent review has discussed such questions as how to decide on the number of compartments to model and what the time-scale of modeling should be, as well as various technical aspects of pharmacokinetic modeling. The quantitative aspects of plasma and tissue binding were also considered. The authors quoted above have elaborated a number of important concepts such as the distinction between local uptake which is limited by drug availability to the region (flow-limited conditions), and local uptake which is determined by the membrane transport properties of the target cells (membrane-limited conditions). The pharmacokinetic consequences of these different assumptions have been worked out in considerable detail, as have those resulting from different degrees of water- and lipid-solubility (Dedrick et al., 1970). Furthermore, the distinction has been made between linear and nonlinear concentration dependence of drug-target binding, and the 'linear binding' approximation has been used in studies of MTX pharmacokinetics (Bischoff, 1975; Bischoff et al., 1970, 1971; Dedrick et al., 1970, 1973b; Zaharko et al., 1971). Dedrick et al. (1975) have extended the MTX model to include such intracellular biochemical events as synthesis of new dihydrofolate reductase molecules. This model will be further considered in Section 7. The pharmacokinetics of AraC presents a very different situation from MTX, in
250
R. C. JACKSONand K. R. HARRAP I
PLASMA
OL-OG
!,
I
/
LIVER
/
~Billllryr.,~retlL
/
G.I. TRACT
¢ ¢ eGutAbeoeptlon Gut Lumen
Ok KIDNEY
1
Urine
OM
FIG.3. Distributionof methotrexatein bodycompartments. Qc, QL, Qx, QM, plasma flowrates to gut, liver, kidney, and muscles, respectively; T, nominal residence time in bile subcompartments; rl,2,3, drug transport rate in bile; cL2.3.4,drug concentration in gut lumen compartments.(FromBischoffet al., 1971.)(Reproduced with permission of the copyright owner.) that (a) AraC is rapidly metabolized by deamination to the inactive arabinosyluracil (AraU), and (b) AraC must be phosphorylated to the triphosphate level to form AraCTP before it can exert its cytotoxic action. Dedrick et al. (1972, 1973a) have presented a pharmacokinetic model for the detailed distribution of AraC and AraU in man. This model employs known enzyme kinetic parameters to describe deamination of AraC to AraU; it does not, however, include the fact that AraC must be converted by phosphorylation to AraCTP. These studies have brought to light a number of interesting aspects of the pharmacology of AraC. Using kinetic parameters for deaminases which had been measured in vitro, the model successfully predicted plasma concentrations of AraC and AraU following injections of AraC; in addition, the important role of 'lean tissue' in acting as a reservoir of AraC and releasing it to the blood was demonstrated. Deamination in human liver was shown to be so rapid that a large fraction of the AraC was inactivated in a single passage through this tissue; the observed kinetics were thus more strongly dependent upon blood flow rate than on the detailed enzyme kinetic properties in the various tissues (Dedrick et al., 1972). Good agreement was obtained between theory and experimental results when models incorporating enzyme parameters for mice, monkeys, dogs and humans were compared (Dedrick et al., 1973a). Lincoln and Freireich (1974) discussed a computer simulation of AraC kinetics; using an interactive program, with display of computer graphics, they simulated the distribution of this drug in various body compartments (i.e. blood, bone marrow, gut and testis) in man. This work made use of a biological modeling program developed at RAND Corporation (Groner and Clark, 1971). Morrison et al. (1975) have recently described a model for pharmacokinetic analysis of AraC which includes some of the intracellular reactions involved in the production and breakdown of AraCTP. This model will be discussed in Section 7. Leme et al.
Computer models of anticancer drug interaction
251
(1975) devised a compartmental pharmacokinetic model for the disposition of moderate- and high-dose MTX in man. When tritiated MTX was administered in large doses, two labeled compartments appeared in the plasma. One of these was MTX, and the other, which was not exchangeable with MTX, was suggested by these authors to be the result of biliary excretion ~of MTX, degradation of the drug by the intestinal flora and subsequent resorption of these degradation products. It was found that a model with nonlinear inter-compartmental transfer gave better agreement with observed results than a simple linear compartmental model. Recent developments in pharmacokinetic modeling have extended the BischoffDedrick model to newer antitumor agents. Harris and Gross (1975) presented a preliminary pharmacokinetic model for adriamycin. This model has ten compartments, including heart and bone marrow, the dose limiting sites of host toxicity for adriamycin. The major site of drug clearance was considered to be the liver, with renal clearance being negligible. A constant S0 per cent of the adriamycin present in the plasma was assumed to be protein-bound. Tissue-plasma concentration ratios were measured experimentally for rabbit tissues at three time points following intravenous drug administration, and the mean values, assumed to be time-independent, were used in the simulation. The model was used to predict adriamycin concentrations in the ten tissues at times of up to 50hr after intravenous administration of adriamycin to rabbits. The agreement between predicted and measured tissue levels was good. The model was then used to predict human plasma levels, which again gave good agreement with measured values. Allen and Creaven (1975) have described a human pharmacokinetic model for isophosphamide. This agent, like the closely related drug cyclophosphamide, must be metabolically activated. This model is discussed in Section 4. Other recent developments in the techniques of pharmacokinetic modeling allow the systems to accommodate multiple unequal drug doses (Niebergall et al., 1974). A further extension of such models allows nonuniform time spacing of the doses (Howell, 1975). The preceding models have been 'deterministic', i.e. they assume that the rate constants and other system parameters are known with precision, and are constant within a particular compartment or subcompartment of the system. Real systems, in fact, are so complex that certain of their properties appear to exhibit random behavior. Some authors, consequently, have suggested that a more realistic treatment of these systems would be to use stochastic, rather than deterministic, assumptions at least in parts of the models. In such systems, the rate constants or other system parameters are treated as random variables. Aspects of stochastic compartmental modeling were discussed by Matis and Hartley (1971), Thakur et al. (1973) and Purdue (1974). Padgett and Tsokos (1972) have discussed the theory of a stochastic model for chemotherapy for a two-organ biological system, and have described computer simulations for a one-organ model (Padgett and Tsokos, 1970). Chuang and Lloyd (1974) have discussed the theory of two-compartment and three-compartment models. These authors conclude: 'A mean of the concentration process is obtained which in effect is a best-fit mean concentration-time ( C - T ) curve for the measured compartment of the group of individuals under observation. It is not possible to predict the C - T curve for any particular individual in a population from the stochastic analysis, but it is possible to predict with some degree of confidence the expected range for tl~e C - T curve. In addition, expected stochastic C - T curves for the unobserved or unobservable compartments may be predicted'. McInnis et al. (1975) described a two-compartment stochastic model of the disposition of the anticancer drug 5-(3,3dimethyl-l-triazeno)imidazole-4-carboxamide (NSC-45388); the predicted disappearance of NSC-45388 from plasma in a dog was compared with experimental data. All of the pharmacokinetic models discussed above have been attempts to describe the compartmental and temporal changes in concentration of a single drug, or of one drug and its metabolites. In general, drug-drug interactions have not been considered, although it is possible that some of these models could be modified to take such interactions into account. A few theoretical approaches to the pharmacokinetics
252
R.C. JACKSONand K. R. HARRAP
of drug combinations have appeared. Buell et al. (1970) have shown that the concepts of control theory may be applied to the computation of multidrug dosage schedules. The problem of drug administration was treated as a discrete-time optimal control problem, and methods for solving this problem by dynamic programing were shown. Preliminary pharmacokinetic studies of drug-drug interactions have appeared, although these do not refer to antitumor drugs (Prescott, 1973). Rowland and Matin (1973) discuss a number of drug combinations for which a careful kinetic analysis has clarified the site of drug--drug interaction. Thus the plasma half-life of the oral hypoglycemic agent tolbutamide is considerably lengthened by sulfaphenazole; this is accomplished by blockage by the latter agent of the hepatic microsomal oxidation of the tolbutamide. These authors also discussed the interaction between the anticoagulant, warfarin, and phenylbutazone. Warfarin is bound to a large extent to plasma protein, but it is, apparently, the unbound fraction which exerts the pharmacological effect. Phenylbutazone displaces bound warfarin and thus augments the anticoagulant effect. In addition, by selective inhibition, it alters the isomeric composition of the warfarin (administered as a racemic mixture) in plasma, hastening the metabolism of the less active isomer, and delaying that of the active form. Such studies have not been carried out with antitumor drugs. However, in view of the importance of such drug interactions, and of the relatively great advances in antitumor drug pharmacokinetics in recent years, it may be hoped that such studies will be conducted in the very near future. 4. MODELS INVOLVING DRUG METABOLISM AT A SITE DISTANT FROM THE TUMOR It is convenient to distinguish between (a) drugs which require no metabolic activation, such as MTX, (b) drugs which are activated in the target cell, such as AraC, and (c) drugs which are activated at some site distant from the tumor, usually the host liver, such as cyclophosphamide. The first class of drugs is the most amenable to pharmacokinetic modeling, as the various successful models of MTX pharmacokinetics have clearly shown. The second class to be modeled in detail requires a knowledge of the intracellular reactions involved in drug activation. Many of the drugs in this category are purine and pyrimidine nucleoside analogs, which exert cytotoxicity at the nucleotide level. The nucleotide forms of these drugs are unable to effectively cross the plasma membrane, and must therefore be administered in forms which are activated in the cell in which they are to act. Models involving anticancer agents of this type are discussed in Sections 5 and 7. Drugs which must be metabolically activated in the liver, or in some other body site, present certain complications to experimental study. For example, these agents cannot be tested simply against tumor cells in culture; this problem has been circumvented in the case of cyclophosphamide by adding isolated hepatic microsomes to the cultures. These complications extend to mathematical modeling. A satisfactory model of a drug of this kind must simulate the disposition of both the original compound, and the metabolically active form in the body, as well as the activation process itself. Such models, therefore, inevitably combine both pharmacokinetic and intracellular biochemical features. What appears to be the first model of this kind to describe the behavior of an anticancer drug, cyclophosphamide, was presented by Cohen et al. (1971). Biotransformation (activation) was assumed to take place within the large central compartment, 'since the relatively large volume of distribution is likely to encompass the liver where oxidation of cyclophosphamide occurs'. Metabolic activation was not modeled in detail, but was treated as first-order with respect to the concentration of cyclophosphamide. This would approximate the situation if the concentration was always well below the Km of the activation process. The model was used to estimate parameters for a total of eleven doses given to seven patients. The authors noted that barbiturates enhanced the activity of the microsomal drug metabolizing enzymes of the liver. They did not explicitly model this activation, but they noted that one of their
Computer modelsof anticancerdrug interaction
253
patients who had received phenobarbital appeared to have metabolized a larger fraction of the administered dose than the other patients. A more recent and detailed model of this kind is that of Allen and Creaven (1975). This study modeled the pharmacokinetics of isophosphamide (NSC-109724), an analog of cyclophosphamide in which one of the /3-chloroethyl groups is replaced by hydrogen, and is transferred to the 3-position of the oxazaphosphorine ring. Isophosphamide is active against some tumors which are resistant to cyclophosphamide, and also differs from cyclophosphamide in that the host toxicity is directed principally at the kidney, rather than the bone-marrow (Allen and Creaven, 1975). Like cyclophosphamide, isophosphamide requires activation by the hepatic endoplasmic reticulum. The model of Allen and Creaven is primarily a pharmacokinetic model, but includes two parameters (K,, and Vm~) to describe the drug activation (Fig. 4). This model indicated that the excretion rate constant for drug metabolites was larger than for the unchanged drug. Furthermore, the activated drug appeared to distribute through a much smaller volume than unchanged drug, indicating limited permeability into many parts of the body. The authors report that the model is currently being used to predict plasma levels of isophosphamide after repeat doses of the drug, and also to assess the influence of other drugs on the pharmacokinetics of isophosphamide, though these studies are not yet published. Perhaps the most complete model of drug-drug interactions involving drug metabolism is the study of Rowland and Matin (1973), already mentioned in Section 3. In their detailed model of the interaction of sulfaphenazole with tolbutamide they successfully simulated the effect of sulfaphenazole in increasing the plasma half-life of tolbutamide from its normal value of approximately 7 hr to about 44 hr. This situation differs from that with cyclophosphamide in that tolbutamide is itself active, and the oxidized drug formed in the hepatic endoplasmic reticulum is inactive. This interesting model will not be discussed further here, since it was not concerned with antitumor drugs. Nevertheless, it is likely that similar studies of this kind, on the pharmacokinetics of drug-drug interactions at the level of hepatic activation and inactivation of anticancer drugs, will appear before too long. •
x----O
DOSE
@
FIG. 4. Multicompartmentpharmacokineticmodel for isophosphamide.X~ and X4 are central compartments;X3 and X6, urine; X2 and X~, peripheral compartmentsfor isophosphamideand metabolite(s), respectively; V.... k~e and kE~ are compartmentsfor metabolism,where E is the substrate-enzymecomplex; kx is the transfer rate constant (hr-~) for irreversible biologic alkylation or drug decomposition; kii (i = 1, 2, 3 . . . . n, j = 1, 2, 3 . . . . . m ) are transfer rate constants (hr-~) for drug transport into compartment j from compartment i; V~ and V~, are volumes of the central compartmentfor isophosphamideand metabolite(s)respectively.(From Allen and Creaven, 1975.)
5. INTRACELLULAR BIOCHEMICAL INTERACTIONS In distinguishing between interactions at the cellular level, and those occurring at a site distant from the tumor, as discussed in the previous section, we are following the distinction made by Grindey e t a l . (1975) between 'intracellular metabolic interactions' and interactions at the level of 'remote sy,stemic biotransformation'. A few rather detailed computer models of intracellular biochemical interactions of anticancer drugs
R. C. JACKSONand K. R. HARRAP
254
have been published. Study of these interactions is facilitated by the ability to conduct experiments with cell cultures, where extraneous variables such as extracellular drug concentration, and nutritional status may be much more rigorously controlled than in animals. Thus, the predictions of models of cellular biochemistry may be subjected to detailed testing. Mathematical models at the intracellular level generally are collections of rate equations describing the interconversion of substrates and cofactors in enzyme-catalyzed reactions. Such models are thus inherently nonlinear. Detailed computer modeling of biochemical control mechanisms in multienzyme pathways was first done by Garfinkel, Hess," Chance and their collaborators (Chance et al., 1960; Garfinkel et al., 1961; Garfinkel and Hess, 1964). These workers have described a detailed model of the reactions of glycolysis in ascites cells. This model accurately predicted the effects of such perturbations as addition of glucose or of the uncoupling agent dibromophenol on the concentrations of glycolytic intermediates. With models of this kind, which are intended to simulate collections of chemical reactions, it is axiomatic that a detailed knowledge is required of the relevant reactions, their stoichiometry, and the rate equations. In the case of enzyme-catalyzed reactions it is necessary to know the kinetic mechanism of the reaction to write the rate equation. These principles were first applied to the study of the intracellular interactions of antitumor drugs by Werkheiser, Grindey, Moran and Nichol (Werkheiser, 1968, 1971; Grindey et al., 1970, 1975; Nichol et al., 1972; Werkheiser et al., 1973). These papers have presented a number of different models; one of them contained features of both biochemical and pharmacokinetic modeling, and this will be discussed in Section 7. The model of DNA synthesis was an intracellular model; its features are summarized in Fig. 5. The assumptions underlying this simulation were as follows: (i) the system being modeled was the L1210 mouse leukemia cell growing in culture; (ii) the cell was treated as an open system in a steady state, with the medium acting as a constant source of precursors, and the growing DNA chain as a sink; (iii) within this simplified system the transfer of material from the source to the sink was as described in Fig. 5, the rates being determined by enzymes which obeyed unidirectional MichaelisMenten kinetics; (iv) 5-FdUMP, AraCTP and AraATP were treated as competitive inhibitors of their respective target reactions, and IQ-1 (1-formylisoquinoline thiosemicarbazone) and MTX were treated as noncompetitive inhibitors of ribonucleotide reductase and thymidylate synthetase, respectively. The details of the rate equations may be found in the original paper (Werkheiser et al., 1973). Using this model, the
RC
i .....!'" " "~ / •"""
i
/
RU ! ...................
D
:
\
....................................................... i
RA,
RG
:
i
.,,iJ .,~I
/
m, A
I
i
\ ,s
J
Fro. 5. Model of metabolic regulation of DNA synthesis. RC, CDP; RU, UDP; ] ~ , ADP; RG, GDP; C, dCTP; U, dUMP; T, dTTP; A, dATP; F, dGTP; D, DNA; feedback inhibition, dotted lines; feedback activation, dashed lines. (From Grindey et al.. ]975.)
Computer models of anticancer drug interaction
255
Roswell Park group simulated the interactions of binary combinations of M T X , AraC, AraA, IQ-I and FUdR. Some of the results plotted as isobolograms are shown in Fig. 6, which compares the computer predictions with the observed effects of the drug pairs on L1210 cell growth in vitro. A number of important conclusions emerged from the various studies conducted with this model: (i) in general, in pathways with multiple feedback inhibitions,the behavior of the pathway when subjected to multiple inhibition could not be predicted from the sequential, concurrent or complementary relationship of the inhibitors, but was a function of the structure of the system as a whole; (ii)the types of drug interactions were relatively independent of the system parameters, and primarily depended upon the structure of the model. Occasional exceptions were noted to this generalization, and these will be discussed further below; (iii)the model was highly interactive, i.e.the presence of an inhibitor caused changes in the concentrations of all the intermediates. The effects of a perturbation in one region were not localized, but were manifested throughout the entire system; (iv) the feedback regulation led to results which could not have been intuitively foreseen• Inhibition of a pathway could result in increased, decreased, or unchanged concentration of the product of that pathway; and (v) the intensity of interaction between two drugs was related to the degree of inhibition of the system, A drug combination which was mildly synergistic or antagonistic at 50 per cent inhibition became intensely so when the system was inhibited by 90 per cent. The model of D N A synthesis was used to predict the response of deoxyribonucleoside triphosphate pools to antimetabolite chemotherapy. Some of the results were those which would have been intuitively predicted; thus, IQ-I, the ribonucleotide reductase inhibitor, caused decreased levels of all deoxyribonucleoside triphosphates. However, many of the predictions could not have boon made without the model: for example, AraA and AraC were both predicted to lead to increases in dCTP, d A T P and dTTP, and decreases in uridine nucleotides. In some cases the predictions describe transient behavior as well as steady-state solutions. The effect of MTX, for example, was predicted to cause an increase in uridine and cytidine nucleotides, with, initially, a sudden fall in d T r P and a rise in dATP; the feedback properties of the model were 1.0~
IX'.
°1 -
1.0~
1.0~
I",,
0.5 10-1 R~s
~.0
•
~
IX:"
,\J
0.5 ARA-C Cc
1 , 0 ~
1.0
•
~
1.
0.~ ARA-C Cc
1.0
Oo •
go o:s
1".o ~
IQ-1 Rn
_
II." . G'
0.5
i.o
ARA-A Ac
= 0~5 10-1 Rn
~'1.O
el
o'.s
1.o
ARA-C C©
/ G"
O15 ARA-A AC
I~O
C"
0.5 ARA-C CC
1.0
FIG. 6. Isobolograms showing interactions of pairs of dr~]gs on the growth of LI210 cells in vitro (points) compared with the effects on DNA synthesis as predicted by the model (solid lines). (From Werkheiser etal., 1973.)
256
R. C. JACKSON and K. R. HARRAP
such, however, that the dTTP level later returned to the normal value, with a concomitant drop in dATP, which became the limiting nucleotide in the steady state. Werkheiser et al. (1973) conclude: 'No claim is made that this model represents the real system in any of its details, but it seems likely that the principles of its construction are reasonable'. Another model described by the same authors (Grindey et al., 1975) was used for a simple simulation of folate metabolism. This model treated dihydrofolate reductase and thymidylate synthetase as separate reactions, and was therefore able to model more accurately the interactions of MTX and 5-FdUMP. A limitation of this model was that simple equations for competitive inhibition were used; in these instances, the rate equations for tight-binding competitive inhibition would have been more appropriate. This model also allowed for simulation of folate synthesis, and the inhibition of that process; thus, the interactions of sulfonamides with MTX or 5FdUMP could be studied. Mammalian cells, unlike many microorganisms, are unable to synthesize folate; however, the model could be used to simulate uptake of preformed folate from the medium. The model predicted that an inhibitor of dihydrofolate synthesis or uptake would be synergistic with an inhibitor of dihydrofolate reductase. It also predicted that MTX and FUdR would be antagonistic, and the combination of sulfonamides with FUdR would be additive. These predictions are consistent with the observed biological behavior of such combinations. These models of Werkheiser et al. in general represent minimal models of the systems they attempt to simulate. That is, they are an attempt to produce the simplest possible mathematical description of the systems under study which still retains the essential control features of the parent systems. These models are thus far less sophisticated than those of, for example, Garfinkel, in which scores of individual rate processes are simultaneously studied. However, especially where the details of a system are not well understood from the experimental data, simplicity is often desirable in a model. These particular simulations have realized the hopes of their authors in accounting for many of the observed interactions between anticancer drugs, and in suggesting profitable experiments. They conclude: 'Simulations that do not lead to new experiments in the laboratory have little value; there should be a continuing interplay between model and experiment'. The introduction of these particular models represented a major advance in the study of the biochemical pharmacology of antitumor drugs. An intracellular biochemical model of some reactions of folate and purine nucleotide metabolism has been described by Jackson and Harrap (1973). This model is limited in range, in that it only describes a small area of cellular metabolism (Fig. 7), but the part described is given in detail. The individual rate equations for 23 reactions are simultaneously integrated. When kinetic mechanisms had been studied in sufficient detail, the appropriate rate equations were used. Where the kinetic mechanism was unknown, it was treated as rapid equilibrium random. For tight-binding inhibitors such as aminopterin, the Goldstein equation was used. This model can therefore be employed to study the effects of inhibitors on steady-state fluxes and pool sizes, and also examine the transient phase following addition of inhibitors or precursors. As with the models discussed above, it is assumed that the nutrient medium acts as a source of precursors, and the exponentially growing cells act as a sink. The uninhibited system is in open steady state; following perturbation a transient phase occurs which may or may not lead to another steady state. The system (as described in Fig. 7) possessed a considerable degree of inherent stability. For example, if dihydrofolate reductase was partially inhibited, the substrate dihydrofolate accumulated until the increased substrate concentration exactly compensated for the degree of inhibition, after which the steady-state rate remained unchanged. If, however, the dihydrofolate reductase was inhibited severely enough to make it the rate-limiting enzyme in the cycle of dihydrofolate oxidation and reduction, then the system assumed a new steady state, in which the flux through the cycle was diminished. Some of the observed effects were explainable on the basis of existing knowledge of the control of multienzyme systems (Webb, 1963; Reiner, 1969; Westley, 1969). Other effects were
Computer models of anticancer drug interaction
.._JJ
.......
257
A,cA\\
1 FIG. 7. Diagram of reactions simulated by a mathematical model. AICAR, 5-amino4-imidazole carboxamide ribotide; CH:FIt4, NS,10-methenyltetrahydrofolate; CH2Ftt,, N5,10methylenetetrahydrofolate; 5-FdUMP, 5-fluorodeoxyuddylate; FI-I2, dihydrofolate; FH4, tetrahydrofolate; GAR, glycinamide ribonucleotide; 10-HCOFH4. N 10-formyltetrahydrofolate; MTX--methotrexate. The rate processes simulated are as follows: vl, dihydrofolate reductase (EC1.5.1.3): v2, serine hydroxymethyltransferase (EC 2.1.2.1); v3, thymidylate synthetase (EC 2.1.1.4.5); v4, serine hydroxymethyltransferase (back reaction); vs, NS,10-methylenetetrahydrofolate dehydrogenase (EC 1.5.1.5); v6, NS,10-methenyltetrahydrofolate cyciohydrolase (EC 3.5.4.9); v~, NS,10-methenyltetrahydrofolate cyclohydrolase (back reaction); 08, NS,10-methenyltetrahydrofolate:glyeinamide ribonucleotide transformylase (EC 2.1.2.2); vg, Nl04ormyltetrahydrofolate:5-amino-4-imidazole carboxamide ribonucleotide transformylase (EC 2.1.2.3); v~0,carrier-mediated cellular uptake of MTX; v., nonenzymic combination of tetrahydrofolate and formaldehyde; v~2, nonenzymic decomposition of NS,10-methylenetetrahydrofolate; v~3, cellular uptake of MTX by passive diffusion; v~4, efllux of MTX by passive diffusion; vl~, biosynthesis of serine; v,6, serine dehydratase (EC 4.2.1.13); v~0,carrier-mediated efflux of MTX; v2,, biosynthesis of deoxyuridylate; v~, biosynthesis of glycinamide ribonucleotide. Rate processes considered by the model but not shown in Fig. 7 are: v~7,biosynthesis of dihydrofolate reductase; v~8,cellular degradation of dihydrofolate reductase; v,9, cellular degradation of dihydrofolate reductaseMTX complex. (From Jackson and Harrap, 1973.) not predictable f r o m the behavior of earlier, simpler models of sequential or cyclic pathways, and a p p e a r e d to be a consequence of the complex t o p o g r a p h y of the model, which contains three overlapping cycles and m a n y branch points. Initially, simulations were run to estimate the effect of changes in p a r a m e t e r s which could determine resistance to MTX. The model, like real L1210 cells, b e c a m e less sensitive to M T X if m e m b r a n e permeability was altered. If the change in m e m b r a n e permeability was the result of an increase in the Kt value of the transport protein, the reduced permeability could be o v e r c o m e by higher concentrations of extracellular drug, but if the change was in the Vmax for m e m b r a n e transport of drug, the mutation remained effective e v e n at v e r y high extracellular levels of drug. Both types of mutation have b e e n found in cultured cells (Jackson et al., 1975; N i e t h a m m e r and Jackson, 1975). An increase in the cellular dihydrofolate reductase concentration in the model resulted in M T X resistance, as is so often o b s e r v e d experimentally. The model simulated the synthesis and turnover of dihydrofolate reductase, and either an increased rate of synthesis or reduced rate of degradation of e n z y m e increased the steady-state level. Since M T X appears to stabilize dihydrofolate reductase against turnover, it increased the cellular dihydrofolate reductase concentration. Of course, the additional e n z y m e with M T X bound to it was inactive, but the inhibited e n z y m e was in equilibrium (or close to it) with free e n z y m e and free drug, and the model showed that under some circumstances the p r e s e n c e of M T X could actually increase the concentration of free (active) dihydrofolate reductase in the cell. Simulation of the effects of M T X on the
R. C. JACKSON and K. R. HARRAP
258
activity of serine hydroxymethyltransferase and thymidylate synthetase showed that normally inhibition of these enzymes was not a major contributor to the effect of the drug. More recent studies (Jackson and Niethammer, 1977) indicate that in extremely resistant cells, in which the dihydrofolate reductase concentration is increased 300fold or more, free MTX may accumulate to sufficient levels to make thymidylate synthetase inhibition the rate-limiting component of MTX action. When this state is reached, leucovorin no longer completely rescues the cells, but only partially protects them from MTX. This surprising result, predicted first by the computer model, has been confirmed with high dihydrofolate reductase containing human cells in culture; the partial protection by leucovorin is due to its conversion to 5,10-methylenetetrahydrofolate which then competes with MTX for the folate binding site on thymidylate synthetase. Since these mutant lines with huge reductase contents have never been described in vivo, the practical importance of this effect is unknown. Another factor predicted to alter cellular response to MTX was the dissociation constant of the cell's dihydrofolate reductase for the inhibitor. This result was also unexpected, since this dissociation constant is in all cases very low, so that the inhibition is often described as 'stoichiometric'. When the titration curves for such inhibitors are plotted, the early part of the curve corresponds to 1:1 inactivation of enzyme molecules by inhibitor molecules, regardless of the exact Ki value. Thus, it seemed that the exact value of the latter parameter was unimportant, so long as it was sufficiently low to keep inhibition 'stoichiometric'. Simulation showed the fallacy in this reasoning. Because dihydrofolate reductase was present in great excess over thymidylate synthetase, the rate-limiting enzyme in the 'dihydrofolate cycle', it was necessary to inhibit it by 95 per cent or more in order to make it rate-limiting, and thus result in MTX-induced growth inhibition. However, at 95 per cent inhibition, the enzyme-inhibitor titration curves are no longer linear, but instead show a 'toe'. In this region the exact value of the Ki makes a great deal of difference to the observed reaction rate. Recent studies (Jackson et al., 1976) have confirmed in four cultured mammalian cell lines (with similar total activities of dihydrofolate reductase, similar rates of drug transport, similar total folate cofactor pools, and similar dTTP levels) that the sensitivity to MTX correlated inversely with Ki value for the respective dihydrofolate reductases and MTX. Table 1 compares the computer predictions of the inhibitory concentrations of the four cell lines, and their predicted relative resistance, with the observed values. The predictions were based on measured values for the relevant enzyme activities and kinetic parameters, the appropriate rate equations, and the assumption that a cell would divide when its content of both d'ITP and purine nucleotides had doubled. Considering the uncertainties inherent in this kind of prediction, the agreement is excellent. Here again, experiments confirmed the predictions of the model. A simpler model, in which inhibition of dihydrofolate reductase by MTX was treated as stoichiometric, could not, of course, have made this prediction. As well as being used to examine the effect of system parameters on response to single drugs, the model was used to predict drug interactions. The drugs simulated were antifolates, thymidylate synthetase inhibitors, glutamine analogs, and purine TABLE 1. Concentrations of M T X Required to Inhibit Cell Growth by 50 per cent as Predicted by Metabolic Simulation. The Appropriate Kinetic Parameters for Dihydrofolate Reductase were Used, but other Factors were Assumed to be Constant
Cell line
K~ (pM)
LI210 W1-L2 L5178Y Yoshida
5 7 29 190
Predicted Resistance* 1I)5o (riM) (-fold) 14.2 18.3 72.6 409
(1) 1.3 5.1 28.8
Observed Resistance* IDso (riM) (-fold) 5.8 13 56 390
(1) 2.2 9.7 67.2
*Resistance is relative to the predicted or observed value for LI210.
Computer models of anticancer drug interaction
259
antimetabolites. When the combination of MTX + FUdR was simulated, the results showed some interesting similarities with the predictions of the model of Grindey et al. (1975) for the same combination, and an important difference. When the kinetic parameters for L1210 mouse leukemia cells were used in the model, it predicted that the combination would be antagonistic, in agreement with the model of Grindey and his co-workers. However, when the parameters for Yoshida sarcoma cells were substituted, the model predicted that the two agents would be almost additive. The reasons for this difference formed the basis for further study with the model, and for a series of experimental studies. If the dihydrofolate oxidation-reduction cycle is considered alone, it is clear the MTX and 5-FdUMP cannot be additive, since one or another of the two enzymes dihydrofolate reductase and thymidylate synthetase will be rate-limiting, and inhibition of the non-rate-limiting enzyme will have little effect on the flux through the cycle. This was the prediction of the model of Grindey et al., which did not take into account the inhibitory effects of MTX on purine nucleotide synthesis. The model of Jackson and Harrap (1973), which simulated the involvement of folates in the pathway of de novo purine nucleotide biosynthesis made a similar prediction for L1210 leukemia cells. Moreover, it predicted that in the L1210 cell, the growth inhibitory effect of MTX was primarily a consequence of a limitation in the availability of purine nucleotides. In the slower-growing Yoshida sarcoma cell, the MTX effect was primarily the consequence of decreased thymidylate production. The balance between the anti-thymidylate effect of MTX and the antipurine nucleotide effect was primarily determined by the competition for 5,10-methylenetetrahydrofolate b y the enzymes thymidylate synthetase and 5,10-methylenetetrahydrofolate dehydrogenase. Thymidylate synthetase requires substrate amounts of the tetrahydrofolate cofactor, since this material acts not only as a one-carbon group source, but also as a reducing agent. The folate cofactor is released from the synthetase as dihydrofolate, which must then be reduced by dihydrofolate reductase, the target enzyme of MTX. 5,10-Methylenetetrahydrofolate dehydrogenase, on the other hand, leads to purine nucleotide biosynthesis, a process in which the two folate-dependent transformylase enzymes utilize the one-carbon groups bound to the tetrahydrofolate cofactors, and release the cofactors as free tetrahydrofolate, which does not require reduction; thus, purine nucleotide biosynthesis requires catalytic amounts of tetrahydrofolate. A consequence of this phenomenon is that the rate at which tetrahydrofolates are depleted in MTX-treated cells is dependent upon the competition for 5,10-methylenetetrahydrofolate by the two enzymes. In rapidly dividing cells, such as L1210 leukemia cells, the ratio of thymidylate synthetase to 5,10-methylenetetrahydrofolate dehydrogenase is higher than in slower-growing cells (Tattersall et al., 1974; Jackson and Weber, 1976). Thus, the model predicted that the ratio of these two enzymes would determine whether the effect of MTX was exerted primarily on the biosynthesis of purine nucleotides or of thymidylate. Experimental work supported this prediction. Thus, in the LS178Y lymphoma, where the ratio of synthetase to dehydrogenase was high, MTX caused decreases in both dTTP and dGTP levels, and in the presence of added thymidine, MTX-treated cells died of purine nucleotide starvation. On the other hand, in the Yoshida sarcoma, where the ratio of synthetase to dehydrogenase was lower, MTX toxicity was partly reversed by thymidine, and only dTTP was decreased among the deoxyribonucleoside triphosphate pools (Tattersall et al., 1974). This discovery had implications for combination chemotherapy, since 5-FdUMP, as an inhibitor of thymidylate synthetase, altered the effective ratio of the enzymes competing for 5,10-methylenetetrahydrofolate, and thus antagonized the antipurine effect of MTX. Thus, the antagonism between the two drugs was most marked in those cells in which the antipurine effect of MTX formed a significant part of its mechanism of action, i.e. in rapidly growing cells. In slowergrowing cells, where MTX functions primarily as an anti-thymidylate agent, the antagonism between MTX and 5-FdUMP is less pronounced. This prediction was recently confirmed in a series of rat hepatoma cell lines in vitro, where MTX and FUdR were strongly antagonistic in the rapidly growing Novikoff tumor, and only
260
R.C. JACKSONand K. R. HARRAP
slightly less than additive in the relatively slowly growing hepatoma 8999S. The model indicated that the interactions between MTX and 5-FU derivatives were more complex than that of solely a 5-FdUMP-antagonizing effect on the antipurine activity of MTX. MTX administration resulted in a very large cellular accumulation of dUMP which competed with 5-FdUMP for binding to thymidylate synthetase. This predicted accumulation of dUMP in the presence of MTX was confirmed by experiment (Tattersall et al., 1973). Thus, not only does FdUMP antagonize the antipurine effect of MTX, but MTX antagonizes the binding of FdUMP to thymidylate synthetase. Furthermore, the sequential sites of action of the two agents leads to the prediction that this combination would be infra-additive. Finally, it has been shown (Santi et al., 1974; Danenberg et al., 1974) that binding of dUMP or FdUMP to thymidylate synthetase is dependent upon the prior binding of 5,10-methylenetetrahydrofolate; when this folate cofactor is depleted, as in the presence of MTX, FdUMP is unable to bind to the synthetase molecule. In this situation the normal substrate dUMP is also unable to bind, and thus the effect would be a neutral one, were it not for the irreversible nature of the FdUMP inhibition. This means that FdUMP permanently reduces the concentration of active enzyme molecules, but requires the presence of the tetrahydrofolate cofactor to do so. MTX, then, interacts with FdUMP in at least four ways, but in all four cases the effect works in the same direction, tending to cause mutual antagonism, to an extent which is most severe in cells containing the enzyme patterns characteristic of rapid proliferation. Drug--drug interactions are even harder to predict when the multiple effects work in opposing directions. The folate and purine model was used to unravel multiple effects of this kind in the interactions between 6-mercaptopurine (6-MP) and MTX or FUdR. 6-MP, after conversion to its ribonucleotide, thioinosinate, inhibits IMP dehydrogenase, the rate-limiting enzyme of guanine nucleotide biosynthesis, and, more weakly, adenylosuccinate synthetase, which is the first enzyme of the pathway leading to the adenine nucleotides. These effects lead to a depression of the ribonucleotide pools, particularly the guanylate pools. GMP and AMP are feedback inhibitors of the first enzyme of the de n o v o biosynthetic pathway of purine nucleotides, PRPP amidotransferase. Thus, a decrease in GMP concentration tends to activate de n o v o purine nucleotide biosynthesis, and antagonize the antipurine effect of MTX. On the other hand, thioinosinate itself is a pseudo, feedback inhibitor of PRPP amidotransferase, so this action opposes the first effect. Which action predominates can readily be determined with a computer model, since instructions can be given to ignore the inhibition of PR PP amidotransferase by thioinosinate. The prediction by the model for the interaction of MTX and 6-MP is shown in isobol form in Fig. 8, together with measured results for L1210 cells in culture; both predicted and observed results correspond to therapeutic summation. Similar results were obtained for the combina1.C
"-...... ~ \
0.~
"%...
0.6
"'"""
x~ 0.4
\ ,.,,,,X, \ "'".i \
0.2
I
0.2
I
0.4
!
0.6
I
0.8
M
1.0
6-MP
FIG. 8. Responseisobologramfor combinationsof methotrexate(MTX)and 6-mercaptopurine (6-MP) acting on in vitro growth of LI210 cells in the presence of 2.5 × 10-SM thymidine (TdR). Solid line, experimental results; dashed line, kinetic prediction. (From Harrap and Jackson, 1975.)
Computer modelsof anticancerdrug interaction
261
tion of FUdR plus 6-MP; approximately additive effects were predicted and found with the cultured L1210 cell system. Figure 9 shows a simulation of a system subjected to inhibition by the combination of azaserine plus thioinosinate (the active form of 6-MP). The dotted line shows the effect of azaserine alone; an initial very pronounced inhibition is seen, then the degree of inhibition decreases as substrates accumulating behind the metabolic block tend to overcome the inhibition. In the third phase of the azaserine effect, the degree of blockade again increases as the slow irreversible phase of the azaserine inhibition of formylglycinamidine synthesis becomes significant. The effect of thioinosinate alone is seen in the solid line. The initial inhibitory effect is almost instantaneous, and it is only subsequently reversed to a slight extent, since the substrates of this reaction, NAD + and IMP both have alternative routes of utilization which result in the metabolism of a greater proportion of these substrates than does IMP dehydrogenase, the principal target enzyme of thioinosinate. Thus, IMP and NAD + do not accumulate, and the system soon reaches a new inhibited steady state. The dashed line shows the effects of azaserine plus thioinosinate in combination. The combined inhibitory effect is not additive, but it is appreciably greater than the inhibition caused by either agent alone. The properties of the model in this instance were very much what would have been predicted for two inhibitors acting sequentially at sites which straddle a divergent metabolic branch point (Harrap and Jackson, 1975). 10C %'..
8C - \ ": • .'" ........ .. ... .-...... ...-'" . . . .
~ ~:..
4C
Azaserine ( 5 . M )
..-'" ... '..
.-,""
ul
TIMP ~
4~)
q~O
~"
~
~-
I(~X)
(2~M)
Azaserine + T I M P
' 120
TIME (sec)
FIG. 9. Simulated effects of a binary combination of a z a s e r i n e and thioinosinic acid (TIMP) on the production of xanthosine m o n o p h o s p h a t e (XMP). (Jackson R. C. (1977) unpublished observations.)
To summarize, the conclusions which may be drawn from the modeling studies of Jackson and Harrap, and those of Werkheiser, Grindey and their colleagues are generally in agreement. Sometimes the interactions of two inhibitors appear to be predictable from the concurrent, sequential, or complementary nature of their sites of action, as was the case for the azaserine plus thioinosinate combination discussed above. In general, however, it cannot be assumed that this will necessarily be the case. Where the metabolic pathways involved contain multiple cycles or feedback inhibitions (and all the pathways where antitumor drugs act are characterized by such features), then the nature of drug interactions cannot be intuitively predicted. Thus, to reason that two or more drugs will be synergistic, additive, or antagonistic on the basis of the relative locations of their sites of action (even assuming single sites of action, which seldom occurs) is futile. Whether a mathematical model will be useful in this situation depends directly upon the extent to which the model reproduces the essential control features of the system under study. There can be no arbitrary rule about how complex a model must be before it becomes a realistic one, but considering the immense complexity of the regulation of the biochemical pathways upon which anticancer drugs act, it seems likely that a model which can accurately predict drug-drug interactions in most cases is likely to be considerably more complex than those in existence at present. What the existing models have done is to embody some of the salient regulatory features of nucleotide metabolism, show that the interactions JPT Vol. 4, No. 2--B
262
R.C. JACKSONand K. R. HARRAP
of drugs which are used to inhibit growth in the relatively simple L1210 cell culture system may be predicted in a reasonable proportion of combinations, and generally establish the feasibility of the approach. If the key regulatory features of a system are not written into a model, the behavior of the model will only fortuitously resemble that of the real system; thus, everything depends upon detailed and accurate quantitative biochemical studies of the target cells, and of the interactions of the drugs and their metabolites with the pathways concerned. The role of the mathematical model is to integrate the results of the analytical biochemical work and show how the pieces may fit together. The final question posed, and still not satisfactorily answered, is how quantitative changes in system parameters can produce qualitative changes in system behavior. The general conclusion drawn from the models discussed in this section was that usually this does not happen. Werkheiser et al. (1973) state: 'The 30 experiments yielded 270 interaction patterns. In many instances the intensity of the interaction was changed by parameter variation, but in only sixteen cases was there a change in the type of interaction'. These authors then quote the earlier conclusion of Garfinkel (1966): 'The behavior of the model is generally not very sensitive to changes in parameters; a drastic change in behavior usually requires the introduction of a new idea'. A similar conclusion was drawn from the detailed model of folate metabolism of Jackson and Harrap (1973) that usually a quantitative change in parameter values produced a quantitative change in the behavior of the model. However, from the regulatory point of view, it may be that the exceptions to this generalization convey more interesting information about the regulation of the system under study. Werkheiser et al. (1973) note: 'In four of these cases it was clear that the quantitative change in parameter value resulted in a qualitative change in the structure of the model'. An exact description of the kinetics which can show such properties is difficult. Clearly, in nature, quantitative changes in the concentrations of compounds can lead to qualitative differences in the systems affected. Examples are the effect of inorganic ion concentrations on the charged/discharged state of a nerve cell, the effect of the environmental lactose concentration on the/3-galactosidase content of E. coli (either uninduced or fully induced), and the fact that antimetabolites below a certain concentration may be tolerated by a tissue culture cell, but above that concentration the cell will die. Also relevant to the present discussion is the fact that malignancy or non-malignancy (clearly a qualitative property) is apparently determined by quantitative differences in the biochemical processes of the cell. These processes have been modeled, but always by introducing an arbitrary discontinuity into the model. For example, if the model is instructed that a cell will die if the concentration of tetrahydrofolate cofactors remains below a certain level for a given time, then the behavior of the model may reflect the properties of the cell, but it does not explain them. Far more interesting is the situation in which the discontinuity arises naturally from the rate equations. A few systems of this kind have been modeled; the study by Prigogine and Babloyantz (1972) on the control of E. coil /3-galactosidase is an example. All of the systems of this type studied so far have certain features in common. (i) There is an autocatalytic process, or a pair of cross-catalytic processes. In the c a s e of the /3-galactosidase system the entry of lactose into the cell is autocatalytic, since after sufficient lactose has entered the cell by passive diffusion to induce the lac operon, a galactoside permease is induced which greatly enhances the transport of lactose through the cell membrane. In the case of purine nucleotide biosynthesis there is a cross-catalytic process, in that GTP is a substrate for adenylosuccinate synthetase, and thus necessary for ATP biosynthesis, and ATP is a substrate for GMP synthetase, and therefore necessary for GTP formation. It is not known whether this cross-catalytic effect is of regulatory significance. (ii) Systems which give rise to discontinuities exhibit the phenomenon of multistability, i.e. they can exist in alternative stable steady states. Figure 10 illustrates this property in a simple model of glycolysis (Reich and Sel'kov, 1974). Figure 10a summarizes the highly simplified model used: v~ represents hexokinase, and product A is glucose
Computer models of anticancer drug interaction
263
6-phosphate, which inhibits its own production; ve represents phosphofructokinase (the intervening reactions between vl and v2 are not modeled, since they are considered to be of secondary importance from the regulatory point of view); v2 is activated by ADP, one of its products; v3 represents the succeeding reactions of glycolysis. Figure 10b plots v2 and v3 as functions of the concentration of B: v3 is considered to show a simple hyperbolic dependence upon B, and v, is not absolutely dependent upon B, but is allosterically activated by it; the solid lines represent v2 at increasing levels of A. Intersections of a solid line and the dotted line correspond to steady states for the subsystem vl, v2, i.e. points where B is produced and utilized at the same rate; note that at some values of A there is one steady-state point, for other values of A there are two, and at others, three. Where there are three steady-state points for a particular value of A, the middle one is always unstable for reasons explained by Reich and Sel'kov (1974). Figure 10c plots the value of the steady-state rate, v(--v2 = v3) as a function of A (solid line). Superimposed are plots of v~ as a function of A (dashed lines) for different amounts of enzyme 1 (hexokinase). At points at which a dashed line intersects the solid line, vt --- v2 - v3 and the system as a whole is in a steady state. Consider the system at point p~: the reaction is close to its maximum. If A is increased, or if E~ is increased, the system will move along the solid line toward Ps, and the steady-state rate, v, will show only a small increase. Now consider the system to be at point pro, when for some reason A is decreased, v2 and v3 can only reach a steady state by a marked reduction in rate, and the system drops to point p2; but now A is being produced (by v0 faster than it is used, so the system moves along the solid line to point p3, where once again vl-- v2 = v3. Now the system has reached a new steady state in which v is much less than before, and this large reduction in v was achieved by a very small reduction in A. Note, too, that the final value of A is greater than the initial value, which demonstrates the difficulties of drawing conclusions concerning changes in rates of processes from the size of precursor pools. Now imagine the system again at point pl, when the level of E~ is slightly decreased; the system in this case can only achieve a steady state by moving to p4, so again a small perturbation has resulted in a large change in v; if El now increases to its original value the system will not return to p~, but will move to p3. A
G6P
FDP ADP
E1
b
,. ,,
•
C
~|(A)
. . . . . .
::
\~
p \ ~.~.p_
~ \\ ~
_ ~ ' ~ 2 3
"~,, ¢\,',~\
v
"'?/' I
I
I
I
I
I
!
i
I
I
(A)
(e)
FIG. 10. Multistability in the glycolytic pathway. 10a. Schematic outline of the model of glycolysis used by Reich and Sei'kov (1974). v, is the rate of hexokinase (HK); v~ that of phosphofructokinase (PFK); v3, the succeeding reactions of glycolysis. Reactions between v, and v2 are not modeled. A represents glucose 6-phosphate (G6P) and B fructose diphosphate (FDP) and adenosine diphosphate (ADP). G6P inhibits its own formation, while v~ is activated by ADP. 10b. Plots of v2 (solid lines) and v3 (dotted line) as functions of B. 10c. Steady state values of v(-- v2 = v3) plotted as a function of A (solid line). Plots of vt as a function of A (dashed lines) are superimposed for different amounts of enzyme 1. (From Reich and Sel'kov, 1974.)
264
R. C. JACKSONand K. R. HARRAP
further increase in E~ will move the system along the solid line to p6, with a slight concomitant increase in v; increasing E~ still further will then cause a sudden large change as the system jumps to ps. Similarly, an increase in A when the system is at P6 results in a jump to P7, and then a decrease in A as the system moves to ps. It is clear that systems of this kind possess the appropriate properties to account for some biological switching processes. It remains to be determined whether this kind of kinetics could be involved in the triggering of cells from a static to a proliferating state, and if so whether antitumor drugs could then be employed to switch the system back again. In this connection, it is of interest that the adenine and guanine nucleotide branches of IMP utilization, as mentioned above, possess a cross-catalytic relationship, which could result in kinetics very similar to the autocatalytic situation just discussed in the glycolytic pathway. It is clearly important to conduct experimentally, detailed kinetic studies of enzymes which might show this kind of behavior, and then use detailed mathematical models to explore the optimal way to selectively switch off the biosynthetic processes of the malignant cells. 6. CYTOKINETIC MODELS The progression of cells through the cell cycle obviously depends upon intracellular biochemical events, which have become increasingly well understood in the last few years. However, it is convenient to separate consideration of models of the cell cycle from the models of metabolic pathways discussed in the last section, because the mathematical models of the cell cycle which have been produced primarily describe the movement of cells from one stage to another in numerical terms, rather than explain the cell cycle in biochemical terms. Cytokinetics, like pharmacodynamics, is a subject which is inevitably dependent upon a mathematical treatment, and whose complexity suggested at an early stage the development of mathematical models of some considerable sophistication. It is not the purpose here to review all of the uses of mathematical models of cell kinetics; instead, this discussion will concentrate on computer models which have been used to simulate drug effects. Some early models of cell growth were those of Barret (1966), Brockwell et al. (1970), Brockwell and Trucco (1970), Bronk et al. (1968), Bronk (1969), Franke (1970), Hahn, (1966, 1967), Hahn and Kallman (1967), Hirsch and Eagelberg (1965, 1966), Kember (1969), Takahashi (1968), Tibaux et al. (1969) and Trucco and Brockwell (1968). A review of some earlier applications of computer techniques in the study of the kinetics of tumor growth has been presented by Aroesty et al. (1973). An early model which made use of a large body of experimental data to study drug effects on the cell cycle was that of Shackney (1970), whose model was formulated to answer the following questions: (i) How can one predict whether an effective drug can be made a curative drug if dosage and schedule of administration are modified? (ii) Given that a schedule is potentially curative, how can one predict the optimal dosage and the optimal number of courses? (iii) Given that a particular drug may be curative in one experimental tumor system, how can one predict its behavior in another tumor system? (iv) Given that a drug is curative in animals when given according to particular schedules, how can one predict its behavior in human cancer, and how can one devise schedules that will be curative in man? Shackney attempted to approach questions of this kind by simulating the response of the LI210 mouse leukemia to AraC. His assumptions were such as to embrace as far as possible the known kinetic properties of the L1210 tumor. It was assumed that multiplying cells progressed in an orderly fashion through G~, S, G2 and M phases of the cycle and then re-entered G~. A reversible state of nonproliferation (Go) was also assumed. Growth was treated as asynchronous and it was assumed that a single viable cell was sufficient to produce a lethal tumor. AraC was treated as an S phase specific agent, and it was assumed that a given dose of drug would destroy the reproductive integrity of a constant fraction, not a constant number, of the susceptible population.
Computer models of anticancer drug interaction
265
The model made provision for nonproliferating pools of cells at each stage of the cell cycle, and for natural cell death. A distribution of cell-cycle times was also assumed. The constants for the model were then derived from certain standard treatment schedules. Experiments with the L1210 leukemia were simulfited at doses, schedules and implant sizes other than those used for standardization. The model's predictive capability was tested by simulating experiments without prior knowledge of results. Some results of such simulations are shown in Fig. 11. In both cases fair agreement was found between predicted and observed increases in life span. The author also discusses the limitations of his model. First, his treatment of slowly proliferating cells as nonproliferating, and the assumption that nonproliferating cells were completely drug resistant introduced errors, especially when large inocula, or noncurative drug doses, were simulated. The model made simplifying assumptions about the shape of the AraC dose-response curve, and did not deal with differences in drug-distribution characteristics at different anatomic sites in the host animal. Also, the model made no provision for host death from drug toxicity, and consequently at very high drug doses, predicted results better than those actually observed. Furthermore, the model did not provide for the emergence, during the course of therapy, of drug-resistant clones of tumor cells. Despite these limitations it was clear that a model of this kind was able to coordinate data concerning the kinetic properties of tumor cells, and allow that information to be used to optimize drug schedules. lO~
¢n .J ..i
+d I I
o
|
~//
:E
-I Z
lo3
1°:o
I
I t I !
t'l
/
I I
/t
/
/
I
Ii
I
/I
V
I/
I
II V
1 i I i
,
v i/I
v
I
I
I¢, I
7 9 11 1 15 17 19 21 23 25 27 TIME
0
I
4
I
I
8
I
I
12
/
I/
|
V
t i t/
'
I
16
~
/
I
20
~
J
24
I
I
30
I
I
32
(DAYS)
FIG. 11. Simulation of arabinosylcytosine (AraC) chemotherapy schedules. Treatment commenced at a tumor cell population of 107 cells. Solid lines represent the expected total number of cells, broken lines the expected number of viable cells. Schedule (a): 106 cells inoculated i.v. on day 0; treatment, 15 mg/kg of AraC every 3 hr x 8 on days 2, 6, 10, and 14; predicted host lifespan 26.6 days; predicted increase in lifespan 351 per cent; observed increase in lifespan: 240, 341,420 and 450 per cent in four experiments, respectively. Schedule (b): 104 cells inoculated i.v. on day 0; treatment,-15 mg/kg of AraC every 3 hr x 8 on days 5, 9, 13, and 17; predicted host lifespan 29.1 days; predicted increase in lifespan 234 per cent; observed increase in lifespan 281 per cent. (From Shackney, 1970.)
Jusko (1971, 1973) described models for the simulation of the effects of phasespecific and non-phase-specific drugs, which incorporated cytokinetic information. These models consider pharmacodynamic aspects of drug action as well, and will be discussed in the next section. Himmelstein and Bischoff (1973a,b) and Bischoff (1973) also described models for the study of phase-specific and non-phase-specific drugs which consider both pharmacodynamic and cytokinetic properties of the system (also discussed in Section 7). Roti Roti and Okada (1973) have described a mathematical model of the cell cycle of LS178Y mouse lymphoma cells. Their treatment makes the assumption that it is possible to divide the cell cycle into small time segments such that any given cell may
266
R. C. JACKSONand K. R. HARRAP
be uniquely classified into one of the states on the basis of cell age (i.e. time after division). Then, given the starting distribution of cells in the cell cycle and the growth conditions, the model will predict the distribution of cells at different times in the cell cycle. The authors of this model used it to simulate the effects of colcemide, thymidine, and colcemide plus thymidine. Donaghey and Drcwinko (1975) have described a flexible computer simulation program which they claim may be used for a variety of cytokinetic studies. The program was evaluated using data from cultured human lymphoma cells. In this model S phase was divided into three equal segments (early, mid, and late S); the length of stay for each cell in each segment of the cycle was independent, and normally distributed. G2 phase was subdivided into two equal segments. The length of stay of each cell in both segments was triangularly distributed and skewed to the right, because it was thought that a cell would need a minimum time to complete the biochemical events required for mitosis, but after that the exit time from G: would depend upon events which could be treated as random. The length of time in mitosis was similarly skewed. The proliferation rate was set at 1.6 (i.e. every 100 cells dividing were assumed to produce 160 daughter cells); this assumption was based upon the experimental observation that the probability of division was 0.82 for these cells. Ten per cent of the cells emerging from division were assumed to enter a Go compartment; since the model was only used to simulate short periods of time, the subsequent fate of these Go cells was not studied. The G1 phase was divided into an early period and a late period. The early GI period was skewed to the left, on the assumption that cells would proceed through a number of randomly distributed events until a state was reached in which cell division became likely, and that when this point was reached the cells proceeded fairly rapidly into the second stage of G1, where the cell was considered to be accumulating precursors necessary for it to enter the DNA synthetic phase. This model was used to predict the distribution of cells through the cell cycle following a 24-hr thymidine block. The model predicted that under these conditions 92 per cent of the cells would accumulate in S phase (the observed value was 85-90 per cent). These authors state their intention of using this model to study the effects of anticancer drugs. They conclude: 'We realize that the model has several deficiencies and that many assumptions are arbitrary. These shortcomings result from the paucity of knowledge about many fundamental biological principles governing the behavior of cell populations'. However, this CELLSIM program is sufficiently flexible that new information may be easily incorporated into models based upon it. Blumenson (1973, 1975) has produced a comprehensive model of the human granulopoietic system, and has applied this mathematical model to the study of cancer chemotherapeutic agents, in particular the relationship between drug effects on bonemarrow cells and the ensuing pattern of leukopenia in the peripheral blood. This author suggests that the origin, maturation, release from the bone-marrow, migration into the tissue spaces, and eventual death of granulocytes are now sufficiently well understood that detailed models are likely to give meaningful information about the relationships between changes in the various compartments of granulocytes and their precursors. These various compartments interact by a number of feedback controls. For example, during infection there is an increased demand by the tissue spaces for granulocytes from the peripheral blood, which in turn, affects proliferative controls in the marrow. The model is a deterministic, rather than a stochastic one; for example, the cell kinetics are described in sufficient detail to permit one to distinguish between the different phases of the cell cycle, and it is assumed that, once triggered into cycle, the cells pass through the different phases in an ordered and obligatory pattern. The cells in the marrow are classified into 130 one-hour compartments of increasing maturity. In the steady state, the size of each such compartment is equal to the steady-state hourly granulocyte turnover rate (GTR); the cells in the oldest (130th) marrow compartment are released into the blood each hour. Thereafter, all cells in the total blood granulocyte pool (TBGP) have an equal chance of leaving for the tissue spaces during a 1-hr period. This probability depends upon the needs of the tissue
'
Computer modelsof anticancerdrug interaction
267
spaces, and may vary from hour to hour. Not only do cells in the 130th marrow compartment always enter the peripheral blood after 1 hr in that compartment, but younger cells also may do so according to need, which is determined by the hourly blood granulocyte c o u n t . The loss of maturing cells from the marrow in turn stimulates stem cell activity. Stem cell block, and its consequences, may be simulated. Certain parts of the modeling of stem cell cycling are an extension of the earlier models of Lajtha et al. (1962) and King-Smith and Morley (1970). The faithfulness of Blumenson's model was tested by simulating the effects of 5-FU, considered to be administered as an intravenous drip for a period of 0.5-3 hr. It was assumed that a certain fraction of cells in S phase were killed during the treatment, this fraction depending upon dose and time of administration. Not only are stem cells killed by the drug, but also myeloblasts, promyelocytes and myelocytes in S phase. This leads to a decrease in total marrow cellularity which increases the probability of uncommitted stem cells differentiating in the granulocyte track. The predicted peripheral white blood cell counts following the simulated 5-FU treatment showed the kind of highly fluctuating pattern actually observed in treated patients. This model may thus be used to deduce the effects of 5-FU on the bone-marrow by projecting back from results found in the peripheral blood, and improve the accuracy of monitoring chemotherapy without the need for frequent marrow biopsies. Other suggested applications for this model are in the study of bacterial infection, to simulate bone-marrow transplants, and to study the breakdown of the feedback mechanisms for granulocyte production and release, and thus examine the implications of the various theories relating the loss of such control mechanisms with the etiology of chronic granulocytic leukemia. Although this model has not apparently been used to study drug interactions, the structure of the model would probably permit studies of this kind. Blumenson's model is of particular interest and importance as an example of a system which studies the effects of antineoplastic drugs, not on the tumor itself, but on a proliferative organ of the host. Woo et al. (1975) have developed a cell-cycle stage-specific multicompartrnental model to investigate the effects of drugs on tumor growth. The model can simultaneously consider the following types of effect: (a) non-cycle-specific cell kill and cycle stage-specific cell kill and (b) progression delay of cells through the cell-cycle phases, which may bring about an accumulation of cells in the various phases of the cycle. The model was an elaboration of an earlier one by Dombernowsky and Hartmann (1972). A schematic diagram of this model is shown in Fig. 12. This model was used to simulate the effects of BCNU, AraC and MTX on the cell-cycle and proliferation kinetics of the L1210 leukemia. BCNU kills cells in all stages of the cell cycle at equal rates and nonproliferating cells at slower rates, and in addition elongates cell-cycle time by increasing transit time in S phase. AraC causes some progression delay through all phases, especially in G1 phase but only kills cells in the S-phase. By comparison of predicted fractional survival curves, obtained using a range of parameters in the model, with experimental values, it was possible to obtain the best estimates for the rates of cell kill, rates of inhibition, and elongation of transit times for each cell compartment, at a number of dosage schedules for these two agents. In the case of MTX, the model was used to examine experimental data of Skipper et al. (1970) obtained at two different concentrations. Simulations in which a large fraction of cells were inhibited (i.e. transferred to a nonproliferative 'quiescent' compartment) and a certain dose-dependent fraction of cells in S phase were killed, gave close agreement with experimental results. These authors discuss how models of this kind may be used to optimize drug scheduling. For example, under conditions in which a drug causes cell kill in S phase, but also delays progression from Gl to S, the first dose of drug will kill the cells already in S phase. Considerably later, a wave of the delayed cells will then enter S, and if this time can be predicted accurately, a second administration of drug can be timed to kill this cohort .of delayed cells. It is suggested that this model can provide reliable predictions using far less detailed experimental information than that required by the model of Shackney.
R. C. JACKSONand K. R. HARRAP
268 -
-',
.
Sphase .
.
.
.
.
rLS--
I;
I L-. ii ~1
; T ;
.
.
.
I '
.
:
.
.
,,o
I :
L.._ ;..~_] 1 L_I 3 ~
.
'~
r
MITOSIS .
.
.
.
.
.
.
.
:
L
"
FIG. 12. Schematicdiagramof the relationshipbetweendrug-tumorinteractionsillustrated by the cell cycle stage-specific model. Each block enclosed by broken lines represents three possible modes of drug action within each cell cycle phase;/3* denotes the drug-inducedcell kill rate, 3' the fractionof cells accumulatingin the drug-inducedquiescentcompartment;transit time in each phase may also be elongated. (From Woo et al., 1975.) Chuang and Lloyd (1975) and Chuang (1975) have presented a mathematical analysis of the tumor model originally presented by Skipper (1969) in which the tumor cells are treated as existing in three compartments, proliferating, temporarily nonproliferating (resting) and permanently nonproliferating. The proliferating cells pass through the four cell-cycle phases, (i.e. G~, S, G2 and M). Chuang and Lloyd (1975) have developed a treatment of such a system in which the effects of antitumor agents are expressed in terms of loss functions for each proliferating (or potentially proliferative) compartment, these functions depending upon the drug concentration at the tumor site. The simulations based upon this approach were used to suggest a 'best' schedule, which did not, however, consider host toxicity. The model was tested by simulating the effects of AraC in a variety of dose schedules. An example of one such simulation is illustrated in Fig. 13. These simulations were made possible by parameter estimations based on data indicating the percentage of labeled mitoses, using a statistical procedure which calculated the theoretical fraction of labeled mitoses from the experimental values. It is shown that a model of this type can be used to distinguish between various hypotheses of drug action and inter-compartmental relationships of cells, by simulating the theoretical situations and determining which system best describes the observed behavior. Bender arid Dedrick (1975) have used mathematical modeling to investigate cytokinetic aspects of resistance to anticancer drugs. Their model contains the usual four proliferating phases of the cell cycle, a viable but nonproliferating Go phase, and a nonproliferating nonviable compartment. It is assumed that spontaneous cell loss may occur from all of these compartments. These authors used their model to calculate theoretical fractions of cell kill for an S phase specific agent for tumors with 16- and 24-hr cycle times, using repeat drug doses spaced at 8- and 16-hr intervals. The factors which may influence possible cell 'recruitment' into cycle are discussed. It is concluded that cytokinetic factors which influence responsiveness to drugs are tumor size and location, rate of cell proliferation, distribution of cell-cycle time and S phase transit time in the cell population, variation in the growth fraction, and the degree of synchrony of the tumor cell population. 7. M U L T I L E V E L MODELS Since the understanding of the interactions between a drug, a tumor, and a host animal requires study at several levels of biological organization, an obvious ex-
269
Computer models of anticancer drug interaction I06
z 0
'
'
'
'
'
'
'
'
'
'
÷ ' G I - 'C E L L I S x D o
iO 5
S-CELLS G2-CELL M-CELLS
S
I.< "J
:E bJ
,04 o
,. ¢,n
tu b-
10 2
e X td
I-
I0 I
U, I
i0 O
,
,
,
,
,
~
i
i
i
i
i
i
.
.
,
.
q
.
~
.
.
~
.
i
i
i
i
o) ,oO ~ E° -
o
,o-,~ io-2t
0
;
'~
3
;,
~
~ TIME,
r
~
~
,o
.
I~
13
,4
,5
DAYS
FIG. 13. C o n c e n t r a t i o n o f a z a b i n o s y l c y t o s i n e ( A r a C ) a n d e x p e c t e d cell p o p u l a t i o n s in p r o l i f e r a t i n g phases of the L1210 leukemia. Animals recieved an inoculum of l0 s cells and were
treated on a schedule of 15 mg/kg/dose every 3 hr x 8 on days 2, 6, l0 and 14. (From Chuang, 1975.)
tension of the types of mathematical model discussed in the preceding sections was a model containing features that embraced two or more such levels. An early example of this kind of system was the model of MTX action described by Werkheiser (1971). Fourteen parameters were used to describe the action of MTX on a tissue and its environment. They were: (a) rate of urinary excretion of drug; (b) tissue permeability; (c) generation time of cells and duration of S phase; (d) rate of dihydrofolate reductase (DHFR) synthesis; (e) rate of decay of free and MTX-bound DHFR; (f) Km and Vmx for D H F R and thymidylate synthetase; (g) total folate concentration; (h) normal D H F R concentration; and (i) critical FH4 concentration. This last parameter is based on the assumption that if the cellular tetrahydrofolate (FH4) concentration falls below a certain value then the cell will die. This assumption is not necessarily biochemically warranted, and is restrictive, in that cells can probably tolerate much lower than the normal availability of reduced folates if a supply of thymidine and purine is available; nevertheless, since the nature of antifolate-induced cell death (as distinct from mere growth inhibition) is not well understood, some kind of arbitrary condition of this kind is hard to avoid. It will be seen from the list of parameters given that some of them describe intracellular biochemical functions (i.e. enzyme kinetic properties, and rates of enzyme synthesis and turnover), some describe cytokinetic qualities of the target cell (i.e. generation time and duration of S phase) and some are concerned with the pharmacodynamics of the system. Thus, the behavior of MTX was concerned with three different levels of biological organization, though in all cases in a highly simplified way. The model as thus formulated resembled the properties of real systems in a number of respects. MTX toxicity, at sufficient doses, was cumulative; the total dose for lethality was 70-times the total when divided into four daily doses; and these properties were in agreement with experimental observations, as was the steepness of the MTX dose-response curve. This model predicted that tumors with low membrane permeability to MTX would be most sensitive to large single doses, while highly permeable tumors were most sensitive to small repeated doses; this prediction was also confirmed experimentally. An example of a pair of predictions is shown in Fig. 14. It was concluded that this model, simple though it was, was sufficiently realistic to have predicting power.
R. C. JACKSON and K. R. HARRAP
270
400.MOLAR ON DAY 0 10
~
700~MOLAR ON DAY 0
EE MTX (~M)
. FeEE F" 2
F~ COFACZORS,~-"-" (.M)
.k .... CELL COUNT
CELL COUNT
103 lO 1
DAYS
DAYS
FIG. 14. Time course of the effects of a single course of methotrexate (MTX) on a simulated tissue without a diffusion barrier. FH~, dihydrofolate; FH4, tetrahydrofolate. (From Werkheiser, 1971.)Reproduced by permission of the New York Academyof Sciences. A more recent model, along similar multilevel principles, is the very detailed description by Morrison et al. (1975) of the disposition of AraC and its metabolites. A pharmacokinetic compartmental model is used, based on that of Bischoff et al. (1971), with the augmentation by an additional compartment to represent the spleen, in the portal circulation, since this organ is believed to make an important contribution to AraC phosphorylation. The pharmacokinetic flow scheme used is illustrated in Fig. 15a. For each compartment shown mass-balance equations were drawn up which described the organ concentrations of AraC and its metabolites in terms of drug-flow in and out of each compartment. Where this model went beyond a purely pharmacokinetic scheme, however, was in the inclusion of a considerable number of intracellular biochemical processes concerned with the activation and inactivation of AraC (Fig. 15b). AraC is converted to its 5'-monophosphate by deoxycytidine kinase, which is inhibited by dCTP and AraCTP; deoxycytidine, as a competing substrate, also inhibits AraC phosphorylation. The enzymes deoxycytidylate kinase and nucleoside diphosphokinase further phosphorylate AraCMP to the triphosphate level, and these reactions are opposed by corresponding nucleotidases. The deamination of AraC and AraCMP to AraU and AraUMP by cytidine deaminase and deoxycytidylate deaminase, respectively, is also modeled, as is the interconversion of AraU and AraUMP. Apart from effects of dCTP and AraCTP on deoxycytidine kinase, no other feedback inhibitions are considered, although there is provision for competition for cytidine deaminase by the competing substrate, deoxycytidine. This model was then used to simulate events following the administration of single doses of AraC to mice. Of particular importance was the intracellular concentration of AraCTP, the active form of the drug. This was shown to persist at cytoxic concentrations long after most of the AraC was cleared from the blood, explaining the fact that a single dose of AraC will kill cells that were not in S phase (as measured by the incorporation of tritiated thymidine) during the period of measurable plasma AraC. The fact that AraCTP follows a very different time-course from AraC is obviously of great significance in understanding the action of AraC. Higher nucleotidase activities, of course, decrease the steady-state concentration of AraCTP and the length of time that it persists at cytotoxic levels. The livers of DBA and C3H strains of mice show a very large difference in deaminase activities, and the effects of this enzymatic divergence were simulated, resulting in major differences in the predicted pharmacokinetics of AraC.
271
Computer models of anticancer drug interaction
! I
~.E,S
r O.,NE ~
LEAN
.......i
I.
r
(o)
dTMP I dGMP J dAMP
membrane oxtmce~dsr JMJ'g,~oJhdw fluid fluid
+
d~/,
. . . . . . . . . . . . .
.® .e~ I
"~,,11~
~
•
- dCTP
® : ~
I
a r a - C , i b - l a r l i - C 4*- I l r l I - C M P .4. a r a - C D P .Q- I m - . C T P . . . . .
®
(i)
•®
~
.NA
® m-u:-;m-u~- m-UMP :- ,,.-UOP = a,.-UrP I
,
0
(b) FIG. 15. Pharmacokinetic simulation of the disposition of arabinosylcytosine (AraC) and its metabolites. (a) Flow scheme of the distribution of AraC to various tissue compartments. (b) Metabolism of AraC. AraU, arabinosyluracil; araCMP, arabinosylcytosine monophosphate; araCDP, arabinosyi cytosine diphosphate; araCTP, arabinosylcytosine tripbosphate; araUMP, arabinosyl uracil monophosphate; araUDP, arabinosyl uracil diphosphate; araUTP, arabinosyl uracil triphosphate; dTMP, thymidine monophosphate; dGMP, deoxyguanosine monophosphate; dAMP, deoxyadenosine monophosphate; dCTP, deoxycytosine triphosphate; dCR, deoxycytidine. (From Morrison et aL, 1975.)
The authors state that future developments of this very interesting model will include more biochemical details, such as further phosphorylation of AraUMP; more body compartments, such as the peritoneal cavity; scale-up for human application; and inclusion of cytokinetic information. It is to be hoped that this approach may be extended to the study of other antineoplastic drugs and to drug-drug interactions. Werkheiser's model (1971) was an attempt to combine simple descriptions of pharmacodynamic, cytokinetic and biochemical properties of a tumor-host system in a single model. The simulations of Morrison et al. ~(1975) included more biochemical information and a great deal more pharmacodynamic information, but did not, at least in the version described, include cytokinetic detail. The remaining work to be
272
R.C. JACKSONand K. R. HARRAP
reviewed in this section is concerned with mathematical models combining cytokinetic and pharmacodynamic aspects. Jusko (1973) described a pharmacodynamic model for cell-cycle-specific chemotherapeutic agents which was a development of his earlier work on the pharmacokinetics of non-phase-specific agents (Jusko, 1971). In the cell-cycle-specific model the cells are considered to exist in either of two interconvertible compartments Cs, the drug-sensitive cells, and Cr, the drug-resistant cells; interconversion of cells from the two compartments was described by first-order rate constants. The nonproliferating cell compartment was susceptible to cell loss or differentiation, and the proliferating compartment could be depleted by drug-induced cell death. Drug-receptor interaction was considered to be bimolecular and irreversible, and assumed to bind a negligible proportion of the total drug dose. The sensitive (Cs) cell compartment was capable of cell proliferation, represented by a single parameter; the model could represent steady state, exponential, or Gompertzian behavior in the absence of drug. The pharmacokinetic aspects of the system were represented either by a one-compartment open model or a two-compartment open model. The model was tested by using it to fit data obtained by Bruce and co-workers with the effects of vinblastine on lymphoma and hematopoietic cells. The effects of repeated doses of AraC were also simulated. As well as predicting cell kill, the model could be used to estimate the survival time of drug-treated animals. Himmelstein and Bischoff (1973a, b) have developed a pharmacokinetic/cytokinetic model which simulates cell-cycle events in considerable detail; the cytokinetic treatment was based upon a model by Rubinow (1968). An important difference from the model of Jusko is that a continuous distribution of cellular maturity is considered, rather than two discrete fractions. The authors consider the possibilities that drugs kill cells or render them incapable of cell division, and do not describe the situation of drug-related growth inhibition by prolongation of generation times. The effects of direct action drugs (i.e. agents which kill cells regardless of their stage in the cell cycle) are first treated on the assumption that cell growth in the absence of drugs would be exponential. Drug-induced cell death is computed from the loss function, which relates cell loss to drug concentration. Various types of loss function are described; for example, first-order, or two-parameter saturable, or a threshold type relationship, where drug is only effective above a certain concentration, or combinations of these three possibilities. On the basis of these simplifying assumptions, an explicit relationship was derived between the cell growth parameters, and the loss function for a non-phase-specific (direct action) drug, and the cell number at a given time. Then the effect of cycle-specific drugs is described, now with the assumption that the loss function varies at different phases of the cell cycle. The drug concentration at a given time is treated by conventional pharmacokinetic expressions. The model of Himmelstein and Bischoff was then tested by simulating the interaction of AraC with mouse leukemia L1210 cells, although curiously, making the assumption that the drug was not phase-specific; this treatment was justified on the grounds that the L1210 cell spends 90 per cent of its time in the S phase, at least when the total cell number is small. Data of Skipper et al. (1964) were considered. The pharmacokinetic parameters were based on measurements of AraC clearance made by Skipper et al. (1967), although as shown by Chou et al. (1975), the kinetics of the active drug form, AraCTP, are quite different from those of AraC. The model of Himmelstein and Bischoff was also modified to describe the results of 12hr continuous infusion of AraC. Despite some of the rather drastic simplifications made during the course of construction of this model, it did appear capable of giving fairly accurate simulations of the results of various AraC treatment schedules with the L1210 system. These authors state 'a comprehensive model for a given system is composed of suitable growth kinetic, drug-cell interaction, and pharmacokinetic models'. Their own studies represent a preliminary step toward this goal of a complete, multilevel description of drug-host-tumor interactions. Lincoln et al. (1974) describe a plan for a wide-ranging 'leukemia simulator', which would not only consider drug clearance and distribution, and the cytokinetics of
Computer models of anticancer drug interaction
273
normal and transformed cells in the marrow, but also gastrointestinal susceptibility, kidney function, liver function, and evaluate risks from present or probable infections, risks of hemorrhage, etc. This grand synthesis represents, at present, a goal rather than an accomplised fact, and the authors describe some of the work in progress towards this objective, illustrating their model's capabilities by simulating the effects of AraC on granulocyte and platelet counts. Finally, one model has been presented which describes the effects of combined chemotherapy and radiotherapy (Andrews, 1969). This preliminary treatment of combined modality therapy will probably form a basis for more detailed theoretical explanations of how to optimize cancer treatment using radiation in combination with drugs. 8. T E C H N I C A L ASPECTS OF COMPUTER MODELING Since this area of research represents the common ground between biology and mathematics, it tends to be pursued by scientists whose primary background lies in one or another of these disciplines, but who move with less confidence in the other. Thus, there are, in the literature, a number of mathematical models of considerable technical sophistication, but which have never been applied to any real biological problem; conversely, there are attempts by biologists to describe pharmacological systems in a quantitative way which are hampered by use of inadequate mathematical techniques. This state of affairs is beginning to change as more pharmacologists become proficient in computer programing and numerical analysis, and as investigators from the traditionally strongly mathematical field of physiological control theory have come to take an interest in quantitative pharmacology. Existing mathematical models of pharmacological systems fall into a number of groups which reflect the different underlying philosophical approaches. 8.1. ABSTRACT AND REALISTIC MODELS
Weinberg and Berkus (1971) review a number of simulations written for abstract systems which resembled living cells in their general organization. These simulations were not models of any particular system, and were not designed so that the results of the simulation could be checked against experimental data. These models may be used, for example, to consider the various possible modes of organization which could be used by living matter to solve particular problems of control. In contrast, are realistic models which attempt to describe a particular biological system in which the aim is to describe the experimentally observed behavior of the system as closely as possible. The simulation of E. coil by Weinberg and Berkus (1971) falls into the latter category, as do most of the simulations discussed in the present review. 8.2. STOCHASTIC AND DETERMINISTIC MODELS
Chemical kinetics, the science from which enzyme kinetics and, more remotely, cytokinetics and pharmacokinetics are derived, is a deterministic science; given a set of rate constants and a set of reactant concentrations, the rate equation may be used to compute the reaction rate. Because the behavior of bulk matter is considered, the probabilistic aspects of the motions of individual molecules have no influence. The work of recent years has shown that processes obeying the laws of chemical kinetics, and the law of diffusion, are capable of giving rise to discontinuities in space and time; thus, in principle, given sufficient knowledge of enzyme kinetics (including the kinetics of multienzyme aggregates and of membrane-bound enzymes), of transport processes, and of the dynamics of circulation, excretion, etc., it should be possible to construct a strictly deterministic model of a pharmacological system which will provide a description to any required degree of detail. In practice, for a model of a very complex system, it is impossible, on a computer of present-day size, to directly simulate every process which might have some trivial influence on the system. This
274
R. C, JACKSONand K. R. HARRAP
problem has been encountered in the field of meteorology, in the construction of weather-predicting models of the atmosphere, where it is known as the 'butterfly-wing effect'; this name reflects the fact that, in a. system of sufficiently low stability, the flutter of a butterfly's wing could set in motion (in theory, at least) a chain of events which could culminate in a hurricane a thousand miles away. The meteorologist cannot simulate every butterfly, and thus must accept the fact that complex systems are subject to small perturbations which must be treated as random. Obviously, exactly analogous considerations apply to pharmacology, and the decision as to whether a particular process is to be treated as random or deterministic depends upon how central the process under consideration is to the questions being asked by the model. For example, if it is observed experimentally that a certain dose of drug, administered for a certain time, kills 40 per cent of tumor cells, then, if a single cell is considered, it might be reasonable to treat the survival or otherwise of that cell as a stochastic process, with a probability of death of 0.4. However, if it was known that the cell population was mixed, and that 60 per cent of the cells had plasma membranes impermeable to the drug, it might be more reasonable to decide the fate of a single cell on a deterministic basis, according to the membrane permeability of the drug. Thus, it is often convenient to incorporate both deterministic and random events into a model, the latter being used to simulate either processes which represent small perturbations rather peripheral to the real purpose of the model, or sometimes processes which cannot be simulated deterministically because of lack of information about what controls them. Technical aspects of the simulation of stochastic processes by MonteCarlo methods are discussed by Sheppard (1969). 8.3. THE SCOPE OF A SIMULATION
Simulations of biological systems vary from broad, shallow descriptions intended to describe the essential features of a wide-ranging process, such as the effects of drugs on progression through the cell cycle, to descriptions in depth (e.g. at the molecular level) of small subsystems (e.g. a dozen enzymes with their substrates and cofactors). A simulation may, of course, attempt to be both wide-ranging and deep; eventually, either programing time, expense, or computer memory availability will become a limiting factor. The decision of what ground should be covered by a simulation, and in how much detail, must be based upon experience of the real system under study, with continual regard to the questions which the model is being constructed to answer. Too much unnecessary detail makes a model expensive and hard to manage, and too little detail makes an unrealistic model, possibly so unrealistic as to have no predictive value. Garfinkel (1969), considering this question, suggests that a useful compromise may be to simulate the subsystem under study in as much detail as possible, with the remainder of the system (the 'environment') described in less detail. 8.4. LINEAR AND NONLINEAR MODELS From the very beginning of the study of enzyme kinetics it was clear that linear (first-order) kinetics were only observed in rather special circumstances, such as when a substrate concentration is very low in comparison to the enzyme Km value, and this science has therefore stressed nonlinear models from the outset. In some qualitative types of biochemical models, 'linearized' transformations of approximate rate equations have been used in a way such that the model retains the qualitative properties of the prototype system (e.g. Reich and Serkov, 1974), but in general, models of multienzyme systems inevitably imply nonlinear equations. The situation has been rather different in pharmacokinetics, where many processes of drug uptake, distribution, and elimination can be described with reasonable quantitative accuracy by first-order kinetics. Linear models have also been used in cytokinetics. Nevertheless, even in pharmacokinetics, it is clear that a detailed description of many processes entails a nonlinear model; for example, the active transport processes responsible for cellular uptake of numerous antitumor drugs may often be described by a Michaelis-
Computer models of anticancer drug interaction
275
Menten equation. An exception to the general rule that biochemical systems cannot be studied rigorously using linear models is provided by tracer experiments. These may be considered as first-order perturbations, and thus permit the study of complex nonlinear systems as a set of simpler, linear subsystems (Berman, (1969). The general rule that simulations in biochemical pharmacology require nonlinear models obviously does not mean that particular components of these models may not show first-order kinetics; this is frequently the case, for example, with the non-carriermediated membrane transport of drugs, or non-enzymic reactions, such as the spontaneous breakdown of 5,10-methylene tetrahydrofolate to tetrahydrofolate plus formaldehyde. Thus, models of pharmacological systems are likely to consist of sets of nonlinear differential equations. In some simple cases, the differential equations may be amenable to transformation into algebraic expressions, but this will not generally be the case, and for the most part the system will have no explicit solution. The usual approach, therefore, is to solve the differential equations numerically, and for this purpose second- or fourth-order Runge-Kutta or predictor-corrector methods seem to have been most used. The only serious problem encountered is that of 'stiffness'. This property, so called because it was encountered first in the study of vibrations of stiff springs, results in a system of equations in which one component has a time constant which may differ by orders of magnitude from those of the rest of the system. In this situation, although the equations can be solved by the usual methods, doing so entails extreme inefficiency. The problem was first reported in biochemical systems by Garfinkel and Hess (1964); these authors were able to alleviate the problem to some degree by altering the molar enzyme concentration (upward) and the turnover number (downward) so that the simulated enzyme activity did not change. For general methods of solution of sets of stiff differential equations, the reader is referred to articles by Butcher (1964), Cooper (1969) and Gear (1971). 8.5. TEMPORAL HIERARCHIES IN BIOLOGICAL MODELS
Any biological model of broad scope may have to simulate processes which take place over very different time scales. Typically, enzyme-substrate interactions take place in fractions of a second, or milliseconds; intermediary metabolite equilibrations may take seconds or minutes; biosynthesis and degradation of enzymes in mammalian cells often involve halfqives of several hours; and the cell division cycle takes between 9 hr and many days. Clearly, the simultaneous simulation of all of these processes, as would be required for a complete understanding of drug-cell interactions, is likely to pose great technical difficulties. This problem has been the subject of a theoretical discussion by Park (1974) who suggests methods for handling problems of this kind. The present authors are not aware of any applications of this approach to a realistic biological model, but it seems to be potentially of great value. The problem of a wide range of time scales was encountered during the construction of the model of folate metabolism in the L1210 mouse leukemia cell (Jackson and Harrap, 1973). Substrate and cofactor interconversions required a few seconds or minutes, MTX transport across the cell membrane typically required from 15 min to several hours before a steady state was reached, and enzyme turnover times ranged between 18 and 90 hr. The approach to the problem was as follows: (i) enzymesubstrate and enzyme-inhibitor complexes were assumed to be in or near the Briggs-Haldane steady state, so that steady-state rate equations could be used, with the exception of the complex between dihydrofolate reductase and MTX, which dissociates with a t~12of about 1 hr, and (ii) the slow processes (i.e. cell division, drug clearance, drug uptake through the cell membrane, and synthesis and breakdown of enzyme molecules) were simulated directly. The rapid processes of substrate and cofactor interconversion were assumed to be close to the steady state throughout. Thus, the cellular inhibitor concentration and enzyme levels at a given time would be calculated by direct simulation, and the substrate and cofactor levels corresponding to
276
R. C. JACKSON and K. R. HAP.RAP
the appropriate enzyme activities would be obtained by interpolation in an array of stored values which had previously been shown by direct simulation to correspond to those enzyme activities. This part of the process was iterative, since, of course, actual enzyme activities would depend upon levels of substrates, feedback inhibitors, etc. The chief disadvantages of this 'two-level' method are, first, the large amount of computer time needed to assemble the table values of steady-state substrate and cofactor concentrations (some of which may never be used) and, second, the large amount of memory required to store these tables of steady-state values when the second stage (slow reactions) of the simulation is running. 8.6. DIGITAL AND ANALOG COMPUTERS Once a model has been constructed, that is, the problem has been stated completely in mathematical terms, the actual solution of the resulting equations is purely a matter of technical convenience. Now that many biologists are conversant with the programing of digital computers, it is likely that most workers who wish to use mathematical modeling as an occasional problem-solving tool will use digital computers. Analog computers have enthusiastic proponents, chiefly on the grounds of their relative simplicity of operation and low initial cost. However, it is probably fair to say that proficiency in analog computing requires more specialized expertise than participation in a large, time-shared digital computer system, which, from the users' point of view, can be treated as a problem-solving black box. Hybrid computers are gaining ground slowly but may be much more widely used in the future; their application to biological problems was discussed by Girling (1969). 8.7. SIMULATION LANGUAGES The digital computer programer has a choice between writing his program entirely from the beginning in a general-purpose language, usually FORTRAN,BASIC,or ALGOL,or of using one of the special-purpose simulation languages. The choice probably depends upon how much the completed simulation is likely to be used. A program completely written especially for a particular function is likely to be more efficient in terms of computer time and memory, but will no doubt require a greater investment of programing time, than one which uses a~ special simulation compiler. The latter generally require a lot of core memory, but they seem to be very suitable for this kind of work; for example, Werkheiser et al. (1973) used MIMIC. 9. CONCLUSIONS: FUTURE USES OF COMPUTER MODELS OF DRUG INTERACTION It is clear that although mathematical modeling has become an established technique in pharmacology, there are still very few instances of the method being applied to the study of drug-drug interactions. Many of the models described in this article have been concerned primarily with reaching an understanding in depth of the action of a single agent on a particular system, or at most with the effects of a small number of agents used singly. The few models described which have studied drug-drug interactions, dealing mainly with interactions at the level of cellular metabolism, however, show the potential of simulation, both as an investigative method in furthering the understanding of drug-host, drug-tumor and drug--drug interactions; as a predictive tool, in suggesting probable positive and negative interactions of various agents; and as a means of optimizing chemotherapeutic regimens, to tailor them to the biochemical characteristics of particular individual tumors. Several of the models discussed here, although not yet used to study drug-drug interactions in the systems they simulate, seem to have the potential to be used for that purpose, and in a few cases their authors have stated their intentions of pursuing such studies. Since the main purpose of these models is to quantitate the properties of systems too complex to be understood intuitively, it is reasonable to believe that the greatest
Computer models of anticancer drug interaction
277
contribution to understanding the complexities of cancer chemotherapy, in particular combination chemotherapy, will be made by those models which are sufficiently broad and deep to represent accurately all the essential features of the interactions between a tumor, a host, and several drugs. This will require multilevel models; that is, models containing biochemical, pharmacodynamic and cytokinetic features. How would such models be most effectively used? Consider that a new antitumor drug has been introduced, and its efficacy as a single agent has been demonstrated. When the time comes to consider possible ways of combining the new drug with other agents, a good mathematical model could make a valuable contribution. The number of possible combinations of, say, three or four drugs, at a number of different dose ranges, is very large, and to test all possible combinations experimentally would be prohibitively time-consuming and expensive. To simulate a large number of such combinations with a computer model would be relatively rapid and inexpensive. Those combinations predicted to be most promising could then be tested experimentally. Thus, the computer model is not a substitute for experiment, but a guide which should help to increase the productivity and cost-effectiveness of limited experimental resources. A second example of how computer models may help to improve chemotherapeutic effectiveness in ways which might otherwise be impossible, is in the optimization of chemotherapy tailored to the biochemical properties of particular tumors. If certain biochemical facts (i.e. enzyme activities, nucleotide pool sizes, etc.) are established for a biopsy specimen, how may every item of available information be most effectively exploited? By incorporating the available information into a model, various drug regimens could be tested against a simulated system which represents the properties of the tumor as accurately as possible, and the 'most probably effective' regimen selected for actual use. Thus, it is felt that mathematical simulation is likely to become an established biological and medical technique, not as a substitute for human intelligence, experimental skill, or clinical judgement, but as a tool which has its place, along with all of the others, in increasing the range and effectiveness of scientific medicine in the increasingly complex field of antitumor chemotherapy. AcknowledgementsmTbe authors thank the following publishers for permission to reproduce various materials cited in the text: Academic Press; American Pharmaceutical Society; Elsevier; National Cancer Institute; New York Academy of Sciences; and Pergamon Press. This work was supported by grants from USPHS (CA-18129), the Medical Research Council and the Cancer Research Campaign.
REFERENCES ALLEN, L. M. and CREAVEN, P. J. (1975) Human pharmacokinetic model for isophosphamide. Cancer Chemother. Rep. 59: 877-882. ANDREWS,J. R. (1969) Combined cancer radiotherapy and chemotherapy: the relevance of cell population kinetics and pharmacodynamics. Cancer Chemother. Rep. 53: 313-324. AROEffrY,J., LINCOLN,T., SnAPmo, N. and Boco^, G. (1973) Tumor growth and chemotherapy: mathematical methods, computer simulations, and experimental foundations. Math. Biosci. 17: 243-300. BAtmE'rr, J. C. (1966) A mathematical model of the mitotic cycle and its application to the interpretation of percentage labeled mitoses data. J. natn. Cancer Inst. 37: 443-450. BELLMAN, R., JACQUEZ,J. A. and KALABA,R. (1960) Some mathematical aspects of chemotherapy: I. One-organ models. Bull. Math. Biophys. 22: 181-198. BENDEa, R. A. and DEDmCK, R. L. (1975) Cytokinetic aspects of clinical drug resistance. Cancer Chemother. Rep. 59: 80.~809. BERMAN,M. (1969) Kinetic modeling in physiology. FEBS Lett. 2: $56--$57. BISCHOFF, K. B. (1973) Pharmacokinetics and cancer chemotherapy. J. Pharmacokinet. Biopharmac. 1: 465--48O. BISCnOFF, K. B. (1975) Some fundamental considerations of the applications of pharmacoidnetics to cancer chemotherapy. Cancer Chemother. Rep. 59: 777-793. BlSCHOFF, K. B., DEDRICK, R. L. and ZAHARKO, D. S. (1970) Preliminary model for methotrexate pharmacoidnetics. J. pharmac. Sci. 59: 149--154. BlSCHOFF, K. B., DEDmCK, R. L., ZAHARKO,D. S. and LONGSTRETH,J. A. (1971) Methotrexate pharmacokinetics. 3". pharmac. Sci. 60: 1128-1130. BLUMENSON, L. E. (1973) A comprehensive modeling procedure for the human granulopoietic system: over-all view and summary of data. Blood 42: 303-313. JFT VoL
4, No. 2.--C
278
R . C . JACKSON and K. R. HARRAP
BLUMENSON, L. E. (1975) A comprehensive modeling procedure for the human granulopoietic system: detailed description and application to cancer chemotherapy. Math. Biosci. 26: 217-239. BROCKWELL, P. J. and TRucco, E, (1970) On the decomposition by generations of the PLM-function. J. theoret. Biol. 26: 149-179. BROCKWELL, P. J., MACLAREN, M. D. and TRUCCO, E. (1970) Monte Carlo simulation of PLM-curves and collection functions. Bull. Math. Biophys. 32: 429-443. BRONK, B. V. (1969) On radioactive labeling of proliferating cells: the graph of labeled mitosis. J. theoret. Biol. 22: 468-492. BRONK, B. V., DIENES, G. J. and PASKIN, A. (1968) The stochastic theory of cell proliferation. Biophys. J. 8: 1353-1398. BUELL, J., JELUFFE, R., KALAaA, R. and SRIDHAR, R. (1970) Modern control theory and optimal drug regimens. II. Combination therapy. Math. Biosci. 6: 67-74. BUTCHER, J. C. (1964) Implicit Runge-Kutta processes. Math. Comp. lg: 50-64. CHANCE, B., GAR~NKEL, D., HIC~XNS,J. and HESS, B. (1960) Metabolic control mechanisms. V. A solution for the equations representing interaction between glycolysis and respiration in ascites turnout cells. J. biol. Chem. 235: 2426-2439. CHOU, T. -C., HUTCmSON, D. J., SCHMID,F. A. and PmLIPS, F. S. (1975) Metabolism and selective effects of 1-/3-D-arabinofuranosylcytosine in LI210 and host tissues in vivo. Cancer Res. 35: 225-236. CHUANG, S. -N. (1975) Mathematical models for cancer chemotherapy; pharmacokinetic and cell kinetic considerations. Cancer Chemother. Rep. $9: 827-842. CHUANt, S. -N. and LLOYD, H. H. (1974) Analysis and identification of stochastic compartment models in pharmacokinetics: implication for cancer chemotherapy. Math. Biosci. 22: 57-74. CHUA~t~, S. -N. and LLOYD, H. H. (1975) Mathematical analysis of cancer chemotherapy. Bull. math. Biol. 37: 147-160. COHEN, J. L., JAO, J. Y. and JUSKO, W. J. (1971) Pharmacokinetics of cyclophosphamide in man. Br. J. Pharmac. 43: 677-680. COOPER, (3. J. (1969) The numerical solution of stiff differential equations. F E B S Lett. 2: $22-$29. DANENBERG, P. V., LANGENBACH,R. J. and HEIDELBERGER, C. (1974) Stuctures of reversible and irreversible complexes of thymidylate synthetase and fluorinated pyrimidine nucleotides. Biochemistry 13: 926-933. DEDRICK, R. L , BISCHOFF,K. B. and ZAHAR~O, D. S. (1970) Interspecies correlation of plasma concentration history of methotrexate. Cancer Chemother. Rep. 54: 95-101. DEDI~CK, R. L., FORRESTER, D. D. and HO, D. H. W. (1972) In vitro-in vivo correlation of drug metabolism-deamination of l-~8-D-arabinofuranosylcytosine. Biochem. Pharmac. 21: 1-16. DEDRICK, R. L., FORRESTER, D. D., CANNON, J. N., EL DAREER, S. M. and MELLETT, L. B. (1973a). Pharmacokinetics of l-/3-D-arabinofuranosylcytosine (Ara-C) deamination in several species. Biochem. Pharmac. 22: 2405-2417. DEDRICK, R. L., ZAHARKO,D. S. and LUTZ, R. J. (1973b) Transport and binding of methotrexate in vivo. J. pharmac. Sci. 62: 882-890. DEDRICK, R. L., ZAnARKO, D. S., BENDER, R. A., BLEYER, W. A. and LUTZ, R. J. (1975) Pharmacokinetic considerations on resistance to anticancer drugs. Cancer Chemother. Rep. 59: 795-804. DOMBERNOWSKY,P. and HARTMANN,N. R. (1972) Analysis of variations in the cell population kinetics with tumor age in the LI210 ascites tumor. Cancer Res. 32: 2452-2458. DOSAtHEY, C. E. and DREWINKO,B. (1975) A comphter simulation program for the study of cellular growth kinetics and its application to the analysis of human lymphoma cells in vitro. Comp. Biomed. Res. 8: 118-128.
FRANKE, E. K. (1970) A mathematical model of synchronized periodic growth of cell populations. J. theoret. Biol. 26: 373-382. GARFINKEL, D. (1966) A simulation study of the metabolism and compartmentation in brain of glutamate, aspartate, the Krebs Cycle and related metabolites. J. biol. Chem. 241: 3918-3929. GAR~NKEL, D. (1969) Construction of biochemical computer models. F E B S Lett. 2: $9-S13. GAR~NKEL, D. and HESS, B. (1964) Metabolic control mechanisms. VII. A detailed computer model of the glycolytic pathway in ascites cells. J. biol. Chem. 239: 971-983. GARFINKEL, D., RUTLEDGE,J. D. and HIOGINS,J. J. (1961) Simulation and analysis of biochemical systems. I. Representation of chemical kinetics. Communs Ass. comput. Mach. 4: 559-562. GEAR, C. W. (1971) DIFSUB for solution of ordinary differential equations [D2]. Communs Ass. comput. Mach. 14: 185-190. GIRLIr4C, B. (1969) An introduction to hybrid computers and their application to optimization problems. F E B S Lett. 2: $58-$62. GR1NDEY, G. B., WER~:HEISER, W. C., MORAN, R. G. and NICHOL, C. A. (1970) Anomalous inhibition of growth by combinations of inhibitors of DNA biosynthesis. Fed. Proc. 29: 609. GmSDEY, G. B., MOP.AS, R. G. and W~.RKHEISF-~,W. C. (1975) Approaches to the rational combination of antimetabolites for cancer chemotherapy. In Drug Design, V01. 5, pp. 169-249, Academic Press, New York. GROSER, G. and CLARK, R. (1971) BIOMOD: An interactive graphics system for analysis through simulation. Proc. 1971 I E E E Conference on Decision and Control, Miami Beach, Florida, pp. 197-201. HAns, G. M. (1966) State vector description of the proliferation of mammalian cells in tissue culture. I. Exponential growth. Biophys. J. 6: 275-290. HAns, G. M. (1967) Cellular kinetics, cell cycles and cell killing. Biophysik. 4: 1-14. HAns, G. M. and K ~ N , R. F. (1967) State vector description of the proliferation of mammalian cells in tissue culture. II. Effects of single and multiple doses of ionizing radiations. Radiat. Res. 30: 702-713. HARRAP, K. R. and JACKSON,R. C. (1975) Enzyme kinetics and combination chemotherapy: An appraisal of current concepts. Adv. Enz. Regul. 13: 77-96. HARRIS, P. A. and GRoss, J. F. (1975) Preliminary pharmacokinetic model for adriamycin. Cancer Chemother. Rep. 59: 819-825.
Computer models of anticancer drug interaction
279
HIMMELSTEIN, K. J. and BISCHOFF, K. B. (1973a) Mathematical representations of cancer chemotherapy effects. J. Pharmacokinet. Biopharmac. 1: 51-68. HIMI~LSTEIN, K. J. and BISCHOFF,K. B. (1973b) Models of AraC chemotherapy of LI210 leukaemia in mice. J. Pharmacokinet. Biopharmac. 1: 69--81. HmscH, H. R. and ENOS,BERG, J. (1965) An analog of ceil-culture growth. J. theoret. Biol. 9: 303-312. HIRSCH, H. R. and ENOELBERG, J. (1966) Decay of cell synchronization: solutions of the cell-gTowth equation. Bull. Math. Biophys. 28: 391--409. HOWELL, J. R. (1975) Mathematical formulation for nonuniform multiple dosing. J. pharmac, Sci. 64: 46 A, A.~6. JACKSON, R. C. and H A I ~ , K. R. (1973) Studies with a mathematical model of folate metabolism. Archs Biochem. Biophys. 158: 827-841. JACKSON, R. C. and NIETHAMMER,D. (1977) Acquired methotrexate resistance in lymphoblasts resulting from altered kinetic properties of dihydrofolate reductase. Eur. J. Cancer 13: 567-575. JACKSON, R. C. and WEBER, G. (1976) Enzyme pattern directed chemotherapy. The effects of combinations of methotrexate, 5-fluorodeoxyuridine and thymidine on rat hepatoma cells in vitro. Biochem. Pharmac. 25: 2613-2618. JACKSON, R. C., NIETHAMMER,D. and HUENNEKENS,F. M. (1975) Enzymic and transport mechanisms of amethopterin resistance in LI210 mouse leukemia cells. Cancer Biochem. Biophys. 1: 151-155. JACKSON, R. C., HART, L. I. and HARRAP, K. R. (1976) Intrinsic resistance to methotrexate of cultured mammalian cells in relation to the inhibition of their dihydrofolate reductases. Cancer Res. 36 1991-1997. JACQUEZ, J. A. (19721 Compartmental Analysis Biology and Medicine: Kinetics of Distribution of Tracer-Labelled Materials. Elsevier, Amsterdam. JACQUEZ,J. A., BELLMAN,R. and KALABA,R. (1960) Some mathematical aspects of chemotherapy. IL The distribution of a drug in the body. Bull. Math. Biophys. 22: 309-322. JUSKO, W. J. (1971) Pharmacodynamics of chemotherapeutic effects. Dose-time-response relationships for phase-nonspecific agents. J. pharmac. Sci. 60: 892-895. JUSKO, W. J. (1973) A pharmacodynamic model for cell-cycle-specific chemotherapeutic agents. J. Pharmacokinet. Biopharmac. 1: 175-200. KEMBER, N. F. (1969) Growing bones on the computer--some pitfalls of a computer simulation of the effects of radiation on bone growth. Cell Tissue Kinet. 2:11-20. KINO-SMITH, E. A. and MORLEY, A. (19701 Computer simulation of granulopoiesis: Normal and impaired granulopoiesis, Blood 36: 254--262. KROGER-THmMER, E. (1966a) Formal theory of drug dosage regimens. J. theoret. Biol. 13: 212-235. KROGER-Th'IEMF~, E. (1966b) Pharmacokinetics and dose-concentration relationships Proc. 3rd International Pharmacological Meeting, Sap Paulo, Brazil, pp. 63-113. K~OGER-THIEMER,E. (1969) Formal theory of drug dosage regimens. IL The exact plateau effect. J. theoret. Biol. 23: 169-190. LArrHA, L. G., OLrVER, R. and GU~EY, C. W. (19621 Kinetic model of a bone-marrow stem-cell population. Br. J. Haemat. 8: 442--460. LEME, P. R., CREAVEN, P. J., ALLEN, L. M. and B E ~ , M. (1975) Kinetic model for the disposition and metabolism of moderate and high-dose methotrexate in man. Cancer Chemother. Rep. 59:811-817. LINCOLN, T. and F~mEICH, E. J. (1974) Pharmacokinetic simulation: A future means for better control of cancer chemotherapy. Rec. Res. Cancer Res. 49: 103--107. LINCOLN, T. L., AROESTY, J., MEIER, G. and GROSS, J. F. (1974) Computer simulation in the service of chemotherapy. Biomed. 28: 9-16. MATIS, J. H. and HARTLEY, H. O, (1971) Stochastic compartmental analysis: Model and least squares estimation from time series data. Biometrics 27: 77-102. M c I ~ I S , B., KAPADIA,A., EL-ASFoUIU, S. and Lop, T. L. (1975) Stochastic compartmental modeling of the disposition of 5-(3,3-dimethyl-l-triazeno)imidazole-4-carboxamide. Cancer Chemother. Rep. 59: 843-845. MOP,PaSON, P. F., LINCOLN, T. L. and AROESTY, J. (1975) Disposition of cytosine arabinoside and its metabolites: A pharmacokinetic simulation. Cancer Chemother. Rep. $9: 861-876. NICHOL, C. A., GPaNDEY, G. B., MORAN,R. G. and WERKHEISER,W. C. (1972) Combinations of drugs which interact as inhibitors of DNA biosynthesis. Adv. Enz. Regul. 10: 63-80. NIEBERGALL, P. J., SUGITA, E. T. and SCHN~,aE, R. L. (1974) Calculation of plasma level versus time profiles for variable dosing regimens. J. pharmac. Sci. 63: 100-105. NIETHAMMER,D. and JACKSON,R. C. (1975) Changes of molecular properties associated with the development of resistance against methotrexate in human lymphoblastoid cells. Fur. J. Cancer 11: 845-854. PhDGEa'r, W. J. and TSOKOS, C. P. (1970) A stochastic model for chemotherapy: Computer simulation. Math. Biosci. 9: 119-133. PADOETT, W. J. and TSOKOS,C. P. (1972) A stochastic model for chemotherapy: two-organ systems. Int. J. Bio-Med. Comp. 3: 29--41. PARK, D. J. M. (1974) The hierarchical structure of metabolic networks and the construction of efficient metabolic simulators. J. theoret. Biol. 46: 31-74. PRESCOTT, L. F. (1973) Pharmacokinetic drug interactions: Mechanisms of toxicity. Excerpta Medica Int. Congr. Series No. 288. 14: 49-58. PRIC,OGINE, I. and BABLOYANTZ,A. (1972) Dissipative structures and their application to biology. F E B S Syrup. 25: 3-18. PURDUE, P. (1974) Stochastic theory of compartments. Bull. Math. Biol. 36: 305-309. REICH, J. G. and SEL'KOV, E. E. (19741 Mathematical analysis of metabolic networks. F E B S Lett. 40: S119-S127. REICHARD, P. (1972) Control of deoxyribonucleotide synthesis in vitro and in vivo. Adv. Enz. Regul. 10: 3-15. REINER, J. M. (1969) Behavior of Enzyme Systems (2nd edn) Van Nostrand, New York.
280
R . C . JACKSON and K. R. HARRAP
ROT! ROTI, J. L. and OKADA,S. (1973) A mathematical model of the cell cycle of LSi78Y. Cell Tissue Kinet. 6:111-124. ROWLAND, M. and MATIN, S. B. (1973) Kinetics of drug-drug interactions. J. Pharmacokinet. Biopharmac. 1: 553-567. RUBINOW, S. I. (1968) A maturity-time representation for cell populations, Biophys. J. 8: 1055-1073. SANTI, D. V., MCHENRY, C. S. and SOMMER,H. (1974) Mechanism of interaction of thymidylate synthetase with 5-fluorodeoxyuridylate. Biochemistry 13: 471-481. SHACKNEY, S. E. (1970) A computer model for tumor growth and chemotherapy, and its application to Ll210 leukemia treated with cytosine arabinoside. Cancer Chemother. Rep. 54: 399-429. SHEPPARD, C. W. (1969) Computer simulation of stochastic processes through model-sampling (Monte Carlo) techniques. F E B S Lett. 2: S14-$21. SKIPPER, H. E. (1969) Improvement of the model systems. Cancer Res. 29: 2329-2333. SKIPPER, H. E., SCHABEL, F. M., JR. and WILCOX, W. S. (1964) Experimental evaluation of potential anticancer agents. XIII. On the criteria and kinetics associated with 'curability' of experimental leukemias. Cancer Chemother. Rep. 35:1-111. SKIPPER, H. E., SCHABEL, F. M., JR. and WILCOX, W. S. (1967) Experimental evaluation of potential anticancer agents. XXI. Scheduling of arabinosylcytosine to take advantage of its S-phase specificity against leukemia cells. Cancer Chemother. Rep. 51: 125-165. SKIPPER, H. E., SCHABEL, F. M., JR., MELLETT, L. B., MONTGOMERY,J. A., WILKOFF, L. J., LLOYD, H. H. and BROCKMAN,R. W. (1970) Implications of biochemical, cytokinetic, pharmacologic and toxicologic relationships in the design of optimal therapeutic schedules. Cancer Chemother. Rep. 54:431-450. TAKAHASHI, M. (1968) Theoretical basis for cell cycle analysis. II. Further studies on labelled mitosis wave method. J. theoret. Biol. 18: 195-209. TATTERSALL, M. H. N., JACKSON, R. C., CONNORS, T. A. and HARRAP, K. R. (1973) Combination chemotherapy: The interaction of methotrexate and 5-fluorouracil. Eur. J. Cancer 9: 733-739. TATI'ERSALL, M. H. N., JACKSON, R. C., JACKSON,S. T. M. and HARRAP, K. R. (1974) Factors determining cell sensitivity to methotrexate: Studies of folate and deoxyribonucleoside triphosphate pools in five mammalian cell lines. Eur. J. Cancer 10: 819-826. TEORELL, T. (1937) Kinetics of distribution of substances administered to the body. I. The extravascular modes of administration. Archs int. Pharmacodyn. ThOr. 57: 205-225. THAKUR, A. K., RESCIGNO, A. and SCHAFER, D. E. (1973) On the stochastic theory of compartments. II. Multi-compartment systems. Bull. math. Biol. 35: 263-271. TIBAUX, G., FIRKET, H. and HOPPER, A. F. (1969) A simple computer simulation model for growing cell populations: Analysis of synchronized HeLa cell cultures. Cell Tissue Kinet. 2: 333-345. TRUCCO, E. and BROCKWELL, P. J. (1968) Percentage labeled mitoses curves in exponentially growing cell populations. J. theoret. Biol. 20: 321-337. WEBa, J. L. (1963) Enzyme and Metabolic lnhibitors. Vol. 1. Academic Press, New York. WEINBERG, R. and BERKUS, M. (1971) Computer simulation of a living cell. Part I. Int. J. Bio-Med. Comp. 2: 95-120. WERKHEISER, W. C. (1968) A mathematical model of cellular parameters relevant to chemotherapy with methotrexate (MTX). Proc. Am. Assoc. Cancer Res. 9:77 (1968). WERKHEISER, W. C. (1971) Mathematical simulation in chemotherapy. Ann. N. Y. Acad. Sci. 186: 343-358. WERKHEISER, W. C., GR/NDEY, G. B., MORAN, R. G. and NICHOL, C. A. (1973) Mathematical simulation of the interaction of drugs that inhibit deoxyribon'ucleic acid biosynthesis, Mol. Pharmac. 9: 320-329. WESTLEY, J. (1969) Enzymic Catalysis. Harper & Row, New York. Woo, K. B., BRENKUS,L. B. and Wnc, K. M. (1975) Analysis of the effects of antitumor drugs on cell cycle kinetics. Cancer Chemother. Rep. 59: 847--860. ZAHARKO, D. S., DEDR/CK, R. L., BISCHOFF, K. B., LONGSTRETH, J. A. and OLIVER/O, g. T. (1971) Methotrexate tissue distribution: Prediction by a mathematical model. J. natn. Cancer Inst. 46: 775-784.
Printed in Great Britain
Specialist Subject Editors: ALANC. SARTORELLI,WILLIAMA. CREASEYand JOSEPHR. BERTINO
COMPUTER
MODELS OF ANTICANCER INTERACTION
DRUG
R. C. JACKSON Laboratory for Experimental Oncology, Indiana University School of Medicine, Indianapolis, Indiana 46202, USA
and K. R. HARRAP Department of Biochemical Pharmacology, Institute of Cancer Research, Sutton, Surrey, UK
1. TYPES OF SYSTEM WHICH HAVE BEEN STUDIED USING COMPUTER MODELS Anticancer drugs, administered to a tumor-bearing host, may interact at a number of levels of biological organization. This report will consider four kinds of interaction. Firstly, p h a r m a c o d y n a m i c interactions: this category of interaction refers to mutual effects on the spatial disposition of two or more agents. Thus, one drug may alter the intestinal absorption, tissue distribution, transport across cellular membranes, or renal excretion of other agents. Secondly, inter-tissue metabolic interactions may arise, e.g. some drugs are only active after being chemically modified by hepatic microsomal enzymes, and such drugs would be antagonized by agents which inhibited those microsomal enzymes; conversely if one drug inhibits the conversion of another to an inactive product, then potentiation may occur. Thirdly, intracellular biochemical interactions may be expected. This category differs from the second kind of interaction, in that the latter concerns interactions involving drug activation or inactivation within organs of the host, and these effects modify the response to a drug or drug metabolite at a distant site; by intracellular biochemical interactions is meant the effect that two or more drugs may have on each other within the same cell. Such interactions fall broadly into two kinds: (a) an influence of one agent on activation or inactivation of another, e.g. an adenosine analog may be active in the nucleotide form, but because the nucleotide cannot cross the plasma membrane, it must be administered as the nucleoside; in this case another agent which is an adenosine kinase inhibitor may well be antagonistic for the first agent. (b) The second kind of biochemical interaction occurs when two or more compounds lead to opposing changes in biochemical pathways; for example, arabinosylcytosine (AraC) acts in many cells as the triphosphate form, which inhibits DNA polymerase, and this inhibition is competitively opposed by deoxycytidine triphosphate (dCTP); thus if another drug leads to elevated nuclear dCTP levels, it may serve to antagonize the cytotoxic action of AraC. Fourthly, drugs may interact in relation to their cyte~'inetic effects. If, for example, a particular compound is active specifically against S-phase cells and another agent blocks progression from G~ to S, then the second agent will antagonize the action of the first compound. The progression through the various stages of the cell cycle is, of course, determined by intracellular biochemical events; however, it is convenient to distinguish this category of drug interaction from the third category in that the dimension of time is here a critical feature. Mathematical models have been written to describe drug behavior at all of these levels, and many of these models are considered below. The models, in general, were written with the intention of understanding better the pharmacology of a specific drug in a particular experimental system. Models of the effects, on tumor and/or host, of JPT Vol. 4, No. 2--A
245
246
R. C. JACKSONand K. R. HARRAP
two or more drugs, which may include, or predict, interactions of each drug on the effects caused by the other drug(s), have not been widely used. Clearly other kinds of interactions between antitumor drugs are possible. For example, if one agent is immunogenic and another is immunosuppressive, they will interact at the level of the immune system. Such interactions are excluded from this discussion, not because they are unimportant, but because they have not been described in sufficient quantitative detail to form the basis of mathematical models. 2. THE NEED FOR MATHEMATICAL MODELS OF DRUG INTERACTION The need for these models reflects the fact that there are certain questions that cannot readily be answered without them. The reasons for this follow from the nature of antitumor drug-target interactions. (i) Antitumor drug-target interactions are complex. The biochemical make up of a cell involves many thousands of components. An anticancer drug will typically react at several points, and each such reaction may lead to a cascade of consequential effects. Consider, for example, the antimetabolites 5-fluorouracil (FU) and 5fluorodeoxyuridine (FUdR) (Fig. 1). These are generally considered to act as antineoplastic agents primarily after conversion to the level of the deoxynucleotide, FdUMP, which is a potent inhibitor of thymidylate synthetase. Figure 1 shows that there are at least three possible routes from FU (or FUdR) to FdUMP. Each of the fluorinated intermediates produced will probably influence to some extent the corresponding reaction with the normal uracil or thymine intermediate; these perturbations, in turn, will alter the rates of other reactions. (ii) Biological systems are highly interactive. Consider, for example, the highly simplified diagram of methotrexate (MTX)-FdUMP interactions (Fig. 2). Thymidylate synthetase (TS) is inhibited by FdUMP, but the inhibition is dependent upon methy-
TPi
I FiTP I
I F
P I
[ FUMP
l_
~1
I
FUR
I
I
Pip~~ R-I-
FUdR I
~dR_l)i
I-P
catabolism FIG. 1. Metabolism of 5-fluorouracil and its nucleoside and nucleotide derivatives. FU, 5-fluorouracil; dR-I-P, deoxyribose-l-phosphate; Pi, inorganic phosphate; FUdR. 5fluorouracildeoxyriboside; FUR, 5-fluorouridine; FUMP, 5-fluorouridine monophosphate; FUDP, 5-fluorouridin¢ diphosphate; FdUDP, 5-fluorodeoxyuridinediphosphate; FdUMP, 5-fluorodeoxyuridine monophosphate; FUTP, 5-fluorouridinetriphosphate; FCTP, 5-fluorocytidine triphosphate.
247
Computer models of anticancer drug interaction MTX I I I
,
MTX"/
~
I~lrines
CH2FHI'~ 4
dUMP r
FIt~. 2. The interaction of MTX and FdUMP. FH2, dihydrofolate: FH4, tetrahydrofolate: CH2FH4, N-S,10-methylenetetrahydrofolate; dUMP, deoxyuridine monophosphate; dTMP, thymidine monophosphate; FdUMP, 5-fluorodeoxyuridine monophosphate; MTX, metho-
trexate. lenetetrahydrofolate (CH2FH4), since FdUMP binds to the TS--CH2FH4 complex. The level of CH2FH4 depends upon the activity of dihydrofolate reductase, which is inhibited by MTX. Thus, MTX can modify the efficacy of FdUMP. However, the effectiveness of MTX, in turn, will depend upon the concentration of the competing substrate, FH2, which is produced by thymidylate synthetase. Thus, the presence of MTX influences the effectiveness of FdUMP, and vice versa. Such interactions are particularly complex where feedback regulation is involved; consider the combination of AraCTP with FdUMP. The latter drug leads to lower levels of cellular dTI'P, which would be expected to decrease the activation of ribonucleotide reductase for purine ribonucleoside diphosphate reduction, thus establishing a deficiency of substrates for DNA polymerase (Reichard, 1972). Now dCTP competes With AraCTP for binding to DNA polymerase, and as a competing substrate would be expected to accumulate behind the AraCTP block; this increased dCTP may then activate deoxycytidylate deaminase, leading to an increased concentration of dUMP. This action would reduce the inhibition by FdUMP, and tend to restore the dTTP level. What actually happens to the nucleotide concentrations in situations like this, of course, must be determined for the particular system in question by experiment. When the result is known, the question arises, can this particular result be explained by what is known about the biochemistry of these pathways? Do the known control mechanisms account for the observed behavior? Where the system under observation is highly interactive, human intuition is unable to keep track of the tangle of interrelationships, and a mathematical model becomes most desirable. A further reason why human intuition or simple non-mathematical models become inadequate is that beyond a certain level of complexity biochemical systems may exhibit the phenomenon of 'multistability', i.e. multiple alternative stable steady states may exist, and small changes in the system parameters may lead to abrupt changes in the levels of ~i particular enzyme or substrate; thus, a small quantitative change in concentration of a compound may trigger a qualitative change in behavior of the system. Prigogine and Babloyantz (1972), for example, discuss the induction of/3-galactosidase in E. coll. This 'is an all-or-none phenomenon'. Bacterial populations grown at low concentrations of inducer in a fixed medium are composed of individuals which are in either of two possible steady states; maximally induced, or noninduced. The bacteria can be shifted from one stable state to the other by transitory changes in the environment... the observed experimental c u r v e s . . , show a sudden jump in the quantity of the enzyme, while the inducer concentration changes relatively slowly'. These authors explain the properties of this system in terms of the interactions between the galactoside substrate, the/3-galactosidase enzyme, the galactoside permease, and the lac operon repressor protein. In such a system, where not only the quantitative aspects of control, but also the qualitative behavior, is affected by small differences in
248
R. C. JACKSONand K. R. HARRAP
concentration of nutrient or metabolite, a mathematical model is necessary to analyze and predict the system behavior. (iii) A further reason why mathematical models of drug interaction are desirable arises from the combinatorial nature of drug treatments. Assume that for a particular tumor twenty agents are known which have some cytotoxic effect. Then the number of combinations of four or fewer drugs which could be selected is (20!/16! + 20!/17! + 20!/18!+20!/19!)= 123,520. This makes the assumption that the order of administration is important, which in many instances is known to be the case. If we allow for different dose-levels and varying time intervals between administration of the successive drugs, it is possible to devise tens of millions of protocols of four or fewer drugs from a choice of twenty useful agents. Since, in practice, only a minute proportion of these protocols can be tested, even in experimental systems, it is interesting to consider the logical processes by which the present combination chemotherapy protocols have arisen. For the most part, these have been statistical; if drug A gave a 40 per cent remission rate for a particular tumor, and drug B gave a 30 per cent remission rate, and if the modes of action were completely independent, then A plus B would be expected to give a remission rate of [40+ (100- 40) x (30/100)] per cent = 58 per cent. This reasoning only applies if the modes of action of the two compounds are completely independent, and in view of the highly interactive nature of cellular biochemistry, this is unlikely to be the case. The likelihood of negative interactions can to some extent be minimized by deliberately choosing agents whose mechanisms of action are apparently unrelated, and perhaps for this reason some successful combinations have been designed consisting of, for example, an alkylating agent, an antimetabolite, and a mitotic spindle inhibitor. Other rules of thumb have been suggested, such as the principle of avoiding cumulative toxicity; thus, if drug A is particularly myelotoxic, drugs used in combination with it should preferably exert their limiting host toxicity in a tissue other than the bone marrow. In this way it may be hoped that the drug effects will be additive toward the tumor, but not toward any one host organ or tissue. Thus, it may be claimed that for the most part presently used clinical combination chemotherapy regimens have been devised primarily on statistical grounds using general principles based on pharmacological and clinical commonsense and intuition. Nevertheless, for the reasons stated above, when a system passes a certain level of complexity, human i0tuition, or even explicit non-mathematical models, become tools of uncertain reliability; it is at this point that quantitative models become desirable. Such models, of course, require detailed accurate data concerning the system to be modeled if their predictions are to be of any value. (iv) The remaining reasons for using mathematical models of drug interaction are that sometimes there are, for technical reasons, no experimental means of answering a question, and a computer simulation may in such cases be the next best thing to a real experiment. In this way we can attempt to answer questions about imaginary situations. Consider an example: the antimetabolite MTX, although a potent inhibitor of dihydrofolate reductase, also inhibits thymidylate synthetase, albeit more weakly, in every system so far examined. Some authorities have considered that this secondary effect is an important component of the drug's cytotoxicity, and others have disagreed. If two tumor cell lines existed which were identical in every respect, except that one had a thymidylate synthetase which was susceptible to MTX inhibition, and the other had a MTX-resistant thymidylate synthetase, then the question could be answered directly by considering the comparative effect of the drug on the two tumor lines. Since such an experimental system does not exist this approach is not possible. However, if a mathematical model of the system is written, containing detailed information concerning the kinetic properties of the relevant enzymes, including their inhibition by MTX, the following question can be asked, 'how does removal of the inhibitory effect of MTX on thymidylate synthetase affect the cytotoxicity predicted by the model?'. This kind of question, although imaginary, is not trivial, since answers to questions of this kind provide information which may be used to help in the design of new inhibitors.
Computer models of anticancer drug interaction
249
Biochemistry, as usually practiced, is largely an analytical subject. The various components of cells are purified, and their properties in isolation are studied in great detail. This analytical approach is necessary if we are to understand the nature of drug-receptor interactions, However, the tools of biochemistry have allowed less progress with synthetic aspects. Sometimes it has been possible to combine purified cellular components into functional reconstituted systems, but this 'synthetic biochemistry' is extremely difficult and is still far behind analytical biochemistry. This is an area where computer techniques can make a contribution by examining the possible ways in which the separate components of a system may fit together and interact. 3. PHARMACOKINETIC (PHARMACODYNAMIC) MODELING OF ANTITUMOR DRUGS The invention of the modern compartmental approach to the science of pharmacokinetics is often credited to Teorell (1937). This subject, dealing with the spatial distribution of drugs in different compartments of the body, and the temporal changes in these concentrations resulting from drug uptake, movement between compartments, and elimination from the body, necessarily relies to a great extent on mathematical models, and makes extensive use of computers. For drugs which are water soluble, and slow acting, and where complex nonlinear transport processes can be ignored, it is often found that a plot of log concentration vs time approximates a straight line. The mathematics of systems of this kind were developed by KriigerThiemer (1966a, b). Such models have been used to predict optimal dosage regimens for a number of agents (Kriiger-Thiemer, 1969). Important aspects of pharmacokinetic processes were developed by Bellman et ai. (1960) and Jacquez et al. (1960) who modeled concentration gradients in capillaries. The mathematical techniques used in compartmental analysis are the subject of a detailed monograph by Jacquez (1972). An early model of MTX action, incorporating some pharmacokinetic assumptions, was written by Werkheiser (1971). This model will be discussed in Section 7. The modern approach to detailed pharmacokinetic modeling of anticancer drugs was pioneered chiefly by Bischoff, Dedrick, Zaharko and their co-workers (Bischoff, 1975; Bischoff et al., 1970, 1971; Dedrick et al., 1970, 1973b; Zaharko et al., 1971). These studies simulated the distribution and excretion of MTX, an appropriate choice for such studies because it is not metabolized in man when administered at conventional doses. The compartmental model of Bischof et al. (1971) considered MTX distribution between blood plasma, liver, kidney, muscle and intestine, secretion into bile, absorption into gut, and losses into urine and feces (Fig. 3). Bischoff (1975) in a recent review has discussed such questions as how to decide on the number of compartments to model and what the time-scale of modeling should be, as well as various technical aspects of pharmacokinetic modeling. The quantitative aspects of plasma and tissue binding were also considered. The authors quoted above have elaborated a number of important concepts such as the distinction between local uptake which is limited by drug availability to the region (flow-limited conditions), and local uptake which is determined by the membrane transport properties of the target cells (membrane-limited conditions). The pharmacokinetic consequences of these different assumptions have been worked out in considerable detail, as have those resulting from different degrees of water- and lipid-solubility (Dedrick et al., 1970). Furthermore, the distinction has been made between linear and nonlinear concentration dependence of drug-target binding, and the 'linear binding' approximation has been used in studies of MTX pharmacokinetics (Bischoff, 1975; Bischoff et al., 1970, 1971; Dedrick et al., 1970, 1973b; Zaharko et al., 1971). Dedrick et al. (1975) have extended the MTX model to include such intracellular biochemical events as synthesis of new dihydrofolate reductase molecules. This model will be further considered in Section 7. The pharmacokinetics of AraC presents a very different situation from MTX, in
250
R. C. JACKSONand K. R. HARRAP I
PLASMA
OL-OG
!,
I
/
LIVER
/
~Billllryr.,~retlL
/
G.I. TRACT
¢ ¢ eGutAbeoeptlon Gut Lumen
Ok KIDNEY
1
Urine
OM
FIG.3. Distributionof methotrexatein bodycompartments. Qc, QL, Qx, QM, plasma flowrates to gut, liver, kidney, and muscles, respectively; T, nominal residence time in bile subcompartments; rl,2,3, drug transport rate in bile; cL2.3.4,drug concentration in gut lumen compartments.(FromBischoffet al., 1971.)(Reproduced with permission of the copyright owner.) that (a) AraC is rapidly metabolized by deamination to the inactive arabinosyluracil (AraU), and (b) AraC must be phosphorylated to the triphosphate level to form AraCTP before it can exert its cytotoxic action. Dedrick et al. (1972, 1973a) have presented a pharmacokinetic model for the detailed distribution of AraC and AraU in man. This model employs known enzyme kinetic parameters to describe deamination of AraC to AraU; it does not, however, include the fact that AraC must be converted by phosphorylation to AraCTP. These studies have brought to light a number of interesting aspects of the pharmacology of AraC. Using kinetic parameters for deaminases which had been measured in vitro, the model successfully predicted plasma concentrations of AraC and AraU following injections of AraC; in addition, the important role of 'lean tissue' in acting as a reservoir of AraC and releasing it to the blood was demonstrated. Deamination in human liver was shown to be so rapid that a large fraction of the AraC was inactivated in a single passage through this tissue; the observed kinetics were thus more strongly dependent upon blood flow rate than on the detailed enzyme kinetic properties in the various tissues (Dedrick et al., 1972). Good agreement was obtained between theory and experimental results when models incorporating enzyme parameters for mice, monkeys, dogs and humans were compared (Dedrick et al., 1973a). Lincoln and Freireich (1974) discussed a computer simulation of AraC kinetics; using an interactive program, with display of computer graphics, they simulated the distribution of this drug in various body compartments (i.e. blood, bone marrow, gut and testis) in man. This work made use of a biological modeling program developed at RAND Corporation (Groner and Clark, 1971). Morrison et al. (1975) have recently described a model for pharmacokinetic analysis of AraC which includes some of the intracellular reactions involved in the production and breakdown of AraCTP. This model will be discussed in Section 7. Leme et al.
Computer models of anticancer drug interaction
251
(1975) devised a compartmental pharmacokinetic model for the disposition of moderate- and high-dose MTX in man. When tritiated MTX was administered in large doses, two labeled compartments appeared in the plasma. One of these was MTX, and the other, which was not exchangeable with MTX, was suggested by these authors to be the result of biliary excretion ~of MTX, degradation of the drug by the intestinal flora and subsequent resorption of these degradation products. It was found that a model with nonlinear inter-compartmental transfer gave better agreement with observed results than a simple linear compartmental model. Recent developments in pharmacokinetic modeling have extended the BischoffDedrick model to newer antitumor agents. Harris and Gross (1975) presented a preliminary pharmacokinetic model for adriamycin. This model has ten compartments, including heart and bone marrow, the dose limiting sites of host toxicity for adriamycin. The major site of drug clearance was considered to be the liver, with renal clearance being negligible. A constant S0 per cent of the adriamycin present in the plasma was assumed to be protein-bound. Tissue-plasma concentration ratios were measured experimentally for rabbit tissues at three time points following intravenous drug administration, and the mean values, assumed to be time-independent, were used in the simulation. The model was used to predict adriamycin concentrations in the ten tissues at times of up to 50hr after intravenous administration of adriamycin to rabbits. The agreement between predicted and measured tissue levels was good. The model was then used to predict human plasma levels, which again gave good agreement with measured values. Allen and Creaven (1975) have described a human pharmacokinetic model for isophosphamide. This agent, like the closely related drug cyclophosphamide, must be metabolically activated. This model is discussed in Section 4. Other recent developments in the techniques of pharmacokinetic modeling allow the systems to accommodate multiple unequal drug doses (Niebergall et al., 1974). A further extension of such models allows nonuniform time spacing of the doses (Howell, 1975). The preceding models have been 'deterministic', i.e. they assume that the rate constants and other system parameters are known with precision, and are constant within a particular compartment or subcompartment of the system. Real systems, in fact, are so complex that certain of their properties appear to exhibit random behavior. Some authors, consequently, have suggested that a more realistic treatment of these systems would be to use stochastic, rather than deterministic, assumptions at least in parts of the models. In such systems, the rate constants or other system parameters are treated as random variables. Aspects of stochastic compartmental modeling were discussed by Matis and Hartley (1971), Thakur et al. (1973) and Purdue (1974). Padgett and Tsokos (1972) have discussed the theory of a stochastic model for chemotherapy for a two-organ biological system, and have described computer simulations for a one-organ model (Padgett and Tsokos, 1970). Chuang and Lloyd (1974) have discussed the theory of two-compartment and three-compartment models. These authors conclude: 'A mean of the concentration process is obtained which in effect is a best-fit mean concentration-time ( C - T ) curve for the measured compartment of the group of individuals under observation. It is not possible to predict the C - T curve for any particular individual in a population from the stochastic analysis, but it is possible to predict with some degree of confidence the expected range for tl~e C - T curve. In addition, expected stochastic C - T curves for the unobserved or unobservable compartments may be predicted'. McInnis et al. (1975) described a two-compartment stochastic model of the disposition of the anticancer drug 5-(3,3dimethyl-l-triazeno)imidazole-4-carboxamide (NSC-45388); the predicted disappearance of NSC-45388 from plasma in a dog was compared with experimental data. All of the pharmacokinetic models discussed above have been attempts to describe the compartmental and temporal changes in concentration of a single drug, or of one drug and its metabolites. In general, drug-drug interactions have not been considered, although it is possible that some of these models could be modified to take such interactions into account. A few theoretical approaches to the pharmacokinetics
252
R.C. JACKSONand K. R. HARRAP
of drug combinations have appeared. Buell et al. (1970) have shown that the concepts of control theory may be applied to the computation of multidrug dosage schedules. The problem of drug administration was treated as a discrete-time optimal control problem, and methods for solving this problem by dynamic programing were shown. Preliminary pharmacokinetic studies of drug-drug interactions have appeared, although these do not refer to antitumor drugs (Prescott, 1973). Rowland and Matin (1973) discuss a number of drug combinations for which a careful kinetic analysis has clarified the site of drug--drug interaction. Thus the plasma half-life of the oral hypoglycemic agent tolbutamide is considerably lengthened by sulfaphenazole; this is accomplished by blockage by the latter agent of the hepatic microsomal oxidation of the tolbutamide. These authors also discussed the interaction between the anticoagulant, warfarin, and phenylbutazone. Warfarin is bound to a large extent to plasma protein, but it is, apparently, the unbound fraction which exerts the pharmacological effect. Phenylbutazone displaces bound warfarin and thus augments the anticoagulant effect. In addition, by selective inhibition, it alters the isomeric composition of the warfarin (administered as a racemic mixture) in plasma, hastening the metabolism of the less active isomer, and delaying that of the active form. Such studies have not been carried out with antitumor drugs. However, in view of the importance of such drug interactions, and of the relatively great advances in antitumor drug pharmacokinetics in recent years, it may be hoped that such studies will be conducted in the very near future. 4. MODELS INVOLVING DRUG METABOLISM AT A SITE DISTANT FROM THE TUMOR It is convenient to distinguish between (a) drugs which require no metabolic activation, such as MTX, (b) drugs which are activated in the target cell, such as AraC, and (c) drugs which are activated at some site distant from the tumor, usually the host liver, such as cyclophosphamide. The first class of drugs is the most amenable to pharmacokinetic modeling, as the various successful models of MTX pharmacokinetics have clearly shown. The second class to be modeled in detail requires a knowledge of the intracellular reactions involved in drug activation. Many of the drugs in this category are purine and pyrimidine nucleoside analogs, which exert cytotoxicity at the nucleotide level. The nucleotide forms of these drugs are unable to effectively cross the plasma membrane, and must therefore be administered in forms which are activated in the cell in which they are to act. Models involving anticancer agents of this type are discussed in Sections 5 and 7. Drugs which must be metabolically activated in the liver, or in some other body site, present certain complications to experimental study. For example, these agents cannot be tested simply against tumor cells in culture; this problem has been circumvented in the case of cyclophosphamide by adding isolated hepatic microsomes to the cultures. These complications extend to mathematical modeling. A satisfactory model of a drug of this kind must simulate the disposition of both the original compound, and the metabolically active form in the body, as well as the activation process itself. Such models, therefore, inevitably combine both pharmacokinetic and intracellular biochemical features. What appears to be the first model of this kind to describe the behavior of an anticancer drug, cyclophosphamide, was presented by Cohen et al. (1971). Biotransformation (activation) was assumed to take place within the large central compartment, 'since the relatively large volume of distribution is likely to encompass the liver where oxidation of cyclophosphamide occurs'. Metabolic activation was not modeled in detail, but was treated as first-order with respect to the concentration of cyclophosphamide. This would approximate the situation if the concentration was always well below the Km of the activation process. The model was used to estimate parameters for a total of eleven doses given to seven patients. The authors noted that barbiturates enhanced the activity of the microsomal drug metabolizing enzymes of the liver. They did not explicitly model this activation, but they noted that one of their
Computer modelsof anticancerdrug interaction
253
patients who had received phenobarbital appeared to have metabolized a larger fraction of the administered dose than the other patients. A more recent and detailed model of this kind is that of Allen and Creaven (1975). This study modeled the pharmacokinetics of isophosphamide (NSC-109724), an analog of cyclophosphamide in which one of the /3-chloroethyl groups is replaced by hydrogen, and is transferred to the 3-position of the oxazaphosphorine ring. Isophosphamide is active against some tumors which are resistant to cyclophosphamide, and also differs from cyclophosphamide in that the host toxicity is directed principally at the kidney, rather than the bone-marrow (Allen and Creaven, 1975). Like cyclophosphamide, isophosphamide requires activation by the hepatic endoplasmic reticulum. The model of Allen and Creaven is primarily a pharmacokinetic model, but includes two parameters (K,, and Vm~) to describe the drug activation (Fig. 4). This model indicated that the excretion rate constant for drug metabolites was larger than for the unchanged drug. Furthermore, the activated drug appeared to distribute through a much smaller volume than unchanged drug, indicating limited permeability into many parts of the body. The authors report that the model is currently being used to predict plasma levels of isophosphamide after repeat doses of the drug, and also to assess the influence of other drugs on the pharmacokinetics of isophosphamide, though these studies are not yet published. Perhaps the most complete model of drug-drug interactions involving drug metabolism is the study of Rowland and Matin (1973), already mentioned in Section 3. In their detailed model of the interaction of sulfaphenazole with tolbutamide they successfully simulated the effect of sulfaphenazole in increasing the plasma half-life of tolbutamide from its normal value of approximately 7 hr to about 44 hr. This situation differs from that with cyclophosphamide in that tolbutamide is itself active, and the oxidized drug formed in the hepatic endoplasmic reticulum is inactive. This interesting model will not be discussed further here, since it was not concerned with antitumor drugs. Nevertheless, it is likely that similar studies of this kind, on the pharmacokinetics of drug-drug interactions at the level of hepatic activation and inactivation of anticancer drugs, will appear before too long. •
x----O
DOSE
@
FIG. 4. Multicompartmentpharmacokineticmodel for isophosphamide.X~ and X4 are central compartments;X3 and X6, urine; X2 and X~, peripheral compartmentsfor isophosphamideand metabolite(s), respectively; V.... k~e and kE~ are compartmentsfor metabolism,where E is the substrate-enzymecomplex; kx is the transfer rate constant (hr-~) for irreversible biologic alkylation or drug decomposition; kii (i = 1, 2, 3 . . . . n, j = 1, 2, 3 . . . . . m ) are transfer rate constants (hr-~) for drug transport into compartment j from compartment i; V~ and V~, are volumes of the central compartmentfor isophosphamideand metabolite(s)respectively.(From Allen and Creaven, 1975.)
5. INTRACELLULAR BIOCHEMICAL INTERACTIONS In distinguishing between interactions at the cellular level, and those occurring at a site distant from the tumor, as discussed in the previous section, we are following the distinction made by Grindey e t a l . (1975) between 'intracellular metabolic interactions' and interactions at the level of 'remote sy,stemic biotransformation'. A few rather detailed computer models of intracellular biochemical interactions of anticancer drugs
R. C. JACKSONand K. R. HARRAP
254
have been published. Study of these interactions is facilitated by the ability to conduct experiments with cell cultures, where extraneous variables such as extracellular drug concentration, and nutritional status may be much more rigorously controlled than in animals. Thus, the predictions of models of cellular biochemistry may be subjected to detailed testing. Mathematical models at the intracellular level generally are collections of rate equations describing the interconversion of substrates and cofactors in enzyme-catalyzed reactions. Such models are thus inherently nonlinear. Detailed computer modeling of biochemical control mechanisms in multienzyme pathways was first done by Garfinkel, Hess," Chance and their collaborators (Chance et al., 1960; Garfinkel et al., 1961; Garfinkel and Hess, 1964). These workers have described a detailed model of the reactions of glycolysis in ascites cells. This model accurately predicted the effects of such perturbations as addition of glucose or of the uncoupling agent dibromophenol on the concentrations of glycolytic intermediates. With models of this kind, which are intended to simulate collections of chemical reactions, it is axiomatic that a detailed knowledge is required of the relevant reactions, their stoichiometry, and the rate equations. In the case of enzyme-catalyzed reactions it is necessary to know the kinetic mechanism of the reaction to write the rate equation. These principles were first applied to the study of the intracellular interactions of antitumor drugs by Werkheiser, Grindey, Moran and Nichol (Werkheiser, 1968, 1971; Grindey et al., 1970, 1975; Nichol et al., 1972; Werkheiser et al., 1973). These papers have presented a number of different models; one of them contained features of both biochemical and pharmacokinetic modeling, and this will be discussed in Section 7. The model of DNA synthesis was an intracellular model; its features are summarized in Fig. 5. The assumptions underlying this simulation were as follows: (i) the system being modeled was the L1210 mouse leukemia cell growing in culture; (ii) the cell was treated as an open system in a steady state, with the medium acting as a constant source of precursors, and the growing DNA chain as a sink; (iii) within this simplified system the transfer of material from the source to the sink was as described in Fig. 5, the rates being determined by enzymes which obeyed unidirectional MichaelisMenten kinetics; (iv) 5-FdUMP, AraCTP and AraATP were treated as competitive inhibitors of their respective target reactions, and IQ-1 (1-formylisoquinoline thiosemicarbazone) and MTX were treated as noncompetitive inhibitors of ribonucleotide reductase and thymidylate synthetase, respectively. The details of the rate equations may be found in the original paper (Werkheiser et al., 1973). Using this model, the
RC
i .....!'" " "~ / •"""
i
/
RU ! ...................
D
:
\
....................................................... i
RA,
RG
:
i
.,,iJ .,~I
/
m, A
I
i
\ ,s
J
Fro. 5. Model of metabolic regulation of DNA synthesis. RC, CDP; RU, UDP; ] ~ , ADP; RG, GDP; C, dCTP; U, dUMP; T, dTTP; A, dATP; F, dGTP; D, DNA; feedback inhibition, dotted lines; feedback activation, dashed lines. (From Grindey et al.. ]975.)
Computer models of anticancer drug interaction
255
Roswell Park group simulated the interactions of binary combinations of M T X , AraC, AraA, IQ-I and FUdR. Some of the results plotted as isobolograms are shown in Fig. 6, which compares the computer predictions with the observed effects of the drug pairs on L1210 cell growth in vitro. A number of important conclusions emerged from the various studies conducted with this model: (i) in general, in pathways with multiple feedback inhibitions,the behavior of the pathway when subjected to multiple inhibition could not be predicted from the sequential, concurrent or complementary relationship of the inhibitors, but was a function of the structure of the system as a whole; (ii)the types of drug interactions were relatively independent of the system parameters, and primarily depended upon the structure of the model. Occasional exceptions were noted to this generalization, and these will be discussed further below; (iii)the model was highly interactive, i.e.the presence of an inhibitor caused changes in the concentrations of all the intermediates. The effects of a perturbation in one region were not localized, but were manifested throughout the entire system; (iv) the feedback regulation led to results which could not have been intuitively foreseen• Inhibition of a pathway could result in increased, decreased, or unchanged concentration of the product of that pathway; and (v) the intensity of interaction between two drugs was related to the degree of inhibition of the system, A drug combination which was mildly synergistic or antagonistic at 50 per cent inhibition became intensely so when the system was inhibited by 90 per cent. The model of D N A synthesis was used to predict the response of deoxyribonucleoside triphosphate pools to antimetabolite chemotherapy. Some of the results were those which would have been intuitively predicted; thus, IQ-I, the ribonucleotide reductase inhibitor, caused decreased levels of all deoxyribonucleoside triphosphates. However, many of the predictions could not have boon made without the model: for example, AraA and AraC were both predicted to lead to increases in dCTP, d A T P and dTTP, and decreases in uridine nucleotides. In some cases the predictions describe transient behavior as well as steady-state solutions. The effect of MTX, for example, was predicted to cause an increase in uridine and cytidine nucleotides, with, initially, a sudden fall in d T r P and a rise in dATP; the feedback properties of the model were 1.0~
IX'.
°1 -
1.0~
1.0~
I",,
0.5 10-1 R~s
~.0
•
~
IX:"
,\J
0.5 ARA-C Cc
1 , 0 ~
1.0
•
~
1.
0.~ ARA-C Cc
1.0
Oo •
go o:s
1".o ~
IQ-1 Rn
_
II." . G'
0.5
i.o
ARA-A Ac
= 0~5 10-1 Rn
~'1.O
el
o'.s
1.o
ARA-C C©
/ G"
O15 ARA-A AC
I~O
C"
0.5 ARA-C CC
1.0
FIG. 6. Isobolograms showing interactions of pairs of dr~]gs on the growth of LI210 cells in vitro (points) compared with the effects on DNA synthesis as predicted by the model (solid lines). (From Werkheiser etal., 1973.)
256
R. C. JACKSON and K. R. HARRAP
such, however, that the dTTP level later returned to the normal value, with a concomitant drop in dATP, which became the limiting nucleotide in the steady state. Werkheiser et al. (1973) conclude: 'No claim is made that this model represents the real system in any of its details, but it seems likely that the principles of its construction are reasonable'. Another model described by the same authors (Grindey et al., 1975) was used for a simple simulation of folate metabolism. This model treated dihydrofolate reductase and thymidylate synthetase as separate reactions, and was therefore able to model more accurately the interactions of MTX and 5-FdUMP. A limitation of this model was that simple equations for competitive inhibition were used; in these instances, the rate equations for tight-binding competitive inhibition would have been more appropriate. This model also allowed for simulation of folate synthesis, and the inhibition of that process; thus, the interactions of sulfonamides with MTX or 5FdUMP could be studied. Mammalian cells, unlike many microorganisms, are unable to synthesize folate; however, the model could be used to simulate uptake of preformed folate from the medium. The model predicted that an inhibitor of dihydrofolate synthesis or uptake would be synergistic with an inhibitor of dihydrofolate reductase. It also predicted that MTX and FUdR would be antagonistic, and the combination of sulfonamides with FUdR would be additive. These predictions are consistent with the observed biological behavior of such combinations. These models of Werkheiser et al. in general represent minimal models of the systems they attempt to simulate. That is, they are an attempt to produce the simplest possible mathematical description of the systems under study which still retains the essential control features of the parent systems. These models are thus far less sophisticated than those of, for example, Garfinkel, in which scores of individual rate processes are simultaneously studied. However, especially where the details of a system are not well understood from the experimental data, simplicity is often desirable in a model. These particular simulations have realized the hopes of their authors in accounting for many of the observed interactions between anticancer drugs, and in suggesting profitable experiments. They conclude: 'Simulations that do not lead to new experiments in the laboratory have little value; there should be a continuing interplay between model and experiment'. The introduction of these particular models represented a major advance in the study of the biochemical pharmacology of antitumor drugs. An intracellular biochemical model of some reactions of folate and purine nucleotide metabolism has been described by Jackson and Harrap (1973). This model is limited in range, in that it only describes a small area of cellular metabolism (Fig. 7), but the part described is given in detail. The individual rate equations for 23 reactions are simultaneously integrated. When kinetic mechanisms had been studied in sufficient detail, the appropriate rate equations were used. Where the kinetic mechanism was unknown, it was treated as rapid equilibrium random. For tight-binding inhibitors such as aminopterin, the Goldstein equation was used. This model can therefore be employed to study the effects of inhibitors on steady-state fluxes and pool sizes, and also examine the transient phase following addition of inhibitors or precursors. As with the models discussed above, it is assumed that the nutrient medium acts as a source of precursors, and the exponentially growing cells act as a sink. The uninhibited system is in open steady state; following perturbation a transient phase occurs which may or may not lead to another steady state. The system (as described in Fig. 7) possessed a considerable degree of inherent stability. For example, if dihydrofolate reductase was partially inhibited, the substrate dihydrofolate accumulated until the increased substrate concentration exactly compensated for the degree of inhibition, after which the steady-state rate remained unchanged. If, however, the dihydrofolate reductase was inhibited severely enough to make it the rate-limiting enzyme in the cycle of dihydrofolate oxidation and reduction, then the system assumed a new steady state, in which the flux through the cycle was diminished. Some of the observed effects were explainable on the basis of existing knowledge of the control of multienzyme systems (Webb, 1963; Reiner, 1969; Westley, 1969). Other effects were
Computer models of anticancer drug interaction
.._JJ
.......
257
A,cA\\
1 FIG. 7. Diagram of reactions simulated by a mathematical model. AICAR, 5-amino4-imidazole carboxamide ribotide; CH:FIt4, NS,10-methenyltetrahydrofolate; CH2Ftt,, N5,10methylenetetrahydrofolate; 5-FdUMP, 5-fluorodeoxyuddylate; FI-I2, dihydrofolate; FH4, tetrahydrofolate; GAR, glycinamide ribonucleotide; 10-HCOFH4. N 10-formyltetrahydrofolate; MTX--methotrexate. The rate processes simulated are as follows: vl, dihydrofolate reductase (EC1.5.1.3): v2, serine hydroxymethyltransferase (EC 2.1.2.1); v3, thymidylate synthetase (EC 2.1.1.4.5); v4, serine hydroxymethyltransferase (back reaction); vs, NS,10-methylenetetrahydrofolate dehydrogenase (EC 1.5.1.5); v6, NS,10-methenyltetrahydrofolate cyciohydrolase (EC 3.5.4.9); v~, NS,10-methenyltetrahydrofolate cyclohydrolase (back reaction); 08, NS,10-methenyltetrahydrofolate:glyeinamide ribonucleotide transformylase (EC 2.1.2.2); vg, Nl04ormyltetrahydrofolate:5-amino-4-imidazole carboxamide ribonucleotide transformylase (EC 2.1.2.3); v~0,carrier-mediated cellular uptake of MTX; v., nonenzymic combination of tetrahydrofolate and formaldehyde; v~2, nonenzymic decomposition of NS,10-methylenetetrahydrofolate; v~3, cellular uptake of MTX by passive diffusion; v~4, efllux of MTX by passive diffusion; vl~, biosynthesis of serine; v,6, serine dehydratase (EC 4.2.1.13); v~0,carrier-mediated efflux of MTX; v2,, biosynthesis of deoxyuridylate; v~, biosynthesis of glycinamide ribonucleotide. Rate processes considered by the model but not shown in Fig. 7 are: v~7,biosynthesis of dihydrofolate reductase; v~8,cellular degradation of dihydrofolate reductase; v,9, cellular degradation of dihydrofolate reductaseMTX complex. (From Jackson and Harrap, 1973.) not predictable f r o m the behavior of earlier, simpler models of sequential or cyclic pathways, and a p p e a r e d to be a consequence of the complex t o p o g r a p h y of the model, which contains three overlapping cycles and m a n y branch points. Initially, simulations were run to estimate the effect of changes in p a r a m e t e r s which could determine resistance to MTX. The model, like real L1210 cells, b e c a m e less sensitive to M T X if m e m b r a n e permeability was altered. If the change in m e m b r a n e permeability was the result of an increase in the Kt value of the transport protein, the reduced permeability could be o v e r c o m e by higher concentrations of extracellular drug, but if the change was in the Vmax for m e m b r a n e transport of drug, the mutation remained effective e v e n at v e r y high extracellular levels of drug. Both types of mutation have b e e n found in cultured cells (Jackson et al., 1975; N i e t h a m m e r and Jackson, 1975). An increase in the cellular dihydrofolate reductase concentration in the model resulted in M T X resistance, as is so often o b s e r v e d experimentally. The model simulated the synthesis and turnover of dihydrofolate reductase, and either an increased rate of synthesis or reduced rate of degradation of e n z y m e increased the steady-state level. Since M T X appears to stabilize dihydrofolate reductase against turnover, it increased the cellular dihydrofolate reductase concentration. Of course, the additional e n z y m e with M T X bound to it was inactive, but the inhibited e n z y m e was in equilibrium (or close to it) with free e n z y m e and free drug, and the model showed that under some circumstances the p r e s e n c e of M T X could actually increase the concentration of free (active) dihydrofolate reductase in the cell. Simulation of the effects of M T X on the
R. C. JACKSON and K. R. HARRAP
258
activity of serine hydroxymethyltransferase and thymidylate synthetase showed that normally inhibition of these enzymes was not a major contributor to the effect of the drug. More recent studies (Jackson and Niethammer, 1977) indicate that in extremely resistant cells, in which the dihydrofolate reductase concentration is increased 300fold or more, free MTX may accumulate to sufficient levels to make thymidylate synthetase inhibition the rate-limiting component of MTX action. When this state is reached, leucovorin no longer completely rescues the cells, but only partially protects them from MTX. This surprising result, predicted first by the computer model, has been confirmed with high dihydrofolate reductase containing human cells in culture; the partial protection by leucovorin is due to its conversion to 5,10-methylenetetrahydrofolate which then competes with MTX for the folate binding site on thymidylate synthetase. Since these mutant lines with huge reductase contents have never been described in vivo, the practical importance of this effect is unknown. Another factor predicted to alter cellular response to MTX was the dissociation constant of the cell's dihydrofolate reductase for the inhibitor. This result was also unexpected, since this dissociation constant is in all cases very low, so that the inhibition is often described as 'stoichiometric'. When the titration curves for such inhibitors are plotted, the early part of the curve corresponds to 1:1 inactivation of enzyme molecules by inhibitor molecules, regardless of the exact Ki value. Thus, it seemed that the exact value of the latter parameter was unimportant, so long as it was sufficiently low to keep inhibition 'stoichiometric'. Simulation showed the fallacy in this reasoning. Because dihydrofolate reductase was present in great excess over thymidylate synthetase, the rate-limiting enzyme in the 'dihydrofolate cycle', it was necessary to inhibit it by 95 per cent or more in order to make it rate-limiting, and thus result in MTX-induced growth inhibition. However, at 95 per cent inhibition, the enzyme-inhibitor titration curves are no longer linear, but instead show a 'toe'. In this region the exact value of the Ki makes a great deal of difference to the observed reaction rate. Recent studies (Jackson et al., 1976) have confirmed in four cultured mammalian cell lines (with similar total activities of dihydrofolate reductase, similar rates of drug transport, similar total folate cofactor pools, and similar dTTP levels) that the sensitivity to MTX correlated inversely with Ki value for the respective dihydrofolate reductases and MTX. Table 1 compares the computer predictions of the inhibitory concentrations of the four cell lines, and their predicted relative resistance, with the observed values. The predictions were based on measured values for the relevant enzyme activities and kinetic parameters, the appropriate rate equations, and the assumption that a cell would divide when its content of both d'ITP and purine nucleotides had doubled. Considering the uncertainties inherent in this kind of prediction, the agreement is excellent. Here again, experiments confirmed the predictions of the model. A simpler model, in which inhibition of dihydrofolate reductase by MTX was treated as stoichiometric, could not, of course, have made this prediction. As well as being used to examine the effect of system parameters on response to single drugs, the model was used to predict drug interactions. The drugs simulated were antifolates, thymidylate synthetase inhibitors, glutamine analogs, and purine TABLE 1. Concentrations of M T X Required to Inhibit Cell Growth by 50 per cent as Predicted by Metabolic Simulation. The Appropriate Kinetic Parameters for Dihydrofolate Reductase were Used, but other Factors were Assumed to be Constant
Cell line
K~ (pM)
LI210 W1-L2 L5178Y Yoshida
5 7 29 190
Predicted Resistance* 1I)5o (riM) (-fold) 14.2 18.3 72.6 409
(1) 1.3 5.1 28.8
Observed Resistance* IDso (riM) (-fold) 5.8 13 56 390
(1) 2.2 9.7 67.2
*Resistance is relative to the predicted or observed value for LI210.
Computer models of anticancer drug interaction
259
antimetabolites. When the combination of MTX + FUdR was simulated, the results showed some interesting similarities with the predictions of the model of Grindey et al. (1975) for the same combination, and an important difference. When the kinetic parameters for L1210 mouse leukemia cells were used in the model, it predicted that the combination would be antagonistic, in agreement with the model of Grindey and his co-workers. However, when the parameters for Yoshida sarcoma cells were substituted, the model predicted that the two agents would be almost additive. The reasons for this difference formed the basis for further study with the model, and for a series of experimental studies. If the dihydrofolate oxidation-reduction cycle is considered alone, it is clear the MTX and 5-FdUMP cannot be additive, since one or another of the two enzymes dihydrofolate reductase and thymidylate synthetase will be rate-limiting, and inhibition of the non-rate-limiting enzyme will have little effect on the flux through the cycle. This was the prediction of the model of Grindey et al., which did not take into account the inhibitory effects of MTX on purine nucleotide synthesis. The model of Jackson and Harrap (1973), which simulated the involvement of folates in the pathway of de novo purine nucleotide biosynthesis made a similar prediction for L1210 leukemia cells. Moreover, it predicted that in the L1210 cell, the growth inhibitory effect of MTX was primarily a consequence of a limitation in the availability of purine nucleotides. In the slower-growing Yoshida sarcoma cell, the MTX effect was primarily the consequence of decreased thymidylate production. The balance between the anti-thymidylate effect of MTX and the antipurine nucleotide effect was primarily determined by the competition for 5,10-methylenetetrahydrofolate b y the enzymes thymidylate synthetase and 5,10-methylenetetrahydrofolate dehydrogenase. Thymidylate synthetase requires substrate amounts of the tetrahydrofolate cofactor, since this material acts not only as a one-carbon group source, but also as a reducing agent. The folate cofactor is released from the synthetase as dihydrofolate, which must then be reduced by dihydrofolate reductase, the target enzyme of MTX. 5,10-Methylenetetrahydrofolate dehydrogenase, on the other hand, leads to purine nucleotide biosynthesis, a process in which the two folate-dependent transformylase enzymes utilize the one-carbon groups bound to the tetrahydrofolate cofactors, and release the cofactors as free tetrahydrofolate, which does not require reduction; thus, purine nucleotide biosynthesis requires catalytic amounts of tetrahydrofolate. A consequence of this phenomenon is that the rate at which tetrahydrofolates are depleted in MTX-treated cells is dependent upon the competition for 5,10-methylenetetrahydrofolate by the two enzymes. In rapidly dividing cells, such as L1210 leukemia cells, the ratio of thymidylate synthetase to 5,10-methylenetetrahydrofolate dehydrogenase is higher than in slower-growing cells (Tattersall et al., 1974; Jackson and Weber, 1976). Thus, the model predicted that the ratio of these two enzymes would determine whether the effect of MTX was exerted primarily on the biosynthesis of purine nucleotides or of thymidylate. Experimental work supported this prediction. Thus, in the LS178Y lymphoma, where the ratio of synthetase to dehydrogenase was high, MTX caused decreases in both dTTP and dGTP levels, and in the presence of added thymidine, MTX-treated cells died of purine nucleotide starvation. On the other hand, in the Yoshida sarcoma, where the ratio of synthetase to dehydrogenase was lower, MTX toxicity was partly reversed by thymidine, and only dTTP was decreased among the deoxyribonucleoside triphosphate pools (Tattersall et al., 1974). This discovery had implications for combination chemotherapy, since 5-FdUMP, as an inhibitor of thymidylate synthetase, altered the effective ratio of the enzymes competing for 5,10-methylenetetrahydrofolate, and thus antagonized the antipurine effect of MTX. Thus, the antagonism between the two drugs was most marked in those cells in which the antipurine effect of MTX formed a significant part of its mechanism of action, i.e. in rapidly growing cells. In slowergrowing cells, where MTX functions primarily as an anti-thymidylate agent, the antagonism between MTX and 5-FdUMP is less pronounced. This prediction was recently confirmed in a series of rat hepatoma cell lines in vitro, where MTX and FUdR were strongly antagonistic in the rapidly growing Novikoff tumor, and only
260
R.C. JACKSONand K. R. HARRAP
slightly less than additive in the relatively slowly growing hepatoma 8999S. The model indicated that the interactions between MTX and 5-FU derivatives were more complex than that of solely a 5-FdUMP-antagonizing effect on the antipurine activity of MTX. MTX administration resulted in a very large cellular accumulation of dUMP which competed with 5-FdUMP for binding to thymidylate synthetase. This predicted accumulation of dUMP in the presence of MTX was confirmed by experiment (Tattersall et al., 1973). Thus, not only does FdUMP antagonize the antipurine effect of MTX, but MTX antagonizes the binding of FdUMP to thymidylate synthetase. Furthermore, the sequential sites of action of the two agents leads to the prediction that this combination would be infra-additive. Finally, it has been shown (Santi et al., 1974; Danenberg et al., 1974) that binding of dUMP or FdUMP to thymidylate synthetase is dependent upon the prior binding of 5,10-methylenetetrahydrofolate; when this folate cofactor is depleted, as in the presence of MTX, FdUMP is unable to bind to the synthetase molecule. In this situation the normal substrate dUMP is also unable to bind, and thus the effect would be a neutral one, were it not for the irreversible nature of the FdUMP inhibition. This means that FdUMP permanently reduces the concentration of active enzyme molecules, but requires the presence of the tetrahydrofolate cofactor to do so. MTX, then, interacts with FdUMP in at least four ways, but in all four cases the effect works in the same direction, tending to cause mutual antagonism, to an extent which is most severe in cells containing the enzyme patterns characteristic of rapid proliferation. Drug--drug interactions are even harder to predict when the multiple effects work in opposing directions. The folate and purine model was used to unravel multiple effects of this kind in the interactions between 6-mercaptopurine (6-MP) and MTX or FUdR. 6-MP, after conversion to its ribonucleotide, thioinosinate, inhibits IMP dehydrogenase, the rate-limiting enzyme of guanine nucleotide biosynthesis, and, more weakly, adenylosuccinate synthetase, which is the first enzyme of the pathway leading to the adenine nucleotides. These effects lead to a depression of the ribonucleotide pools, particularly the guanylate pools. GMP and AMP are feedback inhibitors of the first enzyme of the de n o v o biosynthetic pathway of purine nucleotides, PRPP amidotransferase. Thus, a decrease in GMP concentration tends to activate de n o v o purine nucleotide biosynthesis, and antagonize the antipurine effect of MTX. On the other hand, thioinosinate itself is a pseudo, feedback inhibitor of PRPP amidotransferase, so this action opposes the first effect. Which action predominates can readily be determined with a computer model, since instructions can be given to ignore the inhibition of PR PP amidotransferase by thioinosinate. The prediction by the model for the interaction of MTX and 6-MP is shown in isobol form in Fig. 8, together with measured results for L1210 cells in culture; both predicted and observed results correspond to therapeutic summation. Similar results were obtained for the combina1.C
"-...... ~ \
0.~
"%...
0.6
"'"""
x~ 0.4
\ ,.,,,,X, \ "'".i \
0.2
I
0.2
I
0.4
!
0.6
I
0.8
M
1.0
6-MP
FIG. 8. Responseisobologramfor combinationsof methotrexate(MTX)and 6-mercaptopurine (6-MP) acting on in vitro growth of LI210 cells in the presence of 2.5 × 10-SM thymidine (TdR). Solid line, experimental results; dashed line, kinetic prediction. (From Harrap and Jackson, 1975.)
Computer modelsof anticancerdrug interaction
261
tion of FUdR plus 6-MP; approximately additive effects were predicted and found with the cultured L1210 cell system. Figure 9 shows a simulation of a system subjected to inhibition by the combination of azaserine plus thioinosinate (the active form of 6-MP). The dotted line shows the effect of azaserine alone; an initial very pronounced inhibition is seen, then the degree of inhibition decreases as substrates accumulating behind the metabolic block tend to overcome the inhibition. In the third phase of the azaserine effect, the degree of blockade again increases as the slow irreversible phase of the azaserine inhibition of formylglycinamidine synthesis becomes significant. The effect of thioinosinate alone is seen in the solid line. The initial inhibitory effect is almost instantaneous, and it is only subsequently reversed to a slight extent, since the substrates of this reaction, NAD + and IMP both have alternative routes of utilization which result in the metabolism of a greater proportion of these substrates than does IMP dehydrogenase, the principal target enzyme of thioinosinate. Thus, IMP and NAD + do not accumulate, and the system soon reaches a new inhibited steady state. The dashed line shows the effects of azaserine plus thioinosinate in combination. The combined inhibitory effect is not additive, but it is appreciably greater than the inhibition caused by either agent alone. The properties of the model in this instance were very much what would have been predicted for two inhibitors acting sequentially at sites which straddle a divergent metabolic branch point (Harrap and Jackson, 1975). 10C %'..
8C - \ ": • .'" ........ .. ... .-...... ...-'" . . . .
~ ~:..
4C
Azaserine ( 5 . M )
..-'" ... '..
.-,""
ul
TIMP ~
4~)
q~O
~"
~
~-
I(~X)
(2~M)
Azaserine + T I M P
' 120
TIME (sec)
FIG. 9. Simulated effects of a binary combination of a z a s e r i n e and thioinosinic acid (TIMP) on the production of xanthosine m o n o p h o s p h a t e (XMP). (Jackson R. C. (1977) unpublished observations.)
To summarize, the conclusions which may be drawn from the modeling studies of Jackson and Harrap, and those of Werkheiser, Grindey and their colleagues are generally in agreement. Sometimes the interactions of two inhibitors appear to be predictable from the concurrent, sequential, or complementary nature of their sites of action, as was the case for the azaserine plus thioinosinate combination discussed above. In general, however, it cannot be assumed that this will necessarily be the case. Where the metabolic pathways involved contain multiple cycles or feedback inhibitions (and all the pathways where antitumor drugs act are characterized by such features), then the nature of drug interactions cannot be intuitively predicted. Thus, to reason that two or more drugs will be synergistic, additive, or antagonistic on the basis of the relative locations of their sites of action (even assuming single sites of action, which seldom occurs) is futile. Whether a mathematical model will be useful in this situation depends directly upon the extent to which the model reproduces the essential control features of the system under study. There can be no arbitrary rule about how complex a model must be before it becomes a realistic one, but considering the immense complexity of the regulation of the biochemical pathways upon which anticancer drugs act, it seems likely that a model which can accurately predict drug-drug interactions in most cases is likely to be considerably more complex than those in existence at present. What the existing models have done is to embody some of the salient regulatory features of nucleotide metabolism, show that the interactions JPT Vol. 4, No. 2--B
262
R.C. JACKSONand K. R. HARRAP
of drugs which are used to inhibit growth in the relatively simple L1210 cell culture system may be predicted in a reasonable proportion of combinations, and generally establish the feasibility of the approach. If the key regulatory features of a system are not written into a model, the behavior of the model will only fortuitously resemble that of the real system; thus, everything depends upon detailed and accurate quantitative biochemical studies of the target cells, and of the interactions of the drugs and their metabolites with the pathways concerned. The role of the mathematical model is to integrate the results of the analytical biochemical work and show how the pieces may fit together. The final question posed, and still not satisfactorily answered, is how quantitative changes in system parameters can produce qualitative changes in system behavior. The general conclusion drawn from the models discussed in this section was that usually this does not happen. Werkheiser et al. (1973) state: 'The 30 experiments yielded 270 interaction patterns. In many instances the intensity of the interaction was changed by parameter variation, but in only sixteen cases was there a change in the type of interaction'. These authors then quote the earlier conclusion of Garfinkel (1966): 'The behavior of the model is generally not very sensitive to changes in parameters; a drastic change in behavior usually requires the introduction of a new idea'. A similar conclusion was drawn from the detailed model of folate metabolism of Jackson and Harrap (1973) that usually a quantitative change in parameter values produced a quantitative change in the behavior of the model. However, from the regulatory point of view, it may be that the exceptions to this generalization convey more interesting information about the regulation of the system under study. Werkheiser et al. (1973) note: 'In four of these cases it was clear that the quantitative change in parameter value resulted in a qualitative change in the structure of the model'. An exact description of the kinetics which can show such properties is difficult. Clearly, in nature, quantitative changes in the concentrations of compounds can lead to qualitative differences in the systems affected. Examples are the effect of inorganic ion concentrations on the charged/discharged state of a nerve cell, the effect of the environmental lactose concentration on the/3-galactosidase content of E. coli (either uninduced or fully induced), and the fact that antimetabolites below a certain concentration may be tolerated by a tissue culture cell, but above that concentration the cell will die. Also relevant to the present discussion is the fact that malignancy or non-malignancy (clearly a qualitative property) is apparently determined by quantitative differences in the biochemical processes of the cell. These processes have been modeled, but always by introducing an arbitrary discontinuity into the model. For example, if the model is instructed that a cell will die if the concentration of tetrahydrofolate cofactors remains below a certain level for a given time, then the behavior of the model may reflect the properties of the cell, but it does not explain them. Far more interesting is the situation in which the discontinuity arises naturally from the rate equations. A few systems of this kind have been modeled; the study by Prigogine and Babloyantz (1972) on the control of E. coil /3-galactosidase is an example. All of the systems of this type studied so far have certain features in common. (i) There is an autocatalytic process, or a pair of cross-catalytic processes. In the c a s e of the /3-galactosidase system the entry of lactose into the cell is autocatalytic, since after sufficient lactose has entered the cell by passive diffusion to induce the lac operon, a galactoside permease is induced which greatly enhances the transport of lactose through the cell membrane. In the case of purine nucleotide biosynthesis there is a cross-catalytic process, in that GTP is a substrate for adenylosuccinate synthetase, and thus necessary for ATP biosynthesis, and ATP is a substrate for GMP synthetase, and therefore necessary for GTP formation. It is not known whether this cross-catalytic effect is of regulatory significance. (ii) Systems which give rise to discontinuities exhibit the phenomenon of multistability, i.e. they can exist in alternative stable steady states. Figure 10 illustrates this property in a simple model of glycolysis (Reich and Sel'kov, 1974). Figure 10a summarizes the highly simplified model used: v~ represents hexokinase, and product A is glucose
Computer models of anticancer drug interaction
263
6-phosphate, which inhibits its own production; ve represents phosphofructokinase (the intervening reactions between vl and v2 are not modeled, since they are considered to be of secondary importance from the regulatory point of view); v2 is activated by ADP, one of its products; v3 represents the succeeding reactions of glycolysis. Figure 10b plots v2 and v3 as functions of the concentration of B: v3 is considered to show a simple hyperbolic dependence upon B, and v, is not absolutely dependent upon B, but is allosterically activated by it; the solid lines represent v2 at increasing levels of A. Intersections of a solid line and the dotted line correspond to steady states for the subsystem vl, v2, i.e. points where B is produced and utilized at the same rate; note that at some values of A there is one steady-state point, for other values of A there are two, and at others, three. Where there are three steady-state points for a particular value of A, the middle one is always unstable for reasons explained by Reich and Sel'kov (1974). Figure 10c plots the value of the steady-state rate, v(--v2 = v3) as a function of A (solid line). Superimposed are plots of v~ as a function of A (dashed lines) for different amounts of enzyme 1 (hexokinase). At points at which a dashed line intersects the solid line, vt --- v2 - v3 and the system as a whole is in a steady state. Consider the system at point p~: the reaction is close to its maximum. If A is increased, or if E~ is increased, the system will move along the solid line toward Ps, and the steady-state rate, v, will show only a small increase. Now consider the system to be at point pro, when for some reason A is decreased, v2 and v3 can only reach a steady state by a marked reduction in rate, and the system drops to point p2; but now A is being produced (by v0 faster than it is used, so the system moves along the solid line to point p3, where once again vl-- v2 = v3. Now the system has reached a new steady state in which v is much less than before, and this large reduction in v was achieved by a very small reduction in A. Note, too, that the final value of A is greater than the initial value, which demonstrates the difficulties of drawing conclusions concerning changes in rates of processes from the size of precursor pools. Now imagine the system again at point pl, when the level of E~ is slightly decreased; the system in this case can only achieve a steady state by moving to p4, so again a small perturbation has resulted in a large change in v; if El now increases to its original value the system will not return to p~, but will move to p3. A
G6P
FDP ADP
E1
b
,. ,,
•
C
~|(A)
. . . . . .
::
\~
p \ ~.~.p_
~ \\ ~
_ ~ ' ~ 2 3
"~,, ¢\,',~\
v
"'?/' I
I
I
I
I
I
!
i
I
I
(A)
(e)
FIG. 10. Multistability in the glycolytic pathway. 10a. Schematic outline of the model of glycolysis used by Reich and Sei'kov (1974). v, is the rate of hexokinase (HK); v~ that of phosphofructokinase (PFK); v3, the succeeding reactions of glycolysis. Reactions between v, and v2 are not modeled. A represents glucose 6-phosphate (G6P) and B fructose diphosphate (FDP) and adenosine diphosphate (ADP). G6P inhibits its own formation, while v~ is activated by ADP. 10b. Plots of v2 (solid lines) and v3 (dotted line) as functions of B. 10c. Steady state values of v(-- v2 = v3) plotted as a function of A (solid line). Plots of vt as a function of A (dashed lines) are superimposed for different amounts of enzyme 1. (From Reich and Sel'kov, 1974.)
264
R. C. JACKSONand K. R. HARRAP
further increase in E~ will move the system along the solid line to p6, with a slight concomitant increase in v; increasing E~ still further will then cause a sudden large change as the system jumps to ps. Similarly, an increase in A when the system is at P6 results in a jump to P7, and then a decrease in A as the system moves to ps. It is clear that systems of this kind possess the appropriate properties to account for some biological switching processes. It remains to be determined whether this kind of kinetics could be involved in the triggering of cells from a static to a proliferating state, and if so whether antitumor drugs could then be employed to switch the system back again. In this connection, it is of interest that the adenine and guanine nucleotide branches of IMP utilization, as mentioned above, possess a cross-catalytic relationship, which could result in kinetics very similar to the autocatalytic situation just discussed in the glycolytic pathway. It is clearly important to conduct experimentally, detailed kinetic studies of enzymes which might show this kind of behavior, and then use detailed mathematical models to explore the optimal way to selectively switch off the biosynthetic processes of the malignant cells. 6. CYTOKINETIC MODELS The progression of cells through the cell cycle obviously depends upon intracellular biochemical events, which have become increasingly well understood in the last few years. However, it is convenient to separate consideration of models of the cell cycle from the models of metabolic pathways discussed in the last section, because the mathematical models of the cell cycle which have been produced primarily describe the movement of cells from one stage to another in numerical terms, rather than explain the cell cycle in biochemical terms. Cytokinetics, like pharmacodynamics, is a subject which is inevitably dependent upon a mathematical treatment, and whose complexity suggested at an early stage the development of mathematical models of some considerable sophistication. It is not the purpose here to review all of the uses of mathematical models of cell kinetics; instead, this discussion will concentrate on computer models which have been used to simulate drug effects. Some early models of cell growth were those of Barret (1966), Brockwell et al. (1970), Brockwell and Trucco (1970), Bronk et al. (1968), Bronk (1969), Franke (1970), Hahn, (1966, 1967), Hahn and Kallman (1967), Hirsch and Eagelberg (1965, 1966), Kember (1969), Takahashi (1968), Tibaux et al. (1969) and Trucco and Brockwell (1968). A review of some earlier applications of computer techniques in the study of the kinetics of tumor growth has been presented by Aroesty et al. (1973). An early model which made use of a large body of experimental data to study drug effects on the cell cycle was that of Shackney (1970), whose model was formulated to answer the following questions: (i) How can one predict whether an effective drug can be made a curative drug if dosage and schedule of administration are modified? (ii) Given that a schedule is potentially curative, how can one predict the optimal dosage and the optimal number of courses? (iii) Given that a particular drug may be curative in one experimental tumor system, how can one predict its behavior in another tumor system? (iv) Given that a drug is curative in animals when given according to particular schedules, how can one predict its behavior in human cancer, and how can one devise schedules that will be curative in man? Shackney attempted to approach questions of this kind by simulating the response of the LI210 mouse leukemia to AraC. His assumptions were such as to embrace as far as possible the known kinetic properties of the L1210 tumor. It was assumed that multiplying cells progressed in an orderly fashion through G~, S, G2 and M phases of the cycle and then re-entered G~. A reversible state of nonproliferation (Go) was also assumed. Growth was treated as asynchronous and it was assumed that a single viable cell was sufficient to produce a lethal tumor. AraC was treated as an S phase specific agent, and it was assumed that a given dose of drug would destroy the reproductive integrity of a constant fraction, not a constant number, of the susceptible population.
Computer models of anticancer drug interaction
265
The model made provision for nonproliferating pools of cells at each stage of the cell cycle, and for natural cell death. A distribution of cell-cycle times was also assumed. The constants for the model were then derived from certain standard treatment schedules. Experiments with the L1210 leukemia were simulfited at doses, schedules and implant sizes other than those used for standardization. The model's predictive capability was tested by simulating experiments without prior knowledge of results. Some results of such simulations are shown in Fig. 11. In both cases fair agreement was found between predicted and observed increases in life span. The author also discusses the limitations of his model. First, his treatment of slowly proliferating cells as nonproliferating, and the assumption that nonproliferating cells were completely drug resistant introduced errors, especially when large inocula, or noncurative drug doses, were simulated. The model made simplifying assumptions about the shape of the AraC dose-response curve, and did not deal with differences in drug-distribution characteristics at different anatomic sites in the host animal. Also, the model made no provision for host death from drug toxicity, and consequently at very high drug doses, predicted results better than those actually observed. Furthermore, the model did not provide for the emergence, during the course of therapy, of drug-resistant clones of tumor cells. Despite these limitations it was clear that a model of this kind was able to coordinate data concerning the kinetic properties of tumor cells, and allow that information to be used to optimize drug schedules. lO~
¢n .J ..i
+d I I
o
|
~//
:E
-I Z
lo3
1°:o
I
I t I !
t'l
/
I I
/t
/
/
I
Ii
I
/I
V
I/
I
II V
1 i I i
,
v i/I
v
I
I
I¢, I
7 9 11 1 15 17 19 21 23 25 27 TIME
0
I
4
I
I
8
I
I
12
/
I/
|
V
t i t/
'
I
16
~
/
I
20
~
J
24
I
I
30
I
I
32
(DAYS)
FIG. 11. Simulation of arabinosylcytosine (AraC) chemotherapy schedules. Treatment commenced at a tumor cell population of 107 cells. Solid lines represent the expected total number of cells, broken lines the expected number of viable cells. Schedule (a): 106 cells inoculated i.v. on day 0; treatment, 15 mg/kg of AraC every 3 hr x 8 on days 2, 6, 10, and 14; predicted host lifespan 26.6 days; predicted increase in lifespan 351 per cent; observed increase in lifespan: 240, 341,420 and 450 per cent in four experiments, respectively. Schedule (b): 104 cells inoculated i.v. on day 0; treatment,-15 mg/kg of AraC every 3 hr x 8 on days 5, 9, 13, and 17; predicted host lifespan 29.1 days; predicted increase in lifespan 234 per cent; observed increase in lifespan 281 per cent. (From Shackney, 1970.)
Jusko (1971, 1973) described models for the simulation of the effects of phasespecific and non-phase-specific drugs, which incorporated cytokinetic information. These models consider pharmacodynamic aspects of drug action as well, and will be discussed in the next section. Himmelstein and Bischoff (1973a,b) and Bischoff (1973) also described models for the study of phase-specific and non-phase-specific drugs which consider both pharmacodynamic and cytokinetic properties of the system (also discussed in Section 7). Roti Roti and Okada (1973) have described a mathematical model of the cell cycle of LS178Y mouse lymphoma cells. Their treatment makes the assumption that it is possible to divide the cell cycle into small time segments such that any given cell may
266
R. C. JACKSONand K. R. HARRAP
be uniquely classified into one of the states on the basis of cell age (i.e. time after division). Then, given the starting distribution of cells in the cell cycle and the growth conditions, the model will predict the distribution of cells at different times in the cell cycle. The authors of this model used it to simulate the effects of colcemide, thymidine, and colcemide plus thymidine. Donaghey and Drcwinko (1975) have described a flexible computer simulation program which they claim may be used for a variety of cytokinetic studies. The program was evaluated using data from cultured human lymphoma cells. In this model S phase was divided into three equal segments (early, mid, and late S); the length of stay for each cell in each segment of the cycle was independent, and normally distributed. G2 phase was subdivided into two equal segments. The length of stay of each cell in both segments was triangularly distributed and skewed to the right, because it was thought that a cell would need a minimum time to complete the biochemical events required for mitosis, but after that the exit time from G: would depend upon events which could be treated as random. The length of time in mitosis was similarly skewed. The proliferation rate was set at 1.6 (i.e. every 100 cells dividing were assumed to produce 160 daughter cells); this assumption was based upon the experimental observation that the probability of division was 0.82 for these cells. Ten per cent of the cells emerging from division were assumed to enter a Go compartment; since the model was only used to simulate short periods of time, the subsequent fate of these Go cells was not studied. The G1 phase was divided into an early period and a late period. The early GI period was skewed to the left, on the assumption that cells would proceed through a number of randomly distributed events until a state was reached in which cell division became likely, and that when this point was reached the cells proceeded fairly rapidly into the second stage of G1, where the cell was considered to be accumulating precursors necessary for it to enter the DNA synthetic phase. This model was used to predict the distribution of cells through the cell cycle following a 24-hr thymidine block. The model predicted that under these conditions 92 per cent of the cells would accumulate in S phase (the observed value was 85-90 per cent). These authors state their intention of using this model to study the effects of anticancer drugs. They conclude: 'We realize that the model has several deficiencies and that many assumptions are arbitrary. These shortcomings result from the paucity of knowledge about many fundamental biological principles governing the behavior of cell populations'. However, this CELLSIM program is sufficiently flexible that new information may be easily incorporated into models based upon it. Blumenson (1973, 1975) has produced a comprehensive model of the human granulopoietic system, and has applied this mathematical model to the study of cancer chemotherapeutic agents, in particular the relationship between drug effects on bonemarrow cells and the ensuing pattern of leukopenia in the peripheral blood. This author suggests that the origin, maturation, release from the bone-marrow, migration into the tissue spaces, and eventual death of granulocytes are now sufficiently well understood that detailed models are likely to give meaningful information about the relationships between changes in the various compartments of granulocytes and their precursors. These various compartments interact by a number of feedback controls. For example, during infection there is an increased demand by the tissue spaces for granulocytes from the peripheral blood, which in turn, affects proliferative controls in the marrow. The model is a deterministic, rather than a stochastic one; for example, the cell kinetics are described in sufficient detail to permit one to distinguish between the different phases of the cell cycle, and it is assumed that, once triggered into cycle, the cells pass through the different phases in an ordered and obligatory pattern. The cells in the marrow are classified into 130 one-hour compartments of increasing maturity. In the steady state, the size of each such compartment is equal to the steady-state hourly granulocyte turnover rate (GTR); the cells in the oldest (130th) marrow compartment are released into the blood each hour. Thereafter, all cells in the total blood granulocyte pool (TBGP) have an equal chance of leaving for the tissue spaces during a 1-hr period. This probability depends upon the needs of the tissue
'
Computer modelsof anticancerdrug interaction
267
spaces, and may vary from hour to hour. Not only do cells in the 130th marrow compartment always enter the peripheral blood after 1 hr in that compartment, but younger cells also may do so according to need, which is determined by the hourly blood granulocyte c o u n t . The loss of maturing cells from the marrow in turn stimulates stem cell activity. Stem cell block, and its consequences, may be simulated. Certain parts of the modeling of stem cell cycling are an extension of the earlier models of Lajtha et al. (1962) and King-Smith and Morley (1970). The faithfulness of Blumenson's model was tested by simulating the effects of 5-FU, considered to be administered as an intravenous drip for a period of 0.5-3 hr. It was assumed that a certain fraction of cells in S phase were killed during the treatment, this fraction depending upon dose and time of administration. Not only are stem cells killed by the drug, but also myeloblasts, promyelocytes and myelocytes in S phase. This leads to a decrease in total marrow cellularity which increases the probability of uncommitted stem cells differentiating in the granulocyte track. The predicted peripheral white blood cell counts following the simulated 5-FU treatment showed the kind of highly fluctuating pattern actually observed in treated patients. This model may thus be used to deduce the effects of 5-FU on the bone-marrow by projecting back from results found in the peripheral blood, and improve the accuracy of monitoring chemotherapy without the need for frequent marrow biopsies. Other suggested applications for this model are in the study of bacterial infection, to simulate bone-marrow transplants, and to study the breakdown of the feedback mechanisms for granulocyte production and release, and thus examine the implications of the various theories relating the loss of such control mechanisms with the etiology of chronic granulocytic leukemia. Although this model has not apparently been used to study drug interactions, the structure of the model would probably permit studies of this kind. Blumenson's model is of particular interest and importance as an example of a system which studies the effects of antineoplastic drugs, not on the tumor itself, but on a proliferative organ of the host. Woo et al. (1975) have developed a cell-cycle stage-specific multicompartrnental model to investigate the effects of drugs on tumor growth. The model can simultaneously consider the following types of effect: (a) non-cycle-specific cell kill and cycle stage-specific cell kill and (b) progression delay of cells through the cell-cycle phases, which may bring about an accumulation of cells in the various phases of the cycle. The model was an elaboration of an earlier one by Dombernowsky and Hartmann (1972). A schematic diagram of this model is shown in Fig. 12. This model was used to simulate the effects of BCNU, AraC and MTX on the cell-cycle and proliferation kinetics of the L1210 leukemia. BCNU kills cells in all stages of the cell cycle at equal rates and nonproliferating cells at slower rates, and in addition elongates cell-cycle time by increasing transit time in S phase. AraC causes some progression delay through all phases, especially in G1 phase but only kills cells in the S-phase. By comparison of predicted fractional survival curves, obtained using a range of parameters in the model, with experimental values, it was possible to obtain the best estimates for the rates of cell kill, rates of inhibition, and elongation of transit times for each cell compartment, at a number of dosage schedules for these two agents. In the case of MTX, the model was used to examine experimental data of Skipper et al. (1970) obtained at two different concentrations. Simulations in which a large fraction of cells were inhibited (i.e. transferred to a nonproliferative 'quiescent' compartment) and a certain dose-dependent fraction of cells in S phase were killed, gave close agreement with experimental results. These authors discuss how models of this kind may be used to optimize drug scheduling. For example, under conditions in which a drug causes cell kill in S phase, but also delays progression from Gl to S, the first dose of drug will kill the cells already in S phase. Considerably later, a wave of the delayed cells will then enter S, and if this time can be predicted accurately, a second administration of drug can be timed to kill this cohort .of delayed cells. It is suggested that this model can provide reliable predictions using far less detailed experimental information than that required by the model of Shackney.
R. C. JACKSONand K. R. HARRAP
268 -
-',
.
Sphase .
.
.
.
.
rLS--
I;
I L-. ii ~1
; T ;
.
.
.
I '
.
:
.
.
,,o
I :
L.._ ;..~_] 1 L_I 3 ~
.
'~
r
MITOSIS .
.
.
.
.
.
.
.
:
L
"
FIG. 12. Schematicdiagramof the relationshipbetweendrug-tumorinteractionsillustrated by the cell cycle stage-specific model. Each block enclosed by broken lines represents three possible modes of drug action within each cell cycle phase;/3* denotes the drug-inducedcell kill rate, 3' the fractionof cells accumulatingin the drug-inducedquiescentcompartment;transit time in each phase may also be elongated. (From Woo et al., 1975.) Chuang and Lloyd (1975) and Chuang (1975) have presented a mathematical analysis of the tumor model originally presented by Skipper (1969) in which the tumor cells are treated as existing in three compartments, proliferating, temporarily nonproliferating (resting) and permanently nonproliferating. The proliferating cells pass through the four cell-cycle phases, (i.e. G~, S, G2 and M). Chuang and Lloyd (1975) have developed a treatment of such a system in which the effects of antitumor agents are expressed in terms of loss functions for each proliferating (or potentially proliferative) compartment, these functions depending upon the drug concentration at the tumor site. The simulations based upon this approach were used to suggest a 'best' schedule, which did not, however, consider host toxicity. The model was tested by simulating the effects of AraC in a variety of dose schedules. An example of one such simulation is illustrated in Fig. 13. These simulations were made possible by parameter estimations based on data indicating the percentage of labeled mitoses, using a statistical procedure which calculated the theoretical fraction of labeled mitoses from the experimental values. It is shown that a model of this type can be used to distinguish between various hypotheses of drug action and inter-compartmental relationships of cells, by simulating the theoretical situations and determining which system best describes the observed behavior. Bender arid Dedrick (1975) have used mathematical modeling to investigate cytokinetic aspects of resistance to anticancer drugs. Their model contains the usual four proliferating phases of the cell cycle, a viable but nonproliferating Go phase, and a nonproliferating nonviable compartment. It is assumed that spontaneous cell loss may occur from all of these compartments. These authors used their model to calculate theoretical fractions of cell kill for an S phase specific agent for tumors with 16- and 24-hr cycle times, using repeat drug doses spaced at 8- and 16-hr intervals. The factors which may influence possible cell 'recruitment' into cycle are discussed. It is concluded that cytokinetic factors which influence responsiveness to drugs are tumor size and location, rate of cell proliferation, distribution of cell-cycle time and S phase transit time in the cell population, variation in the growth fraction, and the degree of synchrony of the tumor cell population. 7. M U L T I L E V E L MODELS Since the understanding of the interactions between a drug, a tumor, and a host animal requires study at several levels of biological organization, an obvious ex-
269
Computer models of anticancer drug interaction I06
z 0
'
'
'
'
'
'
'
'
'
'
÷ ' G I - 'C E L L I S x D o
iO 5
S-CELLS G2-CELL M-CELLS
S
I.< "J
:E bJ
,04 o
,. ¢,n
tu b-
10 2
e X td
I-
I0 I
U, I
i0 O
,
,
,
,
,
~
i
i
i
i
i
i
.
.
,
.
q
.
~
.
.
~
.
i
i
i
i
o) ,oO ~ E° -
o
,o-,~ io-2t
0
;
'~
3
;,
~
~ TIME,
r
~
~
,o
.
I~
13
,4
,5
DAYS
FIG. 13. C o n c e n t r a t i o n o f a z a b i n o s y l c y t o s i n e ( A r a C ) a n d e x p e c t e d cell p o p u l a t i o n s in p r o l i f e r a t i n g phases of the L1210 leukemia. Animals recieved an inoculum of l0 s cells and were
treated on a schedule of 15 mg/kg/dose every 3 hr x 8 on days 2, 6, l0 and 14. (From Chuang, 1975.)
tension of the types of mathematical model discussed in the preceding sections was a model containing features that embraced two or more such levels. An early example of this kind of system was the model of MTX action described by Werkheiser (1971). Fourteen parameters were used to describe the action of MTX on a tissue and its environment. They were: (a) rate of urinary excretion of drug; (b) tissue permeability; (c) generation time of cells and duration of S phase; (d) rate of dihydrofolate reductase (DHFR) synthesis; (e) rate of decay of free and MTX-bound DHFR; (f) Km and Vmx for D H F R and thymidylate synthetase; (g) total folate concentration; (h) normal D H F R concentration; and (i) critical FH4 concentration. This last parameter is based on the assumption that if the cellular tetrahydrofolate (FH4) concentration falls below a certain value then the cell will die. This assumption is not necessarily biochemically warranted, and is restrictive, in that cells can probably tolerate much lower than the normal availability of reduced folates if a supply of thymidine and purine is available; nevertheless, since the nature of antifolate-induced cell death (as distinct from mere growth inhibition) is not well understood, some kind of arbitrary condition of this kind is hard to avoid. It will be seen from the list of parameters given that some of them describe intracellular biochemical functions (i.e. enzyme kinetic properties, and rates of enzyme synthesis and turnover), some describe cytokinetic qualities of the target cell (i.e. generation time and duration of S phase) and some are concerned with the pharmacodynamics of the system. Thus, the behavior of MTX was concerned with three different levels of biological organization, though in all cases in a highly simplified way. The model as thus formulated resembled the properties of real systems in a number of respects. MTX toxicity, at sufficient doses, was cumulative; the total dose for lethality was 70-times the total when divided into four daily doses; and these properties were in agreement with experimental observations, as was the steepness of the MTX dose-response curve. This model predicted that tumors with low membrane permeability to MTX would be most sensitive to large single doses, while highly permeable tumors were most sensitive to small repeated doses; this prediction was also confirmed experimentally. An example of a pair of predictions is shown in Fig. 14. It was concluded that this model, simple though it was, was sufficiently realistic to have predicting power.
R. C. JACKSON and K. R. HARRAP
270
400.MOLAR ON DAY 0 10
~
700~MOLAR ON DAY 0
EE MTX (~M)
. FeEE F" 2
F~ COFACZORS,~-"-" (.M)
.k .... CELL COUNT
CELL COUNT
103 lO 1
DAYS
DAYS
FIG. 14. Time course of the effects of a single course of methotrexate (MTX) on a simulated tissue without a diffusion barrier. FH~, dihydrofolate; FH4, tetrahydrofolate. (From Werkheiser, 1971.)Reproduced by permission of the New York Academyof Sciences. A more recent model, along similar multilevel principles, is the very detailed description by Morrison et al. (1975) of the disposition of AraC and its metabolites. A pharmacokinetic compartmental model is used, based on that of Bischoff et al. (1971), with the augmentation by an additional compartment to represent the spleen, in the portal circulation, since this organ is believed to make an important contribution to AraC phosphorylation. The pharmacokinetic flow scheme used is illustrated in Fig. 15a. For each compartment shown mass-balance equations were drawn up which described the organ concentrations of AraC and its metabolites in terms of drug-flow in and out of each compartment. Where this model went beyond a purely pharmacokinetic scheme, however, was in the inclusion of a considerable number of intracellular biochemical processes concerned with the activation and inactivation of AraC (Fig. 15b). AraC is converted to its 5'-monophosphate by deoxycytidine kinase, which is inhibited by dCTP and AraCTP; deoxycytidine, as a competing substrate, also inhibits AraC phosphorylation. The enzymes deoxycytidylate kinase and nucleoside diphosphokinase further phosphorylate AraCMP to the triphosphate level, and these reactions are opposed by corresponding nucleotidases. The deamination of AraC and AraCMP to AraU and AraUMP by cytidine deaminase and deoxycytidylate deaminase, respectively, is also modeled, as is the interconversion of AraU and AraUMP. Apart from effects of dCTP and AraCTP on deoxycytidine kinase, no other feedback inhibitions are considered, although there is provision for competition for cytidine deaminase by the competing substrate, deoxycytidine. This model was then used to simulate events following the administration of single doses of AraC to mice. Of particular importance was the intracellular concentration of AraCTP, the active form of the drug. This was shown to persist at cytoxic concentrations long after most of the AraC was cleared from the blood, explaining the fact that a single dose of AraC will kill cells that were not in S phase (as measured by the incorporation of tritiated thymidine) during the period of measurable plasma AraC. The fact that AraCTP follows a very different time-course from AraC is obviously of great significance in understanding the action of AraC. Higher nucleotidase activities, of course, decrease the steady-state concentration of AraCTP and the length of time that it persists at cytotoxic levels. The livers of DBA and C3H strains of mice show a very large difference in deaminase activities, and the effects of this enzymatic divergence were simulated, resulting in major differences in the predicted pharmacokinetics of AraC.
271
Computer models of anticancer drug interaction
! I
~.E,S
r O.,NE ~
LEAN
.......i
I.
r
(o)
dTMP I dGMP J dAMP
membrane oxtmce~dsr JMJ'g,~oJhdw fluid fluid
+
d~/,
. . . . . . . . . . . . .
.® .e~ I
"~,,11~
~
•
- dCTP
® : ~
I
a r a - C , i b - l a r l i - C 4*- I l r l I - C M P .4. a r a - C D P .Q- I m - . C T P . . . . .
®
(i)
•®
~
.NA
® m-u:-;m-u~- m-UMP :- ,,.-UOP = a,.-UrP I
,
0
(b) FIG. 15. Pharmacokinetic simulation of the disposition of arabinosylcytosine (AraC) and its metabolites. (a) Flow scheme of the distribution of AraC to various tissue compartments. (b) Metabolism of AraC. AraU, arabinosyluracil; araCMP, arabinosylcytosine monophosphate; araCDP, arabinosyi cytosine diphosphate; araCTP, arabinosylcytosine tripbosphate; araUMP, arabinosyl uracil monophosphate; araUDP, arabinosyl uracil diphosphate; araUTP, arabinosyl uracil triphosphate; dTMP, thymidine monophosphate; dGMP, deoxyguanosine monophosphate; dAMP, deoxyadenosine monophosphate; dCTP, deoxycytosine triphosphate; dCR, deoxycytidine. (From Morrison et aL, 1975.)
The authors state that future developments of this very interesting model will include more biochemical details, such as further phosphorylation of AraUMP; more body compartments, such as the peritoneal cavity; scale-up for human application; and inclusion of cytokinetic information. It is to be hoped that this approach may be extended to the study of other antineoplastic drugs and to drug-drug interactions. Werkheiser's model (1971) was an attempt to combine simple descriptions of pharmacodynamic, cytokinetic and biochemical properties of a tumor-host system in a single model. The simulations of Morrison et al. ~(1975) included more biochemical information and a great deal more pharmacodynamic information, but did not, at least in the version described, include cytokinetic detail. The remaining work to be
272
R.C. JACKSONand K. R. HARRAP
reviewed in this section is concerned with mathematical models combining cytokinetic and pharmacodynamic aspects. Jusko (1973) described a pharmacodynamic model for cell-cycle-specific chemotherapeutic agents which was a development of his earlier work on the pharmacokinetics of non-phase-specific agents (Jusko, 1971). In the cell-cycle-specific model the cells are considered to exist in either of two interconvertible compartments Cs, the drug-sensitive cells, and Cr, the drug-resistant cells; interconversion of cells from the two compartments was described by first-order rate constants. The nonproliferating cell compartment was susceptible to cell loss or differentiation, and the proliferating compartment could be depleted by drug-induced cell death. Drug-receptor interaction was considered to be bimolecular and irreversible, and assumed to bind a negligible proportion of the total drug dose. The sensitive (Cs) cell compartment was capable of cell proliferation, represented by a single parameter; the model could represent steady state, exponential, or Gompertzian behavior in the absence of drug. The pharmacokinetic aspects of the system were represented either by a one-compartment open model or a two-compartment open model. The model was tested by using it to fit data obtained by Bruce and co-workers with the effects of vinblastine on lymphoma and hematopoietic cells. The effects of repeated doses of AraC were also simulated. As well as predicting cell kill, the model could be used to estimate the survival time of drug-treated animals. Himmelstein and Bischoff (1973a, b) have developed a pharmacokinetic/cytokinetic model which simulates cell-cycle events in considerable detail; the cytokinetic treatment was based upon a model by Rubinow (1968). An important difference from the model of Jusko is that a continuous distribution of cellular maturity is considered, rather than two discrete fractions. The authors consider the possibilities that drugs kill cells or render them incapable of cell division, and do not describe the situation of drug-related growth inhibition by prolongation of generation times. The effects of direct action drugs (i.e. agents which kill cells regardless of their stage in the cell cycle) are first treated on the assumption that cell growth in the absence of drugs would be exponential. Drug-induced cell death is computed from the loss function, which relates cell loss to drug concentration. Various types of loss function are described; for example, first-order, or two-parameter saturable, or a threshold type relationship, where drug is only effective above a certain concentration, or combinations of these three possibilities. On the basis of these simplifying assumptions, an explicit relationship was derived between the cell growth parameters, and the loss function for a non-phase-specific (direct action) drug, and the cell number at a given time. Then the effect of cycle-specific drugs is described, now with the assumption that the loss function varies at different phases of the cell cycle. The drug concentration at a given time is treated by conventional pharmacokinetic expressions. The model of Himmelstein and Bischoff was then tested by simulating the interaction of AraC with mouse leukemia L1210 cells, although curiously, making the assumption that the drug was not phase-specific; this treatment was justified on the grounds that the L1210 cell spends 90 per cent of its time in the S phase, at least when the total cell number is small. Data of Skipper et al. (1964) were considered. The pharmacokinetic parameters were based on measurements of AraC clearance made by Skipper et al. (1967), although as shown by Chou et al. (1975), the kinetics of the active drug form, AraCTP, are quite different from those of AraC. The model of Himmelstein and Bischoff was also modified to describe the results of 12hr continuous infusion of AraC. Despite some of the rather drastic simplifications made during the course of construction of this model, it did appear capable of giving fairly accurate simulations of the results of various AraC treatment schedules with the L1210 system. These authors state 'a comprehensive model for a given system is composed of suitable growth kinetic, drug-cell interaction, and pharmacokinetic models'. Their own studies represent a preliminary step toward this goal of a complete, multilevel description of drug-host-tumor interactions. Lincoln et al. (1974) describe a plan for a wide-ranging 'leukemia simulator', which would not only consider drug clearance and distribution, and the cytokinetics of
Computer models of anticancer drug interaction
273
normal and transformed cells in the marrow, but also gastrointestinal susceptibility, kidney function, liver function, and evaluate risks from present or probable infections, risks of hemorrhage, etc. This grand synthesis represents, at present, a goal rather than an accomplised fact, and the authors describe some of the work in progress towards this objective, illustrating their model's capabilities by simulating the effects of AraC on granulocyte and platelet counts. Finally, one model has been presented which describes the effects of combined chemotherapy and radiotherapy (Andrews, 1969). This preliminary treatment of combined modality therapy will probably form a basis for more detailed theoretical explanations of how to optimize cancer treatment using radiation in combination with drugs. 8. T E C H N I C A L ASPECTS OF COMPUTER MODELING Since this area of research represents the common ground between biology and mathematics, it tends to be pursued by scientists whose primary background lies in one or another of these disciplines, but who move with less confidence in the other. Thus, there are, in the literature, a number of mathematical models of considerable technical sophistication, but which have never been applied to any real biological problem; conversely, there are attempts by biologists to describe pharmacological systems in a quantitative way which are hampered by use of inadequate mathematical techniques. This state of affairs is beginning to change as more pharmacologists become proficient in computer programing and numerical analysis, and as investigators from the traditionally strongly mathematical field of physiological control theory have come to take an interest in quantitative pharmacology. Existing mathematical models of pharmacological systems fall into a number of groups which reflect the different underlying philosophical approaches. 8.1. ABSTRACT AND REALISTIC MODELS
Weinberg and Berkus (1971) review a number of simulations written for abstract systems which resembled living cells in their general organization. These simulations were not models of any particular system, and were not designed so that the results of the simulation could be checked against experimental data. These models may be used, for example, to consider the various possible modes of organization which could be used by living matter to solve particular problems of control. In contrast, are realistic models which attempt to describe a particular biological system in which the aim is to describe the experimentally observed behavior of the system as closely as possible. The simulation of E. coil by Weinberg and Berkus (1971) falls into the latter category, as do most of the simulations discussed in the present review. 8.2. STOCHASTIC AND DETERMINISTIC MODELS
Chemical kinetics, the science from which enzyme kinetics and, more remotely, cytokinetics and pharmacokinetics are derived, is a deterministic science; given a set of rate constants and a set of reactant concentrations, the rate equation may be used to compute the reaction rate. Because the behavior of bulk matter is considered, the probabilistic aspects of the motions of individual molecules have no influence. The work of recent years has shown that processes obeying the laws of chemical kinetics, and the law of diffusion, are capable of giving rise to discontinuities in space and time; thus, in principle, given sufficient knowledge of enzyme kinetics (including the kinetics of multienzyme aggregates and of membrane-bound enzymes), of transport processes, and of the dynamics of circulation, excretion, etc., it should be possible to construct a strictly deterministic model of a pharmacological system which will provide a description to any required degree of detail. In practice, for a model of a very complex system, it is impossible, on a computer of present-day size, to directly simulate every process which might have some trivial influence on the system. This
274
R. C, JACKSONand K. R. HARRAP
problem has been encountered in the field of meteorology, in the construction of weather-predicting models of the atmosphere, where it is known as the 'butterfly-wing effect'; this name reflects the fact that, in a. system of sufficiently low stability, the flutter of a butterfly's wing could set in motion (in theory, at least) a chain of events which could culminate in a hurricane a thousand miles away. The meteorologist cannot simulate every butterfly, and thus must accept the fact that complex systems are subject to small perturbations which must be treated as random. Obviously, exactly analogous considerations apply to pharmacology, and the decision as to whether a particular process is to be treated as random or deterministic depends upon how central the process under consideration is to the questions being asked by the model. For example, if it is observed experimentally that a certain dose of drug, administered for a certain time, kills 40 per cent of tumor cells, then, if a single cell is considered, it might be reasonable to treat the survival or otherwise of that cell as a stochastic process, with a probability of death of 0.4. However, if it was known that the cell population was mixed, and that 60 per cent of the cells had plasma membranes impermeable to the drug, it might be more reasonable to decide the fate of a single cell on a deterministic basis, according to the membrane permeability of the drug. Thus, it is often convenient to incorporate both deterministic and random events into a model, the latter being used to simulate either processes which represent small perturbations rather peripheral to the real purpose of the model, or sometimes processes which cannot be simulated deterministically because of lack of information about what controls them. Technical aspects of the simulation of stochastic processes by MonteCarlo methods are discussed by Sheppard (1969). 8.3. THE SCOPE OF A SIMULATION
Simulations of biological systems vary from broad, shallow descriptions intended to describe the essential features of a wide-ranging process, such as the effects of drugs on progression through the cell cycle, to descriptions in depth (e.g. at the molecular level) of small subsystems (e.g. a dozen enzymes with their substrates and cofactors). A simulation may, of course, attempt to be both wide-ranging and deep; eventually, either programing time, expense, or computer memory availability will become a limiting factor. The decision of what ground should be covered by a simulation, and in how much detail, must be based upon experience of the real system under study, with continual regard to the questions which the model is being constructed to answer. Too much unnecessary detail makes a model expensive and hard to manage, and too little detail makes an unrealistic model, possibly so unrealistic as to have no predictive value. Garfinkel (1969), considering this question, suggests that a useful compromise may be to simulate the subsystem under study in as much detail as possible, with the remainder of the system (the 'environment') described in less detail. 8.4. LINEAR AND NONLINEAR MODELS From the very beginning of the study of enzyme kinetics it was clear that linear (first-order) kinetics were only observed in rather special circumstances, such as when a substrate concentration is very low in comparison to the enzyme Km value, and this science has therefore stressed nonlinear models from the outset. In some qualitative types of biochemical models, 'linearized' transformations of approximate rate equations have been used in a way such that the model retains the qualitative properties of the prototype system (e.g. Reich and Serkov, 1974), but in general, models of multienzyme systems inevitably imply nonlinear equations. The situation has been rather different in pharmacokinetics, where many processes of drug uptake, distribution, and elimination can be described with reasonable quantitative accuracy by first-order kinetics. Linear models have also been used in cytokinetics. Nevertheless, even in pharmacokinetics, it is clear that a detailed description of many processes entails a nonlinear model; for example, the active transport processes responsible for cellular uptake of numerous antitumor drugs may often be described by a Michaelis-
Computer models of anticancer drug interaction
275
Menten equation. An exception to the general rule that biochemical systems cannot be studied rigorously using linear models is provided by tracer experiments. These may be considered as first-order perturbations, and thus permit the study of complex nonlinear systems as a set of simpler, linear subsystems (Berman, (1969). The general rule that simulations in biochemical pharmacology require nonlinear models obviously does not mean that particular components of these models may not show first-order kinetics; this is frequently the case, for example, with the non-carriermediated membrane transport of drugs, or non-enzymic reactions, such as the spontaneous breakdown of 5,10-methylene tetrahydrofolate to tetrahydrofolate plus formaldehyde. Thus, models of pharmacological systems are likely to consist of sets of nonlinear differential equations. In some simple cases, the differential equations may be amenable to transformation into algebraic expressions, but this will not generally be the case, and for the most part the system will have no explicit solution. The usual approach, therefore, is to solve the differential equations numerically, and for this purpose second- or fourth-order Runge-Kutta or predictor-corrector methods seem to have been most used. The only serious problem encountered is that of 'stiffness'. This property, so called because it was encountered first in the study of vibrations of stiff springs, results in a system of equations in which one component has a time constant which may differ by orders of magnitude from those of the rest of the system. In this situation, although the equations can be solved by the usual methods, doing so entails extreme inefficiency. The problem was first reported in biochemical systems by Garfinkel and Hess (1964); these authors were able to alleviate the problem to some degree by altering the molar enzyme concentration (upward) and the turnover number (downward) so that the simulated enzyme activity did not change. For general methods of solution of sets of stiff differential equations, the reader is referred to articles by Butcher (1964), Cooper (1969) and Gear (1971). 8.5. TEMPORAL HIERARCHIES IN BIOLOGICAL MODELS
Any biological model of broad scope may have to simulate processes which take place over very different time scales. Typically, enzyme-substrate interactions take place in fractions of a second, or milliseconds; intermediary metabolite equilibrations may take seconds or minutes; biosynthesis and degradation of enzymes in mammalian cells often involve halfqives of several hours; and the cell division cycle takes between 9 hr and many days. Clearly, the simultaneous simulation of all of these processes, as would be required for a complete understanding of drug-cell interactions, is likely to pose great technical difficulties. This problem has been the subject of a theoretical discussion by Park (1974) who suggests methods for handling problems of this kind. The present authors are not aware of any applications of this approach to a realistic biological model, but it seems to be potentially of great value. The problem of a wide range of time scales was encountered during the construction of the model of folate metabolism in the L1210 mouse leukemia cell (Jackson and Harrap, 1973). Substrate and cofactor interconversions required a few seconds or minutes, MTX transport across the cell membrane typically required from 15 min to several hours before a steady state was reached, and enzyme turnover times ranged between 18 and 90 hr. The approach to the problem was as follows: (i) enzymesubstrate and enzyme-inhibitor complexes were assumed to be in or near the Briggs-Haldane steady state, so that steady-state rate equations could be used, with the exception of the complex between dihydrofolate reductase and MTX, which dissociates with a t~12of about 1 hr, and (ii) the slow processes (i.e. cell division, drug clearance, drug uptake through the cell membrane, and synthesis and breakdown of enzyme molecules) were simulated directly. The rapid processes of substrate and cofactor interconversion were assumed to be close to the steady state throughout. Thus, the cellular inhibitor concentration and enzyme levels at a given time would be calculated by direct simulation, and the substrate and cofactor levels corresponding to
276
R. C. JACKSON and K. R. HAP.RAP
the appropriate enzyme activities would be obtained by interpolation in an array of stored values which had previously been shown by direct simulation to correspond to those enzyme activities. This part of the process was iterative, since, of course, actual enzyme activities would depend upon levels of substrates, feedback inhibitors, etc. The chief disadvantages of this 'two-level' method are, first, the large amount of computer time needed to assemble the table values of steady-state substrate and cofactor concentrations (some of which may never be used) and, second, the large amount of memory required to store these tables of steady-state values when the second stage (slow reactions) of the simulation is running. 8.6. DIGITAL AND ANALOG COMPUTERS Once a model has been constructed, that is, the problem has been stated completely in mathematical terms, the actual solution of the resulting equations is purely a matter of technical convenience. Now that many biologists are conversant with the programing of digital computers, it is likely that most workers who wish to use mathematical modeling as an occasional problem-solving tool will use digital computers. Analog computers have enthusiastic proponents, chiefly on the grounds of their relative simplicity of operation and low initial cost. However, it is probably fair to say that proficiency in analog computing requires more specialized expertise than participation in a large, time-shared digital computer system, which, from the users' point of view, can be treated as a problem-solving black box. Hybrid computers are gaining ground slowly but may be much more widely used in the future; their application to biological problems was discussed by Girling (1969). 8.7. SIMULATION LANGUAGES The digital computer programer has a choice between writing his program entirely from the beginning in a general-purpose language, usually FORTRAN,BASIC,or ALGOL,or of using one of the special-purpose simulation languages. The choice probably depends upon how much the completed simulation is likely to be used. A program completely written especially for a particular function is likely to be more efficient in terms of computer time and memory, but will no doubt require a greater investment of programing time, than one which uses a~ special simulation compiler. The latter generally require a lot of core memory, but they seem to be very suitable for this kind of work; for example, Werkheiser et al. (1973) used MIMIC. 9. CONCLUSIONS: FUTURE USES OF COMPUTER MODELS OF DRUG INTERACTION It is clear that although mathematical modeling has become an established technique in pharmacology, there are still very few instances of the method being applied to the study of drug-drug interactions. Many of the models described in this article have been concerned primarily with reaching an understanding in depth of the action of a single agent on a particular system, or at most with the effects of a small number of agents used singly. The few models described which have studied drug-drug interactions, dealing mainly with interactions at the level of cellular metabolism, however, show the potential of simulation, both as an investigative method in furthering the understanding of drug-host, drug-tumor and drug--drug interactions; as a predictive tool, in suggesting probable positive and negative interactions of various agents; and as a means of optimizing chemotherapeutic regimens, to tailor them to the biochemical characteristics of particular individual tumors. Several of the models discussed here, although not yet used to study drug-drug interactions in the systems they simulate, seem to have the potential to be used for that purpose, and in a few cases their authors have stated their intentions of pursuing such studies. Since the main purpose of these models is to quantitate the properties of systems too complex to be understood intuitively, it is reasonable to believe that the greatest
Computer models of anticancer drug interaction
277
contribution to understanding the complexities of cancer chemotherapy, in particular combination chemotherapy, will be made by those models which are sufficiently broad and deep to represent accurately all the essential features of the interactions between a tumor, a host, and several drugs. This will require multilevel models; that is, models containing biochemical, pharmacodynamic and cytokinetic features. How would such models be most effectively used? Consider that a new antitumor drug has been introduced, and its efficacy as a single agent has been demonstrated. When the time comes to consider possible ways of combining the new drug with other agents, a good mathematical model could make a valuable contribution. The number of possible combinations of, say, three or four drugs, at a number of different dose ranges, is very large, and to test all possible combinations experimentally would be prohibitively time-consuming and expensive. To simulate a large number of such combinations with a computer model would be relatively rapid and inexpensive. Those combinations predicted to be most promising could then be tested experimentally. Thus, the computer model is not a substitute for experiment, but a guide which should help to increase the productivity and cost-effectiveness of limited experimental resources. A second example of how computer models may help to improve chemotherapeutic effectiveness in ways which might otherwise be impossible, is in the optimization of chemotherapy tailored to the biochemical properties of particular tumors. If certain biochemical facts (i.e. enzyme activities, nucleotide pool sizes, etc.) are established for a biopsy specimen, how may every item of available information be most effectively exploited? By incorporating the available information into a model, various drug regimens could be tested against a simulated system which represents the properties of the tumor as accurately as possible, and the 'most probably effective' regimen selected for actual use. Thus, it is felt that mathematical simulation is likely to become an established biological and medical technique, not as a substitute for human intelligence, experimental skill, or clinical judgement, but as a tool which has its place, along with all of the others, in increasing the range and effectiveness of scientific medicine in the increasingly complex field of antitumor chemotherapy. AcknowledgementsmTbe authors thank the following publishers for permission to reproduce various materials cited in the text: Academic Press; American Pharmaceutical Society; Elsevier; National Cancer Institute; New York Academy of Sciences; and Pergamon Press. This work was supported by grants from USPHS (CA-18129), the Medical Research Council and the Cancer Research Campaign.
REFERENCES ALLEN, L. M. and CREAVEN, P. J. (1975) Human pharmacokinetic model for isophosphamide. Cancer Chemother. Rep. 59: 877-882. ANDREWS,J. R. (1969) Combined cancer radiotherapy and chemotherapy: the relevance of cell population kinetics and pharmacodynamics. Cancer Chemother. Rep. 53: 313-324. AROEffrY,J., LINCOLN,T., SnAPmo, N. and Boco^, G. (1973) Tumor growth and chemotherapy: mathematical methods, computer simulations, and experimental foundations. Math. Biosci. 17: 243-300. BAtmE'rr, J. C. (1966) A mathematical model of the mitotic cycle and its application to the interpretation of percentage labeled mitoses data. J. natn. Cancer Inst. 37: 443-450. BELLMAN, R., JACQUEZ,J. A. and KALABA,R. (1960) Some mathematical aspects of chemotherapy: I. One-organ models. Bull. Math. Biophys. 22: 181-198. BENDEa, R. A. and DEDmCK, R. L. (1975) Cytokinetic aspects of clinical drug resistance. Cancer Chemother. Rep. 59: 80.~809. BERMAN,M. (1969) Kinetic modeling in physiology. FEBS Lett. 2: $56--$57. BISCHOFF, K. B. (1973) Pharmacokinetics and cancer chemotherapy. J. Pharmacokinet. Biopharmac. 1: 465--48O. BISCnOFF, K. B. (1975) Some fundamental considerations of the applications of pharmacoidnetics to cancer chemotherapy. Cancer Chemother. Rep. 59: 777-793. BlSCHOFF, K. B., DEDRICK, R. L. and ZAHARKO, D. S. (1970) Preliminary model for methotrexate pharmacoidnetics. J. pharmac. Sci. 59: 149--154. BlSCHOFF, K. B., DEDmCK, R. L., ZAHARKO,D. S. and LONGSTRETH,J. A. (1971) Methotrexate pharmacokinetics. 3". pharmac. Sci. 60: 1128-1130. BLUMENSON, L. E. (1973) A comprehensive modeling procedure for the human granulopoietic system: over-all view and summary of data. Blood 42: 303-313. JFT VoL
4, No. 2.--C
278
R . C . JACKSON and K. R. HARRAP
BLUMENSON, L. E. (1975) A comprehensive modeling procedure for the human granulopoietic system: detailed description and application to cancer chemotherapy. Math. Biosci. 26: 217-239. BROCKWELL, P. J. and TRucco, E, (1970) On the decomposition by generations of the PLM-function. J. theoret. Biol. 26: 149-179. BROCKWELL, P. J., MACLAREN, M. D. and TRUCCO, E. (1970) Monte Carlo simulation of PLM-curves and collection functions. Bull. Math. Biophys. 32: 429-443. BRONK, B. V. (1969) On radioactive labeling of proliferating cells: the graph of labeled mitosis. J. theoret. Biol. 22: 468-492. BRONK, B. V., DIENES, G. J. and PASKIN, A. (1968) The stochastic theory of cell proliferation. Biophys. J. 8: 1353-1398. BUELL, J., JELUFFE, R., KALAaA, R. and SRIDHAR, R. (1970) Modern control theory and optimal drug regimens. II. Combination therapy. Math. Biosci. 6: 67-74. BUTCHER, J. C. (1964) Implicit Runge-Kutta processes. Math. Comp. lg: 50-64. CHANCE, B., GAR~NKEL, D., HIC~XNS,J. and HESS, B. (1960) Metabolic control mechanisms. V. A solution for the equations representing interaction between glycolysis and respiration in ascites turnout cells. J. biol. Chem. 235: 2426-2439. CHOU, T. -C., HUTCmSON, D. J., SCHMID,F. A. and PmLIPS, F. S. (1975) Metabolism and selective effects of 1-/3-D-arabinofuranosylcytosine in LI210 and host tissues in vivo. Cancer Res. 35: 225-236. CHUANG, S. -N. (1975) Mathematical models for cancer chemotherapy; pharmacokinetic and cell kinetic considerations. Cancer Chemother. Rep. $9: 827-842. CHUANt, S. -N. and LLOYD, H. H. (1974) Analysis and identification of stochastic compartment models in pharmacokinetics: implication for cancer chemotherapy. Math. Biosci. 22: 57-74. CHUA~t~, S. -N. and LLOYD, H. H. (1975) Mathematical analysis of cancer chemotherapy. Bull. math. Biol. 37: 147-160. COHEN, J. L., JAO, J. Y. and JUSKO, W. J. (1971) Pharmacokinetics of cyclophosphamide in man. Br. J. Pharmac. 43: 677-680. COOPER, (3. J. (1969) The numerical solution of stiff differential equations. F E B S Lett. 2: $22-$29. DANENBERG, P. V., LANGENBACH,R. J. and HEIDELBERGER, C. (1974) Stuctures of reversible and irreversible complexes of thymidylate synthetase and fluorinated pyrimidine nucleotides. Biochemistry 13: 926-933. DEDRICK, R. L , BISCHOFF,K. B. and ZAHAR~O, D. S. (1970) Interspecies correlation of plasma concentration history of methotrexate. Cancer Chemother. Rep. 54: 95-101. DEDI~CK, R. L., FORRESTER, D. D. and HO, D. H. W. (1972) In vitro-in vivo correlation of drug metabolism-deamination of l-~8-D-arabinofuranosylcytosine. Biochem. Pharmac. 21: 1-16. DEDRICK, R. L., FORRESTER, D. D., CANNON, J. N., EL DAREER, S. M. and MELLETT, L. B. (1973a). Pharmacokinetics of l-/3-D-arabinofuranosylcytosine (Ara-C) deamination in several species. Biochem. Pharmac. 22: 2405-2417. DEDRICK, R. L., ZAHARKO,D. S. and LUTZ, R. J. (1973b) Transport and binding of methotrexate in vivo. J. pharmac. Sci. 62: 882-890. DEDRICK, R. L., ZAnARKO, D. S., BENDER, R. A., BLEYER, W. A. and LUTZ, R. J. (1975) Pharmacokinetic considerations on resistance to anticancer drugs. Cancer Chemother. Rep. 59: 795-804. DOMBERNOWSKY,P. and HARTMANN,N. R. (1972) Analysis of variations in the cell population kinetics with tumor age in the LI210 ascites tumor. Cancer Res. 32: 2452-2458. DOSAtHEY, C. E. and DREWINKO,B. (1975) A comphter simulation program for the study of cellular growth kinetics and its application to the analysis of human lymphoma cells in vitro. Comp. Biomed. Res. 8: 118-128.
FRANKE, E. K. (1970) A mathematical model of synchronized periodic growth of cell populations. J. theoret. Biol. 26: 373-382. GARFINKEL, D. (1966) A simulation study of the metabolism and compartmentation in brain of glutamate, aspartate, the Krebs Cycle and related metabolites. J. biol. Chem. 241: 3918-3929. GAR~NKEL, D. (1969) Construction of biochemical computer models. F E B S Lett. 2: $9-S13. GAR~NKEL, D. and HESS, B. (1964) Metabolic control mechanisms. VII. A detailed computer model of the glycolytic pathway in ascites cells. J. biol. Chem. 239: 971-983. GARFINKEL, D., RUTLEDGE,J. D. and HIOGINS,J. J. (1961) Simulation and analysis of biochemical systems. I. Representation of chemical kinetics. Communs Ass. comput. Mach. 4: 559-562. GEAR, C. W. (1971) DIFSUB for solution of ordinary differential equations [D2]. Communs Ass. comput. Mach. 14: 185-190. GIRLIr4C, B. (1969) An introduction to hybrid computers and their application to optimization problems. F E B S Lett. 2: $58-$62. GR1NDEY, G. B., WER~:HEISER, W. C., MORAN, R. G. and NICHOL, C. A. (1970) Anomalous inhibition of growth by combinations of inhibitors of DNA biosynthesis. Fed. Proc. 29: 609. GmSDEY, G. B., MOP.AS, R. G. and W~.RKHEISF-~,W. C. (1975) Approaches to the rational combination of antimetabolites for cancer chemotherapy. In Drug Design, V01. 5, pp. 169-249, Academic Press, New York. GROSER, G. and CLARK, R. (1971) BIOMOD: An interactive graphics system for analysis through simulation. Proc. 1971 I E E E Conference on Decision and Control, Miami Beach, Florida, pp. 197-201. HAns, G. M. (1966) State vector description of the proliferation of mammalian cells in tissue culture. I. Exponential growth. Biophys. J. 6: 275-290. HAns, G. M. (1967) Cellular kinetics, cell cycles and cell killing. Biophysik. 4: 1-14. HAns, G. M. and K ~ N , R. F. (1967) State vector description of the proliferation of mammalian cells in tissue culture. II. Effects of single and multiple doses of ionizing radiations. Radiat. Res. 30: 702-713. HARRAP, K. R. and JACKSON,R. C. (1975) Enzyme kinetics and combination chemotherapy: An appraisal of current concepts. Adv. Enz. Regul. 13: 77-96. HARRIS, P. A. and GRoss, J. F. (1975) Preliminary pharmacokinetic model for adriamycin. Cancer Chemother. Rep. 59: 819-825.
Computer models of anticancer drug interaction
279
HIMMELSTEIN, K. J. and BISCHOFF, K. B. (1973a) Mathematical representations of cancer chemotherapy effects. J. Pharmacokinet. Biopharmac. 1: 51-68. HIMI~LSTEIN, K. J. and BISCHOFF,K. B. (1973b) Models of AraC chemotherapy of LI210 leukaemia in mice. J. Pharmacokinet. Biopharmac. 1: 69--81. HmscH, H. R. and ENOS,BERG, J. (1965) An analog of ceil-culture growth. J. theoret. Biol. 9: 303-312. HIRSCH, H. R. and ENOELBERG, J. (1966) Decay of cell synchronization: solutions of the cell-gTowth equation. Bull. Math. Biophys. 28: 391--409. HOWELL, J. R. (1975) Mathematical formulation for nonuniform multiple dosing. J. pharmac, Sci. 64: 46 A, A.~6. JACKSON, R. C. and H A I ~ , K. R. (1973) Studies with a mathematical model of folate metabolism. Archs Biochem. Biophys. 158: 827-841. JACKSON, R. C. and NIETHAMMER,D. (1977) Acquired methotrexate resistance in lymphoblasts resulting from altered kinetic properties of dihydrofolate reductase. Eur. J. Cancer 13: 567-575. JACKSON, R. C. and WEBER, G. (1976) Enzyme pattern directed chemotherapy. The effects of combinations of methotrexate, 5-fluorodeoxyuridine and thymidine on rat hepatoma cells in vitro. Biochem. Pharmac. 25: 2613-2618. JACKSON, R. C., NIETHAMMER,D. and HUENNEKENS,F. M. (1975) Enzymic and transport mechanisms of amethopterin resistance in LI210 mouse leukemia cells. Cancer Biochem. Biophys. 1: 151-155. JACKSON, R. C., HART, L. I. and HARRAP, K. R. (1976) Intrinsic resistance to methotrexate of cultured mammalian cells in relation to the inhibition of their dihydrofolate reductases. Cancer Res. 36 1991-1997. JACQUEZ, J. A. (19721 Compartmental Analysis Biology and Medicine: Kinetics of Distribution of Tracer-Labelled Materials. Elsevier, Amsterdam. JACQUEZ,J. A., BELLMAN,R. and KALABA,R. (1960) Some mathematical aspects of chemotherapy. IL The distribution of a drug in the body. Bull. Math. Biophys. 22: 309-322. JUSKO, W. J. (1971) Pharmacodynamics of chemotherapeutic effects. Dose-time-response relationships for phase-nonspecific agents. J. pharmac. Sci. 60: 892-895. JUSKO, W. J. (1973) A pharmacodynamic model for cell-cycle-specific chemotherapeutic agents. J. Pharmacokinet. Biopharmac. 1: 175-200. KEMBER, N. F. (1969) Growing bones on the computer--some pitfalls of a computer simulation of the effects of radiation on bone growth. Cell Tissue Kinet. 2:11-20. KINO-SMITH, E. A. and MORLEY, A. (19701 Computer simulation of granulopoiesis: Normal and impaired granulopoiesis, Blood 36: 254--262. KROGER-THmMER, E. (1966a) Formal theory of drug dosage regimens. J. theoret. Biol. 13: 212-235. KROGER-Th'IEMF~, E. (1966b) Pharmacokinetics and dose-concentration relationships Proc. 3rd International Pharmacological Meeting, Sap Paulo, Brazil, pp. 63-113. K~OGER-THIEMER,E. (1969) Formal theory of drug dosage regimens. IL The exact plateau effect. J. theoret. Biol. 23: 169-190. LArrHA, L. G., OLrVER, R. and GU~EY, C. W. (19621 Kinetic model of a bone-marrow stem-cell population. Br. J. Haemat. 8: 442--460. LEME, P. R., CREAVEN, P. J., ALLEN, L. M. and B E ~ , M. (1975) Kinetic model for the disposition and metabolism of moderate and high-dose methotrexate in man. Cancer Chemother. Rep. 59:811-817. LINCOLN, T. and F~mEICH, E. J. (1974) Pharmacokinetic simulation: A future means for better control of cancer chemotherapy. Rec. Res. Cancer Res. 49: 103--107. LINCOLN, T. L., AROESTY, J., MEIER, G. and GROSS, J. F. (1974) Computer simulation in the service of chemotherapy. Biomed. 28: 9-16. MATIS, J. H. and HARTLEY, H. O, (1971) Stochastic compartmental analysis: Model and least squares estimation from time series data. Biometrics 27: 77-102. M c I ~ I S , B., KAPADIA,A., EL-ASFoUIU, S. and Lop, T. L. (1975) Stochastic compartmental modeling of the disposition of 5-(3,3-dimethyl-l-triazeno)imidazole-4-carboxamide. Cancer Chemother. Rep. 59: 843-845. MOP,PaSON, P. F., LINCOLN, T. L. and AROESTY, J. (1975) Disposition of cytosine arabinoside and its metabolites: A pharmacokinetic simulation. Cancer Chemother. Rep. $9: 861-876. NICHOL, C. A., GPaNDEY, G. B., MORAN,R. G. and WERKHEISER,W. C. (1972) Combinations of drugs which interact as inhibitors of DNA biosynthesis. Adv. Enz. Regul. 10: 63-80. NIEBERGALL, P. J., SUGITA, E. T. and SCHN~,aE, R. L. (1974) Calculation of plasma level versus time profiles for variable dosing regimens. J. pharmac. Sci. 63: 100-105. NIETHAMMER,D. and JACKSON,R. C. (1975) Changes of molecular properties associated with the development of resistance against methotrexate in human lymphoblastoid cells. Fur. J. Cancer 11: 845-854. PhDGEa'r, W. J. and TSOKOS, C. P. (1970) A stochastic model for chemotherapy: Computer simulation. Math. Biosci. 9: 119-133. PADOETT, W. J. and TSOKOS,C. P. (1972) A stochastic model for chemotherapy: two-organ systems. Int. J. Bio-Med. Comp. 3: 29--41. PARK, D. J. M. (1974) The hierarchical structure of metabolic networks and the construction of efficient metabolic simulators. J. theoret. Biol. 46: 31-74. PRESCOTT, L. F. (1973) Pharmacokinetic drug interactions: Mechanisms of toxicity. Excerpta Medica Int. Congr. Series No. 288. 14: 49-58. PRIC,OGINE, I. and BABLOYANTZ,A. (1972) Dissipative structures and their application to biology. F E B S Syrup. 25: 3-18. PURDUE, P. (1974) Stochastic theory of compartments. Bull. Math. Biol. 36: 305-309. REICH, J. G. and SEL'KOV, E. E. (19741 Mathematical analysis of metabolic networks. F E B S Lett. 40: S119-S127. REICHARD, P. (1972) Control of deoxyribonucleotide synthesis in vitro and in vivo. Adv. Enz. Regul. 10: 3-15. REINER, J. M. (1969) Behavior of Enzyme Systems (2nd edn) Van Nostrand, New York.
280
R . C . JACKSON and K. R. HARRAP
ROT! ROTI, J. L. and OKADA,S. (1973) A mathematical model of the cell cycle of LSi78Y. Cell Tissue Kinet. 6:111-124. ROWLAND, M. and MATIN, S. B. (1973) Kinetics of drug-drug interactions. J. Pharmacokinet. Biopharmac. 1: 553-567. RUBINOW, S. I. (1968) A maturity-time representation for cell populations, Biophys. J. 8: 1055-1073. SANTI, D. V., MCHENRY, C. S. and SOMMER,H. (1974) Mechanism of interaction of thymidylate synthetase with 5-fluorodeoxyuridylate. Biochemistry 13: 471-481. SHACKNEY, S. E. (1970) A computer model for tumor growth and chemotherapy, and its application to Ll210 leukemia treated with cytosine arabinoside. Cancer Chemother. Rep. 54: 399-429. SHEPPARD, C. W. (1969) Computer simulation of stochastic processes through model-sampling (Monte Carlo) techniques. F E B S Lett. 2: S14-$21. SKIPPER, H. E. (1969) Improvement of the model systems. Cancer Res. 29: 2329-2333. SKIPPER, H. E., SCHABEL, F. M., JR. and WILCOX, W. S. (1964) Experimental evaluation of potential anticancer agents. XIII. On the criteria and kinetics associated with 'curability' of experimental leukemias. Cancer Chemother. Rep. 35:1-111. SKIPPER, H. E., SCHABEL, F. M., JR. and WILCOX, W. S. (1967) Experimental evaluation of potential anticancer agents. XXI. Scheduling of arabinosylcytosine to take advantage of its S-phase specificity against leukemia cells. Cancer Chemother. Rep. 51: 125-165. SKIPPER, H. E., SCHABEL, F. M., JR., MELLETT, L. B., MONTGOMERY,J. A., WILKOFF, L. J., LLOYD, H. H. and BROCKMAN,R. W. (1970) Implications of biochemical, cytokinetic, pharmacologic and toxicologic relationships in the design of optimal therapeutic schedules. Cancer Chemother. Rep. 54:431-450. TAKAHASHI, M. (1968) Theoretical basis for cell cycle analysis. II. Further studies on labelled mitosis wave method. J. theoret. Biol. 18: 195-209. TATTERSALL, M. H. N., JACKSON, R. C., CONNORS, T. A. and HARRAP, K. R. (1973) Combination chemotherapy: The interaction of methotrexate and 5-fluorouracil. Eur. J. Cancer 9: 733-739. TATI'ERSALL, M. H. N., JACKSON, R. C., JACKSON,S. T. M. and HARRAP, K. R. (1974) Factors determining cell sensitivity to methotrexate: Studies of folate and deoxyribonucleoside triphosphate pools in five mammalian cell lines. Eur. J. Cancer 10: 819-826. TEORELL, T. (1937) Kinetics of distribution of substances administered to the body. I. The extravascular modes of administration. Archs int. Pharmacodyn. ThOr. 57: 205-225. THAKUR, A. K., RESCIGNO, A. and SCHAFER, D. E. (1973) On the stochastic theory of compartments. II. Multi-compartment systems. Bull. math. Biol. 35: 263-271. TIBAUX, G., FIRKET, H. and HOPPER, A. F. (1969) A simple computer simulation model for growing cell populations: Analysis of synchronized HeLa cell cultures. Cell Tissue Kinet. 2: 333-345. TRUCCO, E. and BROCKWELL, P. J. (1968) Percentage labeled mitoses curves in exponentially growing cell populations. J. theoret. Biol. 20: 321-337. WEBa, J. L. (1963) Enzyme and Metabolic lnhibitors. Vol. 1. Academic Press, New York. WEINBERG, R. and BERKUS, M. (1971) Computer simulation of a living cell. Part I. Int. J. Bio-Med. Comp. 2: 95-120. WERKHEISER, W. C. (1968) A mathematical model of cellular parameters relevant to chemotherapy with methotrexate (MTX). Proc. Am. Assoc. Cancer Res. 9:77 (1968). WERKHEISER, W. C. (1971) Mathematical simulation in chemotherapy. Ann. N. Y. Acad. Sci. 186: 343-358. WERKHEISER, W. C., GR/NDEY, G. B., MORAN, R. G. and NICHOL, C. A. (1973) Mathematical simulation of the interaction of drugs that inhibit deoxyribon'ucleic acid biosynthesis, Mol. Pharmac. 9: 320-329. WESTLEY, J. (1969) Enzymic Catalysis. Harper & Row, New York. Woo, K. B., BRENKUS,L. B. and Wnc, K. M. (1975) Analysis of the effects of antitumor drugs on cell cycle kinetics. Cancer Chemother. Rep. 59: 847--860. ZAHARKO, D. S., DEDR/CK, R. L., BISCHOFF, K. B., LONGSTRETH, J. A. and OLIVER/O, g. T. (1971) Methotrexate tissue distribution: Prediction by a mathematical model. J. natn. Cancer Inst. 46: 775-784.