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 on local concentration of oxygen and nutrients. Although physiologic resistance due to diffusion gradients of these substances and drugs 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 resistant phenotypes of individual cells. We conclude that this integrative method, 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. [Cancer Res 2009;69(10):4484–92]
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.
The three-dimensional 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 − (ϕV + ϕD) and assumed to move freely] flowing through the ECM, 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; refs. 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 ref. 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.
The tumor is a mixture of cells, interstitial fluid, and 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).
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 ref. 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.
Cells are composed entirely of water, which is a reasonable first approximation in terms of volume fraction (see discussion in ref. 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 refs. 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).
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 taken up by the tumor cells.
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 three-dimensional tumor tissue, which is assumed to be fairly compact (noninvasive/infiltrative).
Introduction
We implement a novel quantitative approach that links growth and regression of a tumor mass to the underlying phenotype to study the effect of drug and nutrient delivery as mediators of physiologic 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 (2–5). Both hypoxia and hypoglycemia contribute to physiologic resistance through various mechanisms, including induction of oxidative stress and a decrease in the number of proliferating cells (6–12). 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 nanometer 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 refs. 5, 14, 15, references therein). Doxorubicin cellular pharmacodynamics has been modeled; for example, ref. 16 presented a model providing good fits to in vitro cytotoxicity data. Drug transport was modeled in spheroids versus monolayers (17). A model capable of predicting intracellular doxorubicin accumulation that matched experimental observations was described in ref. 18. Different drug kinetic effects in vitro were compared in ref. 19, showing that a single drug infusion could be more effective than repeated short applications. Models using multiscale approaches [i.e., linking events at subcellular, cellular, and tumor scales (e.g., ref. 20); studying vascularized tumor treatment (e.g., ref. 21); and simulating nanoparticle effects (e.g., ref. 3)] have also been developed. Existing mathematical models are often limited to radially symmetrical tumor representations and not fully constrained through experimentally set parameters.
Here, we use a multiscale computational model, extending a previous formulation of tumor growth founded in cancer biology (22–25), to enable more rigorous quantification of diffusion effects on tumor drug response. This model can represent nonsymmetrical solid tumor morphologies three-dimensionally, 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 physiologic 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 done under conditions with no gradients (i.e., with cells grown as one-dimensional monolayers) and then use these values to calculate cell survival in three-dimensional tumor geometry. These values are compared with experiments in which cells are grown as in vitro tumor spheroids, representing a three-dimensional tumor environment with diffusion gradients. This approach allows us (a) to fully constrain the computational model, using experimentally obtained parameter values, and (b) to validate the hypothesized functional relationships by comparing the computed three-dimensional 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 three-dimensional 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 (∼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 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; ref. 27). We conclude that computational modeling tightly integrated with tumor biology extends the information that can be learned from pharmacokinetic experiments (28–30), 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 wild-type (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) and 0.02% (w/v) EDTA (Sigma) for 4 min at 37°C. Complete medium was added to this mixture to inhibit enzyme activity, and the cell suspension was passed through a 24-gauge needle six times to ensure a single-cell population. A 10-mL High Aspect Ratio Vessel (HARV; Synthecon) 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 h and were maintained in RCCS culture for a total of 30 d. During this time, medium 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 μmol/L pimonidazole (Chemicon) in complete medium for 2 h at 37°C before fixation as above. Proteinase K digestion was used for antigen retrieval before primary antibody incubation. The mouse IgG primary antibody was used at 1:600 dilution and incubated with 4-μm sections for immunolocalization of hypoxia through detection of pimonidazole protein adducts. Immunohistochemical staining for GLUT-1 was done with rabbit polyclonal antibody against the COOH-terminal portion (Abcam). The NHE-1 antibody is a rabbit polyclonal antibody (Santa Cruz). Immunohistochemistry was done with the Discovery XT Automated Staining platform (Ventana Medical Systems, Inc.). Deparaffinization and antigen retrieval of tissue sections were done online. Antigen retrieval was done on the Ventana instrument with a borate buffer (pH 8) for 40 min at 95°C. Once the tissue was conditioned in this way, primary antibody was applied manually at 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 medium 1640 (Life Technologies Invitrogen) supplemented with 3% fetal bovine serum (Life Technologies Invitrogen), 2 mmol/L l-glutamine, and 1% penicillin/streptomycin in humidified 7.5% CO2 at 37°C. Cells for monolayer culture were seeded at 20,000 per well in 24-well Costar 3527 cell culture plates (Corning). Cells for spheroids were seeded at 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 d and then exposed to doxorubicin (Bristol-Myers Squibb) concentrations ranging from 0 to 16,384 nmol/L in 4× nanomolar increments (0, 4, 16, 64, etc.) for 96 h, representing at 256 nmol/L a typical area under the curve (AUC)12
In pharmacokinetics, the AUC is a calculation to evaluate the total exposure of the body to a given drug. In a graph plotting plasma drug concentration versus time, the area under this curve is the AUC.
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 Fig. S1; refs. 33–36). Mass conservation equations describe growth (proliferation as a function of total number of cycling cells) and death from the cytotoxic effects of the drug (apoptosis as a function of a rate constant dependent on 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 on 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 smaller 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, intracellular 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/s (37) and uptake rate may be as large as 1 × 10−3 s−1 (38), whereas for oxygen, the corresponding values are ∼1 × 10−5 cm2/s (39) and 1 × 10−1 s−1 (40). Diffusion penetration lengths thus calculated (L ≈ 100 μm) are consistent with our experimental observations (immunohistochemistry) 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 doxorubicin 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 ∼90 μm for a 2/3 drop across the spheroid viable region (∼100 μm). This is consistent with in vitro doxorubicin penetration into Chinese hamster lung cell spheroids at 24 h (43) reporting an average 2/3 drop of external concentration across the viable region. Other studies have shown doxorubicin penetration in vivo (mouse model) decreasing by half within 40 to 50 μm of blood vessels by 3 h (4), as well as penetration in humans (biopsies) decreasing by half within 50 μm of blood vessels after 2 h of a bolus and within 60 to 80 μm after 24 h in breast cancer tissue (2). By fitting the solution of the unsteady Eq. 3 to these data, we could indirectly estimate average cellular uptake rates η of doxorubicin by breast cancer cells in three-dimensional 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 λM/λM,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
Cytotoxicity of MCF-7 cells to doxorubicin 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 doxorubicin (10, 11) through a decreased expression of topoisomerase II in MCF-7 cells (12).
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 histologic data under conditions with no drug, as described previously (45), 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 ∼0.8 mm, with viable region of thickness ∼0.1 mm (see also ref. 46). The latter is consistent with the estimated cell substrate penetration length L calculated above. The model calibration based on these quantities consistently yielded λN/λN = 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.
Results
The hypothesized functional relationships of the computational model link growth and regression of the tumor mass to the underlying phenotype as described in Eqs. 1–3 (Quick Guide). The model is fully constrained (Materials and Methods) through experimentally based parameters (Fig. 1). To quantify physiologic resistance in three-dimensional tumor tissue, diffusion gradients are incorporated in the functional relationships and the model is calibrated to compute these gradients. Accordingly, immunohistochemical measurements were done on drug-sensitive tumor spheroids cultured in vitro (Fig. 2) to confirm that in our experiments (a) cell substrate concentration decreases across the radius of the spheroid, (b) apoptosis correlates with hypoxia, and (c) 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, up-regulation of Na+/H+ transporters detected by increasing positivity for NHE-1 is observed, showing a gradient of pH. Na+/H+ transporters are up-regulated in response to acidosis. Up-regulation of GLUT-1 transporters (adaptation to hypoxic/hypoglycemic stress) is shown in the increasing positivity for GLUT-1 (48) in perinecrotic regions, indicating gradients of glucose and oxygen. The frequency of apoptotic cells increases toward the center as shown by increasing positivity of terminal deoxyribonucleotidyl transferase–mediated dUTP nick end labeling (TUNEL) staining. End-stage apoptotic bodies are concentrated at the perinecrotic necrotic regions, with most cells in the center being apoptotic or necrotic. Cell proliferation is limited to the viable rim, as shown by Ki-67 nuclear staining.
Computational model parameters. Functional relationships linking tumor mass growth and regression to the underlying phenotype are based on these parameters, with values derived from experimental observations (obtained previously and in this study). (*) Value shown is an example, corresponding to drug-sensitive MCF-7 WT cells at 256 nmol/L doxorubicin (Dox).
Computational model parameters. Functional relationships linking tumor mass growth and regression to the underlying phenotype are based on these parameters, with values derived from experimental observations (obtained previously and in this study). (*) Value shown is an example, corresponding to drug-sensitive MCF-7 WT cells at 256 nmol/L doxorubicin (Dox).
Diffusion gradients are incorporated in the model functional relationships to simulate the experimental conditions. Clockwise from top left: Immunohistochemical staining of tumor spheroids showing increasing positivity for pimonidazole protein adducts (darker brown) in the direction of the center, indicating a decrease in oxygen across the viable rim of ∼100 μm (v.r., viable rim; n, necrotic region); NHE-1 staining showing up-regulation of Na+/H+ transporters (darker brown) toward the necrotic region in response to a decrease in pH; increasing positivity for GLUT-1 staining (brown) in the perinecrotic region showing cellular adaptation to hypoglycemic and hypoxic stress due to decreasing concentration of glucose and oxygen across the spheroid radius; increasing positivity of TUNEL staining (dark brown) toward the center indicating correlation of apoptosis with hypoxia; Ki-67 nuclear staining (dark gray) showing cell proliferation correlating with substrate availability. Bar, 50 μm.
Diffusion gradients are incorporated in the model functional relationships to simulate the experimental conditions. Clockwise from top left: Immunohistochemical staining of tumor spheroids showing increasing positivity for pimonidazole protein adducts (darker brown) in the direction of the center, indicating a decrease in oxygen across the viable rim of ∼100 μm (v.r., viable rim; n, necrotic region); NHE-1 staining showing up-regulation of Na+/H+ transporters (darker brown) toward the necrotic region in response to a decrease in pH; increasing positivity for GLUT-1 staining (brown) in the perinecrotic region showing cellular adaptation to hypoglycemic and hypoxic stress due to decreasing concentration of glucose and oxygen across the spheroid radius; increasing positivity of TUNEL staining (dark brown) toward the center indicating correlation of apoptosis with hypoxia; Ki-67 nuclear staining (dark gray) showing cell proliferation correlating with substrate availability. Bar, 50 μm.
According to the computational model, the solution of the reaction-diffusion, Eq. 3, yields substrate gradients (Fig. 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., ref. 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 〈σ/σ∞〉 ≈ 0.7 in the living tumor tissue. This is in agreement with measurements of glucose (41) and oxygen (e.g., ref. 42) in spheroid viable regions from previously reported in vitro measurements (Materials and Methods).
Computational model computes the gradients in three-dimensional (3-D) tissue. Cross section of a tumor spheroid in the computational model (A; dotted line, tumor boundary; dashed line, perinecrotic area) showing calculated diffusion gradients (arrow) of cell substrate concentration σ/σ∞ (as fraction of external concentration in the medium) using diffusion penetration length parameter L ≈ 100 μm. Necrosis becomes significant at substrate level σ/σ∞ < 0.5. B, doxorubicin concentration d/d∞ in tumor tissue (normalized with concentration in the medium) versus distance from the tumor/medium interface, estimated for the model to be consistent with data in refs. 2, 43 and others (Materials and Methods). By solving the unsteady Eq. 3, the drug gradient across the viable region can be calculated at various times: bottom curve, 2 h; middle curve, 24 h; topmost curve, steady state, showing ∼2/3 drop from the tumor edge to the perinecrotic region.
Computational model computes the gradients in three-dimensional (3-D) tissue. Cross section of a tumor spheroid in the computational model (A; dotted line, tumor boundary; dashed line, perinecrotic area) showing calculated diffusion gradients (arrow) of cell substrate concentration σ/σ∞ (as fraction of external concentration in the medium) using diffusion penetration length parameter L ≈ 100 μm. Necrosis becomes significant at substrate level σ/σ∞ < 0.5. B, doxorubicin concentration d/d∞ in tumor tissue (normalized with concentration in the medium) versus distance from the tumor/medium interface, estimated for the model to be consistent with data in refs. 2, 43 and others (Materials and Methods). By solving the unsteady Eq. 3, the drug gradient across the viable region can be calculated at various times: bottom curve, 2 h; middle curve, 24 h; topmost curve, steady state, showing ∼2/3 drop from the tumor edge to the perinecrotic region.
Similarly, the model computes the drug gradients that develop within the viable tumor spheroid tissue. The calculated drug concentration profile reaches steady state in ∼24 hours and decays by 2/3 across the viable region to yield an average drug concentration of 50% that of outside the spheroid (Fig. 3B). Note that the drug penetration length is comparable to the penetration length of substrates (Fig. 2). Because 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 also 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) versus monolayers was lower by at least 50% for both cell lines at drug concentrations d > 64 nmol/L representing clinically relevant dosages14
In vitro exposure to 256 nmol/L for 96 h yields an AUC of ∼2.5 × 10−5 mol/L h, whereas for a typical patient doxorubicin plasma levels after a bolus injection represent an exposure of ∼3 × 10−5mol/L h.
Cell proliferation and apoptosis. A, cell proliferation (ratio of proliferation-to-viability counts; n = 3) in vitro was lower by at least 50% under conditions with diffusion gradients (spheroids) versus monolayers for both cell lines at drug concentrations d > 64 nmol/L representing clinically relevant dosages. B, drug-induced cell apoptosis death rate parameter λA′(σ∞, d) (inverse time) versus doxorubicin concentrations d calculated (Materials and Methods) from measured in vitro monolayer viability counts (n = 3); glucose concentration σ∞ = 2.0 g/L.
Cell proliferation and apoptosis. A, cell proliferation (ratio of proliferation-to-viability counts; n = 3) in vitro was lower by at least 50% under conditions with diffusion gradients (spheroids) versus monolayers for both cell lines at drug concentrations d > 64 nmol/L representing clinically relevant dosages. B, drug-induced cell apoptosis death rate parameter λA′(σ∞, d) (inverse time) versus doxorubicin concentrations d calculated (Materials and Methods) from measured in vitro monolayer viability counts (n = 3); glucose concentration σ∞ = 2.0 g/L.
The functional relationships of the computational model (Eqs. 1–3) linking growth and regression of the tumor mass to the underlying phenotype are validated by quantifying the physiologic resistance introduced by the diffusion gradients. Cell viabilities predicted by the model agree well with the ones observed experimentally in three-dimensional spheroids (Fig. 5). In contrast, monolayer cell viability data are a poor prognosticator of drug resistance in three-dimensional 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 (Fig. 5A), whereas for drug-resistant (MCF-7 40F) spheroids, the corresponding increase was 280% (Fig. 5B). The median dose15
Median dose is the drug concentration required to achieve 50% cell kill.
Validation of the hypothesized functional relationships. This is done by the computational model quantifying the physiologic resistance introduced by the diffusion gradients. Cell viabilities as fraction of control N/NC versus doxorubicin concentrations d; glucose concentration σ∞ = 2.0 g/L; time T = 96 h of drug exposure; n = 3. A, MCF-7 WT drug-sensitive cells. B, MCF-7 40F drug-resistant cells. White columns, three-dimensional in vitro tumor spheroids (dark gray columns, in vitro monolayer data reported for comparison); light gray columns, predictions of the computational model (error bars reflect variation in the apoptosis rate parameter—see Materials and Methods).
Validation of the hypothesized functional relationships. This is done by the computational model quantifying the physiologic resistance introduced by the diffusion gradients. Cell viabilities as fraction of control N/NC versus doxorubicin concentrations d; glucose concentration σ∞ = 2.0 g/L; time T = 96 h of drug exposure; n = 3. A, MCF-7 WT drug-sensitive cells. B, MCF-7 40F drug-resistant cells. White columns, three-dimensional in vitro tumor spheroids (dark gray columns, in vitro monolayer data reported for comparison); light gray columns, predictions of the computational model (error bars reflect variation in the apoptosis rate parameter—see Materials and Methods).
We noted that drug-resistant spheroids maintained more compact, nearly spherical shapes (Fig. 6A) that simply decreased in size as cells were killed, whereas drug-sensitive spheroids formed irregular, looser shapes (Fig. 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).
Correlation of tumor morphologic stability with drug resistance. A, compact, nearly spherical morphology of a drug-resistant MCF-7 40F tumor spheroid in vitro is contrasted with the looser, irregular morphology of a drug-sensitive MCF-7 WT spheroid (B). In the model, stable (C) and unstable (D) morphologies depend on variation in the parameter values for cell adhesion forces 45. Bar, 100 μm.
Correlation of tumor morphologic stability with drug resistance. A, compact, nearly spherical morphology of a drug-resistant MCF-7 40F tumor spheroid in vitro is contrasted with the looser, irregular morphology of a drug-sensitive MCF-7 WT spheroid (B). In the model, stable (C) and unstable (D) morphologies depend on variation in the parameter values for cell adhesion forces 45. Bar, 100 μm.
Discussion
We have used an integrative method that tightly couples computational modeling with biological data to quantify physiologic resistance in three-dimensional 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 chemoprotective effects as cell phenotypic properties and cell-cell and cell-extracellular matrix (ECM) 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 economical 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 (compare Fig. 6). How this relates to cell adhesion molecule expression (e.g., integrin and 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 cell adhesion molecules. We have previously investigated the effect of the tumor microenvironment on tissue morphology (45, 50), suggesting that marginally stable environmental conditions could directly affect morphogenesis and present an additional challenge to therapy (45, 48) 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, noninvasive 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 effect (5). Synergism may be due to increased drug binding to ECM in tumor areas proximal to the drug source, whereas 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, perinecrotic areas. Thus, inducing a large necrotic fraction during treatment may paradoxically further select for drug resistance.
The class of prediction modeling presented in this article, based on an assimilation of complex processes with a fully three-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 toward 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.
Disclosure of Potential Conflicts of Interest
No potential conflicts of interest were disclosed.
Acknowledgments
Grant support: National Cancer Institute, NIH grants U56CA113004 and R01 CA093650 and National Science Foundation, The Cullen Trust of Health Care, NSF-DMS 0818104, National Cancer Institute, Department of Defense.
The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.
We acknowledge John Sinek and Steven Wise (Mathematics, University of California-Irvine), Sandeep Sanga (Biomedical-Engineering, University of Texas-Austin), Paul Macklin (School of Health Information Sciences, University of Texas-Houston), Ernest Han (Obstetrics/Gynecology, University of California-Irvine), Hoa Nguyen (Medicine, University of California-Irvine) for helpful comments and discussions, and the reviewers for their valuable input.