PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Cancer Res. Author manuscript; available in PMC 2010 May 15.
Published in final edited form as:
PMCID: PMC2720602
NIHMSID: NIHMS109814

Prediction of drug response in breast cancer using integrative experimental/computational modeling

Abstract

Nearly 30% of women with early stage breast cancer develop recurrent disease attributed to resistance to systemic therapy. Prevailing models of chemotherapy failure describe three resistant phenotypes: cells with alterations in transmembrane drug transport, increased detoxification and repair pathways, and alterations leading to failure of apoptosis. Proliferative activity correlates with tumor sensitivity. Cell cycle status, controlling proliferation, depends upon local concentration of oxygen and nutrients. Although physiological resistance due to diffusion gradients of these substances and drug is a recognized phenomenon, it has been difficult to quantify its role with any accuracy that can be exploited clinically. We implement a mathematical model of tumor drug response that hypothesizes specific functional relationships linking tumor growth and regression to the underlying phenotype. The model incorporates the effects of local drug, oxygen and nutrient concentrations within the three-dimensional tumor volume, and includes the experimentally observed individual cells’ resistant phenotypes. By extracting mathematical model parameter values for drug and nutrient delivery from monolayer (one-dimensional) experiments and using the functional relationships to compute drug delivery in MCF-7 spheroid (three-dimensional) experiments, we use the model to quantify the diffusion barrier effect, which alone can result in poor response to chemotherapy both from diminished drug delivery and from lack of nutrients required to maintain proliferative conditions. We conclude that this integrative methodology tightly coupling computational modeling with biological data enhances the value of knowledge gained from current pharmacokinetic measurements, and, further, that such an approach could predict resistance based on specific tumor properties and thus improve treatment outcome.

Keywords: breast cancer, drug response, mathematical model

INTRODUCTION

We implement a novel quantitative approach that links growth and regression of a tumor mass to the underlying phenotype to study the impact of drug and nutrient delivery as mediators of physiological resistance. It is well known that inefficient vascularization may prevent optimal transport of oxygen, nutrients, and therapeutics to cancer cells in solid tumors [1]. As a result, the drug agent must diffuse through the tumor volume to reach the entire tumor cell population, and there is mounting evidence to indicate that these diffusion gradients may significantly limit drug access [25]. Both hypoxia and hypoglycemia contribute to physiological resistance through various mechanisms, including induction of oxidative stress and a decrease in the number of proliferating cells [612]. The myriad of stresses can lead to selection of cells that resist apoptotic conditions, thus adding to pathologic resistance in tumors [13]. This evidence strongly suggests that the diffusion process alone can lead to the evolution of drug resistance in tumor cells that exceeds predictions based on individual cell phenotype [5]. It has not been easy to quantify the resistance effects of diffusion gradients with any accuracy that can be exploited in a clinical setting. The different physical scales in a tumor spanning from the nano- to the centimeter scale present a complex system that to be better understood could benefit from appropriate mathematical models and computer simulations in addition to laboratory and clinical observations.

In particular, biocomputational modeling of tumor drug response has endeavored in the last two decades to address this need. Space limitations preclude a full description (see [5,14,15], references therein). Dox cellular pharmacodynamics has been modeled; e.g., [16] presented a model providing good fits to in-vitro cytotoxicity data. Drug transport was modeled in spheroids vs. monolayers [17]. A model capable of predicting intracellular Dox accumulation that matched experimental observations was described in [18]. Different drug kinetics effects in-vitro were compared in [19], showing that a single drug infusion could be more effective than repeated short applications. Models employing multi-scale approaches, i.e., linking events at sub-cellular, cellular, and tumor scales (e.g., [20]), studying vascularized tumor treatment (e.g., [21]), and simulating nanoparticle effects (e.g., [3]) have also been developed. Existing mathematical models are often limited to radially-symmetric tumor representations and not fully constrained through experimentally-set parameters.

Here, we employ a multi-scale computational model, extending a previous formulation of tumor growth founded in cancer biology [2225], to enable more rigorous quantification of diffusion effects on tumor drug response. This model can represent non-symmetric solid tumor morphologies in 3-D, thus providing the capability to capture the physical complexity and heterogeneity of the cancer microenvironment. More significantly, we fully constrain the model through functional relationships with parameters set from experiments. We hypothesize the simplest relationships that would at the same time be biologically founded and which could be calibrated by the experimental data. These relationships link tumor mass growth and regression to the underlying phenotype. They provide the mathematical basis for describing cell mitosis, apoptosis and necrosis modulated by diffusion gradients of oxygen, nutrients, and cytotoxic drugs, and enable quantification of the physiological resistance introduced by these gradients. Input parameters include the diffusion coefficients for these substances, and the rate constants for proliferation, apoptosis, and necrosis. We measure parameter values from independent experiments performed under conditions with no gradients, i.e., with cells grown as 1-D monolayers, and then use these values to calculate cell survival in 3-D tumor geometry. These values are compared with experiments in which cells are grown as in-vitro tumor spheroids, representing a 3-D tumor environment with diffusion gradients. This approach allows us 1) to fully constrain the computational model, using experimentally-obtained parameter values, and 2) to validate the hypothesized functional relationships by comparing the computed 3-D tumor viability with the spheroid tumor growth experiments. By quantifying the link between tumor growth and regression and the underlying phenotype, the work presented here provides a quantitative tool to study tumor drug response and treatment.

