1 1 Introduction

While drug resistance to solid tumors is often a consequence of genetic factors, such as upregulation of anti-apoptotic proteins or overexpression of efflux mechanisms, factors at coarser physiological scales may also profoundly influence tumor therapeutic response [50]. A tumor is a heterogeneous three-dimensional composite of fibrous and connective tissues, stromal components, vasculature, and multiple clones of cancer cells. Atop this intrinsic heterogeneity is layered the anatomical and functional irregularity of tumoral vasculature, characterized by erratic flow, collapsed vessels, diminished oxygen tension, and a large mean tissue-to-vessel distance [5,29,30,34, 54]. As a consequence, the tumor microenvironment is highly variable, marked by gradients of nutrient and oxygen, resulting in regions of hypoxia, acidity, and necrosis, and heterogeneous proliferation. In order for an anticancer agent to be effective, it must extravasate, diffuse through lesion tissue, and be transported into cells, where it must bind to its target and induce cell apoptosis or mitotic inhibition. The tumor environment is not conducive to these processes [50]. The vessel bed’s blood flow and spatial distribution hinder uniform extravasation, calling into question the capability of drug molecules to adequately distribute throughout tissue. Recent experiments in vitro [38,73,81] and in vivo [40,57] demonstrate limited drug penetration through tumors, especially highly protein-bound molecules like doxorubicin and paclitaxel. Once a drug molecule has traversed lesion tissue from its point of extravasation and is presented to a cancer cell, the path from extracellular space to intracellular target is fraught with difficulties ranging from protonation due to the acidic environment, which renders anthracyclines incapable of traversing membrane, to intracellular removal by drug efflux pumps, to cellular processes that effect DNA repair and drug clearance [3,15,31,65,71,79]. In addition to pharmacokinetics, drug pharmacodynamics is equally impaired. Significant hypoxia and hypoglycemia throughout may induce cell quiescence, reducing the efficacy of cell-cycle chemotherapeutic agents like doxorubicin and cisplatin [18,19]. Hypoglycemia causes the glucose-regulated stress response detrimental to the action of topoisomerase II-directed drugs like doxorubicin [41,49,64,68].

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

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

The past two decades have witnessed explosive growth in the mathematical and computational modeling of vascular and avascular tumors [1,4,7–9,13,36,43,44,56, 82]. As the basic elements of tumoral growth models have matured, specialized treatments of angiogenesis and flow [2,46,47,69], drug delivery and response [33,40,52,66,80], and effects of the tumoral microenvironment [23,45] have been developed. Some simulation and analysis has provided great insight using one-dimensional modeling (employing cylindrical or spherical symmetry) without discrete vasculature. More recently, powerful numerical methods have been developed so as to simulate multidimensional complex morphological progression and its relation to cell phenotype and the microenvironment involving, for example, nutrient and biomechanical tissue response [23,45,82].

In the present paper we examine therapeutic efficacy of two common drugs, cisplatin and doxorubicin, in relation to phenotypic and microenvironmental conditions. Model parameters govern extracellular drug/tissue diffusivity; cellular uptake, efflux, and metabolism; cellular density; and the effect of nutrient heterogeneity on drug action. Simulations are performed in a two-dimensional (non-symmetric) setting employing discrete vasculature, which enables the incorporation of morphological and topological influence on drug and nutrient distributions. The effect of these distributions on therapeutic efficacy is of special interest. Sinek et al. [66] had earlier performed a similar investigation; however, the pharmacokinetics and pharmacodynamics (PKPD) component was rudimentary, assuming one homogenous lesion compartment and not based upon experimentally acquired parameter values. In the present work we implement an extensive multi-compartment PKPD component whose parameter values are calibrated via published experimental data. This enables a comparison of the tissue and cell-level drug dynamics of the two drugs, and facilitates the generation of hypotheses to explain their in vivo characteristics. We ask that the reader consider that if doxorubicin and cisplatin were discovered only today, the simulations herein could be seen as providing insight into their anticipated in vivo performance, potentially streamlining and reducing costs of development. Indeed, the methodology presented herein could, with additional development, be applied to both established and nascent drugs to the end of refining clinical trials and assisting in clinical therapeutic strategy to improve patient comfort and survival.

2 2 Mathematical model and parameters

2.1 2.1 Model description

