Predicting Drug Pharmacokinetics and Effect in Vascularized Tumors Using Computer Simulation

John P. Sinek¹, Sandeep Sanga², Xiaoming Zheng¹, and Vittorio Cristini^1,2

Affiliations

Departments of Mathematics¹, Biomedical Engineering², and Medicine—Hematology/Oncology³, University of California, Irvine, USA

Financial Support

National Cancer Institute, National Science Foundation

Corresponding Author and Mailing Address

Vittorio Cristini

Department of Biomedical Engineering, 3120 Natural Sciences II

University of California

Irvine, CA 92697-2715

Phone: (949) 824-9132

Fax: (949) 824-1727

E-mail: cristini@math.uci.edu

Running Title: Predicting Drug PKPD in vivo Using Simulation

Key Words: In Silico, Simulation, Pharmacokinetics, Prediction, Therapy

Abstract

In silico modeling is an increasingly powerful tool in probing processes driving cancerous behavior and, in particular, tumor response to therapeutics. It has the capability of integrating individually well-understood phenomena across a wide range of scales, such as the vascular extravasation of oxygen and drug, interstitial diffusion, and cells’ response to local concentrations in a spatially multi-dimensional setting, thereby giving rise to behavior that is more than the sum of its inputs. As such, the model can be used as an in silico laboratory, investigating tumor therapeutic response under almost any conceivable set of conditions. The results can suggest hypotheses to be corroborated by further in vitro and in vivo experimentation, the convergence of results shedding light on cancer phenomena and treatment.

In this paper we investigate the pharmacokinetics and effect of doxorubicin and cisplatin in three simulated two-dimensional vascularized tumors, taking into account especially vascular and morphological heterogeneity. To do this we construct a multi-compartment PKPD model calibrated from published experimental data and simulate two-hour bolus administrations followed by 18-hour drug washout. Our results show that, although doxorubicin is well noted for its penetration difficulties, the AUC (concentration-time product) of DNA-bound drug continues to become more homogenous long after bolus cessation due to retention in tissue removed from vessels. This leads to a more homogenous final AUC and therefore more homogenous cell inhibition than would otherwise be expected. This is a known and experimentally verified result. We further predict that, under model assumptions, the effect of hypoxia and hypoglycemia (“nutrient effect”) does not appreciably alter this result, although it increases the amount of drug required for 50 percent total inhibition (IC₅₀) in vivo. Cisplatin maintains homogenous drug distribution throughout therapy, and so results in very homogenous AUC and cell inhibition. We then demonstrate that the improvement obtained by therapies that increase doxorubicin penetration is contingent upon tumor characteristics, such as cell density and Pgp efflux. For dense tumors, the improvement is great, while it is negligible for those expressing Pgp. We offer an explanation for this phenomenon. Again, the nutrient effect does not significantly alter these results. Finally, we investigate permeabilization therapy with respect to cisplatin. Results predict that the improvement obtained in monolayer carries over unaffected in vivo, thus implying that assay assisted chemotherapy selection might be especially efficacious for this drug. In summary, the utility of in silico investigation is demonstrated using doxorubicin and cisplatin, with some results accurately post-dicted as a partial validation.

Introduction

While drug resistance to solid tumors is often considered to be a consequence of genetic factors, such as over-expression of anti-apoptotic proteins or production of drug efflux pumps, factors residing at coarser physiological scales may have profound influence on tumor therapeutic response. Lesion mass is a heterogeneous three-dimensional composite of fibrous and connective tissues, stromal components, vasculature, and multiple clones of cancer cells. Atop this intrinsic heterogeneity is layered the anatomical and functional irregularity of tumoral vasculature, characterized by erratic flow, collapsed vessels, diminished oxygen tension, and a large mean tissue-to-vessel distance (fff, SSS, TTT, UUU, VVV). As a consequence, the tumor microenvironment is highly variable, marked by gradients of nutrient (oxygen and glucose), regions of hypoxia, acidity, and necrosis, and heterogeneous proliferation. In order for an anticancer agent to work, it must extravasate, diffuse through lesion, and then be transported into cells, where it must bind to its target and effect cell apoptosis or mitotic inhibition. Clearly, the tumor environment is not conducive to these processes. The vessel bed’s blood flow and spatial distribution hinder uniform extravasation, calling into question the capability of drug molecules to adequately distribute throughout tissue. In vitro experiments by Tannock and others underscore this concern, demonstrating limited penetration of drugs through lesion tissue, especially highly protein-bound molecules like doxorubicin and paclitaxel (ggg, OOO, mmmm). Once a drug molecule has traversed lesion from its point of extravasation and is presented to a cell, the path from extracellular space to intracellular target is fraught with difficulties ranging from protonation due to the acidic environment, which renders anthracyclines incapable of traversing membrane, to intracellular removal by drug efflux pumps, to nuclear processes that effect DNA repair and drug clearance (REFS). In addition to pharmacokinetics, drug pharmacodynamics is equally impaired. Significant hypoxia and hypoglycemia throughout lesion tend to induce quiescence, and there is evidence that under these conditions, the efficacy of chemotherapeutic agents can be reduced (PPP, rrrr). Quiescence may shield cells directly from the action of cycle phase-specific drugs like doxorubicin and cisplatin, and hypoglycemia can induce changes in topoisomerase II detrimental to the functioning of doxorubicin (nnnn, oooo, pppp).

The heterogeneity and three-dimensionality of the tumoral environment presents a challenge to drug assessment, both during development and in the clinic. Whereas a particular drug may show marked activity against a particular cancer line in vitro, its potency may vanish or become far less reliable in vivo. This is evidenced by the differential between positive predictive accuracy of in vitro-assisted therapy selection (around 70%) and negative predictive accuracy (around 90%), a situation not remarkably changed over the years (ccc, eee). Supraoptimal delivery of drug to cultured cells eliminates the gauntlet of biobarriers in vivo described above, precluding the variability they induce. A drug that consistently works in vitro can therefore be expected to only sometimes work in vivo. Unraveling the myriad interactions of therapeutic determinants within the complex three-dimensional tumoral environment is evidently difficult, resulting in high costs of drug development and patient suffering.