Although in-vitro spheroids are a gross simplification of the complex in-vivo condition, they allow capturing the effect of 3-D tissue architecture in generating diffusion gradients of drug and cell substrates under controlled conditions without the complicating effects of blood flow [1]. Spheroids develop a layer of viable cells surrounding a necrotic core [26]. The thickness of this layer (ca. 100 µm), maintained by substrate diffusion gradients, “mimics” the viable tissue thickness supported by diffusion from surrounding blood vessels in-vivo [26]. This provides a more controllable experimental model than in-vivo, which is necessary in order to test and calibrate the computational model parameters. In turn, computer simulations of the model enable formulating and testing the functional relationships that quantify the dependence of drug resistance on diffusion gradients and particular tissue characteristics (e.g., tissue compactness as a function of cell packing density [27]). We conclude that computational modeling tightly integrated with tumor biology extends the information that can be learned from pharmacokinetic experiments [2830], offering a promising possibility of ultimately quantifying and predicting treatment response from individual tumor characteristics.

MATERIALS AND METHODS

Doxorubicin was used as the cytotoxic drug in the experiments. Experiments used MCF-7 WT breast cancer cell lines to represent drug-sensitive tumor cells and MCF-7 40F breast cancer cell lines to represent drug-resistant tumor cells. Details of experimental methods are given below.

Observation of cell substrate gradients

Confluent MCF-7 WT (breast cancer) cell monolayers were incubated with 0.25% (w/v) trypsin (Sigma, UK) and 0.02% (w/v) EDTA (Sigma, UK) for 4 minutes at 37°C. Complete media was added to this mixture to inhibit enzyme activity and the cell suspension passed through a 24-gauge needle 6 times to ensure a single cell population. A 10 mL High Aspect Ratio Vessel (HARV; Synthecon, USA) was seeded at 2×106 cells/mL and rotated at 15 rpm at 37°C in a humidified atmosphere of 5% CO2 in air. Aggregates formed within 6 hours and were maintained in RCCS culture for a total of 30 days. During this time, media in the HARVs was replenished 50:50 every other day. Spheroids were fixed in formalin, embedded in paraffin and 5 µm sections cut for immunohistochemistry. To observe hypoxia, unfixed spheroids were incubated in 200 µM pimonidazole (Chemicon, UK) in complete media for 2 hours at 37°C prior to fixation as above. Proteinase K digestion was used for antigen retrieval prior to 1° Ab incubation. The mouse IgG 1° Ab was used at 1:600 dilution and incubated with 4 µm sections for immunolocalization of hypoxia through detection of pimonidazole protein adducts. Immunohistochemical (IHC) staining for GLUT-1 was performed using rabbit polyclonal antibody against the C-terminal portion (Abcam). The NHE-1 antibody is a rabbit polyclonal antibody (Santa Cruz). Immunohistochemistry was performed using the Discovery XT Automated Staining platform (Ventana Medical Systems, Inc, Tucson, Arizona). Deparaffinization and antigen retrieval of tissue sections was performed online. Antigen retrieval was performed on the Ventana instrument with a borate buffer (pH 8) 40 minutes at 95°C. Once the tissue was conditioned in this way, primary antibody was applied manually at a 1:800 dilution. (60 min. incubation at 37°C). Primary antibodies were visualized using VMSI validated detection and counterstaining reagents. Images were captured using an Olympus BX50 camera with an RT SPOT (Diagnostic Instruments, Inc), and standardized for light intensity.

Drug response

