Theory and Experimental Validation of a Spatio-temporal Model of Chemotherapy Transport to Enhance Tumor Cell Kill

It has been hypothesized that continuously releasing drug molecules into the tumor over an extended period of time may significantly improve the chemotherapeutic efficacy by overcoming physical transport limitations of conventional bolus drug treatment. In this paper, we present a generalized space- and time-dependent mathematical model of drug transport and drug-cell interactions to quantitatively formulate this hypothesis. Model parameters describe: perfusion and tissue architecture (blood volume fraction and blood vessel radius); diffusion penetration distance of drug (i.e., a function of tissue compactness and drug uptake rates by tumor cells); and cell death rates (as function of history of drug uptake). We performed preliminary testing and validation of the mathematical model using in vivo experiments with different drug delivery methods on a breast cancer mouse model. Experimental data demonstrated a 3-fold increase in response using nano-vectored drug vs. free drug delivery, in excellent quantitative agreement with the model predictions. Our model results implicate that therapeutically targeting blood volume fraction, e.g., through vascular normalization, would achieve a better outcome due to enhanced drug delivery. Author Summary Cancer treatment efficacy can be significantly enhanced through the elution of drug from nano-carriers that can temporarily stay in the tumor vasculature. Here we present a relatively simple yet powerful mathematical model that accounts for both spatial and temporal heterogeneities of drug dosing to help explain, examine, and prove this concept. We find that the delivery of systemic chemotherapy through a certain form of nano-carriers would have enhanced tumor kill by a factor of 2 to 4 over the standard therapy that the patients actually received. We also find that targeting blood volume fraction (a parameter of the model) through vascular normalization can achieve more effective drug delivery and tumor kill. More importantly, this model only requires a limited number of parameters which can all be readily assessed from standard clinical diagnostic measurements (e.g., histopathology and CT). This addresses an important challenge in current translational research and justifies further development of the model towards clinical translation.


Introduction
The biological drivers of cancer promote a physical microenvironment that differs significantly from normal tissues [1]. The underpinnings of the dysregulated physical properties of cancer are an active area of investigation. Research efforts in this area are highlighted by modern oncologic concepts of how the physical microenvironment influences tumor behavior (i.e., theories such as Transport Oncophysics and vascular normalization) [2,3], and are supported by experimental data in animal cancer models and cancer patients [4,5]. Mathematical modeling has provided quantitative insights into the understanding of how these physical properties impact therapeutic drug transport, which will eventually help to optimize, predict and provide treatment strategy on an individual patient basis [6][7][8][9][10]. Notable modeling successes with validation against in vitro, in vivo, and/or patient data include mathematical modeling of cancer tumor development [11][12][13][14] and traditional drug delivery in neoadjuvant [15,16], adaptive [17], and standard chemotherapy treatment conditions [18]. We have also studied how mathematical modeling of the mass transport of drugs can mechanistically describe the therapeutic response to chemotherapy [19][20][21][22] and enable an understanding of the drug delivery process in humans [23]. However, to date, these efforts have been limited by the inability to account for spatial and temporal heterogeneity in drug dosing and tumor characteristics.
Chemotherapy drugs need to traverse the vasculature, interstitial space (i.e., stroma and microenvironment), and cancer cell membranes to finally reach intracellular targets. At multiple biological scales, tumors have properties that impede the delivery of chemotherapy, and the poor delivery of drugs prevents the killing of cancer cells [3,24,25]. In solid tumors, the disorganized and leaky vasculature of the tumor microenvironment influences the delivery of and response to systemic chemotherapy [2,[26][27][28]. Furthermore, the stroma not only acts as a physical barrier to drug delivery, but also impairs the function of blood vessels [29,30]. Conceivably, the ability to characterize these physical properties of tumors would enable scientists and clinicians to not only predict response, but also rationally design therapeutics for an individual cancer patient.
Toward this goal, we have developed a theory of chemotherapy response based on the quantitative physical transport properties of cancer cells [20,21] as well as solid tumors [19,22,23]. Our approach describes how a tumor's biophysical properties affect drug delivery using "quantitative pathology," a translational approach combining mathematical modeling and measurements from histopathology of individual patients' tumors to predict treatment outcome [19][20][21][22]31]. In particular, the mechanistic model presented in [19] accurately predicted the fraction of tumor killed (denoted by f kill ) from chemotherapy in patients with colorectal cancer (CRC) metastatic to liver and in patients with glioblastoma. We also demonstrated the feasibility of the model for clinical translation by predicting patient-specific treatment efficacy using computed tomography (CT) scan data [19,23]. The model presented in [20] is a time-dependent model based on first principles of cell biophysics to predict in vitro the nonlinear dose response curves for two types of delivery methods, free drug and targeted nano-carriers. We found that drug delivery mediated by nano-carriers targeting the tumor cells overcomes resistance to free drug because of improved cellular drug uptake rates.
We previously hypothesized that tumor kill would be significantly enhanced through the elution of drug from a nano-carrier that deposits in the tumor [3]. A series of pilot studies have tested this hypothesis using experimental tumor models for delivering drugs via multistage vectors (MSVs) [32][33][34]; MSVs are nested particles (with smaller particles nested inside larger ones), particularly designed to lodge in tumor vasculature to release nanoparticles over an extend period of time (e.g., several weeks). To date, mathematical modeling validated against experimental laboratory results has resulted in insights into fine-tuning nano-carriers for optimized drug delivery through charge and size optimization [35,36], drug release duration [37,38], and targeting molecule surface loading [39]. Work to this end has also resulted in increased understanding of how vasculature structure and the resulting interstitial fluid behavior are involved in distribution of the nano-carriers in the blood vessels in vivo [40]. Here, we develop a mathematical model to quantitatively formulate the hypothesis of drug delivery via loaded nano-carriers. Specifically, we extend our previous time-dependent modeling work [20] by also accounting for spatial dependence in predicting tumor response to systemic agents. This generalized model allows us to consider a variety of treatment strategies, including systemic drug delivery via nano-carriers, and helps to predict the tumor response to different forms of drug delivery methods before the start of treatment. Model predictions are further validated using experiments on a breast cancer mouse model in vivo.