The multiscale tumor growth and angiogenesis simulator developed by Zheng et al. [82] is used to grow the lesions upon which we simulate chemotherapy. This is a nonlinear, continuum scale, two-dimensional growth engine whose accuracy is made possible by an adaptive finite element mesh due to Cristini et al. [11]. Themesh enables multi-scale computation for finely resolving tumor morphology, especially around important areas such as the necrotic/tumor and tumor/host interfaces and around capillary sprouts. Realistic and heterogenous vasculature through which nutrient and drug is provided is generated via the angiogenesis model of Anderson and Chaplain [2]. The simulation field incorporates three phases (see Fig. 2): viable cancerous tissue, normal host tissue, and necrotic debris. The lesion/host interface is demarked by thick black contours, while the microvasculature appears as a web of thin red curves. Dark interior regions are necrotic debris.

Fig. 2
figure 2

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

Full size image

Briefly, nutrient and oxygen are provided through the discrete microvasculature, which is generated in response to angiogenic regulators produced from perinecrotic cells. This results in proliferation and tumor growth. The simple steady-state diffusion equation

$$0 = {k_{\rm{v}}}\left( {1 - n} \right)\delta + {D_n}{\nabla ^2}n - {k_n}n$$

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

Once the tumors are grown, drug administration via the vasculature is simulated by our multi-compartment pharmacokinetics model, based upon earlier work of [15,17, 20,21]. For cisplatin, there are three compartments corresponding to (1) extracellular, (2) cytosolic, and (3)DNA-bound drug. For doxorubicin, there is a fourth compartment corresponding to intracellular organelles, e.g., lysosomes. The system of equations governing transport for both drugs (with different parameter values) is

$$\eqalign{ & {{\dot s}_1} = {k_{\rm{v}}}\left( {{s_{\rm{v}}} - {s_1}} \right)\delta + {D_s}{\nabla ^2}{s_1} - {{k'}_{12}}{s_1} + {{k'}_{21}}({s_2}/{10^6}{V_{\rm{c}}}) \cr & {{\dot s}_2} = {k_{12}}{10^6}{V_{\rm{c}}}{s_1} - {k_{21}}{s_2} + {k_{32}}{s_3} - {k_{23}}{s_2}\left( {1 - {s_3}/{s_{\rm{m}}}} \right) + {k_{42}}{s_4} - {k_{24}}{s_2} \cr & {{\dot s}_3} = {k_{23}}{s_2}\left( {1 - {s_3}/{s_{\rm{m}}}} \right) - {k_{32}}{s_3} - {k_3}{s_3} \cr & {{\dot s}_4} = {k_{24}}{s_2} - {k_{42}}{s_4} \cr & \cr} $$

where s i represents drug concentration in compartment i, k ij represents a transfer rate from compartment i to j, and k i represents a rate of permanent removal from compartment i and the system. s v is intravascular drug concentration during bolus, and s m is a DNA saturation parameter relevant to doxorubicin. V c is the volume of a cell (assumed spherical with diameter 10 µm, yielding V c = 520 fL cell−1) and appears in the first two equations to reconcile the dimensions of s v and s 1 (µM) with the dimensions of all other compartments (fmoles/cell). k v and δ are the same as in Eq. (1). The primed rates appearing in the first equation are related to their unprimed counterparts via k ′ ij = k ij /F where F is the extracellular fraction of whole tissue. Taking a baseline tumor density of ρ = 1.0E9 cells mL−1, a well-known representative value, in combination with the cell volume previously quoted results in

$$F = 1 - \rho {V_{\rm{c}}}\left( {{{10}^{ - 12}}{\rm{mLf}}{{\rm{L}}^{ - 1}}} \right) = 0.48$$

, also a reasonable value. Finally, D s is the diffusivity of the drug through extracellular space.

Both cisplatin and doxorubicin pass through cell membrane according to k 12 (which includes possible pump and transporter activity, as do all other rates). From there, the drugs may efflux according to k 21 or may bind to DNA according to k 23. The kinetics differ from here for the two drugs. Cisplatin may be removed according to the rate k 3, which destroys the functioning of the drug and repairs the DNA [79]. Doxorubicin, however, has an off rate given by k 32, and moreovermay be sequestered and released by lysosomes according to k 24 and k 42 [3,31,59]. Although lysosomal flow to membrane and exocytosis of sequestered drug plays a role in some drug resistant cell lines, we are not necessarily modeling drug resistance via this function, and so assume this process to be negligible in accordance with [16]. On the other hand, we are concerned with the quantity of drug lysosomes can sequester, as this contributes to the cellular uptake of drug, and hence, its penetration characteristics.

The pharmacodynamics model consists of the Hill-type equation along the lines of those employed in [20,21]

$$E = {{N\left( n \right)} \over {1 + {A^{ - 1}}{x^{ - m}}}}$$

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

2.2 2.2 Pharmacokinetics model parameters

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