MCF-7 WT (drug-sensitive) and MCF-7 40F (drug-resistant) cell lines were cultured in RPMI Media 1640 (Life Technologies Invitrogen, Carlsbad, California) supplemented with 3% FBS (Life Technologies Invitrogen, Carlsbad, California), 2 mM L-glutamine, and 1% penicillin/streptomycin in humidified 7.5% CO2 at 37°C. Cells for monolayer culture were seeded 20,000 per well in 24-well Costar 3527 cell culture plates (Corning, New York). Cells for spheroids were seeded 50,000 per well in 24-well Costar 3473 cytophobic plates, and shaken at 100 rpm for 10 min. on day 1. Both the cells in monolayer and spheroids were incubated for 3 days, and then exposed to doxorubicin (Dox) (Bristol-Myers Squibb, Princeton, New Jersey) concentrations ranging from 0 to 16384 nM in 4× nM increments (0, 4, 16, 64, etc.) for 96 hours, representing at 256 nM a typical area-under-the-curve1 in-vivo [29]. Negative controls were seeded and incubated under the same conditions without drug. Three endpoints were concurrently measured as fraction of negative control: proliferation using tritiated thymidine (Amersham, Buckinghamshire, Great Britain) incorporation assay [31], viability using trypan blue exclusion counts, and metabolic activity using XTT (sodium 3'-[1-[(phenylamino)-carbonyl]-3,4-tetrazolium]-bis(4-methoxy-6-nitro)benzene-sulfonic acid hydrate) assay [32]. All experiments were done in triplicate. Photographs were taken with a digital camera through a Zeiss microscope (100× magnification).

Mathematical model of tumor growth and drug response

Briefly, the equations [22,23] are formulations of mathematical models used in engineering to describe phase separation of two partially miscible components (viable and necrotic tumor tissue), diffusion of small molecules (cell substrates and drug), and conservation of mass (Quick-Guide, Supplementary Figure 1, see also [3336]). Mass conservation equations describe growth (proliferation as a function of total number of cycling cells) and death from the drug’s cytotoxic effects (apoptosis as a function of a rate constant dependent upon the unique cell sensitivity and as a function of cycling cells). These are combined with diffusion of small molecules to a reaction-diffusion equation. Rate constants for proliferation and apoptosis are modified by functions that represent their dependence upon cell nutrients and oxygen (proliferation) and drug concentration (death), along with a dependence on spatial diffusion of these substances. By dimensional analysis of the reaction-diffusion equation the diffusion penetration length L = (D/η)1/2, where D is the corresponding diffusion constant and η is a characteristic cellular uptake rate. Large differences in concentration will occur if L is much less than the tumor radius (linear dimension) in the absence of blood vessels, creating a steep gradient. By using L and knowing the spheroid size, one can determine the steepness of these gradients by solving the reaction-diffusion equation. We note that gradients are affected by cell packing density, intra-cellular uptake, and duration of drug exposure, as well as conditions specific to particular experiments.

Estimation of mathematical model parameters

In our model, cell substrates are represented by glucose and oxygen with L estimated from independent measurements taken from the literature. For glucose D ~ 1×10−7 cm2/sec [37] and uptake rate may be as large as 1×10−3 sec−1 [38], while for oxygen the corresponding values are ~ 1×10−5 cm2/sec [39] and 1×10−1 sec−1 [40]. Diffusion penetration lengths thus calculated (L ≈ 100 µm) are consistent with our experimental observations (IHC) and with previous (murine mammary tumor EMT6/Ro spheroid) observations showing glucose concentration decreasing by 65% [41] and oxygen decreasing by 90% [42] across the tumor viable region.

Although gradients of drug into tumor tissue have been observed with a number of different cell types, many of these observations have been qualitative. Thus, published values for diffusion constants and uptake rates can vary considerably. We estimate Dox diffusion length by considering that the distance at which the penetrating concentration will be 50% of the source is ln 2 times the diffusion length – equivalent L here is ca. 90 microns for a 2/3 drop across the spheroid viable region (~ 100 µm). This is consistent with in-vitro Dox penetration into Chinese hamster lung cell spheroids at 24 hours [43] reporting an average 2/3 drop of external concentration across the viable region. Other studies have shown Dox penetration in-vivo (mouse model) decreasing by half within 40–50 µm of blood vessels by 3 hours [4], as well as penetration in humans (biopsies) decreasing by half within 50 µm of blood vessels after 2 hours of a bolus and within 60–80 µm after 24 hours in breast cancer tissue [2]. By fitting the solution of the unsteady Eq. (3) to this data we could indirectly estimate average cellular uptake rates η of Dox by breast cancer cells in 3-D tissue of ~ 1 day−1 , consistent with the independent penetration length estimate L = (D/η)1/2 from the steady-state profiles.

The rate λM (inverse time) measures the change in cell number in a population due to mitosis, normalized by total cell number (control = cells with no drug), and thus was calibrated by matching proliferation data from our experiments, using as initial guess for λMM,C, the in-vitro cell proliferation as fraction of control divided by in-vitro cell viability (total number of cells N in monolayer) as fraction of control, yielding an estimate of cell proliferation ~1 day−1.

Apoptosis rate λA (inverse time) over the period 0 < t < T (total experiment time) was estimated as a function of local concentrations of substrate σ and drug d as

λA(σ,d)=λA'(σ,d)+λA"(σ),    where  λA'=T1log(N(d)/NC)

and λA"=T1log(1+α(1σ/σ)). Death rate λA’ was calculated from measurements of cell viability N as fraction of negative control NC over time in monolayer cell cultures at various drug concentrations. Physiological resistance based on cell cycle status was modelled with λA". Decreased substrate concentration in 3-D results in cell populations that cycle less and therefore have reduced sensitivity to Dox2. In the above formula, enhanced survival was assumed linearly dependent on substrate concentration with a fitting constant calibrated from previous observations [12] showing an average MCF-7 viability increase 2.2 to 4 times for glucose deprivation of up to 50% at drug concentrations similar to our study, i.e., 2.2≤1+α/2≤4.0. In our experiments, glucose concentration in the medium was σ = 2.0 g/L. Note that for σ = σ, the above formula gives λA"=0 and thus λA"=λA'.

Necrosis parameters for diffusion-limited growth λN (rate of necrosis) and σN (substrate limit for cell viability) were calibrated by matching simulation results to in-vitro spheroid growth curves and histological data under conditions with no drug, as described previously [45], i.e., so that the simulated spheroid radius and viable region thickness matched in-vitro observations. These parameters are directly responsible for the steady spheroid size (and extent of necrosis) after a growth period because they regulate the balance of mass growth due to cell proliferation in the viable region and mass loss due to cell disintegration in the necrotic center [45]. Based on our in-vitro measurements of spheroid and necrotic volumes we found a stationary average spheroid radius of ca. 0.8 mm, with viable region of thickness ~0.1 mm (see also [46]). The latter is consistent with the estimated cell substrate penetration length L calculated above. The model calibration based on these quantities consistently yielded:λNM = 0.7 and σN = 0.5.

Higher cell densities were observed for resistant spheroids, which had a more compact morphology [27]. Based on these observations, higher cell adhesion parameter values [22] were used in the model equations to simulate resistant spheroids than for sensitive ones. We used a recently presented calibration procedure [45] to determine the relative values for sensitive and resistant cell cultures. This observation and the resulting different parameters for the two cell phenotypes underscore the need to explicitly model cell adhesion forces when constructing models of cancer growth.

QUICK-GUIDE TO EQUATIONS AND ASSUMPTIONS

ϕV/t+·(uVϕV)=M·(ϕV(δE/δϕV))+SV,ϕD/t+·(uDϕD)=M·(ϕD(δE/δϕD))+SD.
Equation (1)

The 3D in-vitro environment is modeled as a mixture of viable tumor tissue (volume fraction ϕV), dead tumor tissue (ϕD), and interstitial fluid (ϕW, given indirectly by 1−(ϕVD) and assumed to move freely) flowing through the extracellular matrix, which we treat as a porous medium. These partial differential equations are derived from the conservation of mass of these quantities [22,23]. From left to right in each equation, the terms represent: change of volume fraction with respect to time; tissue advection (bulk transport) by the overall mixture as it moves with local velocity u; net tissue diffusion due to the balance of mechanical forces, including cell-cell adhesion (related to mobility constant M [22,23]), and cell-cell adhesion and repulsion (oncotic pressure), modeled using the variational derivative of the adhesion energy E (for specific forms of this energy, see [22]); and net source of tissue due to the balance of mitosis, apoptosis, and necrosis.

In Words

The temporal rate of change in viable and dead tumor tissue at any location within the tumor equals the amount of mass that is pushed, transported, and pulled due to cell motion, adhesion, and tissue pressure, plus the net result of production and destruction of mass due to cell proliferation and death.

Major Assumptions of the Model

The tumor is a mixture of cells, interstitial fluid, and extracellular matrix (ECM), and we treat the motion of interstitial fluid and cells through the ECM as fluid flow in a porous medium. Therefore, no distinction between interstitial fluid hydrostatic pressure and mechanical pressure due to cell-cell interactions is made. Cell velocity is a function of cell mobility and tissue oncotic (solid) pressure (Darcy’s law); cell-cell adhesion is modeled using an energy approach from continuum thermodynamics [22].

SV=(λM(σ,d)λA(σ,d)λN[σNσ])ϕV,SD=(λA(σ,d)+λN[σNσ])ϕVλLϕD.
Equation (2)

These equations specify net sources SV and SD of mass for viable and dead tumor tissue, respectively. This directly links the conservation of mass equations Eq. (1) to the cell phenotype through hypothesized phenomenological functional relationships that include diffusion of cell substrates and drug (of local concentration σ̣ and d, respectively) through tumor interstitium (see [22] and references therein). For the viable tissue (first equation), the terms on the right hand side represent, from left to right, volume fraction gained from cell mitosis (rate λM), and lost to cell apoptosis (λA) and necrosis (λN). For the dead tissue (second equation), the terms on the right hand side represent, from left to right, increase in volume fraction from cell apoptosis and necrosis, and decrease in volume fraction from cell disintegration by lysis (λL). All rates are inverse time.

In Words

Mass of viable tumor tissue increases through cell proliferation and decreases through cell apoptosis and necrosis. Mass of dead tumor tissue increases through cell apoptosis and necrosis, and decreases through cell lysing. These equations provide a means of incorporating the biology of the problem into the physical conservation laws in Eq. (1).

Major Assumptions of the Model

Cells are composed entirely of water, which is a reasonable first approximation in terms of volume fraction (see discussion in [22]). Cell mitosis is proportional to substrates present [47]. As mitosis occurs, an appropriate amount of water from the interstitial fluid is converted into cell mass. Substrate depletion below a level σN leads to necrosis [26,7]. Cell lysis represents a loss of mass as cellular membranes are degraded and the mass is converted into water that is absorbed into the interstitial fluid. The simulated interface between viable and necrotic tissue is forced to be of biologically realistic thickness (10–100 µm) through the Heaviside function H (see [22,23]. Mitosis and apoptosis rates are functions of drug (d) and substrate (σ) concentrations. Note that the assumption that parameters measured for drug sensitivity in monolayer, modulated by diffusion gradients, can adequately represent response in-vivo may not be correct for real tumors; e.g., it disregards drug resistance from growing with an ECM.

η1σ/t=L22σσ
Equation (3)

This is a partial differential equation that describes cell substrate concentration σ across the tumor viable region. The first term on the right hand side models diffusion of the substrates into the tumor spheroid with penetration length L. The second term models substrate uptake by the cells (with rate η). An analogous equation is posed for drug concentration d.

In Words

The temporal rate of change of cell substrates or drug across the tumor viable region equals the net amount that diffuses into the region minus the amount uptaken by the tumor cells.

Major assumptions of the Model

Diffusion of cell substrates and drug into the tumor, combined with uptake in the tumor interior, creates and maintains a gradient of these substances through the 3-D tumor tissue, which is assumed to be fairly compact (non-invasive/infiltrative).

RESULTS

The hypothesized functional relationships of the computational model link growth and regression of the tumor mass to the underlying phenotype as described in Eq.1Eq.3 (Quick-Guide). The model is fully constrained (Methods) through experimentally-based parameters (Figure 1). In order to quantify physiological resistance in 3-D tumor tissue, diffusion gradients are incorporated in the functional relationships and the model is calibrated to compute these gradients. Accordingly, immunohistochemical measurements were performed on drug-sensitive tumor spheroids cultured in-vitro (Figure 2) to confirm that in our experiments 1) cell substrate concentration decreases across the radius of the spheroid, 2) apoptosis correlates with hypoxia, and 3) proliferation correlates with substrate availability [43,47]. Increasing positivity for pimonidazole protein adducts across the viable rim (~ 100 µm) in the direction of the spheroid center indicates a gradient of oxygen. Across the viable region, upregulation of Na+/H+ transporters detected by increasing positivity for NHE-1 is observed, demonstrating a gradient of pH. Na+/H+ transporters are upregulated in response to acidosis. Upregulation of GLUT-1 transporters (adaptation to hypoxic/hypoglycaemic stress) is shown in the increasing positivity for GLUT-1 [48] in peri-necrotic regions, indicating gradients of glucose and oxygen. The frequency of apoptotic cells increases towards the center as shown by increasing positivity of TUNEL3 staining. End stage apoptotic bodies are concentrated at the peri-necrotic necrotic regions with most cells in the center being apoptotic or necrotic. Cell proliferation is limited to the viable rim, as demonstrated by Ki-67 nuclear staining.