Generalized mathematical model
We extend the time-dependent drug-cell interaction model [20] by incorporating spatial dependence to describe perfusion and diffusion heterogeneities. The governing equations for drug concentration σ(x, t) and the volume fraction of tumor cells φ(x, t) are @s @t ¼ Dr 2 s À l u φs; ð1Þ where D is the diffusivity of the drug, λ u the per-volume cellular uptake rate of drug, and λ k the death rate of tumor cells per unit cumulative drug concentration. While drug uptake is known to have a nonlinear dependency on drug concentration [41,42], we assume a linear drug uptake for simplicity. Because drug diffusion time and the plasma half-life of drug are both much shorter than the time scale for cell death (on the order of minutes vs. hours or days), and also because the model will be examined on time scales of days to weeks, rather than minutes, Eq 1 can be solved at the steady state (i.e., @σ / @t ffi 0; note that this is actually a quasi-steady state, meaning that σ(x, t) quickly relaxes to the instant steady state defined by φ(x, t)). Thus, without the time derivative in Eq 1, the solution σ(x, t) is independent of initial conditions; for boundary conditions, we set a drug concentration σ 0 at the blood vessel wall. We further clarify that the steady state solution of σ(x, t) still has spatial dependence, and hence a gradient; it is only steady state with respect to time. The main assumption underlying the use of a steadystate equation with a constant boundary condition was that the entire time curve of several bolus injections over several months of treatment for each patient could be satisfactorily replaced with a constant σ 0 equal to the time averaged drug concentration throughout the entire multi-month treatment regimen. That is, we assume that a drug administered as bolus at a certain dose level has the same effect as the same total amount of drug administered over several months at a constant, smaller dose level; see [19,23] for validation of this assumption (with patient data) on the use of a constant boundary condition. For a cylindrically symmetric domain surrounding a blood vessel (Fig 1), the boundary conditions can be set to where r denotes the radial position from the center of the cylinder, and r b represents the blood vessel radius; the second boundary condition reflects that the far-field drug concentration flattens out. Furthermore, Eq 2 has no spatial derivatives, and thus only requires the initial conditions for φ(x, t), which we set to i.e., a homogeneous initial tumor volume fraction. As detailed below, this generalized model allows us to examine not only successive (conventional) bolus chemotherapy, characterized by a time-varying intravenous drug concentration σ 0 according to a specific dosing and timing regimen, but also drug release through loaded nano-carriers where drugs are released at a nearly constant rate over a certain time interval, approximated here by a constant σ 0 . Illustration of transport-based hypothesis. By diffusion, a blood vessel supplies substrates to the cylindrical tissue volume surrounding the vessel. We hypothesize that at each position inside the tissue, the substrate supply is supported by the closest blood vessel. Thus, the influenced tissue surrounding a vessel can be estimated to be between a cylinder of radius r b / (L BVF 1/2 ) in dimensionless form, and the vessel itself with dimensionless radius r b / L. Theoretically, chemotherapeutic drugs delivered by a blood vessel kill the tissues immediately adjacent to the vessel, leaving some viable tissues on the far end. Here, we propose that through drug-loaded nano-carriers that can accumulate within tumors and continuously release drugs for a longer time (e.g., lasting several cell cycles), the drugs can penetrate further into the surrounding tissue volume and thus achieve a higher tumor killing ratio.

