| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Mathematical Oncology |
1 School of Health Information Sciences and 2 Division of Nanomedicine, and 3 Department of Biomedical Engineering, University of Texas Health Science Center, Houston, Texas; Departments of 4 Anatomic Pathology and 5 Experimental Therapeutics, and 6 Systems Biology, The University of Texas M. D. Anderson Cancer Center, Houston, Texas; 7 Department of Bioengineering, Rice University, Houston, Texas; 8 Department of Biomedical Engineering, The University of Texas, Austin Texas; 9 Division of Hematology/Oncology, Department of Medicine, University of California, Irvine Medical Center, Orange, California; 10 School of Pharmacy, Centre for Biomolecular Sciences, University Park, University of Nottingham, United Kingdom; 11 Departments of Radiology and Integrated Mathematical Oncology, Moffitt Cancer Center, Tampa, Florida
Requests for reprints: Vittorio Cristini, School of Health Information Sciences, University of Texas Health Science Center, 7000 Fannin #600, Houston, TX 77030. Phone: 713-500-3965; Fax: 713-500-3929; E-mail: vittorio.cristini{at}uth.tmc.edu.
 |
Abstract |
|---|
| Major Findings: 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 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 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 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 |
|---|
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 |
|---|
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 x 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 4x nanomolar increments (0, 4, 16, 64, etc.) for 96 h, representing at 256 nmol/L a typical area under the curve (AUC)12 in vivo (29). Negative controls were seeded and incubated under the same conditions without drug. Three end points were concurrently measured as fraction of negative control: proliferation using tritiated thymidine (Amersham) 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 (x100 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 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 x 10–7 cm2/s (37) and uptake rate may be as large as 1 x 10–3 s–1 (38), whereas for oxygen, the corresponding values are 1 x 10–5 cm2/s (39) and 1 x 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 
, where 
and 
. 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. Physiologic resistance based on cell cycle status was modeled with A''. Decreased substrate concentration in three dimensions results in cell populations that cycle less and therefore have reduced sensitivity to doxorubicin.13 In the above formula, enhanced survival was assumed to be 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 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 |
|---|
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.
|
|
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).
|
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 (Fig. 4A
). Further, increased viability for drug-sensitive cells in three dimensions corresponded with 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 (Fig. 4B).
|
20% higher in spheroid versus monolayer than for MCF-WT, implying not only that the three-dimensional morphology promoted the net survival for both cell types but also that the resistant phenotype was further favored by the three-dimensional configuration.
|
|
|
Discussion |
|---|
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 |
|---|
|
Acknowledgments |
|---|
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.
|
Footnotes |
|---|
13 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).
14 In vitro exposure to 256 nmol/L for 96 h yields an AUC of 2.5 x 10–5 mol/L h, whereas for a typical patient doxorubicin plasma levels after a bolus injection represent an exposure of 3 x 10–5mol/L h.
15 Median dose is the drug concentration required to achieve 50% cell kill.
Received 9/25/08. Revised 1/26/09. Accepted 2/25/09.
|
References |
|---|
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |
Top
Abstract
Quick Guide to Equations...
V), dead tumor tissue (