Figure 1
Computational model parameters
Figure 2
Diffusion gradients are incorporated in the model functional relationships to simulate the experimental conditions

According to the computational model, the solution of the reaction-diffusion Eq. (3) yields substrate gradients (Figure 3A) developing across a viable region of width ~ 0.1 mm. The computed width of this region is consistent with previous observations by other investigators (e.g., [26]). Both glucose and oxygen substrate concentrations were calculated in the model to drop by at least 50% across the spheroid viable region, with an average concentration left angle bracketσ/σright angle bracket ≈ 0.7 in the living tumor tissue. This is in agreement with measurements of glucose [41] and oxygen (e.g., [42]) in spheroid viable regions from previously reported in-vitro measurements (Methods).

Figure 3
Computational model computes the gradients in 3-D tissue

Similarly, the model computes the drug gradients that develop within the viable tumor spheroid tissue. The calculated drug concentration profile reaches steady state in ca. 24 hours and decays by 2/3 across the viable region to yield an average drug concentration of 50% that of outside the spheroid (Figure 3B). Note that the drug penetration length is comparable to the penetration length of substrates (Figure 2). Since Doxorubicin is a cell–cycle specific drug and the proportion of cells that are cycling correlates with the substrate availability, the drug would be less effective as the proportion of cycling cells is diminished. Not only is its concentration decreased by the diffusion barrier, but its efficacy is impaired as a result of the substrate gradients. The physiologic resistance predicted by the model based on the hypothesized functional relationships would in principle apply to any cell-cycle specific drug.