Perhaps the crystal ball we are attempting to build is incomplete when made only of glass typically found in the wet-lab; the in silico lab, i.e., computer simulation, might fulfill a key aspect of the lens. A significant capability of in silico experimentation (including simulated assays) is the complete control over and monitoring of the simulated in vivo tumor environment. Moreover, computer modeling can create hypothetical environments and conditions impossible to achieve otherwise, the study of which is nonetheless instrumental in unraveling disease and drug mechanisms. This expansive control, founded upon an adequately mechanistic mathematical basis, could facilitate the discovery of hypotheses as to why certain drugs or therapeutic strategies would or would not be effective, potentially on a patient-by-patient basis. The relative ease and cost-efficiency of performing simulation could furthermore enable a thorough investigation of strategies, revealing the optimal among them. The judicious combination of this burgeoning technology with the capabilities of the wet-lab is an attractive development in both drug discovery and the clinical management of cancer leading to the easing of patients’ burdens.

The 1990’s witnessed explosive growth in the mathematical and computational modeling of vascular and avascular tumors (aaaa, bbbb, cccc, dddd, eeee, ffff, kkk, gggg). It was not long before similar models were employed to investigate lesion response to anticancer drugs (hhhh, jjjj, kkkk, llll). Notable limitations of these models, however, are their one-dimensionality (employing cylindrical or spherical symmetry) and lack of realistic vasculature. Without true multi-dimensionality and discrete vasculature, it is difficult to simulate the heterogeneities of nutrient and drug that profoundly affect therapeutic efficacy. In the present paper we use the two-dimensional (non-symmetric) tumoral growth engine of Zheng et al., (hhh) which employs the discrete vasculature algorithm of Chaplain and Anderson (mmm), to perform chemotherapy simulations using cisplatin and doxorubicin. In particular, we focus on the effects of drug and nutrient distribution heterogeneity and their impact on in vivo therapeutic efficacy. Sinek et al. had earlier performed a similar investigation (nnn); however, the pharmacokinetics and pharmacodynamics (PKPD) component was rudimentary, assuming one homogenous lesion compartment and not based upon externally acquired parameter values. In the present work we implement an extensive multi-compartment PKPD component whose parameter values are calibrated via published experimental data. This enables a comparison of the tissue- and cell-level drug dynamics of the two drugs, and facilitates the generation of hypotheses to explain their in vivo characteristics. We ask that the reader consider that if doxorubicin and cisplatin were discovered only today, the simulations herein could be seen as providing great insight into their in vivo performance, potentially streamlining and reducing costs of development and clinical trials, and of assisting in clinical therapeutic strategy to improve patient comfort and survival.

Materials and Methods

Model Description

The multiscale tumor growth and angiogenesis simulator developed by Zheng et al. (hhh) is used to grow the lesions upon which we simulate chemotherapy (see Figure AUC). The simulation field incorporates three phases: viable cancerous tissue, normal host tissue, and necrotic debris. The lesion/host interface is demarked by thick black contours, while the microvasculature appears as thin red curves. Dark regions are necrotic debris. Briefly, nutrient (oxygen and glucose) is provided through the discrete microvasculature, which is generated in response to angiogenic regulators produced from perinecrotic cells. This results in proliferation and tumor growth. The simple steady-state diffusion equation

(1)

is used to model nutrient delivery and uptake, where n is the local nutrient normalized by the intravascular level, k_v is a measure of vascular porosity (0 is impermeable, ¥ is completely porous), d is the Dirac delta function located along the vasculature, D_n is nutrient diffusivity, and k_n is the local rate of consumption by cells. The characteristically high porosity of tumor vasculature (REF) implies a very high setting of k_v so that, essentially, vasculature provides a constant boundary condition of 1. Experiments given in (rrr) demonstrate that oxygen penetrates approximately 150 mm into in vitro spheroids before falling to about 10 percent of serum level. At this point necrosis ensues. Combining this with a diffusivity D_n of around 60,000 mm²min^-1 (sss, ttt), the nutrient uptake rate is calculated to be k_n = 24 min^-1. Waste resulting from necrotic cell degradation is assumed to be removed via convection towards and through the tumor-host interface as well as via scavenger cell phagocytosis. In regions where nutrient is sufficient to maintain viability, mitosis is assumed to be directly proportional to its concentration, with the proportionality constant dependent upon the average cell cycle time of the malignant population.

Once the tumors are grown, drug administration via the vasculature is simulated by our multi-compartment pharmacokinetics model. For cisplatin, there are three compartments corresponding to (1) extracellular, (2) cytosolic, and (3) DNA-bound drug. For doxorubicin, there is a fourth compartment corresponding to intracellular organelles, e.g., lysosomes. The system of equations governing transport for both drugs is

(2)

where s_i represents drug concentration in compartment i, k_ij represents a transfer rate from compartment i to j, and k_i represents a rate of permanent removal from compartment i. s_v is intravascular drug concentration during bolus, and s_M is a DNA saturation parameter relevant to doxorubicin. V_C is the volume of a cell (assumed spherical with diameter 10 mm, yielding V_C = 524 fL) and appears in the first equation to reconcile the dimensions of s₁ (mM) with the dimensions of all other compartments (fmoles/cell). k_v and d are the same as in Equation 1. The primed rates appearing in the first equation are related to their unprimed counterparts via a constant due to the fact that extracellular volume is only a fraction F of total tumor volume. Taking a baseline tumor density of 1.0E9 cells/ml in combination with the cell volume previously quoted results in an extracellular fraction of 0.48, all values within standard limits (REFS). Finally, D_s is the diffusivity of the drug through interstitial space.