Table 1 A complete summary of baseline pharmacokinetics and pharmacodynamics parameters
Full size table

Cisplatin parameters We begin with cisplatin, setting k 24 and k 42 to 0 since we assume only three compartments, and k 32 to 0 since we assume the repair rate k 3 is the dominant removal rate of DNA-bound drug. k 3 is next obtained as follows. In experiments performed by Sadowitz et al. [60], adducts per million nucleotides on isolated peripheral blood mononuclear cell DNA fall from 75 to 5 and 185 to 40 in 2 h in two different experiments. Thus, assuming the exponential repair model ˙s3 = k 3 s 3, we calculate the repair rate to be about 0.015 min−1. An initial estimate of k 23 is then made as follows. Sadowitz shows that for 7 µM cisplatin, in 2 h peripheral blood mononuclear cells accumulate from about 25 (non-thiol-blocked cells) to 175 (thiolblocked cells) adducts per million nucleotides. Assuming that DNA consists of about 1.25E6 kbp, this converts to from 1.04E-4 to 7.27E-4 fmoles of Pt docked on the DNA (1 atom/adduct). Neglecting the cell membrane and supposing DNA to be exposed directly to the drug, we have the ODE ˙s3 = 7λ23 − k 3 s 3, where λ23 is a clearance parameter (fL min−1). The solution is s 3 = 7(λ}23/k 3)(1 − exp(−k 3 t)). Substituting values of k 3 = 0.015 min−1, t = 120 min, and 1.04E-4 ≤ s 3 ≤ 7.27E0-4 fmole yields 0.27 ≤ λ23 ≤ 1.9 fL min−1. To convert this to a rate we use the relation k 23 = λ23/V c, arriving at 5.19E-4 ≤ k 23 ⩽ 3.65E-3 min−1. The assumption that DNA was exposed directly to the cisplatin solution means that this rate is only a bootstrap approximation and must be refined. We note that the extremely low ratio of adducts per kbp implies that the saturation capacity of DNA with respect to cisplatin is never approached, and so set s m to ∞.

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

Fig. 1 a
figure 1

Data from Troger et al. [76] used to calibrate k 12, k 21, and k 23 for the cisplatin model. Parameters are fit simultaneously to all four curves treated as one set of data. b Cell inhibition fits for Eq. (4) using Levasseur’s [42] data on A2780 ovarian cancer cells exposed in monolaye

Full size image

Doxorubicin parameters Proceeding to doxorubicin we first obtain an acceptable range for k 12 and k 21 from the literature. For a variety of anthracyclines, including doxorubicin, initial estimates of cell membrane permeability P are taken from experiments with SU-4 and SU-4R wildtype and resistant human lymphoma cells [17], from experiments with EHR2 and EHR2/DNR+ wildtype and resistant Ehrlich ascites tumor cells [15,16], and from experiments with MDA-468 breast cancer cells [40]. The range reported is 2.4 ≤ P ≤ 1000 µ min−1. The relation k 12 = PA c/V c, where A c represents the cell membrane area, can then be used to arrive at an initial range of 1.4 ≤ k 12 ≤ 600 min−1, which will be refined later. In the case of passive diffusion, k 21 = k 12. We note that these values are far larger than those obtained for cisplatin previously. More generally, it has been remarked that cell membrane permeability for cisplatin is much lower than for doxorubicin, etoposide, and vinblastine, although all four drugs are thought to enter cells by passive diffusion [35].

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

DNA binding kinetics of the anthracyclines is nontrivial, perhaps requiring multiple steps and demonstrating sequence specificity [58,59]. Bearing this in mind, as an approximation it will suffice to assume non-specific, one-step binding and unbinding according to the chemical reaction A representative value for the binding coefficient in the above equation for doxorubicin is reported as k on = 4.2E8 M−1 min−1 and a value of the unbinding coefficient (identical with k 32) as k off = 1800 min−1 [59]. From k on we calculate a clearance parameter (as with cisplatin) given as λ23 = k on s m = 4.2E8 fL min−1 (being cautious with the scales of our dimensions). k 23 can then be calculated as λ23/V c, given in Table 1.

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

To refine our initial range of k 12 and k 21, we use the modified version of Eqs. (2) to compare our simulated monolayer uptake profiles of total intracellular drug with those of DeGregorio et al. [14] using human Ewing’s sarcoma and rhabdomyosarcoma cells. At 5.40 min−1 both uptake profiles and equilibrium values compare favorably at three test concentrations.