Model normalization
The generalized model (Eqs 1 and 2) is nondimensionalized as where x 0 = x / L and t ' = t / T are the dimensionless space and time coordinates, with T = (λ k λ u φ 0 σ 0 ) −1/2 (a characteristic time of drug-induced apoptotic cycles; as shown previously [20], this corresponds to the short-term apoptosis time, as the drug uptake reaches significant levels). The dimensionless drug concentration is σ ' = σ / σ 0 , and the tumor volume fraction is made dimensionless by φ ' = φ / φ 0 . Similarly, the boundary conditions in Eqs 3 and 4 can be nondimensionalized as with r ' = r / L, while the initial conditions for the tumor volume fraction in Eq 5 become Note that the nondimensionalization leaves no parameters in the governing differential equations (Eqs 6 and 7), and only special boundary conditions will introduce biophysical parameters into the model.
Considering the cylindrically symmetric domain surrounding a blood vessel, to which the boundary conditions Eqs 8 and 9 apply, Eq 9 can be further revised to Thus, the boundary conditions Eqs 8 and 11 introduce two dimensionless parameters in our numerical analysis: r b / L and BVF. Numerically, we use a two-step process to integrate the generalized model over time. First, we use φ ' from the previous time step and solve Eq 6 for the steady-state solution of σ ' to the first order of Δr, the spatial mesh size. Next, we substitute this σ ' into Eq 7 and compute φ ' for the current time step using a backward Euler method [43], which is accurate to the first order of Δt, the numerical time step size.

Formulation of tumor kill from treatment
By integrating the viable tumor volume fraction at each time point over the cylindrical tissue domain surrounding a blood vessel and affected by the drug diffusion, we calculate f kill as the ratio of the killed tumor volume to the total initial tumor volume: as a function of parameters: r b , blood volume fraction (BVF), and L ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi D=ðφ 0 l u Þ p (the effective diffusion penetration length of the drug). As shown in Fig 1, the model domain is comprised of the space between two concentric cylinders. The inner cylinder has a radius r b / L in dimensionless units, representing the blood vessel at the center of the domain. With the hypothesis that the substrate supply for any spot in a tissue is supported by the closest blood vessel, we estimate r b / (L BVF 1/2 ) to be the dimensionless radius of the influenced tissue volume of the vessel, where BVF < 1. The influenced tissue volume refers to a specific region of tissue that relies on this blood vessel for supply of oxygen and other essential chemicals.
Special formulation of f kill for time-course tumor measurements Histology data are not always available for determining values (for parameters BFV, r b , and L) needed for using Eq 12. Thus, we develop an alternative form of f kill as a function of another set of experimental parameters, the values of which can be obtained from in vivo cytotocixity experiments. By substituting Eq 1 into Eq 2, we have Integrating φ in Eq 13 over the total tissue volume V, we obtain the changing rate of the tumor volume: We denote the tumor volume as V T ; by the convergence theorem, Eq 14 can be rewritten as where @V represents the boundaries of the total tissue volume in question here, and @σ/@n is the flux of the drugs across the boundaries. It is safe to hypothesize that the flux of the drugs becomes negligible at the tissue boundaries far away from the blood vessel, and hence the only contribution in the boundary integral we consider is the flux at the boundaries next to the blood vessel. For simplicity, we define We have previously demonstrated that in vitro porous silicon particles achieve a constant release rate of doxorubicin for up to two weeks at neutral pH [44]. Note that F is dependent on nanoparticle size; for the same type of nanoparticle, the bigger the size, the slower is the drug release rate [45]. We thus further hypothesize that the rate of change of flux for the first several days is approximately zero, dF / dt = 0 (i.e., this initial time period is too short for F to change significantly); hence F is constant, and we have: which leads to or equivalently, As defined previously [19], 1−(V T /V T,0 ) is exactly our definition of f kill , i.e., the fraction of tumor volume killed from chemotherapy. Hence, we obtain a new mathematical formula for calculating the amount of f kill through the delivery method of loaded nano-carriers: Note that there is a quadratic increase in f kill with time, which is consistent as previously observed in vitro [20] (cfr. Eq. 4a in that reference).