Both cisplatin and doxorubicin pass through cell membrane according to k₁₂ (which includes possible pump and transporter activity, as do all other rates). From there, the drugs may efflux according to k₂₁ or may bind to DNA according to k₂₃. The kinetics differ from here for the two drugs. Cisplatin may be removed according to the rate k₃, which destroys the functioning of the drug and repairs the DNA. Doxorubicin, however, has an off rate given by k₃₂, and moreover may be sequestered and released by lysosomes according to k₂₄ and k₄₂. Although lysosomal flow to membrane and exocytosis of sequestered drug plays a role in some drug resistant cell lines, we are not necessarily modeling drug resistance via this function, and so assume this process to be negligible, a valid assumption (Demant et al., 1990). On the other hand, we are concerned with the quantity of drug lysosomes can sequester, as this contributes to the cellular uptake of drug, and hence, its penetration characteristics.

The pharmacodynamics model consists of a simple Hill-type equation

(3)

where E is cell inhibition (1 minus surviving fraction), x is DNA-bound drug-time product (area under the curve, or AUC), and A and m are phenomenologically fit parameters. N(n) is a function of nutrient n ranging from 0 to 1 used to mimic the effect of hypoxia and hypoglycemia. Results with doxorubicin show that cells in deeper layers of spheroids do not respond as well to drug as do cells on the surface, even when intracellular drug levels are taken into account (PPP, rrrr). Other experiments demonstrate reduced response in monolayer when cells are forced into quiescence due to reduced oxygen (nnnn). Still others show that hypoglycemia can deplete Topoisomerase II, thus reducing the effect of some anthracyclines (pppp). These results imply that the response of cells to doxorubicin in vivo might correlate to the local nutrient, which we herein refer to as the “nutrient effect.” For our purposes, the exact form of N is not important. For simplicity, we choose N = n^p, where p is a phenomenological parameter derived from the data of (PPP), and equals 0.4. Since in our model n is normalized with respect to the intravascular level, it runs from 0 to 1, and thus so does N. Furthermore, at full nutrient levels, N = 1, and so cell inhibition (3) is maximal. In our simulations, drug pharmacokinetics (Equation 2) is allowed to proceed from bolus initiation to washout 20 hours later. During this time the locally varying DNA-bound AUC is calculated and used to find cell inhibition (Equation 3).

Pharmacokinetics Model Parameters

It will be useful to bear in mind that a generally acceptable theoretical setup for performing experiments to measure compartmental concentrations (and therefore to derive the rate constants we are after) is either a suspension or monolayer in an inexhaustible drug-laden medium corresponding to s₁. Under these conditions, the relevant model consists of the last three equations in Equation 2, with s₁ held constant. We will refer to this model as the modified version of Equation 2. All model parameters and values are summarized in Figure Parameters along with references. These will be referred to as the baseline values, some of which will be adjusted later to simulate different tumor characteristics and therapeutic treatments. We emphasize that parameter values, having been derived from a variety of published experimental data spanning many years and cell types, correspond to a prototypical tumor and cancer cells suitable for the simulations herein, but not necessarily representative of any particular clinical specimen.

We begin with cisplatin, setting k₂₄, k₄₂, and k₄₁ to 0 since we assume only three compartments, and k₃₂ to 0 since we assume the repair rate k₃ is the dominant removal rate of DNA-bound drug. k₃ is next obtained as follows. In experiments performed by Sadowitz et al. (BBB), adducts per million nt on isolated peripheral blood mononuclear cell DNA fell from 75 to 5 and 185 to 40 in two hours in two different experiments. Thus, assuming the exponential repair model we calculate the repair rate to be about 0.015 min^-1. An initial estimate of k₂₃ is then made as follows. Sadowitz shows that for 7 mM cis, in two hours peripheral blood mononuclear cells accumulate from about 25 (non-thiol-blocked cells) to 175 (thiol-blocked cells) adducts per million nt. Assuming that DNA consists of about 1.25 E6 kbp, this converts to from 1.04 E-19 to 7.27 E-19 moles of Pt docked on the DNA (1 atom/adduct). Neglecting the cell membrane and supposing the DNA to be exposed directly to the drug, we have the ODE , where l₂₃ is a clearance parameter (fL-min^-1). The solution of this is . Substituting values of k₃ = 0.015 min^-1, t = 120 min, and 1.04 E-4 £ s₃ £ 7.27 E-4 fmole yields 0.27 £ l₂₃ £ 1.9 fL-min^-1. To convert this to a rate we use the relation k₂₃ = l₂₃/V_C, arriving at 3.82E-4 min^-1. The assumption that DNA was exposed directly to the cis solution means that this rate is only a bootstrap approximation and must be refined. At this point we note that the extremely low ratio of adducts per kbp implies that the saturation capacity of DNA with respect to cisplatin is never approached, and so set s_M to µ.

Figure Parameters: A complete summary of pharmacokinetics and pharmacodynamics parameters along with references, Tumor growth and angiogenesis parameters can be found in (hhh).

Next, we estimate k₁₂ and k₂₁. While doing this we will refine our initial estimate of k₂₃. The whole procedure involves fitting the best curves to Troger et al.’s data (AAA) (Figure Troger). Troger exposed CAL-27 cells in monolayer to four different concentrations of cisplatin and then measured the total intracellular amount of Pt at selected times. This corresponds to s₂+s₃ in our model. Beginning with the previous estimate of k₂₃ and setting s₁ to the concentrations used by Troger, we adjust k₁₂ and k₂₁ in the modified version of Equation 2 until a good fit of Troger's data is obtained. Simultaneously, we adjust k₂₃ to keep the DNA-bound drug true to the results of (BBB) previously discussed. We remark that the disparity between the inward and outward rates derived for cisplatin may be due in part to carrier-mediated transport, e.g., the CTR1 influx transporter.

Figure Troger: Data from (AAA) used to calibrate k₁₂, k₂₁, and k₂₃ for the cisplatin model. Parameters are fit simultaneously to all four curves treated as one set of data; they are not tailored to each curve.