Cell proliferation (ratio of proliferation-to-viability counts) in-vitro under conditions with diffusion gradients (i.e., in spheroids) vs. monolayers was lower by at least 50% for both cell lines at drug concentrations d > 64 nM representing clinically relevant dosages4 (Figure 4A). Further, increased viability for drug-sensitive cells in 3-D corresponded with ca. 50% lower cell metabolic activity (ratio of metabolic-to-viability counts) in-vitro (data not shown). These data indicate that cell quiescence was significant at drug dosages similar to in-vivo conditions. The cell apoptosis parameter of the computational model calculated from the monolayer cell viability data increases with drug concentration (Figure 4B).

Figure 4
Cell proliferation and apoptosis

The functional relationships of the computational model (Eq. 1Eq.3) linking growth and regression of the tumor mass to the underlying phenotype are validated by quantifying the physiological resistance introduced by the diffusion gradients. Cell viabilities predicted by the model agree well with the ones observed experimentally in 3-D spheroids (Figure 5). In contrast, monolayer cell viability data is a poor prognosticator of drug resistance in 3-D tumors. An average survival increase (as fraction of control) of 250% over monolayer was predicted for drug-sensitive (MCF-7 WT) spheroids under various drug concentrations (Figure 5A), while for drug-resistant (MCF-7 40F) spheroids the corresponding increase was 280% (Figure 5B). The median dose5 for MCF-7 40F was ~20% higher in spheroid vs. monolayer than for MCF-WT, implying that not only did the 3-D morphology promote the net survival for both cell types, but also that the resistant phenotype was further favored by the 3-D configuration.

Figure 5
Validation of the hypothesized functional relationships

We noted that drug-resistant spheroids maintained more compact, nearly spherical shapes (Figure 6A) that simply decreased in size as cells were killed, while drug-sensitive spheroids formed irregular, looser shapes (Figure 6B) that fragmented, particularly at higher drug concentrations. The tumor morphology may depend on the competition between heterogeneous cell proliferation caused by diffusion gradients of cell substrates, driving shape instability, and stabilizing mechanical forces such as cell-cell and cell-matrix adhesion, as we quantitatively showed in previous work [45]. Here, we observe a similar instability for drug-sensitive cells, indicating that their adhesion is below the threshold necessary to maintain tumor shape stability. This is confirmed through the computational model, where variation in the parameter values for cell adhesion [22] reproduces these morphologies (Fig. 6C and D).