Drug-tissue diffusivity The last pharmacokinetics parameter values needed are the diffusivities D s of cisplatin and doxorubicin through tumor interstitium. Formolecules of their size (dox M.W. = 544, cis M.W. = 300), diffusivity should be about 30,000 µm2 min−1 [53,70]. However, doxorubicin faces particularly severe barriers due to binding to extracellular constituents such as hyaluronic acid [37,38], and its diffusivity in some tissues has been estimated to be as low as 1000 µm2 min−1 [40], which we take as our baseline value.

2.3 2.3 Pharmacodynamics model parameters

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

3 3 Computational simulations

3.1 3.1 Non-dimensionalization and numerical methods

Non-dimensionalization of Eqs. 2 is via the length- and time-scales

$$ \mathcal{L} = \sqrt {D_s /k_{12} } \mathcal{T} = k_{12}^{ - 1} $$

with the compartmental concentrations transformed as

$$\eqalign{ & {{\bar s}_1} = {s_1}/{s_{\rm{v}}}{\rm{ }}{{\bar s}_2} = {s_2}/1{\rm{E6}}{V_{\rm{c}}}{s_{\rm{v}}}{\rm{ }}{{\bar s}_3} = {s_3}/1{\rm{E}}6{V_{\rm{c}}}{s_{\rm{v}}} \cr & {{\bar s}_4} = {s_4}/1{\rm{E}}6{V_{\rm{c}}}{s_{\rm{v}}}{\rm{ }}{{\bar s}_{\rm{m}}} = {s_{\rm{m}}}/1{\rm{E}}6{V_{\rm{c}}}{s_{\rm{v}}} \cr} $$

resulting in

$$\eqalign{ & {{\dot \bar s}_1} = {{\bar k}_{\rm{v}}}\left( {1 - {{\bar s}_1}} \right)\delta + {{\bar \nabla }^2}{{\bar s}_1} - {{\bar s}_1}/F + {{\bar k}_{21}}{{\bar s}_2}/F \cr & {{\dot \bar s}_2} = {{\bar s}_1} - {{\bar k}_{21}}{{\bar s}_2} + {{\bar k}_{32}}{{\bar s}_3} - {{\bar k}_{23}}{{\bar s}_2}(1 - {{\bar s}_3}/{{\bar s}_{\rm{m}}}) + {{\bar k}_{42}}{{\bar s}_4} - {{\bar k}_{24}}{{\bar s}_{\scriptstyle 2 \hfill \atop \scriptstyle \hfill} } \cr & {{\dot \bar s}_3} = {{\bar k}_{23}}{{\bar s}_2}\left( {1 - {{\bar s}_3}/{{\bar s}_{\rm{m}}}} \right) - {{\bar k}_{32}}{{\bar s}_3} - {{\bar k}_3}{{\bar s}_3} \cr & {{\dot \bar s}_4} = {{\bar k}_{24}}{{\bar s}_2} = {{\bar k}_{42}}\bar s \cr} $$

The numerical methods for tumor growth and angiogenesis have been described in detail in Zheng et al. [82]. For the reaction-diffusion equations (Eqs. 2) we first use Strang splitting. Then the trapezoidal rule is applied to the reaction part, and the Crank-Nicolson scheme, to the diffusion part. For a description of these methods, see Tyson et al. [78] and references therein.

3.2 3.2 In silico experiments

Four in silico experiments are performed in the following manner:

  1. 1.

    Three simulated tumors are grown using Zheng et al.’s model [82]. Each lesion represents one replication of each experiment.

  2. 2.

    The pharmacokineticsmodel (Eqs. 5) is used to deliver drug to the lesions. In each case we hold the intravascular concentration of drug s v constant for 2 h, then set it to zero for eighteen more hours to allow washout.

  3. 3.

    DNA-bound AUC is computed by re-dimensionalizing s 3 and time and integrating using Matlab. The result is then used in the pharmacodynamics model (Eq. 4) to compute cell inhibition. Nutrient (Eq. 1) is relevant when the nutrient effect is employed.

Each lesion is produced based upon the same set of growth and vasculature parameters (see Zheng et al. for a complete description), but randomness in the angiogenesis algorithm and slightly different initial shapes produce different vasculatures and morphologies. It is assumed that a tumor in vivo does not grow or regress appreciably during the 20-h course of the therapy we are attempting to simulate, hence we stop tumor and vascular growth during our in silico therapies. Intravascular concentrations are calibrated in each case to produce a total cellular growth inhibition of 50%. This concentration is referred to as the IC50. The sharp “square wave” delivery of drug is perhaps a caricature of clinical bolus administration, but it allows for consistent analysis and comparison of results.