Determination of reference values for model parameters
For the generalized space-and time-dependent model (Eqs 1 and 2), patient data from multiple time points are usually necessary to determine the model parameters. However, for each patient in our data set obtained from a cohort of 21 patients with CRC metastatic to liver, we had histology data from one time point (after chemotherapy) and contrast CT imaging data from two time points (before and after chemotherapy). Thus, we directly applied Eq. S2 to obtain patient-derived parameter values for BVF, r b , and L (S1 Text). Briefly, least-squares fitting of Eq. S2 was performed using Mathematica routine "NonlinearModelFit" [46] to the f kill and BVF measured from the histopathology (S1 Fig and S1 Text). This resulted in estimates of r b and L parameters which produced the best fit (Fig 2A); the values are consistent with published data [47,48]. Statistically significant P values were obtained (inset). Moreover, using regression analysis, we also determined a linear correlation between contrast CT enhancement (in Hounsfield Units) and measured BVF values (Fig 2B and S1 Text). Standard deviations for CT data reflect 25% variability in physiology and contrast-injection protocols across patients, and was estimated by calculating standard deviation of CT measurements taken within the patient aorta. CT measurements reflect perfusion of tissue, which relies on the volume fraction of blood vessels, hence resulting in P < 0.001. The linear correlation (Fig 2B) was used to inform model Eq. S2 directly from pre-treatment CT data to predict f kill for each patient. The results in Fig 2C show that there is no statistically significant difference between the predictions obtained from the linear correlation and direct measurements of kill from histopathology of the resected tumors. Lastly, the predicted f kill based on CT scans prior to treatment (Fig 2D, open circles) matched the measured f kill (filled circles) after treatment with an average relative error of % 24%. Hence, we used the obtained, verified parameter values (for BVF, r b , and L) as reference values for subsequent simulation studies.