Figure 6
Correlation of tumor morphological stability with drug resistance

DISCUSSION

We have used an integrative methodology that tightly couples computational modeling with biological data to quantify physiological resistance in 3-D tumor tissue in-vitro. The model hypothesizes functional relationships that mechanistically link the growth and regression of a tumor mass to the underlying phenotype. These relationships incorporate the complex interplay between tumor growth, cell phenotype, and diffusion gradients, caused by heterogeneous delivery of oxygen and cell substrates and removal of metabolites, and are thus able to calculate tumor drug response as a predictable process dependent on biophysical laws. The results underscore the importance of a quantitative approach in evaluating chemotherapeutic agents that takes into consideration diffusion in addition to such chemo-protective effects as cell phenotypic properties and cell-cell and cell-ECM (extracellular matrix) adhesion [49]. Although the model was calibrated with in-vitro data for breast cancer, this approach may also apply to drug response of other tissues.

By modeling transport of drug and cell substrates, this work creates a quantitative tool that could in the future be used to predict resistance in patients based on their tumor cell properties, thus improving treatment outcome. In particular, the proposed functional relationships help to quantify resistance in human breast cancer due to local cell substrate depletion. Our integrative experimental/computational approach provides insight into the physical dynamics of solid tumors, and validates the hypothesized functional relationships as adequate to describe the observed phenomena. Although more complex functional forms could conceivably also provide a reasonable match between model and experimental results, quantifying the diffusion effect with a minimal mathematical description is the preferred approach, unless it is known to be wrong, as it facilitates insight into the system and provides for more economic simulations. We explicitly chose not to use complex relationships that would contain more parameters than the number of independent measurements available for calibration, thus resulting in an underdetermined problem.

Cell packing density affects the magnitude of diffusion gradients and may pose a barrier to effective drug penetration [27]. This density can vary between cell lines and tumors and reflect variations in drug resistance (cf. Figure 6). How this relates to cell adhesion molecule (CAM) expression (e.g., integrin, E-cadherin) is unclear. Mechanisms of cell packing may include stronger cell adhesion forces, due to higher E-cadherin expression; this should also limit proliferation, which does not seem to be the case here for the drug-resistant cells at lower drug concentrations. Additionally, cellular stress affects the quantity and strength of CAMs. We have previously investigated the effect of the tumor microenvironment on tissue morphology [50,45], suggesting that marginally stable environmental conditions could directly affect morphogenesis and present an additional challenge to therapy [48,45] in-vivo by increasing tumor cell invasiveness and leading to complex infiltrative morphologies, depending on the magnitude of cell adhesion forces that tend to maintain compact, non-invasive tumors [50]. Diffusion gradients combined with higher cell packing density augment drug resistance synergistically, as observed in our experiments with the higher median dose for the more compact, drug-resistant spheroids, yet it may be difficult to deterministically gauge their combined impact [5]. Synergism may be due to increased drug binding to ECM in tumor areas proximal to the drug source, while substrate and drug penetration to distal areas is significantly reduced due to higher cell packing, thus exacerbating the resistance effect of the diffusion barrier.

Note that the presumption that the necrotic region is not a risk factor for progression may not be necessarily true because selection pressures for resistance and induction of mutation are strongest in hypoxic, peri-necrotic areas. Thus, inducing a large necrotic fraction during treatment may paradoxically further select for drug resistance.

The class of prediction modeling presented in this paper, based on an assimilation of complex processes with a fully 3-dimensional physical environment, offers the capability of complementing current pharmacokinetic measurements. A more expansive integration of theoretical model parameters with biological data could help to move towards prediction of treatment response based on in-vitro and in-vivo tumor information, and determine the correspondence between in-vitro measurements and the in-vivo condition to refine and validate the assumption that this approach can describe in-vivo tumors. Incorporation of patient-specific tumor phenotypic and microenvironmental parameters into the model could enhance clinical strategies and prognosis evaluation. We further envision that these methods will enhance current pharmacokinetic models used in designing and interpreting Phase II clinical trials.

Supplementary Material

Supp. Figure 1

Suppl. Material Text

ACKNOWLEDGEMENTS

We acknowledge John Sinek and Steven Wise (Mathematics, UC-Irvine), Sandeep Sanga (Biomedical-Engineering, UT-Austin), Paul Macklin (SHIS, UT-Houston), Ernest Han (Obstetrics/Gynecology, UC-Irvine), Hoa Nguyen (Medicine, UC-Irvine) for helpful comments and discussions, and the reviewers for their valuable input. Grant support: The Cullen Trust of Health Care, NSF-DMS 0818104, National Cancer Institute, Department of Defense.

Footnotes

1In pharmacokinetics, the area-under-the-curve (AUC) is a calculation to evaluate the body's total exposure to a given drug. In a graph plotting plasma drug concentration vs. time, the area under this curve is the AUC

2Cytotoxicity of MCF-7 cells to Dox may be affected more by glucose deprivation than by hypoxia [6,44]. The glucose-regulated stress response [10] has been correlated with resistance to topoisomerase II-directed chemotherapeutic drugs such as Dox [10,11] through a decreased expression of topoisomerase II in MCF-7 cells [12].

3Terminal deoxynucleotidyl Transferase Biotin-dUTP Nick End Labeling.