Our first set of experiments compares DNA-bound drug AUC distributions of doxorubicin and cisplatin under the baseline conditions in Table 1. We furthermore show the homogenizing effect of doxorubicin retention on final DNA-bound AUC [19]. We next investigate the impact of inhibition heterogeneity on dosing requirements, paying particular attention to the nutrient effect for doxorubicin under baseline conditions and improved penetration by, for example, removing hyaluronic acid [37,38]. In our third set of simulations we more deeply investigate the effect of doxorubicin penetration therapies under three circumstances: baseline tumor density, high tumor density, and baseline tumor density with Pgp efflux activity. These are chosen because they demonstrate a spectrum of possibilities due to their effect on cellular drug uptake. High tumor density increases uptake, while Pgp efflux decreases it. In order to simulate increased penetration, we increase D s for doxorubicin from its baseline value to 5000 µm2 min−1 for a moderate increase, and 30,000 for the maximum increase, thus matching the performance of cisplatin. To simulate high tumor density we increase ρ by 50% to 1.5E9 cells ml−1. This has the effect of lowering the interstitial fraction F to 0.22, which in turn increases k′12 and k′21 while leaving all other rates unchanged. Pgp efflux is simulated by increasing k 21 by a factor of 10, which has the effect of reducing all intracellular compartment concentrations by approximately the same factor. This is consistent with results of [10] that show Pgp activity can reduce intracellular concentrations of daunorubicin (an anthracycline related to doxorubicin) by up to a factor of 100. In our fourth and final set of experimentswe investigate permeabilization therapy with respect to cisplatin, whereby a detergent, such as digitonin, or electropermeabilization is used to increase the permeability of cell membrane [35,72]. We take an extreme case, increasing the rate constants k 12 and k 21 from baseline both by a factor of 100. Note that this does not increase the limiting intracellular or DNA-bound levels of drug attained in simulated monolayer, only the rate at which these come to equilibrium. Thus highly permeabilized, DNA-bound max95 is attained at 3.4 h of exposure; further permeabilization reduces this negligibly. For comparison, max95 is greater than 27 h for unpermeabilized cells. This therapy is simulated under both in vivo baseline and very high cell densities achieved by increasing the baseline density 75% to 1.75E9 cells ml−1. At this density, the interstitial fraction F drops to a mere 0.08. Both of these are further compared to monolayer results to probe the conditions under which in vitro assays can be used to predict clinical efficacy.

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

3.3 3.3 Results

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

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

Second experiment We next investigate the impact of drug and nutrient heterogeneity on cell inhibition distributions and IC50’s. Bolus administrations are simulated for cisplatin using baseline lesions exactly as in Fig. 2. The PD model (Eq. 4) is then used to calculate cell inhibition. For doxorubicin we use baseline lesions as well as lesions in which drug penetration therapy is applied. Experiments for dox are run both with and without the nutrient effect.

A table of average IC50’s and log(IC50/ IC50,mono)’s for these experiments is given in Table 2. Here and throughout this paper “IC50,mono” refers to baseline cells exposed in monolayer and serves as a reference. Note that, as these are simulated monolayer exposures, IC50,mono is deterministic. Figure 3 shows a typical nutrient profile, using Lesion II as an example with an IQR of 36%. This measurement is completely analogous to that used in Fig. 2 except that here it is applied to nutrient distribution and there is no normalization since nutrient levels are bounded absolutely from 0 to 100 percent, the level within the vasculature, itself. Nutrient IQR’s for the other two lesions are within 2% of this value.

Table 2 Mean ± SD of the IC50 and the logs of their ratios with respect to monolayer treatments for experiments to investigate the impact of drug and nutrient heterogeneity
Full size table
Fig. 3
figure 3

Contour plot shows nutrient distribution in Lesion II demonstrating significant heterogeneity. Other lesions are similar

Full size image

At the 5% significance level, one-tailed t tests show that the average log IC50 ratio is not greater than 0 for cisplatin, underscoring the homogeneity of its distribution. In contrast, out of the four experiments performed for doxorubicin from the combinations of nutrient effect and penetration therapy, three indicate that the average log ratios are greater than 0 at the 5% significance level. Within this group of four we can analyze the strength of the nutrient effect. For the baseline lesion, the nutrient effect increases the log IC50 ratio by 0.406 units (a factor of about 2.5). For the lesion with penetration therapy, the increase is 0.274 units (a factor of about 1.9). Paired t tests show that these differences are significant at the 5% level.