Parameter analysis of the generalized model
We performed a set of simulations with the generalized time-and space-dependent model in a cylindrically symmetric domain surrounding a blood vessel. Model parameter values were set as follows (Fig 2A): L = 155.06 μm and r b = 15.83 μm, giving r b /L = 0.102; we also set BVF = 10 −2 as the reference value for the simulations, as measured BVF values ranged approximately from 10 −3 to 10 −1 . We ran the model (Eqs 6 and 7) for 10 drug-induced apoptotic cycles. To examine the impact of each parameter on f kill , we further simulated nine parametervariation combinations, using three r b /L values, i.e., 0.05, 0.1, 0.5, paired with three BVF values, i.e., 0.005, 0.01, 0.05. Fig 3 shows the numerical results of the model. In the presence of a boundary condition σ 0 at the vessel wall (r = r b ), successive cell layers next to the blood vessel die out (Fig 3A) due to enhancement of drug penetration (Fig 3B), in turn leading to accelerated cell kill ( Fig 3C). As cell kill occurs, tumor volume fraction φ decreases, leading to an increase in local drug concentration σ (because dead cells no longer take drugs), and thus accelerating cell kill in the locations further away from the vessel and deep into the tumor. Fig 4A shows the temporal evolution curves of f kill calculated from Eq 12 by varying the parameters r b /L and BVF, representing conditions where drug-loaded nano-carriers are employed, as well as the estimates using our simplified "bolus" model (i.e., Eq. S2). The results indicate that the "bolus" killing ratios are readily achieved by drug-loaded nano-carriers after 1 or 2 cell cycles. To estimate the benefits of drug release by the nano-carriers over a longer period of time, we normalized the f kill curves using their corresponding bolus f kill values ( Fig  4B). The results suggest that we may achieve 2-to 4-fold of the bolus killing ratios if the drug release from nano-carriers administration lasts for 3 or 4 apoptotic cycles. However, for large BVF values (representing highly vascularized tumors), cell killing effects from both methods of delivery are roughly equivalent. This is expected because the majority of the tumor cells are killed within just one or two apoptotic cell cycles; a 50% increase in tumor kill is nevertheless expected from loaded nano-carriers releasing over a longer period of time. This suggests an alternative strategy to improve chemotherapeutic efficacy by promoting or normalizing angiogenesis at the target site before administrating chemotherapy drugs [2,[49][50][51], or by promoting perfusion by other means such as mild hyperthermia [52], both of which would lead to an increase in BVF. We also note that, in our simulations, the domain size (i.e., the radius of the outer cylinder) was determined by r b / (L BVF 1/2 ). Thus, a larger r b / L or a smaller BVF represented a larger tissue volume relying on the blood vessel for drug transport, and hence required a longer time to achieve the same f kill .  Experimental validation Ethics statement. The animal studies were performed in accordance with the guidelines of the Animal Welfare Act and the Guide for the Care and Use of Laboratory Animals following protocols approved by the Institutional Animal Care and Use Committee (IACUC).
Cytotoxicity experiments: Comparison of free and nano-particle-based drug delivery in vivo. A total of 40 six-week-old female BALB/c mice were obtained from Charles River Laboratories and kept in a specific pathogen-free facility at the Houston Methodist Research Institute. On day zero, mice were inoculated with 5 x 10 4 green fluorescent protein (GFP)-labeled murine 4T1 cells in their inguinal mammary fat pad. Breast tumors were allowed to grow for two weeks to reach a volume V T,0 = 100-200 mm 3 . Mice were then randomly placed into four groups (10 mice per group) and administered intravenously beginning on day 14 according to a predefined protocol. The four groups were: i) control: phosphate buffered saline (PBS), administered twice a week; ii) free doxorubicin (3 mg/kg, i.v.), administered twice a week; iii) 1.0 μm porous silicon particle loaded with chemotherapy drug (iNPG/pDox 1.0, 6 mg/kg, i.v.) [53], administered once a week; and iv) 2.6 μm porous silicon particle loaded with chemotherapy drug (iNPG/pDox 2.6, 6 mg/kg, i.v.), administered once a week. Tumor volume was measured for each mouse on days 14, 17, 21, 25, 28, and 31 after tumor cell injection. Mice were sacrificed on day 31 via CO 2 asphyxiation and tumors were removed. For comparison with model predictions (i.e., Eq 20), the tumor volume measurements were normalized across the four treatment groups to the measurements from the PBS control group and to the initial tumor volume, and for each tumor, f kill was calculated as 1 minus the normalized tumor volume.
Model predictions. From the time-evolution of f kill for the three groups of BALB/c mice (Fig 5), it is evident that, after roughly three days of first treatment with rapid growth, f kill remains approximately constant until the end of the experiments (f kill = 0 at the onset of the treatment on day zero). This rapid growth of the fraction of dead cells is consistent with the quadratic time-dependence predicted by model Eq 20. The measured tumor kill from nanovectors is about 0.5, and roughly 3 times that from free drug, in excellent agreement with the model predictions of a 2-4 fold increase in kill depending on the parameter values (cfr. Fig  4B). From the experimental protocol, we know that the total amount of drug released by the particles is FÁt % 1.2Á10 −4 g for a typical mouse weight of 20 g [44,[54][55][56]. We then analyzed the tumor growth curves (S3 Fig) and estimated the approximate (linear) growth rates. For example, the controls were found to grow at an average rate of Λ % 70 mm 3 / day (proliferation only; no death), while the tumors in mice treated with iNPG/pDox 1.0 μm grow at a rate of % 35 mm 3 / day (net outcome of proliferation minus death rates). The latter result produces a net death rate for the iNPG/pDox 1.0 μm treated tumors of Λ k % 35 mm 3 / day, which, since the specific rate of kill (per molecule of drug) is λ k = Λ k / (Ft), gives, together with Eq 20, an estimate for t kill ¼ 2V 0 L k Á f kill % 4 days (corresponding to f kill % 0.5 and V 0 = 130 mm 3 ), in excellent agreement with the observed time to plateau of the cell kill reported in the experiments (Fig 5). Note that subsequent treatments in the protocol are relatively irrelevant, as the amount of drug used is the same whereas the tumors have already grown by one order of magnitude by the time the second treatment is applied (i.e., one week).

Discussion
A fundamental principle of oncologic therapy is that cancer cell kill must be absolute to achieve a cure. For patients with metastatic cancer, conventional systemic chemotherapy alone does not eradicate all of the disease [57]. This therapeutic challenge is evident in clinical practice with both cytotoxic and cytostatic therapies, and our mathematical model of cellular response to chemotherapy highlights a physical mechanism for this clinical observation [19,22]. Indeed, in our original report of this mathematical description of cancer cell kill due to local drug concentrations, we found that complete eradication of the cells was impossible using current methods to infuse drugs intravenously [23]. Here, we have confirmed this finding in a larger dataset of patients with colorectal liver metastases. Notably, we find that the pathological response to chemotherapy is heterogeneous within a given tumor and that the local physical properties of the tumor describe this response. This is significant in understanding therapeutic resistance, suggesting that the physical microenvironment naturally selects cancer cells that reside in areas with poor drug penetration.
This observation of heterogeneous response inspired a concept of improving drug delivery using nanoparticles that accumulate within the tumor, delivering sustained local release of drug to the cancer cells. We simulated this concept by generalizing our mechanistic model and applying it to the same patient group with colorectal liver metastases. This generalized model now accounts for the spatial and temporal heterogeneities of drug dosing, which builds on our prior work [19,20]. From the parameter analysis (Figs 3 and 4), we find that an extended period of treatment achieves a better treatment outcome (i.e., killing more cancer cells). More importantly, we found that the delivery of systemic chemotherapy using these nanoparticles would have enhanced cell kill by a factor of 2 to 4 over the standard therapy that the patients actually received (Fig 4B). However, this strategy may not be feasible in practice due to the constraint of tolerable cytotoxicity to healthy cells; further investigation into tumor-targeted nano-carriers may be necessary to realize the results in patients that our model predicts.
In our model, for the case of cells that stay further away from the blood vessel, the reduction of cell death is only due to lowered local drug concentration. Indeed, these cells may also experience a change in other microenvironmental conditions (such as hypoxia and nutrient deprivation), which can lead to reduced proliferation rates and therefore reduced sensitivity to cytotoxic drugs. To address this issue, we could add oxygen concentration as a parameter to the model, as demonstrated in our prior work [37]. However, in light of the insight we gained from our model for enhanced drug penetration (Fig 3), oxygen may also penetrate deeper into the surrounding tissue because dead cells no longer uptake oxygen. Therefore, our assumption of a fixed cell sensitivity to drugs may remain valid even without taking oxygen into consideration. In our ongoing work, we are validating our model results shown in Fig 3 with in vivo experiments.
Another observation from our model is that a larger BVF results in a higher tumor killing ratio. One may envision extending our generalized model to other therapeutic strategies that aim to improve efficacy through enhanced drug delivery by increasing BVF. These strategies include, but are not limited to, vascular normalization [2], metronomic dosing [58], and stromal targeting [29,30]. Each of these physical sciences-based therapeutic strategies could be optimized through our predictive mathematical model.
As an experimental demonstration of this concept of drug-loaded nano-carriers, we used a porous silicon particle delivery system loaded with chemotherapy in a mouse model of breast cancer. We have found that these particles localize within tumors [59], likely due to the geometry of the particles. As predicted by our model, we measured a 3.5-fold difference in tumor growth rate in the mice tumors treated with drug-loaded nanoparticles compared to the tumors treated with systemic delivery of doxorubicin (Fig 5). This experimental confirmation of the model results provides rationale to translate these findings and concepts to patients.
In the future, additional layers of complexity involving other factors or physiological barriers, such as tumor cell regrowth or proliferation that repopulates the killed region, effect of cell to cell contact, and effect of chemotherapy on the tumor vasculature [25,37,38], can be added to the model in an effort to improve the accuracy of model predictions. Importantly, this model applies to any systemic agent, including immunotherapy. For example, the model may be generalized even more to account for cancer cell kill by the immune system, accounting for the physical barriers to immune cell infiltration into the cancer. Moreover, through the use of non-invasive characterization of transport prior to therapy using diagnostic CT or MRI imaging, as well as the biological characterization of molecular targets for an individual tumor [60], one could optimize both drug delivery and therapeutic selection for a given patient. This biophysical characterization and prediction strategy would complement genetically-based, patient-specific cancer therapy methods by individualizing drug administration regimens.