4256 nM exposure in-vitro for 96 hours yields an area-under-the-curve (AUC) of ~2.5×10−5 Molar.Hr, while for a typical patient Dox plasma levels after a bolus injection represent an exposure of ~3×10−5 Molar.Hr.

5Median dose is the drug concentration required to achieve 50% cell kill.

REFERENCES

1. Trédan O, Galmarini CM, Patel K, Tannock IF. Drug resistance and the solid tumor microenvironment. J Nat Cancer Inst. 2007;99:1441–1454. [PubMed]
2. Lankelma J, Dekker H, Luque RF, et al. Doxorubicin gradients in human breast cancer. Clin Cancer Res. 1999;5:1703–1707. [PubMed]
3. Sinek J, Frieboes HB, Zheng X, Cristini V. Two-dimensional Chemotherapy Simulations Demonstrate Fundamental Transport and Tumor Response Limitations Involving Nanoparticles. Biomed Microdev. 2004;6:297–309. [PubMed]
4. Primeau AJ, Rendon A, Hedley D, Lilge L, Tannock IF. The distribution of the anticancer drug doxorubicin in relation to blood vessels in solid tumors. Clin Cancer Res. 2005;11:8782–8788. [PubMed]
5. Sinek JP, Sanga S, Zheng X, Frieboes HB, Ferrari M, Cristini V. Predicting drug pharmacokinetics and effect in vascularized tumors using computer simulation. J Math Biol. 2009;58:485–510. [PMC free article] [PubMed]
6. Greijer AE, de Jong MC, Scheffer GL, Shvarts A, van Diest PJ, van der Wall E. Hypoxia-induced acidification causes mitoxantrone resistance not mediated by drug transporters in human breast cancer cells. Cellular Oncology. 2005;27:43–49. [PubMed]
7. Spitz DR, Sim JE, Ridnour LA, Galoforo SS, Lee YJ. Glucose deprivation-induced oxidative stress in human tumor cells. Annals NY Acad Sci. 2000;899:349–362. [PubMed]
8. Lee YJ, Galoforo SS, Berns CM, et al. Glucose deprivation-induced cytotoxicity and alterations in mitogen-activated protein kinase activation are mediated by oxidative stress in multidrug-resistant human breast carcinoma cells. J Biol Chem. 1998;273:5294–5299. [PubMed]
9. Brown NS, Bicknell R. Hypoxia and oxidative stress in breast cancer. Oxidative stress: its effects on the growth, metastatic potential and response to therapy of breast cancer. Breast Can Res. 2001;3:323–327. [PMC free article] [PubMed]
10. Li J, Lee AS. Stress induction of GRP78/BiP and its role in cancer. Curr Mol Med. 2006;6:45–54. [PubMed]
11. Tomida A, Tsuruo T. Drug resistance mediated by cellular stress response to the microenvironment of solid tumors. Anti-Cancer Drug Design. 1999;14:169–177. [PubMed]
12. Yun J, Tomida A, Nagata K, Tsuruo T. Glucose-regulated stresses confer resistance to VP-16 in human cancer cells through a decreased expression of DNA topoisomerase II. Oncol Res. 1995;7:583–590. [PubMed]
13. Pusztai L, Hortobagyi GN. High-dose chemotherapy: how resistant is breast cancer? Drug Resist Updat. 1998;1:62–72. [PubMed]
14. Sanga S, Sinek JP, Frieboes HB, Ferrari M, Fruehauf JP, Cristini V. Mathematical modeling of cancer progression and response to chemotherapy. Expert Rev Anticancer Ther. 2006;6:1361–1376. [PubMed]
15. Sanga S, Frieboes HB, Zheng X, Bearer E, Cristini V. Predictive oncology: a review of multidisciplinary, multi-scale in-silico modeling linking phenotype, morphology and growth. Neuroimage. 2007;37:S120–S134. [PMC free article] [PubMed]
16. El-Kareh AW, Secomb TW. Two-mechanism peak concentration model for cellular pharmacodynamics of Doxorubicin. Neoplasia. 2005;7:705–713. [PMC free article] [PubMed]
17. Ward JP, King JR. Mathematical modeling of drug transport in tumour multicell spheroids and monolayer cultures. Math Biosci. 2003;181:177–207. [PubMed]
18. Jackson TL. Intracellular accumulation and mechanism of action of doxorubicin in a spatiotemporal tumor model. J Theor Biol. 2003;220:201–213. [PubMed]
19. Norris ES, King JR, Byrne HM. Modelling the response of spatially structured tumours to chemotherapy: Drug kinetics. Math Comp Model. 2006;43:820–837.
20. Byrne HM, Owen MR, Alarcón T, Murphy J, Maini PK. Modelling the response of vascular tumours to chemotherapy: a multiscale approach. Math Models Meth Appl Sci. 2005;16:1219–1241.
21. Panovska J, Byrne HM, Maini PK. A theoretical study of the response of vascular tumours to different types of chemotherapy. Math Comp Model. 2007;47:560–579.
22. Wise SM, Lowengrub JS, Frieboes HB, Cristini V. Nonlinear simulations of three-dimensional multispecies tumor growth–I. Model and numerical method. J Theor Biol. 2008;253:524–543. [PMC free article] [PubMed]
23. Frieboes HB, Lowengrub JS, Wise S, et al. Computer simulation of glioma growth and morphology. NeuroImage. 2007;37:S59–S70. [PMC free article] [PubMed]
24. Zheng X, Wise S, 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. [PubMed]
25. Cristini V, Lowengrub J, Nie Q. Nonlinear simulation of tumor growth. J Math Biol. 2003;46:191–224. [PubMed]
26. Sutherland RM. Cell and environment interactions in tumor microregions: the multicell spheroid model. Science. 1988;240:177–184. [PubMed]
27. Grantab R, Sivananthan S, Tannock IF. The penetration of anticancer drugs through tumor tissue as a function of cellular adhesion and packing density of tumor cells. Cancer Res. 2006;66:1033–1039. [PubMed]
28. Fruehauf JP, Brem H, Brem S, et al. In vitro drug response and molecular markers associated with drug resistance in malignant gliomas. Clin Cancer Res. 2006;12:4523–4532. [PubMed]
29. Fruehauf JP. In vitro assay-assisted treatment selection for women with breast or ovarian cancer. Endocrine-Related Cancer. 2002;9:171–182. [PubMed]
30. Mehta RS, Bornstein R, Yu IR, et al. Breast cancer survival and in vitro tumor response in the Extreme Drug Resistance Assay. Breast Cancer Res Treat. 2001;66:225–237. [PubMed]
31. Kern DH, Weisenthal LM. Highly specific prediction of antineoplastic drug resistance with an in vitro assay using suprapharmacologic drug doses. J Nat Cancer Inst. 1990;82:582–588. [PubMed]
32. Roehm NW, Rodgers GH, Hatfield SM, Glasebrook AL. An improved colorimetric assay for cell proliferation and viability utilizing the tetrazolium salt XTT. J Immunol Meth. 1991;142:257–265. [PubMed]
33. Kim J, Kang K, Lowengrub J. Conservative multigrid methods for ternary Cahn-Hilliard systems. Comm Math Sci. 2004;2:53–77.
34. Jiang GS, Shu CW. Effcient implementation of weighted ENO schemes. J Comput Phys. 1996;126:202–228.
35. Trottenberg U, Oosterlee C, Schüller A. Multigrid. New York: Academic Press; 2001.
36. Isaacson E, Keller H. Analysis of Numerical Methods. New York: Wiley; 1966.
37. Casciari JJ, Sotirchos SV, Sutherland RM. Glucose diffusivity in multicellular tumor spheroids. Cancer Res. 1988;48:3905–3909. [PubMed]
38. Kallinowski F, Vaupel F, Runkel S, et al. Glucose uptake, lactate release, ketone body turnover, metabolic micromilieu, and pH distributions in human breast cancer xenografts in nude rats. Cancer Res. 1988;48:7264–7272. [PubMed]
39. Nugent LJ, Jain RK. Extravascular diffusion in normal and neoplastic tissues. Cancer Res. 1984;44:238–244. [PubMed]
40. Casciari JJ, Sotirchos SV, Sutherland RM. Variations in tumor cell growth rates and metabolism with oxygen concentration, glucose concentration, and extracellular pH. J Cell Physiol. 1992;151:386–394. [PubMed]
41. Teutsch HF, Goellner A, Mueller-Klieser W. Glucose levels and succinate and lactate dehydrogenase activity in EMT6/Ro tumor spheroids. Eur J Cell Biol. 1995;66:302–307. [PubMed]
42. Acker H, Carlsson J, Mueller-Klieser W, Sutherland RW. Comparative pO2 measurements in cell spheroids cultured with different techniques. Br J Cancer. 1987;56:325–327. [PMC free article] [PubMed]
43. Durand RE. Slow penetration of anthracyclines into spheroids and tumors: a therapeutic advantage? Cancer Chemother Pharmacol. 1990;26:198–204. [PubMed]
44. Kalra R, Jones A-M, Kirk J, Adams GE, Stratford IJ. The effect of hypoxia on acquired drug resistance and response to epidermal growth factor in Chinese hamster lung fibroblasts and human breast-cancer cells in vitro. Int. J Cancer. 1993;54:650–655. [PubMed]
45. Frieboes HB, Zheng X, Sun C-H, Tromberg B, Gatenby R, Cristini V. An integrated computational/experimental model of tumor invasion. Cancer Res. 2006;66:1567–1604. [PubMed]
46. Freyer JP, Sutherland RM. Regulation of growth saturation and development of necrosis in EMT6/Ro multicellular spheroids by the glucose and oxygen supply. Cancer Res. 1986;46:3504–3512. [PubMed]
47. Carlsson J. A proliferation gradient in three-dimensional colonies of cultured human glioma cells. Int J Cancer. 1977;20:129–136. [PubMed]
48. Gatenby RA, Smallbone K, Maini PK, et al. Cellular adaptations to hypoxia and acidosis during somatic evolution of breast cancer. Br J Cancer. 2007;97:646–653. [PMC free article] [PubMed]
49. Morin PJ. Drug resistance and the microenvironment: nature and nurture. Drug Resist Updates. 2003;6:169–172. [PubMed]
50. Cristini V, Frieboes HB, Gatenby R, Caserta S, Ferrari M, Sinek J. Morphologic instability and cancer invasion. Clin Cancer Res. 2005;11:6772–6779. [PubMed]