Cell inhibition distributions closely mirror their AUC distributions, with that of cisplatin being virtually uniform at 50% inhibition throughout. Conversely, doxorubicin displays heterogeneity, increased with the addition of the nutrient effect. Using Lesion II as a representative example for doxorubicin, the upper block of frames in Fig. 4 demonstrates the inhibition distributions for the baseline lesion with and without the nutrient effect. While the broadening of the cumulative probability plot as well as a comparison of the color distribution plots indicate that the nutrient effect increases heterogeneity, inhibition IQR is reduced from 81 to 77% (again, not normalized). The effect of penetration therapy in the lower block of frames is readily apparent. IQR’s, color plots, and probability graphs all indicate more uniform inhibition, ranging moderately from 35 to 65%. Again we see increased heterogeneity in the plots with the addition of the nutrient effect. This time the IQR also reflects the increase. Lesions I and III yield similar results.

Fig. 4
figure 4

(Upper block) Cell growth inhibition profile of Lesion II at baseline settings with and without the nutrient effect after bolus administration depicted in Fig. 2. Probability plot and IQR are now of inhibition distribution and are not normalized with respect to any average. Although the IQR indicates decreased heterogeneity with the nutrient effect, both the color distribution plot and the probability plot indicate increased heterogeneity as is evidenced by the broadening of the curve. (Lower block) The same experiment, except with doxorubicin penetration increased. Now both the plots and IQR show increased heterogeneity. The appropriate IC50 is used in each experiment

Full size image

Third experiment In our third set of simulations, we investigate the effect of therapies designed to improve doxorubicin penetration under several combinations of drug/interstitum diffusivities, cell densities, and drug efflux activities (e.g., Pgp). Figure 5 gives bar graphs of (a) log IC50 ratios and (b) AUC interquartile ranges for three scenarios. The leftmost triplet corresponds to baseline tumor density and no efflux, resulting in a condition of “normal” cellular uptake. The middle triplet corresponds to high density with no efflux, a condition of high uptake. The rightmost corresponds to baseline density with efflux, a condition of low uptake. For the log IC50 ratio (a), in the baseline tumor case there is a change of −0.388 log units in going from no removal of hyaluronic acid to almost complete removal. When density is increased, the change increases to −0.709; however, when Pgp efflux is activated, ANOVA reveals there is no statistical difference, and in fact, the measured change is positive. Results are similar when the nutrient effect is included, with all bars essentially increased by a constant, approximately 0.37. For the AUC interquartile range (b), it is seen that heterogeneity is greatest in the high density case and least in the Pgp efflux case. Within each triplet the heterogeneity decreases with increasing penetration therapy, as expected. The magnitudes of change mirror those for the log IC50 ratios, with the baseline case experiencing a moderate change (from 146 to 13%), the high density case experiencing a dramatic change (from 217 to 52%), and the Pgp efflux case experiencing the least change (from 57 to 4%).

Fig. 5
figure 5

The effect of increasing doxorubicin penetration on a log(IC50/IC50,mono) and b interquartile range shown in three cases: baseline tumor (excepting penetration therapy), high-density tumor, and normal density tumor with Pgp efflux. High density has the effect of increasing drug uptake, while Pgp efflux has the opposite effect. Three replications per bar with standard deviations and results of ANOVA displayed

Full size image

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

Fig. 6
figure 6

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

Full size image

4 4 Discussion

In agreement with experimental observation, our simulations show that heterogeneities of drug, nutrient, and oxygen, caused in part by irregular vasculature and lesion morphology, exist and can significantly impact therapeutic results. Moreover, the sense and magnitude of their influence is not always intuitively obvious. A good example of this is that, despite its well-noted penetration difficulties, doxorubicin performs well clinically. Our simulations show that this may be somewhat explained by its retention in tissue removed from vasculature, causing homogeneity of exposure to increase long after the bolus has been terminated (Fig. 2). This phenomenon has been experimentally verified in [19]with spheroids. Because of this, the resulting cell inhibition distribution is more homogenous than would otherwise be expected. On the other hand, cisplatin maintains an homogenous DNA-bound distribution at all times from bolus initiation to 20 h later, resulting in an extremely uniform cell inhibition distribution. This result, as well as the near equality of its IC50 and IC50,mono demonstrated in Table 2, has also been experimentally verified with spheroids [18,32,37].