Proceeding to doxorubicin we first obtain an acceptable range for k₁₂ and k₂₁ from the literature. For a variety of anthracyclines, including doxorubicin, initial estimates of cell membrane permeability P are taken from Dordal et al.'s (1995) experiments with SU-4 and SU-4R wildtype and resistant human lymphoma cells, from Demant and Friche's (1998) and Demant et al.'s (1990) experiments with EHR2 and EHR2/DNR+ wildtype and resistant Ehrlich ascites tumor cells, and from Lankelma et al.'s (2000) experiments with MDA-468 breast cancer cells. The range reported is 2.4 £ P £ 1000 mm-min^-1. The relation k₁₂ = PA_C/V_C, where A_C represents the cell membrane area, can then be used to arrive at an initial range of 1.4 £ k₁₂ £ 600 min^-1, which will be refined later. In the case of passive diffusion, k₂₁ = k₁₂. We note that these values are far larger than those obtained for cisplatin previously. More generally, it has been remarked that cell membrane permeability for cis is much lower than for doxorubicin, etoposide, and vinblastine, although all four drugs are thought to enter cells by passive diffusion (Jekunen et al., 1993).

We next turn our attention to DNA-binding affinity. Given the great DNA affinity of the anthracyclines, saturability of the DNA must be taken into account, requiring an estimate of s_M. A difficulty arises since there is evidence a typical anthracycline molecule intercalation occludes from 3 to 10 binding sites in a manner that cannot be corrected exactly by a factor (Rizzo et al., 1989; Tarasiuk et al., 1989; McGhee and von Hippel, 1974). Again, to a first approximation, we assume that such a correction can be applied. Demant and Friche (1998) report a DNA binding site concentration of about 5 mM within a cell volume of 1000 fL, yielding 5 fmoles of sites. A low value of 0.7 fmoles is obtained by using our assumed value of 1.25 E6 kbp and Rizzo et al.'s (1989) reported site exclusion parameter of about 3. Tarasiuk et al. (1989) find that the DNA of human lymphocytes is comprised of about 6.0 E6 kbp and that one intercalating molecule of doxorubicin requires 10 base pairs. Thus, Tarasiuk's data implies a factor-corrected quantity of 1 fmole of binding sites, which we take as a representative value of s_M.

DNA binding kinetics of the anthracyclines is nontrivial, perhaps requiring multiple steps and demonstrating sequence specificity (Rizzo et al., 1989; Qu et al., 2001). Bearing this in mind, as an approximation it will suffice to assume non-specific, one-step binding and unbinding according to the chemical reaction

.

A representative value for the binding coefficient in the above equation for doxorubicin is reported as k_on = 4.2 E8 M^-1min^-1 and a value of the unbinding coefficient as k_off = 1800 min^-1 (Rizzo et al., 1989). k_off is identical with k₃₂. From k_on we calculate a clearance parameter (as with cisplatin) given as l₂₃ = k_ons_M = 4.2 E8 fL-min^-1 (being cautious with the scales of our dimensions). k₂₃ can then be calculated as l₂₃/V_C, given in Figure Parameters.

We next turn our attention to the rates k₂₄ and k₄₂ governing lysosomal sequestration. Experiments of Hurwitz et al. (1997) using U-937 myeloid leukemia cells and their dox-resistant variant U-A10 show that the ratio of DNA-bound to lysosomally-sequestered drug is about 3. (Hurwitz uses daunorubicin, an anthracycline related to doxorubicin.) In our modified model equations with all other parameters set as described above, the amount of sequestered drug at equilibrium is dependent only upon the ratio k₂₄/k₄₂. This ratio furthermore does not affect the equilibrium quantity of DNA-bound drug. Arbitrarily selecting k₂₄ = 1 min^-1, we find that the appropriate DNA-bound to lysosomally-sequestered ratio is obtained by setting k₄₂ to 0.007. Considering that lysosomal membrane permeability is quite high (Demant et al., 1990), the lysosomally-bound drug must achieve equilibrium quickly, which can be modified by changing k₂₄ while keeping the ratio k₂₄/k₄₂ constant. We find that increasing k₂₄ by a factor of 10 reduces the time required to achieve 95% of equilibrium value (max₉₅) to about 300 minutes, below which further increases in k₂₄ only reduce this time negligibly. Thus, we conservatively set k₂₄ = 10 and k₄₂ = 0.07.

To refine our initial range of k₁₂ and k₂₁, we use the modified PK model to compare our simulated monolayer uptake profiles of total intracellular drug with those of DeGregorio et al (1984) using human Ewing’s sarcoma and rhabdomyosarcoma cells. At 5.40 min^-1 both equilibrium values and uptake rates compare favorably at three test concentrations.

The last pharmacokinetics parameter values needed are the diffusivities D_s of cisplatin and doxorubicin through tumor interstitium. For molecules of their size (dox M.W. = 544, cis M.W. = 300), diffusivity should be about 30,000 mm²-min^-1 (sss, ttt). However, doxorubicin faces particularly severe barriers due to its binding to extracellular constituents such as hyaluronic acid (mmmm, qqqq), and its diffusivity in some tissues has been estimated to be as low as 1000 mm²-min^-1, which we take as our baseline value (www).

Pharmacodynamics Model Parameters

In order to calibrate the pharmacodynamics model (Equation 3), we use the in vitro data of Levasseur et al. (NNN) with A2780 ovarian cancer cells exposed in monolayer to both doxorubicin and cisplatin over a range of times and concentrations. We assume the previously discussed modified pharmacokinetics model along with the values derived, and simulate Levassuer’s exposures followed by approximately 24 hours of drug washout in drug-free medium (s₁ is set to 0). During this time, DNA-bound AUC is calculated. These data are then used in conjunction with Levasseur’s surviving fraction data to fit the parameters A and m in Equation 3 using Excel (Figure PD). During this process, nutrient is assumed plentiful (n = 1) thus bypassing the nutrient effect.

Figure PD: Cell inhibition fits for Equation 3 using Levasseur’s data on A2780 ovarian cancer cells exposed in monolayer.

In Silico Experiments

We first grow three vascularized two-dimensional in silico lesions, shown in Figure AUC. Each lesion is produced based upon the same set of growth and vasculature control parameters (see (hhh) for a complete description of these), but randomness in the vasculogenesis algorithm and slightly different initial shapes produce different morphologies. We then perform simulations to demonstrate and analyze distributions of DNA-bound drug, DNA-bound AUC (concentration-time product), and cell inhibition resulting from intravascular bolus administrations of cisplatin and doxorubicin. Each experiment is replicated in each of the three lesions. In each case we hold the intravascular concentration of drug (s_V in Eq. 2) constant for two hours, then set it to zero for eighteen more hours to allow washout. Although this sharp “square wave” is perhaps a caricature of clinical bolus administration, it allows for consistent analysis and comparison of results. Intravascular concentrations are calibrated in each case to produce a total cellular growth inhibition of 50 percent (IC₅₀ concentration). It is assumed that a true in vivo tumor does not grow or regress appreciably during the 20 hour course of the therapy we are attempting to simulate, hence we freeze tumor and vascular growth during our in silico therapies

Our first set of experiments compares DNA-bound drug distributions of doxorubicin and cisplatin under baseline conditions (see Figure Parameters) and examines the effect of doxorubicin retention on final DNA-bound AUC and cell inhibition. We next investigate the impact of inhibition heterogeneity on dosing requirements, paying particular attention to the nutrient effect for doxorubicin. In our third set of simulations we show the effect of therapies designed to improve doxorubicin penetration by, for example, removing hyaluronic acid (mmmm, qqqq). We investigate this under three circumstances: baseline tumor density, high tumor density, and baseline tumor density with Pgp efflux activity. These are chosen because they demonstrate a spectrum of possibilities due to their effect on cellular drug uptake. High tumor density increases uptake, while Pgp efflux decreases it. In order to simulate increased penetration, we increase D_s for doxorubicin from its baseline value to 5000 m²-min^-1 for a moderate increase, and 30,000 for the maximum increase, thus matching the performance of cisplatin. To simulate high tumor density we increase r by 50 percent to 1.5 E9 cells-ml^-1 (REF). This has the effect of lowering the interstitial fraction F to 0.22, which in turn increases k’₁₂ and k’₂₁ while leaving all other rates unchanged. Pgp efflux is simulated by increasing k₂₁ by a factor of 10, which has the effect of reducing all intracellular compartment quantities by approximately the same factor. This is consistent with results of (ssss) that show Pgp activity can reduce intracellular concentrations of daunorubicin (an anthracycline related to doxorubicin) by up to a factor of 100. In our fourth and final set of experiments we investigate permeabilization therapy with respect to cisplatin, whereby a detergent, such as digitonin, or electropermeabilization is used to increase the permeability of cell membrane (Tanaka et al., 2001; Jekunen et al, 1993). We take an extreme case, increasing the rate constants k₁₂ and k₂₁ from baseline both by a factor of 100. Note that this does not increase the limiting intracellular or DNA-bound levels of drug attained in simulated monolayer, only the rate at which these come to equilibrium. Thus highly permeabilized, max₉₅ is attained at 3.4 hours of exposure; further permeabilization reduces this negligibly. For comparison, max₉₅ takes longer than 27 hours for unpermeabilized cells. This therapy is simulated under both baseline and very high cell densities achieved by increasing the baseline density 75 percent to 1.75E9 cells-ml^-1. Both of these are further compared to monolayer results to attempt to probe the conditions under which in vitro assays can be used to predict clinical efficacy.

Results

Although all treatments described in this section are duplicated in each of the three in silico tumors, we display only representative plots with appropriate summaries of all data. Note that the nutrient effect is only used where noted.

We begin by examining DNA-bound AUC distributions at various times in the baseline simulated lesions (each lesion corresponding to a column, A, B, or C), shown in Figure AUC. From top to bottom, the times correspond to two hours, eight hours, fourteen hours, and twenty hours post bolus initiation. Levels are normalized relative to the average AUC within viable lesion for comparison of heterogeneity. Although surrounding host tissue cells uptake and bind with drug differently than cancer cells, we make no distinction in these plots; however, analytical results only consider DNA-bound drug within viable lesion. The two left column sequences (Lesions A and B) show doxorubicin AUC, while the rightmost column shows cisplatin. For both Lesions A and B, at 2 hours doxorubicin AUC is seen to be about 3 times the average (dark red) close the the vasculature, and almost 0 (blue) elsewhere. The distribution is only slightly more homogeneous by 8 hours. By 14 hours the heterogeneity has lessened, with the peaks close to the vasculature reaching only about 2.2. Finally, at the conclusion of washout 20 hours after bolus initiation, the distribution has become much more homogeneous, with the peaks only reaching about 1.7 times the average. In contrast, the cisplatin distribution within Lesion C remains extremely homogeneous, right at the average, throughout the entire treatment.

The probability distributions at the bottom, corresponding to AUC at 20 hours post bolus initiation, allow for a more quantitative comparison. The two corresponding to doxorubicin show some heterogeneity relative to cisplatin on the right. Using the leftmost distribution as an example, the average AUC is found to be 6.04 fmole-min. 25 percent of tumor cells receive less than 1.66 fmole-min each, while 25 percent of tumor cells receive more than 9.54 fmole-min. The remaining 50 percent of the tumor cells receive between these two values, a range of 7.88 fmole-min. When normalized with respect to the average and expressed as a percent, this yields 131% (the interquartile range, or IQR), and gives a concise measure of distribution heterogeneity. IQR’s are given at each of the other timepoints as well. All three tumors, despite varied lesion and vasculature morphologies, demonstrate similar results (not all shown). Doxorubicin AUC IQR’s typically lessen from about 250 percent at 2 hours to 150 percent at 20 hours; cisplatin AUC IQR’s drop from about 10 to 2 percent. Interestingly, in the run shown, the heterogeneity for cisplatin increases slightly from 14 to 20 hours. This happens in some of the other cisplatin simulations as well.