While retention in tissue contributes to the performance of doxorubicin in vivo, Table 2 demonstrates that its heterogeneity of distribution contributes to increased serum drug concentrations to match the same cell inhibition in monolayer. In one case, the average amount of drug increases by nearly one log unit. It is reasonable to expect that heterogeneity of nutrient, resulting in hypoxia and hypoglycemia, should compound this problem for doxorubicin. Indeed, this is the case as can be seen by the approximate doubling of the IC50’s (0.482 µM vs. 1.34 and 0.197 vs. 0.371) when the nutrient effect is applied. By graphically and quantitatively showing corresponding cell inhibition distributions. Figure 4 offers further insight into these phenomena. It is easily seen that cell inhibition distributions are as heterogeneous as their corresponding DNA-bound AUC distributions, with areas of lesion removed from vasculature experiencing reduced cell inhibition. An examination of the probability plots in the upper block, corresponding to baseline lesions, shows that a full 24% of viable lesion undergoes no inhibition at all. It is clear from these graphs that penetration therapy greatly decreases heterogeneity of cell inhibition (and commensurately, IC50) as does removal of the nutrient effect. The latter may be clinically feasible through, for example, carbogen breathing or recombinant human Epo (rHuEPO) administration [63,75].

One puzzling behavior is that while both the color and probability plots demonstrate consistently increased heterogeneity brought about by the nutrient effect (as is evidenced by the broadening of the probability curves), the IQR actually decreases in the baseline case from 81 to 77%. This occurs with Lesions I and III as well. A solution to the mystery is obtained by noticing that, absent the nutrient effect, a large portion of the tumor either experiences no inhibition (about 24%, which is distant from the vasculature) or an already heterogeneous inhibition (about 31%). The remaining 45% receives a near homogeneous level of inhibition (the vertical portion of the curve), and this fraction corresponds to tissue close to the vasculature. Thus, the only significant heterogeneity that can be induced by the nutrient effect is within this fraction. Indeed, it is just this part of the curve that broadens in the second probability distribution, indicating greater heterogeneity, as expected. As the IQR is designed to measure heterogeneity somewhat more globally, it misses—in fact, misdiagnoses—the change occurring within this fraction.

In addition to dosing requirements, there is a second and subtler reason to consider heterogeneity of the inhibition distribution when treating clinical tumors. Heterogeneities in microenvironmental conditions have been linked to increased lesion fragmentation and invasiveness [12,25,39,45,55]. While the mechanisms underlying this phenomenon are complex, involving myriad protein signaling events and activities at the cellular level, they may at least partly rely on gross lesion effects, including those caused by drug administration.

Figure 5 shows what might be expected from therapies that increase doxorubicin penetration by, for example, removing hyaluronic acid. As expected, for the baseline tumors, greater homogeneity and level of AUC is achieved, resulting in reductions of IC50. This effect has been experimentally verified using spheroids [37,38]. That it should be more pronounced for high-density in silico tumors and completely absent in the presence of Pgp efflux is intriguing. A potential explanation is availed by simplifying the pharmacokinetics model (Eqs. 2), reducing it to the one-dimensional, one-compartment steady state diffusion equation 0 = D≡2 s − ks with diffusivity D and uptake rate k. In two dimensions, a segment of blood vessel acting as a source next to a section of tissue approximates the one-dimensional case. This equation has one governing parameter, the characteristic diffusion length L = √D/k, and a (non-unique) solution s v exp(−x/L), where x is distance from the source and s v is the constant level of drug in the vasculature. Considering a section of tissue of thickness d next to a vessel and a fixed exposure time, average AUC is proportional to the integral

$$\int\limits_{\rm{o}}^d {{s_{\rm{v}}}\exp \left( { - x/L} \right)dx = {s_{\rm{v}}}L\left( {1 - \exp \left( { - d/L} \right)} \right)} $$

. Let AUC50 be the fixed average AUC required for fifty percent cell inhibition. Then, ignoring the constant of proportionality,

$${\rm{AU}}{{\rm{C}}_{50}} = {\rm{I}}{{\rm{C}}_{50}}L\left( {1 - \exp ( - d/L)} \right)$$

. Increasing the diffusivity to simulate penetration therapy results in an increased characteristic length CL, where C > 1, and hence a new IC50. The ratio of IC50’s is therefore

$${{{\rm{AU}}{{\rm{C}}_{50}}/CL\left( {1 - \exp ( - d/CL)} \right)} \over {{\rm{AU}}{{\rm{C}}_{50}}/L(1 - \exp ( - d/L))}} = {{\left( {1 - \exp ( - d/L)} \right)} \over {C\left( {1 - \exp ( - d/CL)} \right)}}$$

, which approaches 1 as L → ∞, and approaches 1/C as L → 0. Now, increasing cell density has the effect of increasing k, resulting in a smaller L, thus manifesting the differential in IC50’s. Conversely, activating Pgp efflux has the effect of decreasing k, resulting in a larger L, thus nullifying the differential.

That something along these lines is probably happening is demonstrated in Fig. 5b by the generally high IQR’s for the high density case versus the low IQR’s for the Pgp, and the large change in the high density triplet versus the Pgp triplet. Indeed, there is experimental evidence that Pgp and other transporters might facilitate drug transport further away from vasculature [77]. This could be an important point when deciding upon appropriate therapies for tumors exhibiting different characteristics such as efflux mechanisms and relatively high or low densities. Any therapy involves risk. There may be less to gain under certain conditions, advising that the therapy not be performed or perhaps, that concomitant therapy be performed designed to optimize conditions.

The great homogeneity of both cisplatin AUC and cell inhibition demonstrated in Fig. 2 and Table 2 indicate that in vitro assays using this drug may have relatively high positive predictive accuracy. While our simulations do not yield enough resolution to claim that the results herein answer this question, the outcome shown in Fig. 6 is of interest. There is no statistical difference at the 0.01 significance level between improvement obtained by permeabilization therapy in the simulated tumors (even very dense ones) and monolayer. This compares favorably with the findings of [72] inwhich the improvement obtained via cisplatin permeabilization therapy in vitro is strongly reflected by the improvement in isolated lung perfusion in rats.

While it can and should be argued that the simulations herein fail to account for some (many!) critical aspects of tumor growth and drug response (such as clonal heterogeneity, cell phase sensitivity, and signaling pathways) and that parameter settings may in some cases be inexact, it should not be concluded that these shortcomings invalidate characteristics the simulations have revealed. Indeed, we have correctly post-dicted several results: that doxorubicin retention results in a more uniform AUC and cell inhibition than would otherwise be indicated by its penetration difficulties [19]; that Pgp and other transporters might facilitate drug penetration into lesion tissue [77]; that cisplatin achieves a highly uniform AUC, and its cellular in vitro monolayer inhibition can closely match that of in vitro spheroids [18,32,37]; and that improvement due to cisplatin permeabilization therapy in vitro has the potential to predict improvement in vivo [72]. Through the model we have provided evidence that not only do macroscopic environmental conditions, namely, drug and nutrient distribution heterogeneity, potentially greatly impact therapeutic efficacy, but also that the outcome of therapeutic strategies can depend upon them in nonlinear and a priori unpredictable ways. The results of our third set of experiments (doxorubicin penetration therapy) provide an example. In light of this, it would be prudent to devote attention to factors residing at coarser and more global scales than solely the genetic.

One of our broader goals is to demonstrate how increasingly sophisticated in silico technology, driven by mathematical modeling and calibrated with experimental data, can and is being developed to provide an alternate investigative and clinical tool complementary to traditional methods [6,24,61,62,67]. It can well be imagined that were doxorubicin and cisplatin discovered today, the in vivo simulations herein presented could be used to anticipate their lesion- and cellular-scale pharmacokinetics, helping to refine clinical trial design and lower costs. In clinical application, the results could be used to guide therapeutic strategy. For example, any risks associated with doxorubicin penetration therapy could be minimized if it were known that the patient’s tumor were expressing Pgp or otherwise had lowered cellular uptake, according to the results given in Fig. 5. With further development, we anticipate that in silico models and methods similar to those in this paper will become increasingly accurate and useful. Towards that end we are working on incorporating a more veridical model of vasculature along the lines of McDougall and Stephanou, in which blood flow plays a key role in the formation of vasculature [46,69]. We are furthermore developing 3D models, which are proving successful at accurately simulating morphological evolution [23]. One area of special interest is that of liposomal and nanovectored delivery. These have the potential of unprecedented accuracy and specificity of delivery [22]. It would be interesting to adapt the PKPD model herein, paying particular attention to transport in and extravasation from blood vessels [28], to see what predictions can be made.

The power of in vitro experimentation lies in its ease of implementation while remaining in the biological realm. By its very nature, in vitro experimentation attempts to refine and isolate. Yet, much of what happens in vivo is the result of a nonlinear system whose behavior is more than the sum of its parts. The power of in silico simulation lies in its ability to integrate components into a virtual system capable of reproducing such behavior, implicitly taking into account circuits of information flow difficult to explicitly analyze. Accurately calibrated and rigorously validated, such an integrated model could provide a “dry-lab” to be used as a powerful complement to the traditional wet-lab in fundamental research, drug discovery, and the clinic. It could be used to probe scenarios and test hypotheses that are either difficult or impossible to instantiate in the body. Results could then suggest supportive in vitro and in vivo experimentation, the end result being new therapeutic targets or strategies. Simultaneously, weaknesses (or strengths!) of the in silico model could be uncovered and addressed. Computational models have the potential to facilitate an era of great discovery and progress in understanding and treating cancer, and providing new hope to its victims.