Figure AUC: DNA-bound AUC at four times (rows: 2, 8, 14, and 20 hours) post bolus initiation for three two-dimensional simulated tumor lesions (columns). Results are normalized to average lesion AUC at the time taken to enable comparison of distribution heterogeneities. Thick black contours are tumor boundaries. Thin red curves are vasculature. Dark regions are necrotic areas. Each unit represent 200 mm. Bottom probability distributions show final AUC distributions at 20 hours. A concise measure of heterogeneity is given by the inter-quartile range (IQR), depicted in the lower left graph and explained in the text. Although AUC in host tissue is also shown in plots, analysis considers only DNA-bound drug in viable lesion.

We next investigate the inhibition distributions and IC₅₀’s corresponding to the treatments to the baseline lesions depicted in Figure AUC. Particular attention is paid to the nutrient effect for doxorubicin. A table of average IC₅₀’s and log(IC₅₀/ IC_50,mono)’s for the two drugs in the three replicate lesions is given in Figure IC50. Here and throughout this paper IC_50,mono refers to baseline cells exposed in monolayer, while IC₅₀ refers to cells in lesion. Note that, as these are simulated monolayer exposures, IC_50,mono is deterministic. At the 5% significance level, one-tailed t-tests show that the logs IC₅₀ ratios are greater than 0 for doxorubicin while for cisplatin it is not. This underscores the homogeneity of cisplatin. Figure IC50 also shows a typical nutrient profile, using Lesion B as an example with an IQR of 36%. This measurement is completely analogous to those used in Figure AUC, except that here it is applied to the nutrient distribution and there is no normalization since nutrient levels are bounded absolutely from 0 to 100 percent, the level within the vasculature, itself. The nutrient IQR’s for the other two lesions are within 2% of this value. For doxorubicin, average IC₅₀ for the three replicate baseline lesions is increased by 0.406 log units, or by a factor of about 2.5, when the nutrient effect is used. A paired t-test shows that this difference is significant at the 5% level.

Figure IC50: Means ± SD’s of the IC₅₀’s and the logs of their ratios with respect to monolayer treatments for the three replicate baseline lesions. IC_50,mono is the IC₅₀ of baseline cells in monolayer. At the 5% significance level using a one-tailed t-test, the log ratios for doxorubicin exceed 0 with or without the nutrient effect, while cisplatin does not. A paired one-tailed t-test shows that the log IC₅₀ ratio for doxorubicin with the nutrient effect is greater than that without. Contour plot shows nutrient distribution in Lesion B demonstrating significant heterogeneity. Other lesions are similar.

The inhibition distributions for all three baseline lesions in Figure AUC mirror their respective final AUC distributions: those for doxorubicin are heterogeneous, while that for cisplatin is virtually uniform at 50 precent inhibition throughout. Using Lesion B as a representative example for doxorubicin, the upper pair of frames in Figure Cell_Inhibition demonstrates the cell inhibition distribution for the baseline lesion along with its cumulative probability plot. The second pair of frames shows the result of adding the nutrient effect. While the broadening of the probability plot indicates an apparent increase in heterogeneity, the IQR is reduced from 81 to 77 percent (again, not normalized). A calculation and comparison of the variances (not shown) confirms this reduction. To investigate this further, we simulate therapy designed to improve doxorubicin penetration by, for example, removing hyaluronic acid (mmmm, qqqq). With D_s in Equation 2 set to the maximum of 30,000 mm²-min^-1, the third and fourth pairs in Figure Cell_Inhibition now reveal that inhibition heterogeneity has substantially increased due to the nutrient effect. Lesions A and C yield similar results, with inhibition heterogeneity actually decreasing under baseline conditions, but increasing significantly when drug penetration is augmented.

Figure Cell_Inhibition: Top pair shows the cell growth inhibition profile of Lesion B resulting from the treatment depicted in Figure AUC. Probability plot and IQR are now of inhibition distribution and are not normalized with respect to any average. Second pair shows the effect of including nutrient heterogeneity. Although the broadening of the probability plot indicates increased heterogeneity, the IQR indicates decreased. Bottom two pairs repeat the same experiment, except with doxorubicin penetration increased. Now both the plots and IQR’s show dramatically increased heterogeneity.

In our third set of simulations, we investigate the effect of therapies designed to improve doxorubicin penetration under several combinations of drug/interstitum diffusivities, cell densities, and drug efflux activities (e.g., Pgp).Figure Penetration gives bar graphs of log(IC₅₀/IC_50,mono)’s for three scenarios. The leftmost triplet corresponds to baseline tumor density and no efflux, resulting in a condition of “normal” cellular uptake. The middle triplet corresponds to high density with no efflux, a condition of high uptake. The rightmost corresponds to baseline density with efflux, a condition of low uptake. In the baseline tumor case there is a change of -0.388 log units in going from no removal of hyaluronic acid to almost complete removal. When density is increased, the change becomes –0.709; however, when Pgp efflux is activated, ANOVA reveals there is no statistical difference, and in fact, the measured change is positive. Results are similar when the nutrient effect is included, with all bars essentially increased by a constant, approximately 3.7.

In our fourth and final set of simulations we investigate the effect of permeabilization therapy vis-à-vis cisplatin. Figure Splatin shows log(IC_50,perm/IC_50,unperm) for three cases: monolayer, in vivo with baseline cell density, and in vivo with high cell density. Here, the subscripts “perm” and “unperm” denote the application or withholding of permeabilization therapy. Permeabilization results in a decrease of 0.154 log IC₅₀ units for simulated monolayers, i.e., a reduction of IC₅₀ by a factor of 0.7, and is thus effective in vitro. An interesting question is whether this carries over in vivo, i.e., whether a monolayer assay can be used to predict clinical efficacy. Figure Splatin demonstrates this by displaying the log IC₅₀ differences between permeabilization and standard therapy under both normal and high density in vivo conditions. Improvements are comparable to monolayer results, with all three log-differences about –0.14, and no statistical difference for the baseline case at the 5% significance level using a two-tailed t-test.

Figure Penetration: The effect of increasing doxorubicin penetration is shown in three cases: baseline tumor, high-density tumor, and tumor with Pgp efflux. High density has the effect of increasing drug uptake, while Pgp efflux has the opposite effect. Bars are log(IC₅₀/IC_50,mono), where IC_50,mono is the IC₅₀ for baseline cells exposed in monolayer. For reference, its value is 0.175 mM. Results are similar for simulations carried out under the effects of hypoxia and hypoglycemia, with all bars approximately increased by a constant. Three replications per bar with results of ANOVA displayed.

Figure Splatin: Effect of permeabilization therapy with respect to cisplatin is shown in three cases. Bars are of log(IC_50,perm/IC_50,unperm) where IC_50,perm and IC_50,unperm correspond to permeabilized and unpermeabilized conditions. Three replications per bar with results of two-tailed t-tests relative to monolayer displayed. While there is a statistical difference at the 0.05 significance level for the high-density tumor, this disappears at the 0.01 significance level.

Discussion

That doxorubicin should perform well clinically given its well-noted penetration difficulties can be somewhat explained by its retention in tissue removed from vasculature. This is demonstrated in Figure AUC, where it can be seen that homogeneity of exposure increases long after the bolus has been terminated. This phenomenon has been experimentally verified in (PPP) with spheroids. Because of this, the resulting cell inhibition distribution is more homogenous than would otherwise be expected. On the other hand, cisplatin maintains an homogenous DNA-bound distribution at all times from bolus initiation to 20 hours later, resulting in an extremely uniform cell inhibition distribution.

It is reasonable to expect that heterogeneity of nutrient, resulting in hypoxia and hypoglycemia, should compound resulting heterogeneity of cell inhibition for doxorubicin. Interestingly, this proves not to be the case for the baseline lesions, as Figure Cell_Inhibition shows. However, when doxorubicin penetration is increased, resulting in a more homogeneous AUC, the difference becomes manifest. In fact, this difference is maximized when the AUC is most homogeneous, such as when Pgp efflux is activated and drug penetration is maximized. Conversely, in cases where the AUC is rendered very heterogeneous, such as with increased density and no penetration therapy, the inhibition distribution becomes more homogeneous as measured by the IQR when the nutrient effect is employed. This puzzling behavior can be explained by examining the probability distributions in Figure Cell_Inhibition. In the top graph, the bottom and top quartile endpoints yielding the interquartile range occur at nearly 0% and approximately 82% inhibition. This exhausts almost the entire range, leaving little room for the upper quartile to grow, and no room for the lower quartile to drop. Indeed, the broadening of the curve in the second graph actually indicates greater heterogeneity, but this effect is localized around the region where the curve was vertical, forcing the upper quartile, hence IQR, to drop. The situation is different in the third and fourth probability plots, with sufficient room for both of the quartile endpoints to separate. In a sense, the greater heterogeneity of doxorubicin AUC in the first case masks any further potential effect of nutrient heterogeneity.

It is worthwhile to consider the inhibition distribution when treating clinical tumors for at least two reasons. First, increased heterogeneity typically results in higher dosing requirements since the factors that precipitate heterogeneity usually do so by decreasing drug effectiveness in regions. This observation is borne out by the IC₅₀ comparisons given in Figure IC50, where the average IC₅₀ for doxorubicin in the baseline tumors with the nutrient effect is about 2.5 times greater than without. It is further supported by noting that the IC₅₀ for cisplatin, which appears to have a very homogeneous AUC, is not statistically different from its IC_50,mono. A second and subtler reason is that heterogeneities in growth and regression have been linked to increased lesion fragmentation and invasiveness (ooo, ppp, XXX, AAAA). While the mechanisms underlying this phenomenon are complex, involving myriad protein signaling events and activities at the cellular level, they may at least partly rely on gross lesion effects.

The results in Figure Penetration show what might be expected from therapies that increase doxorubicin penetration by, for example, removing hyaluronic acid. As expected, for the baseline tumors, greater homogeneity and level of AUC is achieved, resulting in reductions of IC₅₀. That this effect should be more pronounced for high-density tumors and completely absent in the presence of Pgp efflux is intriguing. A potential explanation is availed by grossly oversimplifying the pharmacokinetics model (Equation 2), reducing it to the one-dimensional steady state diffusion equation with diffusivity D and uptake rate k. In two dimensions, a segment of blood vessel acting as a source next to a section of tissue approximates the one-dimensional case. This equation has one governing parameter, the characteristic diffusion length , and a solution , where x is distance from the source. Considering a section of tissue of thickness d next to a vessel, AUC is proportional to the integral, , and is inversely proportional to the IC₅₀. Increasing L by some factor C to simulate penetration therapy results in a new characteristic length of CL. The ratio of IC₅₀’s is therefore , which for large characteristic lengths approaches 1, and for small lengths approaches 1/C. Now, increasing cell density has approximately the effect of increasing k, resulting in a small L, while activating Pgp efflux has the effect of decreasing k, resulting in a large L, thus explaining the phenomenon in question. This could be an important point when deciding upon appropriate therapies for tumors exhibiting different characteristics such as efflux mechanisms and relatively high or low densities.

With respect to penetration, cisplatin behaves in vivo almost as if it were acting on a monolayer. This is evident by the similar IC₅₀’s in both cases given in Figure IC50 as well as by the homogenous AUC distributions shown in Figure AUC. This would help to explain the fact that the improvement obtained by permeabilization therapy is similar in the simulated tumors (even very dense ones not likely to be encountered) and monolayer, as shown in Figure Splatin. Perhaps for drugs like cisplatin, which demonstrate little penetration difficulties, in vitro assays used to predict the effect of a therapeutic strategy might carry over mostly unchanged to the in vivo case under a variety of circumstances.

While it can and should be argued that the simulations herein fail to account for some (many!) critical aspects of tumor growth (such as clonal heterogeneity), and that parameter settings may in some cases be inexact, it should not be concluded that these shortcomings invalidate the characteristics the simulations have revealed. Indeed, we have post-dicted results regarding doxorubicin retention and drug and inhibition hetereogeneites given in (PPP) (Other refs?) by only using fundamental cell, tissue, and drug parameters from published experiments. Thus, our results are mechanistically founded, arising naturally through the functioning of the model.

One of our broader goals is to demonstrate how increasingly sophisticated in silico technology, driven by mathematical modeling and calibrated with experimental data, can and is being developed to provide an alternate investigative and clinical tool complementary to traditional methods. The power of in vitro experimentation lies in its ease of implementation while remaining in the biological realm. Yet this is obtained at the cost of accuracy. By its very nature, in vitro experimentation attempts to refine and isolate. Yet, much of what happens in vivo is the result of a nonlinear system whose behavior is more than the sum of its parts. The power of in-silico simulation lies in its ability to integrate components into a virtual system capable of reproducing such behavior, implicitly taking into account circuits of information flow difficult to explicitly analyze. Accurately calibrated and rigorously validated, such an integrated model could provide a “dry-lab” to be used as a powerful complement to the traditional wet-lab in fundamental research, drug discovery, and in the clinic. “Help bridge the invitro/in vivo gap”.

It can well be imagined that were doxorubicin and cisplatin discovered today, the in vivo simulations herein presented could be used to anticipate their lesion- and cellular-scale pharmacokinetics, helping to refine clinical trial design and lower costs. In clinical application, the results could be used to guide therapeutic strategy. For example, any risks associated with doxorubicin penetration therapy could be avoided if it were known that the patient’s tumor were expressing Pgp or otherwise had lowered cellular uptake, according to the results given in Figure Penetration. As a research tool, the dry-lab could be used to probe scenarios and test hypotheses that are either difficult or impossible to instantiate in the body. We have made use of this feature in Figure Cell_Inhibition where we isolated the effect of hypoxia and in Figure Splatin where we have produced a tumor of density unrealizable in the wet-lab so as to probe therapeutic limits. Results could then suggest supportive in vitro and in vivo experimentation, the end result being new therapeutic targets or strategies. Simultaneously, weaknesses (or strengths!) of the in-silico model could be uncovered and addressed. Above a certain threshold of reliability, the model could play a significant role in the clinic, using patient-specific histological and genetic information to set simulator parameters, allowing accurate prediction of response to proposed therapies.

Acknowledgments

The authors would like to thank Steve Wise (U.C. Irvine) for assistance with numerical methods and programming, Ardith El-Kareh (U. of Arizona) for insights into the cellular pharmacokinetics and pharmacodynamics of doxorubicin and cisplatin, Hermann Frieboes (U.C. Irvine) for helpful discussions regarding in vitro experimental protocols and drug response, and John Fruehauf (U.C. Irvine) for suggestions and corrections covering a broad spectrum of oncological knowledge.

References

aaa. Black MM, Spear FD. Effects of cancer chemotherapeutic agents on dehydrogenase activity of human cancer tissue in vitro. Am J Clin Pathol 1953;23:218-227.

bbb. Black MM, Spear FD. Further observations on the effects of cancer chemotherapeutic agents on the in vitro dehydrogenase activity of cancer tissue. J Natl Cancer Inst 1954;14:1147-1158.

ccc. Fruehauf JP. In vitro assay-assisted treatment selection for women with breast or ovarian cancer. Endocr Relat Cancer 2002;9:171-182.

ddd. Loizzi V, Chan JK, Osann K, Cappuccini F, DiSaia PJ, Berman ML. Survival outcomes in patients with recurrent ovarian cancer who were treated with chemoresistance assay–guided chemotherapy. Am J Obstet Gynecol 2003;189:1301-1307.

eee. Fruehauf JP, Bosanquet AG. In vitro determination of drug response: a discussion of clinical applications. Principles and Practice of Oncology Updates 1993;7:1-16.

fff. Jain RK. Delivery of molecular medicine to solid tumors: lessons from in vivo imaging of gene expression and function. J Control Release 2001;74:7-25.

ggg. Tannock IF, Lee CM, Tunggal JK, Cowan DSM, Egorin MJ. Limited penetration of anticancer drugs through tumor tissue: A potential cause of resistance of solid tumors to chemotherapy. Clin Cancer Res 2002;8:878-884.

hhh. Zheng X, Wise SM, Cristini V. Nonlinear simulation of tumor necrosis, neo-vascularization, and tissue invasion via an adaptive finite-element/level-set method. Bull Math Biol 2005;67:211-259.

jjj. Plank MJ, Sleeman BD. Lattice and non-lattice models of tumour angiogenesis. Bull Math Biol 2004;66:1785-1819.

kkk. Cristini V, Lowengrub J, Nie Q. Nonlinear simulation of tumor growth. J Math Biol 2003;46:191-224.

lll. McDougall SR, Anderson ARA, Chaplain MAJ, Sherratt JA. Mathematical modelling of flow through vascular networks: Implications for tumour-induced angiogenesis and chemotherapy strategies. Bull Math Biol 2002;64:673-702.

mmm. Anderson ARA, Chaplain MAJ. Continuous and discrete mathematical models of tumor-induced angiogenesis. Bull Math Biol 1998;60:857-900.

nnn. Sinek J, Frieboes H, Zheng X, Cristini V. Two-dimensional chemotherapy simulations demonstrate fundamental transport and tumor response limitations involving nanoparticles. Biomed Microdevices 2004;6:297-309.

ooo. Cristini V, Frieboes HB, Gatenby R, Caserta S, Ferrari M, Sinek J. Morphologic instability and cancer invasion. Clin Cancer Res 2005;11:6772-6779.

ppp. Frieboes H, Zheng X, Sun C-H, Tromberg B, Gatenby R, Cristini V. An integrated computational/experimental model of tumor invasion. Cancer Res 2006;66:1597-1604.