Towards Multiscale Modeling of Nanovectored Delivery of Therapeutics to Cancerous Lesions

H.B. Frieboes¹, J.P. Sinek², S. Sanga¹, F. Gentile⁴, A. Granaldi⁵, P. Decuzzi^4,5, C. Cosentino⁶, F. Amato⁶, M. Ferrari³ and V. Cristini^1,2,*

Departments of Biomedical Engineering¹ and Mathematics², 3120 Natural Sciences II, University of California, Irvine, CA 92697-2715

^*Corresponding author: cristini@math.uci.edu. Present affiliation: School of Health Information Sciences, University of Texas Health Science Center at Houston

³Davis Heart and Lung Research Institute, The Ohio State University, Columbus OH 43210-1002

Center for BioNanotechnology and Engineering for Medicine⁴, University of Magna Græcia, 88100, Catanzaro, ITALY,

Center of Excellence in Computational Mechanics⁵, Politecnico di Bari, 70125, Bari, ITALY

Biomechatronics Lab, Department of Experimental and Clinical Medicine⁶, University of Magna Græcia, 88100, Catanzaro, ITALY

Abstract

If an oncologist could know beforehand the effects, both therapeutic and toxic, of a given drug on a given patient, then not only would Hypocrites’ maxim, “Primum non nocere,” remain inviolate, but also would the patient stand to gain tremendous benefit from receiving the right drug. Although standards of care for many cancers have been established, this alchemic desire is frustrated by a host of systemic, tissue-scale and cellular resistance mechanisms that yield disappointing differentials between treatment predictions and clinical results. Computational models may provide a much-needed bridge between the two, producing highly realistic in silico tumors upon which alternate therapies may be conducted. The power of such models over experimental assays lies in their ability to integrate processes over a multitude of scales, approximating the complex in vivo interplay of phenomena such as heterogeneous vascular delivery of drug and nutrient, diffusion through tissue, heterogeneous lesion growth, apoptosis, necrosis, and cellular uptake, efflux, and target binding. In this paper we present and discuss strategies to build multiscale computational models integrated with experimental models to study drug delivery via nanovectors to cancerous lesions and also discuss their predictive capability and limitations. These models describe events spanning the nanovector transport and drug release at the sub-cellular scale to cancer growth and drug response at the tumor scale. Model components are based on previous work and work in progress, particularly simulations of multi-dimensional tumor growth and chemotherapy that link cellular-level drug release and kinetics to the macro level effects. Here we propose how these components or “modules” can be associated into a higher order model, and describe how determination of parameters from experimental data would ensure correspondence with experimental phenomena, enabling model predictive capability. This multiscale approach represents an important step to further the understanding of cancer and facilitate the design of optimal therapies through nanotechnology.

1. Introduction

Therapeutic molecules are prevented from reaching their target in desired mass fractions by a multitude of barriers along their route from the points of administration. Among barriers are those of exquisite biological nature, such as insufficient specificity of target-ligand complexes, endothelial and epithelial barriers preventing extravasation in general (Jain 1988, Jain 1990,), and specifically the blood-brain-barrier; further immunological obstacles such as sequestration by the reticuloendothelial system (RES) (Kreuter, 1983, Kreuter, 1985), and sensitization phenomena. Potential biophysical barriers include non-specific, physical determinants of targeting; hemodynamics within the lesion microvasculature (Haroon 1999; Jain 1990; Jain 2001b); adverse osmotic pressure, which results in an outward net convective force, ejecting drugs from cancer lesions into the vascular compartment (Jain 2001b, Padera et al., 2004). Yet another barrier to therapeutic efficacy is posed by tumor tissue morphologic instabilities that arise in 3-D if therapeutic molecules, especially those with anti-angiogenic properties directed against cancer lesions. This may result in collective migration of tumor cell clusters and strands from the primary tumor, thus rendering therapeutic intervention substantially more complex (Kunkel et al., 2001, Pennacchietti et al., 2003, Lamszus et al., 2003, Bello et al., 2004, Rubenstein , 2000). These obstacles across different physical scales, from the tumor and blood vessel (10^-2 m) to the molecular (10^-9 m) scale, can be a major component of drug resistance beyond that offered by genetic mechanisms.

Recently, delivery of therapeutic agents via nanovectors such as nanoparticles and liposomes has become possible through the combined efforts of materials, nano, and engineering sciences applied to cancer medicine (Ferrari, 2005; La Van et al., 2003). These devices transport drug to targeted locations and then release it with more precise scheduling and dosage than conventional drug administration, thereby increasing therapeutic specificity and efficacy. Cancer cell recognition can be improved by conjugating nanovectors with ligands, and their physical and chemical properties can be tuned to optimize intravascular delivery and avoid biological barriers such as the reticuloendothelial system. Nanovectors can also transport contrast agents for monitoring in vivo tumor progression. Extensive experimental efforts are currently expended in this field, including surface modification to prolong circulation and ligand-particle conjugation to enhance selectivity (Bazile et al., 1995, Demoy et al., 1999, Gabizon et al., 2004, Illum & Davis, 1983, Illum & Davis, 1984, Matsumura et al., 2004, Moghini & Gray, 1997, Rudt & Muller, 1993, Sapra & Allen, 2003, Soppimath et al., 2001, Tröster et al., 1990, Tröster & Kreuter, 1992, Wang et al., 1995, Wilkins & Myers, 1966, Wilkins, 1967, Wilkins, 1967, Yoo & Park, 2004).

The design of these nanotechnology devices and methodologies, and the assessment of their laboratory and clinical performance have been largely accomplished with ad-hoc methods. A primarily empirical approach for optimizing nanovector therapy may not be sufficient to choose amongst a vast number of strategies or to obtain their full potential. In addition, it may be very challenging to overcome biobarriers that can seriously diminish therapeutic efficacy. Taking a cue from other fields, it has become apparent that only in conjunction with mathematical modeling based on fundamental principles can these designs succeed. To complement the extensive experimental research in this field, specific mathematical models ranging from nanovector performance to tumor growth and drug response have been proposed. The study of tumor biology via physical formulations has produced numerous models, too numerous to mention here, using continuum and cellular automaton approaches. Cancer growth, angiogenesis, metastasis, etc. have all been abstracted to a mathematical level, as reviewed by Araujo and McElwain (2004), Moreira and Deutsch (2002), and Mantzaris et al. (2004). A comprehensive list of mathematical models developed to study drug delivery and release at the nanoparticle level was recently compiled by Feng and Chien (2003). For instance, particle protection from RES uptake as a function of type and density of polyethylene glycol chains grafted to its surface has been predicted by Torchilin et al. (1994). Abnormal hemodynamics and vasculogenesis, which has a key role in nanovector extravasation, has been modeled by Baish et al. (1996), Anderson and Chaplain (1998), and McDougall et al. (2002). Deterministic and stochastic models of receptor-ligand binding due to Bell (1978) and Cozens-Roberts et al. (1990, 1990) have elucidated requirements for efficient particle-target binding and endocytosis.

A recent biocomputational implementation by Zheng et al. (2004) encompassed some of the main physical characteristics of tumor growth and implemented an in silico system that exhibited combined two-dimensional, non-symmetric tumor growth. Angiogenesis was incorporated through the models of Anderson & Chaplain (1998) as extended by McDougall et al. (2002). Sinek et al. (2004) studied nanoparticle mediated drug delivery and tumor response using the simulator of Zheng et al. (2004). Their multi-dimensional and multi-scale simulations demonstrated the potential increased efficacy of nanoparticle-based therapy and its limitations due to biobarriers. Wise et al. (in review) expanded this formulation to enable three-dimensional, non-symmetric tumor growth as well as multi-species representation.

An integrative approach is needed to combine experimental work and mathematical modeling. Here we propose a multidisciplinary effort towards the establishment of a multiscale theoretical model aimed at optimizing bio-barrier avoidance by injected therapeutic moieties, with and without nanoscale delivery vectors. The model encompasses molecular phenomena such as affinity complexation on sparse antigens, and single-file, non-Fickian molecular transport through nano-sized channels. It involves cellular-scale and extra-cellular matrix modeling components, such as architectural disruptions in angiogenic neovasculature, leading to enhanced permeation and retention (EPR) targeting, and ultra-structural, or cell distribution markers of disease. The proposed model further involves larger domains, relating the spatial distribution of targeted delivery, to the formation of complex tissue morphologies within the treated lesion.

In order to more completely capture the nonlinear coupling between the different physical scales of cancer, the multiscale model proposed herein incorporates the most advanced components describing cancer progression and therapy response based on fundamental physical and biological processes that influence nanovectored chemotherapy. The main model component is the Tumor Growth module, which reflects cancer progression from early to advanced stages and response to treatment. Clinically relevant tumors for the most part require Angiogenesis in order to obtain oxygen and nutrients. The intra-vascular and trans-vascular Nutrient Transport has an important role in tumor growth and response. Similarly, the vasculature affects intra-vascular and trans-vascular Nanovector Transport. Once in the vicinity of a tumor, Nanovector Cell Targeting and Specificity are key elements. Drug release mechanisms influence Pharmacokinetic and Pharmacodynamic considerations. Finally, the actual Genotype of tumor cells has important effects on the drug response.

Our fundamental hypothesis is that cancer treatment and response can be quantified via this all-encompassing predictive multiscale model with the long-term goal of optimizing variables for nanotherapy application to a specific patient’s disease. A rigorous quantitative analysis of drug action predicted in this manner will require experimental and clinical observations for model validation and refinement. Initially, values for model parameters originate from several sources. Information about patient’s tumors can be obtained from the clinic, such as biopsies for primary culture in vitro and dynamic contrast enhanced magnetic resonance imaging (DC MRI). Visualization of tumors in patients or through in vivo animal models can provide information about tumor growth, angiogenesis, drug pharmacokinetics, and vascular and, nanovector transport and interaction with the biological target. In vitro monolayer and spheroid models can provide tumor growth, drug pharmacodynamics/-kinetics, genotype information and nanovector adhesive performances. Nanotechnology can provide nanovector targeting and specificity parameter values. This “training phase” of the multiscale model is currently underway; once completed, prediction of treatment response will be obtained by setting groups of parameters to values assigned for specific conditions and running simulations of different therapy protocols. For example, if a patient evaluation indicates that an aggressive tumor strain with limited angiogenic potential is present, parameters for Tumor Growth, Angiogenesis, Nutrient Transport, Nanovector Transport, Pharmacokinetics, and Pharmacodynamics could be set to values that reflect this particular situation. Simulation of treatment through the multiscale model would then provide insight into the optimal protocol to maximize tumor regression for this patient. If actual response in vivo (or in vitro) does not match model prediction, then the model formulation is fine-tuned, i.e., the model “learns” from its results. This iterative methodology, mimicking the process of learning through cycles of research, development, and testing, is summarized in Figure 1.

Figure 1. Methodology for building a multiscale model of nanovectored therapeutics to cancerous lesions. Values of parameters for the model are initially obtained from patient tumor vascular imaging (DC MRI), in vitro models of tumor growth and drug response, in vivo animal models, and basic nanotechnology science. Model predictions are compared to actual outcomes, and model formulation is fine-tuned based on these results. This iterative approach ensures correspondence with in vitro and in vivo phenomena, and enables model predictive capability.

2. Model Components

2.1 Tumor Growth

A one-millimeter-radius tumor contains about 10⁶ cells, enabling at this scale the introduction of continuum models of cancer growth and progression. Nonlinear models, however, are required in order to capture complex tumor growth dynamics and morphologies. Most mathematical modeling of tumor growth has been carried out in cylindrically or spherically symmetric coordinates; the complexity of three-dimensional living tissue ultimately requires fully three-dimensional modeling. The tumor growth model decribed in detail here (Wise et al., in review) builds on previous work by Burton (1966), Greenspan (1972 and 1976), Maggelakis (1990), Byrne and Chaplain (1995 and 1996), Chaplain (1996), Friedman and Reitich (1999), Bellomo and Preziosi, 2000; Breward et al. (2002, 2003), Byrne and King (2003), Cristini et al. (2003) and Zheng et al. (2005). In particular, the latter presented a two-dimensional, finite-element/level-set (sharp-interface) simulator that incorporated the coupled neovascularization and nonlinear growth of solid tumors in order to examine the morphological evolution and stability of solid tumors in silico. The simulator was multiscale, based on a sophisticated unstructured adaptive mesh (Cristini et al., 2001; Anderson et al., 2005), and captured the nonlinear coupling between the growth and angiogenesis components. Wise et al. (in review) reformulated these well-known continuum-scale models, significantly extending the results of Cristini et al. (2003), where a boundary integral method had been used to follow the evolution of two-dimensional, nonnecrotic tumors. Cristini et al. determined that a set of two non-dimensional parameters (related to mitosis rate, apoptosis rate, cell mobility, and cell adhesion) could regulate tumor morphology and growth. However, their model was unable to capture features such as topological change, the onset of necrosis, or the transition from avascular to vascular growth.

A major difficulty with the sharp interface model of Zheng et al. (2005) was that it could incorporate multiple tumor species in a mass-conserving way only at the expense of introducing complicated conditions at the moving tumor boundary. The description of multiple cell species is important because the tumor microenvironment and impaired cell genetic mechanisms can select cells for survival under abnormal conditions (Graeber et al., 1996). Furthermore, the microenvironment of invasive tumors may be characterized by non-sharp boundaries between host and tumor tissues, and between multiple viable cell species (Liotta and Kohn, 2001; Kunkel et al., 2001; Lamszus et al., 2003). For example, it would be challenging to study how an invading tumor affects the host at the cellular scale with a sharp interface model. Wise et al. (in review) seeks to remedy these limitations by proposing a true multispecies numerical model of growth that accounts for non-sharp (diffuse) host/ tumor boundaries and processes such as cell mutation and necrosis.

Please et al. (1998 and 1999) applied multiphase modeling to tumor growth by representing both tumor cells and extracellular fluid as separate continuum phases. Ward and King and Breward et al. (2002) modeled avascular cancer growth as a two-phase description comprised of tumor tissue and dead tissue (extracellular space). This multiphase modeling is needed to capture avascular tumor growth as live tumor cells proliferate into the extracellular space. Breward et al. (2003) extended their avascular model to describe vascular tumor growth, thus incorporating a third phase to describe temporal and spatial distribution of blood vessels. Araujo and McElwain (2005) proposed a multiphase model of tumor growth that included a solid phase representing the extracellular matrix, in order to more accurately represent residual stresses. Along these lines, Wise et al. (in review) uses a multiphase modeling approach by capturing the evolution and interactions between intermixing multiple tumor species, necrotic tissue, host tissue, and interstitial fluid.

The model of Wise et al. (in review) includes a diffuse interface model and employs an energy formulation to accurately represent species intermixing at interfaces, such as that between tumor and host tissues, in order to more accurately capture the heterogeneity of the dynamic tumor environment and its interactions at multiple scales. Three tissue species are included: viable tumor (spatially varying density given by r_V), necrotic (r_N) and host tissue (r_H). For simplicity, a constant normalized density is enforced via r_V + +r_N + r_H = 1. Total tumor tissue is given as r_T = r_V + r_N. Principal quantities characterizing and affecting growth are normalized nutrient n (generically representing oxygen, glucose, and possibly growth factors), normalized pressure p, growth velocity u, and cell adhesion simulated by using an equivalent surface tension at the tumor/host interface (Chaplain, 1996). Although an actual tumor exhibits at least two separate and unequal pressures—one due to interstitial fluid, the other due to cell-cell contact (Jain, 2001, Jain, 1990, Padera et al., 2004)—for simplicity the model makes no distinction between the two.

As nutrient diffusion through tumor interstitium is a fast process compared to tumor growth, it is modeled by a steady-state diffusion-reaction equation with its source along the vasculature (Byrne & Chaplain, 1995, 1997). The tumor is assumed to be saturated with growth factors so that only nutrient availability limits cell proliferation. Therefore the fraction of cycling cells is given as identical to the normalized nutrient n. Regions where nutrient falls below a specified minimum n_N become necrotic, and the mass is assumed to be removed through lysing. Local viable tumor response S_V and necrotic tissue response S_N is given by the equations

(1)

, (2)

where l_M , l_A, l_N, and l_L are mitosis, apoptosis, necrosis, and lysis rates, H(·) is the Heavyside function, and Q(·) is an interface sharpening interpolation. These mass sources are in turn incorporated into the following generalized Cahn-Hilliard equations governing growth:

(3)

. (4)

These equations are analogues of standard reaction-diffusion equations, but with the additional capability of maintaining diffuse interfaces through a dimensionless adhesion energy given by

. (5)

In the above f(r_T ) = r_T ²(1 - r_T)²/4 is a bulk energy density with minima at r_T = 0 and 1 producing the tendency for the separation of tumor tissue (r_T = 1) from host tissue (r_T = 0). The tendency for interfacial mixing is given by , which exacts an increasing penalty as the interface becomes sharp. The influence of this effect is controlled by the coefficient e²/2. The density g enforces the inequalities 0 ≤ r_V ≤ r_T and

0 ≤ r_D ≤ r_T, and has the form

, (6)

where w_g is constant. In Eqs. 3 and 4 M is a positive mobility constant, and is the variational derivative of the adhesion energy E with respect to density r. Finally, since tumor cells and extracellular matrix are treated as comprising a viscous fluid flowing through porous media, the relation between pressure p and tumor cell velocity u is governed by the Darcy-Stokes equation:

,

where m is variable cell mobility, m_* is cell mobility just inside the tumor, g is surface tension, related to cellular adhesion, and k is local total curvature.

This growth component forms the core of the multiscale model proposed herein. It can accept inputs from other components that govern formation of new vasculature, mutations to more malignant cell species, standard chemotherapy, anti-angiogenic therapy, and nano-vectored therapy. Along with angiogenesis, tumor gowth is the most technically and computationally demanding process to render. An efficient nonlinear multigrid/finite-difference algorithm provides for the first time the capability to simulate fully three-dimensional tumor growth, as shown in Figure 2 (top), and with multiple cell species (Wise et al., in review; Frieboes et al., in review (a)). An adaptive mesh version, utilizing a method called Adaptive Full Approximation Scheme (Trottenberg et al., 2001), with the potential to further speed computation, is also in progress (Wise et al., in review).

Figure 2. Top: Simulation of three-dimensional tumor growth through the diffuse interface model of Wise et al. Adapted with permission from Wise et al., (In review). Left: Small tumor at beginning of simulation (time t=0). Right: By using a low value surface tension parameter to represent weak cell-cell adhesion, as is the case with aggressive glioblastomas, the tumor grows into an unstable mass that eventually breaks up (t=100 days). Bottom: Simulation of a two-dimensional capillary network showing migration of proliferating endothelial cells starting from the sprout tips (along y=0.0) towards a source of VEGF emanating from tumor cells (along y = 1.0). Reprinted with permission from Plank and Sleeman, Bull Math Biol 66, 1807 and 1801 (1998). Copyright © Elsevier. Left: The lattice-based random walk model of Anderson & Chaplain produces fractal-looking angiogenesis. Right: The non-lattice-based random walk model of Plank & Sleeman produces more interconnected vessels.

2.2 Angiogenesis

The modeling of abnormal vasculature (Haroon et al., 1999, Jain, 1990, Jain, 2001) and associated hemodynamics (Jain, 1988, Jain, 1990) generated by tumors as they grow and invade host tissue is important because drug delivery through nanovectors is typically done through the bloodstream. Although nanovectors can be assumed to preferentially extravasate from the relatively larger pores of tumor vessels (Hobbs et al., 1998, Yuan et al., 1995), increased tumoral fluid extravasation due to interstitial pressure (Jain, 1990) may present a major obstacle for consistent drug delivery. Previous modeling work has revealed additional insights into how the topology of a tumor’s vasculature affects blood circulation. For instance, irregularities in the vascular geometry could lead to a two-fold increase in vascular resistance relative to resistance measured in a uniform tube with the same mean diameter (Secomb & Hsu, 1996). Excessive vascular compliance and leakiness could cause blood flow to be diverted from the center of the tumor to its periphery (Baish et al., 1997). Although vasculature irregularity could be as detrimental to nanovector drug delivery systems as in the case of free drug administration (Sinek et al., 2004), nanovectors are different from drug molecules in that they could be designed and engineered to specifically overcome or reduce the influence of such biophysical barriers on drug delivery.

Over the past few decades, various models have been developed in order to capture the morphological structure of tumor vasculature. These models typically fall into one or more of three categories: 1) continuum models which describe variables such as endothelial cell density, nutrient concentrations, and growth regulator concentrations as continuous fields using differential equations; 2) mechanical-chemical models that involve both the effect of diffusible chemical species, such as angiogenic regulators, and of mechanical interactions between endothelial cells and the extracellular matrix; and 3) discrete, cellular-automaton-style models that describe the dynamics of individual elements, such as endothelial cells, based on a set of rules (Mantzaris et al. (2004). The latter provide a comprehensive, historical review of these angiogenesis modeling efforts.

In order to represent the physiology of tumor vasculature in a multidimensional simulated environment, we consider the use of hybrid continuum-discrete models of tumor angiogenesis (Anderson & Chaplain, 1998, Sleeman & Wallis, 2002, Plank & Sleeman, 2003, Plank & Sleeman, 2004). In particular, the lattice-based, random-walk model by Anderson and Chaplain (1998) and the non-lattice-based, random-walk model by Plank and Sleeman (2004) are representative of two different types of angiogenesis models capable of generating realistic vascular topology through physiological stimuli such as vascular endothelial growth factor (VEGF), angiopoietin-1 (Ang-1), and fibroblast growth factor-2 (FGF-2). For the sake of simplicity in these and other similar angiogenesis models, these tumor angiogenic factors (TAF) are combined into a single continuum variable. Whereas endothelial cell (EC) motion for models such as Anderson and Chaplain (1998) is typically determined by probabilities dependent on local TAF concentration and correspond to cells standing still or moving in a particular direction along the lattice grid, models like Plank and Sleeman (2004) describe EC migration by velocities, which also depend upon local TAF concentration. The primary difference between the two types of models is whether or not EC migration is geometrically constrained. Both types similarly model EC migration towards TAF via chemotaxis and along the ECM towards gradients of adhesive factors, such as fibronectin, via haptotaxis. In addition, both models assume a respective set of rules to capture branching and looping events (i.e., anastomoses) that are critical for the establishment of blood flow in a newly formed vasculature.

The outcomes of these models are capillary networks showing realistic dendritic structures (Figure 2, bottom) that reproduce the extensive branching observed experimentally just before the capillary network penetrates a tumor (Less et al., 1991, Skinner et al., 1990). Mathematically, the continuum formulation by Anderson and Chaplain (1998) models endothelial cell density e, TAF, c, and fibronectin f as follows:

(7)

In the first equation, the first term on the right hand side represents relatively weak cell diffusion with diffusivity D_e, while the second and third terms represent chemotaxis up the TAF gradient and haptotaxis up the fibronectin gradient, respectively. a_c and a_f can be constant, although it is more realistic to have a_c be a decreasing function of TAF (Anderson & Chaplain, 1998). In the second equation, the first term on the right represents production of fibronectin by endothelial cells, while the second term represents fibronectin uptake, with n_f and h_f being constant. In the last equation, the right term represents uptake of TAF with constant rate h_c. Seeding several small regions of high density “sprouts” along a parental vessel sets initial conditions for endothelial cell density. The parental vessel and perinecrotic rim just within a tumor are assumed to produce initial fibronectin and TAF, respectively. Concentrations are thus assumed to decay with distance from their sources.

2.3 Intra-vascular Flow

Acute regulation of local blood flow through vasculature is achieved largely by contraction or relaxation of smooth muscle in vessel walls. Accommodation to chronic changes involves long-term structural adaptation of vessel diameters and wall thicknesses, in which each vessel segment responds to local mechanical and biochemical stimuli (Pries et al., 1998)). Many studies have shown that vessels adapt specifically to mechanical forces exerted by flowing blood and transmural pressure (Pries et al., 1998, 2000, 2001, 2005, Zakrzewicz et al., 2002). Changes in tissue metabolism, blood flow and intravascular pressure induce stimuli in the form of shear stress, , circumferential stress, , and biochemical signals, resulting in vessel structural adaptations (i.e., changes in diameter and wall thickness) (Zakrzewicz et al., 2002, Pries et al., 2005), that cause heterogeneity in vessel diameter and wall thickness throughout in vivo vasculature beds.

The distribution of oxygen, nutrients, and drug in the tumor microenvironment that govern tumor growth are critically dependent on both the vasculature morphology and blood flowing through it. Specifically, blood flow in vascular networks is highly sensitive to changes in vessel diameter and wall thickness, sincewhere is flow resistance and is luminal diameter (Pries et al., 1998). Heterogeneity in the tumor microenvironment is thought to drive invasive tumor morphologies through hypoxia and nutrient deprivation (Young et al., 1988, 1990, Cairns et al., 2001, Posotovit et al., 2002, Rofstad et al., 2002, Montesana et al., 1991). The abnormality of tumor vasculature due to poor organization, high tortuousity, hyperpermeability, and high leakiness (Jain et al., 2001), is widely believed to cause heterogeneity in the tumor microenvironment. In order to more accurately capture the in vivo nature of tumor vasculature and its influence on the microenvironment, a model must be able to capture the dynamic morphology and hemodynamics of the vasculature bed as well as the transfer of substances from blood into the lesion tissue.

The works of McDougall et al. (2002) and Stephanou et al. (2005) are among the first to consider including blood flow descriptions with a mathematical model of angiogenesis, specifically Anderson and Chaplain (1998). Their follow-up work (Stephanou et al., 2006, McDougall et al., in press) in addition to Alarcon et al. (2003) represent the leading theoretical work in applying models that capture vascular remodeling in addition to blood flow to a vascular description, which is a necessary to derive information regarding the mechanical stresses that primarily drive vascular remodeling. McDougall et al. (2002) extended Anderson and Chaplain’s angiogenesis model to include a blood flow description by modeling the generated vascular networks as a series of straight, rigid, cylindrical capillaries joined at adjacent nodes. Flow was modeled through each cylindrical element by Poiseuille’s Law describing flow-rate as a function of capillary lumen radius, fluid viscosity, capillary length, and pressure drop. A radius R_ij is randomly assigned from a probability distribution to each element (vessel segment) joining nodes i and j on the grid. At each node i there exists a pressure P_i, and through each element joining nodes i and j there exists a flux Q_ij. This flux is assumed to obey Poiseuille’s law with the segment’s inner radius randomly assigned,

(8)

where m is fluid viscosity and L_ij is the length of element ij. Conservation of mass at each node i requires that

, (9)

where j varies over the four adjacent lattice nodes. Given the pressure drop across the parental vessel, which feeds the capillaries going into the tumor, this approach results in an exactly determined system of linear equations. Based on this modified vasculature model, McDougall et al. (2002) identified tumor vasculature as a biobarrier affecting chemotherapy efficacy. Their simulations suggested that the highly interconnected nature of tumor vasculature can reduce the amount of drug delivered to the tumor, and in some cases, completely bypass the entire tumor mass.

Stephanou et al. (2005) furthered the work of McDougall et al. (2002) by developing an algorithm that normalizes simulated vasculature produced by Anderson and Chaplain’s angiogenesis model. They examined how pruning vessels by anti-angiogenic drugs might affect blood flow distribution, and consequently drug delivery to tumors. Stephanou and coworkers recently advanced their prior work to include vascular adaptation effects (Pries et al., 1998, 2000, 2001, 2005, Zakrzewicz et al., 2002) to the angiogenesis model; simulations provided insight into the effect of vascular remodeling on oxygen and drug supply to tumors (Stephanou et al., 2006). Recent work by McDougall et al. (in press) extended Stephanou et al. (2006) to more realistically capture in vivo angiogenesis by simultaneously coupling vascular adaptation with growing vasculature beds and blood flow. The primary advantage of this extension is that vasculature can be adapted dynamically as the individual vessels grow rather than adapt the vasculature after it has been fully developed.

2.4 Trans-vascular and Interstitial Flow

A combination of irregular vasculature form and function, high interstitial fluid pressure, high cell-cell pressure leading to collapse of tumor vessels, and interstitial drug-binding result in highly heterogeneous nutrient and drug distribution (Jain, 1988; Jain, 1990; Jain, 2001; Padera et al., 2004). Tumors thus develop diffusion gradients of oxygen and nutrient in vivo, as has been observed with many cancer cell types in vitro (Acker et al., 1984, Acker et al., 1987, Carlsson & Acker, 1988, Casciari et al., 1988, Franko & Sutherland, 1979, Mueller-Klieser & Sutherland, 1982, Mueller-Klieser, 1987, Sutherland & Durand, 1976, Teutsch et al., 1995). In vitro spheroid models evince cell proliferation gradients, with outer cells having highest mitotic activity and cells near necrotic regions being mostly quiescent (Bredel-Geissler et al., 1992, Carlsson, 1977, Freyer & Sutherland, 1980, Kunz-Schughart, 1999, Wartenberg et al., 1998). Not surprisingly then, studies have suggested that tumor morphology may play a key role in determining drug response (e.g,, Kobayashi et al., 1993). Since rapidly cycling cells have been observed to be more sensitive to anti-proliferative agents like doxorubicin (Tannock, 1994; DeGregorio et al., 1984), those exposed to minimal nutrient could present a more resistant species. Furthermore, drug administration under such conditions may influence gene expression to enhance drug resistance (Sutherland, 1988, Jain, 2001b) and contribute to morphologies that increase invasiveness (Cristini et al., 2005, Frieboes et al., 2006). Morphological barriers were recently studied using a computational model of drug response with parameter values determined from in vitro cytotoxicity data (Frieboes et al., in prep.). Simulation results supported the hypothesis that the aforementioned mechanisms and phenomena do, indeed, significantly contribute to reduced drug efficacy. It is apparent that drug released from nanoparticles and liposomes after extravasation along tumoral vasculature would likewise be affected by tumor morphology.

The net local production S_n of nutrient can be modeled as (Sinek et al., 2004)

, (10)

where n_n is a vasculature transfer function controlling extravasation of nutrient n (normalized with respect to the concentration in the blood), d is a dimensionless line delta function supported at the vasculature, and h_n is the rate of nutrient uptake by tumor cells. The function n_n can take many forms, but usually includes information regarding local normalized nutrient and pressure p. The form , where is a transfer coefficient, has been successfully employed in our simulations. The notation (·)⁺ means max{·, 0} and is used for the blood-to-tumor pressure difference term because if it becomes negative the blood vessel will collapse (Padera et al., 2004). Transport through tumor tissue is modeled by the following quasi-steady state equation (Sinek et al., 2004), since the time scale of nutrient diffusion is much smaller than that of cell proliferation:

(11)

where D_n is nutrient diffusivity through the lesion, assumed constant. A typical contour plot showing a nutrient distribution generated from the above equations is given in Figure 3.

Figure 3. Two-dimensional simulation of nutrient gradients in the tumor microenvironment due to irregular tumor morphology. Left: Tumor is shown by dark lines with its capillary network after simulated chemotherapeutic treatment via nanoparticles. The mass has fragmented as a result of treatment. Right: Nutrient distribution within the lesion. Note strong gradients resulting from vascular heterogeneity and interstitial pressure. Adapted from Biomedical Microdevices, Vol. 6, 2004, p. 307, Sinek et al., Figure 5, © 2004 Kluwer Academic Publishers. With kind permission of Springer Science and Business Media.

2.5 Intra-vascular and Trans-vascular Nanovector Transport

In Section 2.3, the intra-vascular transport of nutrients and molecules was analysed employing a classical one-dimensional (in space) model. In this section a higher order model is considered to analyse the transport of nanovectors within a capillary flow. This leads to define an effective longitudinal diffusivity, which can be used to derive accurate results from lower order models. Should review modelling work in the field and put your into perspective.

Assuming that the nanovector is sufficiently small not to perturb the flow field within the capillary, the local nanovector concentration C(r,z;t), function of the axial distance z along the capillary, the radial distance r (Figure 3), and time t, is governed by a classical convection-diffusion equation as

(12)

where w is the local fluid velocity, and Pe=R_ew_max/D is the Peclet number, with w_max being the maximum fluid velocity on the capillary axis, R_e the capillary radius, and D the molecular diffusivity of the injected nanovector.

By solving the advection-diffusion equation and integrating the local concentration C over the cross-section, the variation of the mean concentration C_m of nanovectors along the channel is derived at different time steps (Figure 4, middle left), or conversely the variation of the mean concentration as a function of time at different cross sections z (Figure 4, middle right). In the case of an impermeable channel C_m can be described as a Gaussian like curve moving downstream and spreading along the channel. How does this solution depend on the velocity profile assumed?

Figure 4. This is wrong: there is no Pouissile flow in capillaries. The velocity profile for a laminar flow in a cylindrical tube - Pouissile flow. Middle Left & Center: The mean concentration distribution along the capillary at different times (Left) and its variation with time at different locations (Center). Parameter values are: capillary radius R_e=10^-5 m, mean velocity at the entrance of the capillary w₀=10^-4 m/sec, D=10^-11 m²/sec. Middle Right: The effective diffusion coefficient D_eff as a function of the nanovector radius a for different values of the parameter R_eU. Bottom Left: The effective diffusion along the tube as a function of the permeability parameter Π. Bottom Right: The mean concentration profile, against z, at different Π.

More interesting then the variation of the concentration C or of the mean concentration C_m is the effective longitudinal diffusivity of the nanovectors along the channel, that is to say how fast a nanovector can move and be transported along a capillary. In the absence of a fluid flow, the dispersion of a passive species is governed by the classical molecular (Brownian) diffusion and for a spherical nanovector of radius a the classical relation of Einstein-Stokes can be used to estimate the diffusion coefficient D in an infinite medium being

D=(k_BT)/(6πηa) (13)

where η is the viscosity of the medium and k_BT the Boltzmann thermal energy. In the presence of a fluid flow, the dynamics of the nanovector is also governed by advection due to the non-uniform velocity field across the capillary. This enhances significantly the longitudinal diffusivity depending on the nanovector size and flow conditions. Already in 1953, Taylor introduced for the transport of a passive species in a circular straight channel an apparent diffusion coefficient D_app, comprising the sole convective contribution, given as

D_app=R_e²U²/(48D) (14)

where D is the molecular diffusion given in (13), and U the fluid mean velocity. Subsequently, Aris (1954) derived the longitudinal effective diffusion coefficient D_eff, considering also the contribution of the molecular diffusion D, leading to

D_eff=D+R_e²U²/(48D) (15)

This seems trivial superposition of 13 and 14. can you explain? Also, what is the definition of effective diffusivity? It can be plugged into a diffusion equation?

Notice that the second term on the right hand side of (15), i.e. the convective contribution to the longitudinal diffusion, is proportional to the second power of the mean fluid velocity and inversely proportional to the molecular diffusivity D. Also, D_eff coincides with D in the case of small molecules, as proteins and nutrients, which are nanometric or even sub-nanometric in size (< 1-2 nm). The Taylor and Aris formula holds in the limit of sufficiently large times after injection (> R_e²/D), or in other words, in the limit of sufficiently long channels (> UR_e²/D), being derived from the asymptotic solution of the general advection-diffusion equation (12). Therefore, it is just important to observe that assuming a capillary radius R_e of the order of 10 mm (this looks more like a microvessel than a capillary) with a mean fluid velocity U of about 100 μm/s (Ganong, 2003), and a brownian diffusion coefficient D of the particle ranging between 10^-11 - 10^-9 m²/s, the Taylor's regime would be fully developed after a critical distance ranging between 10 μm and 1 mm, which is relatively small compared to the usual length of capillaries (really? Can you support this?). Such a critical distance would be even smaller in the case of tumoral capillaries where the mean fluid velocity can be up to two order of magnitude smaller than that of normal capillaries (1 – 10 μm/sec).

In the case of impermeable blood vessels, substituting (13) in (15), it follows

D_eff=(k_BT)/(6πηa) +R_e²U²(6πηa)/(48k_BT) (16)

This relation shows that the effective diffusivity has a biphasic (non-monotonic) behavior with respect to the radius a of the nanovector: the molecular diffusion contribution decreases as a increases and it becomes less and less important for larger and larger nanovectors; the convective contribution increases as a increases and becomes less and less important for smaller and smaller nanovectors. Considering a spherical nanovector in plasma (η=1.8×10^-3 Pa s), at ambient temperature, the effective diffusion coefficient D_eff is plotted in Figure 4 (middle right) as a function of the radius a of the nanovector for different values of the parameter R_eU, ranging between 10^-12 and 10^-10 m²/s. The values chosen for R_eU are physiologically relevant in human capillaries. For a fixed R_eU curve (Figure 4, middle right), the effective diffusivity D_eff decreases starting from small a, reaches a minimum and then increases steadily with a. The minimum in this plot indicates the so-called critical radius a_cr which identifies the nanovector with the minimum effective diffusivity within a capillary of radius R_e and with a mean velocity U. In other words, given the “R_eU capillary”, the nanovector with the corresponding critical radius has the smallest longitudinal diffusivity compared to any other larger or smaller nanovector. Notice that the critical radius a_cr increases as the parameter R_eU and the viscosity η of the fluid reduce.

In Decuzzi et al. (2006), the effective diffusivity and the critical radius a_cr has been derived also in the case of permeable channels. Notice that tumor capillaries are much more permeable than normal capillaries [], and the vascular fluid pressure can not be anymore considered as linear along the channel, as well as the mean velocity U which reduces moving downstream the channel. As a consequence, the effective diffusivity D_eff varies along the channel in a fairly complicated manner as shown in Figure 4 (bottom left) for a nanovector with a radius of 200 nm and for different values of the permeability parameter P. This parameter Π, introduced in Decuzzi et al. (2006), has the form

, (17)

where l is the length of the capillary, w₀ the mean velocity at the entrance of the capillary, L_p is the vascular hydraulic conductivity and p_i is the interstitial pressure. Notice that as the channel permeability increases, the effective diffusion coefficient reduces and in the limit tends to approach the pure molecular diffusion value (dashed horizontal line).

Interestingly, for sufficiently large, but still physiologically relevant Π, the effective diffusivity in the central zone of the capillary drops to the pure molecular diffusion value. This result has to be ascribed to the existence of zones within the channel where, due to the blood vessel permeability, the fluid velocity is null and this consequently leads to a negligible or even zero convective contribution to the longitudinal transport. A critical radius can be defined even in the case of the permeable channels and it is larger than that derived for a non-permeable channel, due to the smaller convective contribution. In particular Decuzzi et al. (2006) have shown that the critical radius increases almost linearly with the vascular hydraulic conductivity L_p.

These results have several implications for the delivery and transport of intra-vascularly injected nanovectors and molecules. The circulatory system is made up of conduits of various size and mechanical characteristics carrying blood from the heart to the tissues and back to the heart again. Moving from the aorta to the capillaries the product between the mean flow velocity and the vessel radius (R_eU) reduces from about 10^-2 m²/s to 10^-12 m²/s. Consequently, convective diffusion dominates in large vessels as aorta, arteries and veins; whereas as the vessels become smaller, the critical radius grows ranging between about 1 nm and 100 nm for impermeable capillaries. Based on this, it can be concluded that:

(i) drug molecules, monoclonal antibodies, nano-spheres, dendrimers with characteristic sizes up to 1-5 nm are transported within the capillary network mostly by molecular diffusion;

(ii) nanovectors with radii whose values range between 10 nm and 100 nm, as liposomes and polymeric particles, would diffuse less effectively within the capillary network and would tend to concentrate only at the entrance of the capillary;

(iii) nanovectors with radii whose value is larger than 100 nm (super critical radius a>a_cr) are transported within the capillary network mainly by convective diffusion and would tend more easily to move downstream along the capillary, even in the case of capillaries with small mean fluid velocities.

To conclude, the variation of C_m along the channel at a fixed time step for different values of the wall permeability is shown in Figure 4 (bottom right): as the wall permeability increases, the C_m curve becomes steeper and moves more slowly downstream.

Can you substantiate finding from experimental data? Perhaps theory not valid for capillary flow?

Need to explain how this can be used as input to blood-to-tissue tranfer parameters in the multiscale code.

This preliminary results support the idea that most of the free molecules or nanovectors entering a permeable capillary would tend to leave the vessel (extravasate) right after the capillary entrance leading to a highly non-uniform drug distribution within the extra-vascular space. How is this connected to the tissue-scale modules?

2.6 Nanovector-Cell Targeting and Specificity

A broad spectrum of nanovector types have been presented in the literature (LaVan et al., 2003; Ferrari, 2005) with different composition and chemico-physical properties which can be used as drug carriers (drug delivery systems) and as contrast agents in medical imaging (imaging systems). Examples are nanospheres (Duncan, 2003) where the pay-load (drug molecules or contrast dies) is dispersed within a polymer matrix; multilayered nano/microcapsules and liposomes (Crommelin and Schreier, 1994) where the pay-load is contained in the internal capsule; nanoporous Si particles (Cohen et al., 2003) where the pay-load binds to the pores surface.

A systemically administered nanovector before reaching its final target and execute its mission has to make its way into the circulatory system and reach the diseased microvasculature. In the case of a “vascular targeting strategy”, the nanovector has to adhere firmly at the tumoral endothelial cells and from there release the therapeutic agent or pay load towards the extravascular space; whereas in the case of a “tumoral microenvironment targeting”, the nanovector has to extravasate crossing the fenestrations in the blood vessels, diffuse through the extracellular matrix (ECM) and eventually bind to the target tumor cell.

2.6.1. Non-Specific Interactions and Initial Nanovector-Biological target approach

How is this related to Section 2.5? The motion of the nanovectors towards the target surface (target cell) is governed by diffusion and convection, as seen in the previous paragraph, and by non-specific forces as gravitation, van der Waals, electrostatic, and steric interactions, as described in Decuzzi et al. (2005). Therefore the “recognition” and approach of the nanovector to the target cell is influenced by the chemicophysical properties of the nanovector, of the medium and of the biological target.

Considering the most vastly treated case in the literature of two parallel planes facing each other in a liquid medium (Israelachvili, 1992):

- the van der Waals interaction energy per unit area is given by

Wvdw (δ)= A/(12πδ²)

(18)

where δ is the separation distance between the two planar surfaces and A is the Hamaker's constant depending on the dielectric properties of the interacting bodies (nanovector and cell) and of the medium;

- the electrostatic interaction energy per unit area at relatively low potential between two similar surfaces is given by

Wel (δ)=(64/κ) kBT ρ∞ ev ec Exp(κδ), ei = tanh[(zv eψi)/(4kBT)];

(19)

where ρ_∞ is the number density of ions in the medium (m^-3), ψ_v and ψ_c are the electrostatic potential at the nanovector (subscript v) and cell (subscript c) surface respectively, e is the elementary charge of an electron (e=1.602×10^-19 C), κ is the inverse of the Debye length (m^-¹), and z_v is the ion valence;

- the steric interaction energy per unit area between surfaces containing adsorbed polymer layers at low surface coverage is given by

Ws=36ΓkBTexp[δ/Rg]

(20)

where R_g is the unperturbed radius of gyration of the adsorbed polymer, Γ is the number of polymer chains per unit area.

Such physical forces are strongly influenced by the geometry of the interacting bodies. A simple formula exists in the case of two spheres, or in the limit of a sphere interacting with a a planar interface, known as the Derjaguin approximation, which relates the interaction force F to the interaction energies defined above depending on the geometry of the system. For a sphere of radius a (nanovector) interacting with a planar substrate (cellular target much smaller than the nanovector), the Derjaguin approximation for the force simply reads as

Fsphere(δ)=2πaWplane(δ)

(21)

When adjusted appropriately, the additive effects of these forces result in an energy-separation curve with a minimum at the equilibrium separation distance. Such adjustments are typically achieved through (i) grafting polyethylene glycol (PEG) or polyethylene oxide (PEO) chains of appropriate density Γ and length R_g to the particle surface; (ii) changing the surface charge of the vector and its bulk material chemical composition; and (iii) changing size and shape. A singular difficulty in the targeted delivery of nanovectors is the interaction with biological substrates other than the final target (tumor cell or tumor endothelium), especially macrophages and specialized cells lining the liver, spleen, bone marrow, and lymphatic tissue (collectively known as the reticuloendothelial system), which rapidly sequester and remove nanovectors from circulation by binding certain proteins to them in a process called opsonization. The challenge in developing effective nanovectors is the selection of optimal surface properties to avoid this biological barrier, yet specific enough to recognize and eventullay bind to the target cells.

2.6.2 Specific Interactions and Firm Nanovector-Biological target adhesion

When the nanovector is in close proximity with the biological substrate (target cell), molecular-specific reactions with the target cell become possible leading eventually to a firm adhesion of the nanovector. In fact the firm adhesion is regulated and mediated mainly by the selective binding between molecules expressed over the target surface (receptors, R) and counter-molecules (ligands, L) conjugated or grafted at the nanovector surface. The specific ligand/receptor interaction is generally described as a reversible bimolecular reaction in which the molecule L comes together to the molecule R to form the non-covalent bond LR, that is

L+R « LR

(22)

with a forward reaction rate k_f (association rate) and a reverse reaction rate k_r (dissociation rate). The ratio K_A=k_f/k_r is the binding affinity constant, and it is well known since the work of Bell (1978) to depend on the force F applied to pull the ligands away from the receptors and tending to open the bond LR (Figure 5, top left). Bell introduced the phenomenological law k_r=k_r0 Exp[cF/(k_BT)], where c is a characteristic length scale of the LR bond and k_r0 is the zero load reverse rate.

Figure 5. Top Left: The interface between a nanovector and the biological target Top Right: y-axis is the right hand side of Eq. 16 and represents rate of bond formation under various separation forces. In all cases, when above the x axis, the number of bonds moves to the right along the curve; otherwise, it moves to the left. Forces represented from top to bottom are none, a subcritical force F_sub, the critical force F_c, and a supercritical force F_sup. Bottom: The probability p_i^st of forming i bonds for different values of the applied force F.

The most elementary formulations of the ligand/receptor bond formation assume one-to-one binding with a surface density of bound receptor/ligand complexes much greater than the density of total available receptors. Under such conditions, bond formation between ligands and receptors can be modeled by the differential equation (Bell, 1978):

(23)

where N_o is the total number of bound receptor/ligand pairs at time t, N_r and N_l are the total number of available receptors and ligands, respectively, and k_f and k_r are the rates of bond formation and distruption. The force F tending to open the LR bond can be the resultant of the non-specific interactions cited above in the case of a tumoral microenvironment targeting or is mainly related to the hydrodynamic force in the case of the vascular targeting, being generally the hydrodynamic force much larger than the non-specific interactions. Under conditions of no force a linear relationship between the rate of bond formation and the number of bonds obtains (see Figure 5, top right). Here, the number of bonds moves to the right when above the x axis, and to the left when below, usually seeking equilibrium at the x intercept. However, when F > 0, the exponential produces a qualitatively different dynamics. Past a certain critical force F_c the rate of bond formation is always negative, and no equilibrium can be reached, i.e., the nanovector-cell complex will separate.

A large variety of receptors are expressed over the cell membranes with different functions, surface density and distribution (Alberts et al., 2002). For instance on leukocytes more than 250 protein species have been identified, half of which are thougth to act as adhesion molecules (Barclay, 1998). Several factors contribute to the strength of the nanovector/cell adhesion including the surface density of the receptors m_r and ligands m_l; the area of the surface of interaction A_c between the cell and the opposing substrate, the binding affinity K_A between ligands and receptors; and the external applied force F. Understanding which of these key parameters is more effective in controlling the adhesive strength in a fairly general case is of great importance. Mathematical models for specific problems are availaible in the open literature as for the peeling of membranes (Evans, 1985, Bell et al. 1984), and for the rolling and adhesion of spherical beads (Hammer and Apte, 1992). Despite this there is still no a general and clear picture of which of the above parameters have the largest effect on the adhesive strength. In addition to this, there is a growing bulk of evidences showing that cell-surface adhesion is rather mediated by clusters of ligand/receptor pairs with few adhesive molecules acting cooperatively and within which the applied external load is shared equally. This has been observed for the case of Selectins, adhesive molecules which tend to be localized in a larger number at the tip of the microvillus in circulating lymphocytes whilst are only minimally present on the flat portions of the cell membrane favouring tethering, rolling and final adhesion of the blood cell to the vascular wall; and in the case of Integrins, adhesive molecules whose clustering favours cell adhesion and migration across the extracellular matrix.

Differently from a one-to-one binding, within a cluster at any given time t each of the bond can be either in an open or closed state and rebinding of initially open bonds is allowed as well as opening of initially closed bonds (Figure 5, top left). With such a scenario, the adhesive interaction between a cell and a nanovector has to be treated as a dynamic and stochastic process: a cell that was initially firmly adherent to a surface could be fully detached after some time, or conversely a cell in proximity of a biomietic surface could firmly adhere given a certain time.

Thus the classical deterministic approach, where the number of ligand/receptor bonds is known exactly at a certain time t, is substituted by a stochastic approach where the probability p_i of having i ligand/receptor bonds at time t is considered. The master equation for the kinetics of small systems (McQuarrie, 1963) is used, which describes the rate of change in probability p_i(t) as given by

(24)

where i ranges from 0 to N_o, being N_o the minimum between the absolute number of ligands N_L (=A_cm_l) or receptors N_R (=A_cm_r). The probabilities p_i are zero for i values outside this range. Therefore, the relation (24) gives rise to a system of N_o linear differential equations with the conservative condition

at each time t.

(25)

The following general constitutive equation (Piper et al. 1998) has been considered for the reverse and reactio

(26)

where the constants χ, b, c and d depend on the ligand/receptor pair, and k_r and k_f are the reaction rates at zero load. The bond length parameter χ can vary significantly depending on the ligand/receptor pair and also on the applied force (Evans et al. 2001). It is assumed to be χ=1 Å. The parameters b, c and d can be derived experimentally using the constitutive relation (26). It is assumed b, c and d to be unity. Notice that for b=1, c=0 the classical Bell (1978) model is recovered with a constant forward reaction rate k_f=k_f0. The values for k_r0 and k_f0 are also dramatically affected by the ligand/receptor pair and the chemical properties of the surrounding fluid, for instance the pH.

The steady state solution is independent of the initial conditions and is readily derived from (24). The steady state prrobability p_i^st is plotted in Figure 5 (bottom) for different applied loads F, ranging between 100 and 200 pN.

For a relatively small force (=100 pN), the probability of having N_o (=10) bonds is 0.614 and it decreases as i reduces being smaller than 10^-3 already for i=6. As the applied force F increases, the probability of having a number of bonds smaller than N_o grows so that at a relatively large force (=200 pN), the probability of having N_o bonds is almost 0, whereas p_o=1.

More important than the probability p_i of having i active bonds is the adhesion probability P_a defined as the probability of having at least one bond active. Therefore the adhesion probability has the mathematical form P_a=1-p_o and can be considered as a measure of the strength of the nanovector-cell adhesion. The variation of the adhesion probability P_a with the applied force F is given in Figure 6 (top left). It is clearly shown that P_a decreases drammatically passing from unity to zero in a narrow interval of forces. This behavior leads to define a critical force F_cr above (below) which the probability of adhesion P_a is smaller (larger) than 0.5. In the present case the 50% probability of adhesion occurs for F»160 pN which is then the critical force: for F<F_cr nanovector-cell adhesion (detachment) is more likely than detachment (adhesion).

Figure 6. Top Left: The probability of adhesion P_a (=1-p_o) as a function of the applied load F. The dashed line is for F_cr, corresponding to a P_a=0.5. Top Right: The probability of adhesion P_a (=1-p_o) as a function of the applied load F and for differente values of the interaction area A_c. Bottom Left: The probability of adhesion P_a (=1-p_o) as a function of the applied load F and for differente values of the surface density m_l of ligands on the nanovector surface. Bottom Right:The probability of adhesion P_a (=1-p_o) as a function of the applied load F and for differente values of the reverse reaction rate k_r of the ligand/receptor pair.

Based on the formulas presented above, it is readily observed that the probability of adhesion P_a is affected by key governing parameters as the interaction area A_c, which is related directly to the characteristi size of the nanovector; the surface density of ligands m_l over the nanovector surface; and the reverse reaction rate k_r, which is a measure of the strength of the ligand/receptor bond. The surface density of the receptors can be modulated only up to a certain extent, whereas the forward reaction rates are much less important than k_r that dominated detachment.

The effect of the interaction area A_c is shown in Figure 6 (top right), where the adhesion probability P_a is plotted as a function of the applied force F for different values of A_c ranging between 0.5 and 2.0 μm². As the interaction area increases, at a given density m_l, the absolute number of ligands and receptors available for binding increases favouring attachment thus leading to larger critical forces. In addition, it can be observed for the data considered a linear increase of the probability of adhesion with interaction area A_c, for sufficiently large A_c; whereas at small A_c the linearity is lost (results not shown). The effect of the surface ligand density m_l is shown in Figure 6 (bottom left), where the adhesion probability P_a is plotted as a function of the applied force F for m_l ranging between 0.5 and 2.0×10¹³ sites/μm².

As the surface ligand density increases, the adhesion probability increases and so does the critical force. This is due to the larger absolute number of ligands and receptors available for binding in the fixed interaction area A_c. However, the increase of P_a with m_l is not linear as tend to saturates for large m_l. Finally the effect of the reverse reaction rate k_r is shown in Figure 6 (bottom right), where the adhesion probability P_a is plotted as a function of the applied force F for k_r ranging between 0.5 and 2.0×10^-4 s^-1. As the reverse reaction rate increases, that is to say as the strenght of the ligand/receptor pair decreases, the adhesion probability decreases and so does the critical force. Similarly as for m_l, the decrease of P_a with k_r is not linear as tend to saturates for large k_r.

From these results (Decuzzi et al., in preparation), it can be concluded that the surface of interaction A_c has the largest influence on the strength of adhesion. Thus to improve the adhesive strength between the nanovector and the target cell is of fundamental importance to increase the area of interaction. The contribution of the density of ligands and the binding avidity among receptors and ligands to the adhesive strength is generally much smaller. This is in agreement with several experimental evidences and theoretical predictions, and it is also interesting to notice that leukocytes employ this same strategy to adhere firmly to the endothelial walls.

2.7 Drug Release Mechanisms

2.7.1 Nanoparticles and Liposomes

Whether nanoparticles and liposomes merely extravasate at the lesion site or bind to target cells and become endocytosed, their payload must eventually be released. The physics of liposomal and nanoparticle drug release has undergone extensive study, with the Higuchi, power law, and Weibull models sometimes used as phenomenological approximations (Sinek et al., in press). The Weibull model, which is a single exponential asymptotically approaching 100% release in time, can be interpreted as a mechanistic model (Kosmidis et al., 2003). Whereas liposomes generally release drug by lysing, nanoparticles do so through surface or bulk erosion (Langer & Peppas, 1983, Zhang et al., 2003). For polymeric particles, an initial rapid release usually occurs due to drug adsorbed on the surface, followed by a more sustained release of incorporated drug (Soppimath et al., 2001, Kreuter, 1994, Illum et al., 1986). This release, in turn, is due to three primary mechanisms: dissolution of drug from the solid phase, diffusion of dissolved drug through the matrix, and erosion of the polymer matrix itsself (Feng & Chien, 2003). Often the Higuchi (1961), power law (Peppas, 1985), and Weibull (1951) models provide adequate phenomenological approximations of these processes. More recent models and experimental results are discussed by Siepmann and Goepferich (2001), Kreuter (1994), and recently reviewed by Feng and Chien (2003) and Frieboes et al. (2006). Nanoparticle release profiles usually show a simple bi-exponential release pattern described by

, (27)

where C_t is the amount of released drug at time t, C_¥ is the total amount of drug in the nanovector, C₁ and a are parameters corresponding to the rapidly released portion, and C₂ and b are parameters corresponding to the slowly released portion (Sinek et al., 2004, Kreuter, 1994, Feng & Chien, 2003). If the release is sustained long enough, the rapid release term becomes negligible, leading to , or a constant release at a rate of C₂b. Depending upon nanoparticle design, release can approach 0^th order for an appreciable time, a desirable trait since this would expose lesion to fairly uniform levels of drug over long periods (Feng & Chien, 2003). Even with this simplification, cellular level drug kinetics and transport is highly non-uniform because of inhomogeneous transport of vectors through and extravasation from tumoral vasculature and because of drug gradients due to cellular uptake and metabolism. The simplified case was assumed in a study by Sinek et al. (2004), where it was shown that, while therapy was improved in terms of mass regression, heterogeneous drug effect contributed to incomplete regression and fragmentation.

2.7.2 Implanted Microfabricated Devices

The predominant means of drug administration are subcutaneous injection and oral delivery, which do not allow precise control of drug delivery rate or target area, resulting in non-homogeneous release that is fast in the beginning, possibly reaching toxic levels, and lower at the end of the treatment period, possibly at ineffective concentrations. Newer methods for controlled drug release range from transdermal patches to implants, bioadhesive systems, and injectable peptide/protein drugs from biodegradable polymers (Breimer, 1999). A sustained release implant based on nanopore membrane technology is an example of such a device, where a nanopore membrane is fitted into a capsule for subcutaneous implantation and exploited as a delivery device for a variety of chemotherapeutic drugs by precisely controlling the membrane pore length, size, and density. This approach exhibits a constant drug release rate over a period of several weeks, thus avoiding the burst effects of other methods. Mathematical modeling of this application can be useful to quantify effective release rates and maintain constant delivery over the required time interval. Molecular diffusion dynamics of a solute across semi-permeable membranes can be described based on Fick’s laws of diffusion (Saatdjian, 2000), keeping in mind that diffusion through silicon nanochannels having sizes comparable to the solute molecular dimensions has been shown to be non-Fickian (Martin et al., 2005). We employ the mathematical model of Cosentino et al. (2005) that describes these experimental results by means of a reasonable interpretation (as we discuss later) and, at the same time, recovers the classical diffusion laws in the unconstrained case.

Most of non-Fickian diffusion cases observed experimentally have been attributed to wall drag effects or single file diffusion (Clark et al., 2000, Gupta et al., 1996, Hahn et al., 1996, Kukla, 1996, Meersmann et al., 2000, Wei et al., 2000). The case analyzed in Martin et al. (2005) differs from these studies because the membrane is made up of silicon and fabricated by photolithographic techniques, and the pores are rectangular and nanometric only in one dimension (other dimensions are in the mm range). The observable macroscopic effect is a prolonged linear release of several molecules, eventually assuming an exponential Fick’s profile. Fick’s first law for a binary mixture is , where is the diffusion coefficient of solute A in solvent B in a reservoir of volume V_T and J_A is mass (or molar) flux (with respect to the mass average velocity). The following assumptions hold in order to obtain a suitable model: (a) The experimental volume, V_T, contains a total mass of drug A; it can be divided into two compartments of volume V₁ (the reservoir) and V₂ (the sink), with respective initial mass concentrations and (); (b) concentrations are homogeneous in each compartment and concentration variation is spatially defined in a thin boundary region of depth L; (c) given a Cartesian reference system (x,y,z), the concentration gradient has zero components along the y and z axes. The aim is to calculate the mass flux of drug through a generic surface of area S that is assumed to be perpendicular to the diffusion path. Taking into account assumptions b) and c), this flux can be approximated as. Denoting the mass concentrations at time t in the two compartments by and , the conservation principle produces (for all t), . It is then possible to compute the governing equation for the concentration: , which yields , where . The corresponding flux is

(28)

Therefore in the free diffusion case the release profile is exponential.

A constrained diffusion model can be described as follows. Experimental results in Martin et al. (2005) showed that release profiles are linear for a certain period, and then switch to a Fickian exponential trend. This observation suggests that the model developed above for the unconstrained case is suitable in some regime. The effect of the membrane can be modeled by means of saturation of the mass flux:

(29)

where is the saturation threshold, and the saturation function has been defined as

. (30)

This assumption is reasonable if each nanochannel is viewed as a bottleneck. Thus, the molecular flux through a channel will remain invariant over a certain concentration level regardless of the number of particles per unit volume in the reservoir compartment. This description coincides with the classical diffusion laws if the threshold value is very large, as in the case of unconstrained diffusion. Finally, the switch between linear and exponential diffusion can be explained by the concentration decreasing in the reservoir, eventually forcing the concentration gradient (and flux) to fall below the threshold value.

Two different cases are therefore possible, depending on the value of concentration:

Constrained diffusion: If , we obtain , with a linear release profile.

Free diffusion: If , we derive the same evolution laws as in the unconstrained case, with an exponential release profile.

This model can be used to fit the experimental data presented in Martin et al. (2005). The saturation threshold parameter value is found by a least square interpolation of the data; interpolations using another model (Peskir, 2003) can also be obtained through a generalized diffusion law on the basis of the van der Waals equation for real gases, obtaining a two-parameter model that provides a more physical interpretation at the price of a greater mathematical complexity. The two models yield similar results, which suggests a novel interpretation of non-Fickian diffusion through nanochannels: it is possible to hypothesize that van der Waals forces are influenced by the presence of the nanochannel, and that this perturbation induces a saturation on the molecular flux when concentration gradients are high enough and the channel height is comparable to molecular size. Model results are compared to experimental data in Figure 7.

Denoting the channel height by h and the molecular hydrodynamic radius by r_s we can identify a linear dependence of with respect to the ratio . Recalling that the flux determines the solute delivered per unit area during the linear regime, the following considerations can be made: a) In order to achieve sustained zero-order kinetics, the nanochannel height has to be comparable with molecular size of solute molecules; b) Once a suitable value of q has been chosen, the corresponding value of can be derived by interpolation using the linear relation given in Figure 7 (bottom right) or through experimental testing (the dependence of on q would need further validation by interpolating with more experimental data points); c) the desired release rate can be achieved by acting on the effective porous area S, i.e. the greater the porous area, the greater the amount of drug released per unit time;

d) the duration of the linear release regime is tuned by varying the amount of drug in the reservoir. Provided that the initial concentration gradient is large enough to yield the flux saturation, the flux reverts to a Fickian behavior only after a certain amount of drug in the reservoir has been delivered. The concentration gradient is thus diminished under a certain concentration threshold.

2.8 Pharmacokinetics and Pharmacodynamics

Although the ideal for nanoparticles and liposomes would be to preferentially and uniformly extravasate through lesion vasculature, and then to uniformly release drug over long periods, such performance has not yet been achieved. We must thus incorporate a pharmacokinetic component in the multiscale model. Pharmacokinetic modeling employs compartment modeling (Holz & Fahr, 2001) to study cellular drug-uptake and intracellular drug interactions and to provide insight into modeling of cellular-scale mechanisms of drug resistance. Dordal et al. (1995) used a standard 3-compartment model to investigate cellular drug uptake and to quantify decreased intracellular sequestration, increased efflux, and decreased membrane permeability as they relate to reduced drug effectiveness. By fitting experimental data to the model, they obtained kinetic parameters for both inward and outward transport, and used them to quantify the relative importance of the cellular mechanisms. Their results indicated that of the three cellular mechanisms modeled, decreased intracellular sequestration in a non-exchangeable compartment is the most significant contributor towards drug resistance (Dordal et al., 1995). Similarly, compartment modeling can be applied to investigate additional components affecting drug delivery such as target repair mechanisms and extracellular drug binding (Sanga et al., in press).

An accounting of tumor-level and cellular pharmacokinetics can be approximated by studying bolus administration of free drug, which by virtue of its variable nature typically yields an initial rapid rise in plasma concentration followed by a slower clearance. Drug diffusion through interstitium, cellular uptake, and target binding (e.g., DNA) must be accurately tracked in order to determine target, as opposed to plasma, AUC. To describe tissue-level and cellular pharmacokinetics in our multiscale model we employ the four compartment model of Sinek et al. (in press), following earlier work by Dordal et al. (1992, 1995), Jackson (2003), and El-Kareh and Secomb (2003). Although the model was formulated with cisplatin and doxorubicin in mind, it is easily extendable. Compartment 1 is extracellular interstitium, Compartment 2 is intracellular cytosol, Compartment 3 is DNA, the target of both drugs, and Compartment 4 is intracellular lysosomes, implicated in the sequestration and removal of doxorubicin (Hurwitz et al., 1997). The system of equations is

(31)

where s_i is drug concentration in Compartment i, s_M is the saturation capacity of the DNA, the k_ij’s are transfer rates from Compartment i to j, n_s is a spatially and temporally variable drug production rate related to quantity of extravasated nanovectors and their drug release characteristics, and D_s is interstitial drug diffusivity. The dimensionless delta function d is supported along the vasculature, close to which nanovectors are assumed to remain after extravasation. The term in the second equation represents drug removal from the cytosol via glutathione conjugation, while in the third represents repair mechanisms that remove bound drug from DNA. The model is flexible and allows for the incorporation of many cellular drug resistance mechanisms through appropriate rate settings. For example, glutathione conjugation of intracellular cisplatin and DNA repair of cisplatin-DNA adducts are reflected in the values of k₂ and k₃. Another well-known mechanism, the removal of intracellular doxorubicin via P-glycoprotein or multidrug resistance protein efflux pumps, can be incorporated by increasing the ratio of k₂₁ to k₁₂.

Simulations of bolus injections have demonstrated penetration differences between doxorubicin and cisplatin that are in agreement with experiment (Tannock, 2001; Tannock et al., 2002; Zheng et al., 2001). Figure 8 shows the DNA-bound drug distribution in a 150 mm thick section of tissue (the approximate diffusional limit of oxygen before necrosis sets in) adjacent to a blood vessel. While cisplatin enters cells relatively slowly, doxorubicin enters rapidly and binds with high affinity to intracellular components, especially DNA and lysosomes (Demant & Friche, 1998, Dordal et al., 1992, Dordal et al., 1995, Erlansson et al., 1992, Jekunen et al., 1993, Paul et al., 1979, Tannock et al., 2002). Thus, doxorubicin drug gradients are particularly strong.

Figure 8. Upper two panels: DNA-bound drug (left, cisplatin, right, doxorubicin) within a 150 micron thick section of tissue adjacent to a blood vessel as depicted by inset. Times indicated are referenced from initiation of a two-hour bolus. Lower two panels: Resulting survival distribution within the same tissue sections at a total survival of 50 percent. Adapted with permission from Sanga et al., Exp. Rev. Anticancer Ther. (in press).

While the study of drug distribution throughout a lesion is important, perhaps even more important is the resulting distribution of cell inhibition, which is what determines tumor regression, and resulting tumor morphology. In vitro and in vivo data suggest that hypoxia due to chemotherapeutic treatment (especially anti-angiogenesis) promotes invasiveness (Pennacchietti et al., 2003; Bello et al., 2004). Furthermore an important mechanism underlying this process involving inhomogeneous growth and regression has been hypothesized and recently studied using the technology presented herein (Cristini et al., 2005). In order to model these and other phenomena of interest, a pharmacodynamics component is necessary to translate nutrient and drug distributions into growth and regression.

Because of the difficulty in capturing the intricacy and abundance of intracellular signalling pathways governing apoptosis, pharmacodynamic models tend to be highly phenomenological, often employing Hill-type equations yielding cell survival as a function of some “damage measure,” such as extracellular concentration-time exposure (AUC) (see Figure 9). Such a model is largely an empirical fit to data, yielding little insight into the mechanics of drug action and cell response. An alternative is a simple exponential kill model. An improvement over both these models is the use of DNA-bound drug as the damage measure, as demonstrated by El-Kareh and Secomb (2003 and 2005). Their investigation was prompted by the observation that models employing extracellular AUC consistently overestimated cytotoxicity in cases of extended exposure to the drugs cisplatin and doxorubicin (Sanga et al., in press). Toxicity would experimentally achieve a plateau above that continued exposure to fresh drug would have no effect. Thus, they hypothesized that it was not the time of exposure per se that correlated with cytotoxicity, but rather the peak level of DNA-bound drug (El-Kareh and Secomb, 2003). By using this measure they showed that for short exposure times, the delay in achieving DNA-bound drug equilibrium could explain increasing cytotoxicity in time. Their model consistently proved to be the best fit even for long exposure in in vitro datasets (Troger et al., 1992) establishing that peak DNA-bound cisplatin is a stronger indicator of cytotoxicity than extracellular or intracellular concentrations. Later, the model was extended to doxorubicin, providing better fits to in vitro cytotoxicity data than previous models (El-Kareh and Secomb, 2005). Lankelma et al. (2003) further applied the concept of “fading memory” to capture the history of cell injury and rejuvenation.

Figure 9. A Hill-type model of in vitro cell survival as a function of extracellular doxorubicin concentration-time (AUC) exposure. Markers represent experimental results of Levasseur et al. (1998) using A2780 human ovarian cancer cells. Curve is produced by the Hill-type equation shown, where S is survival and x is extracellular AUC.

Since drug-induced cell death is fundamentally an exponential process, albeit complicated by the intricacies of the cell cycle and mechanisms designed to repair injury, we opt for the exponential kill model recently proposed by Sinek et al. (in prep.), which recently incorporated a multi-compartmental tissue- and cell-level pharmacokinetics and pharmacodynamics (PKPD) model for cisplatin and doxorubicin based on experimentally derived parameter estimations, thereby establishing a more rigorous platform for analyzing the effectiveness of chemotherapeutic drugs. PKPD modeling is certainly an established field with representative examples given by Panetta (1997), Gardner (2002), Lankelma (2003), and Jackson (2003). While the model developed by Sinek et al. (in prep.) is founded on similar principles as prior work, its novelty lies in its coupling with cancer progression simulators tracking tumor growth and angiogenesis, i.e. Zheng et al. (2005), thus providing a platform for simulating and analyzing chemotherapy applied to vascularized in vivo tumors.

Although the PKPD model of Sinek et al. (in prep.) is designed specifically for cisplatin and doxorubicin, the intention is to present a protocol for adapting the model to simulate therapy with other drugs. The model presumes that drug concentration remains constants along the vasculature during a simulated intravenous bolus administration. The vasculature acts as a source for oxygen, nutrient, and drug throughout the tumor. The concentration of drug decreases due to cellular uptake as it diffuses from the vasculature into the tumor. Through the appropriate adjustment of transfer-rate coefficients between compartments, the model explicitly accounts for tissue- and cell-level biobarriers, and tracks drug penetration from its source (i.e., the vasculature) through the lesion interstitium, cell membranes, and into intracellular organelles and eventual target. The pharmacodynamic effect is based on the amount and time of exposure of DNA-bound drug. The drug undergoes various forms of transport resistance along its path towards its intended target (e.g., DNA), which consequently diminish the drug’s efficacy. By providing precise control over the parameters corresponding to PKPD elements, the model of Sinek et al. (in prep.) can delivery hypothesis-testing capability (Sanga et al., in press). The system of equations modeling PK is capable of tracking the amount of drug both spatially and temporally through multiple compartments based on governing rate-parameters. These can be tuned to specific drug and cancer type, and as well as to cancer grade and cellular subtleties of individual patients. Sinek et al. (in prep.) first derived parameter values through experimental data reported in published literature; in the future, it is intended that a targeted histological, cellular, and genetic analysis of tissue biopsies will provide a protocol for determining model parameter values in a consistent manner.

The pharmacodynamic component is a Hill-type, phenomenological model similar to El-Kareh and Secomb’s work (2003, 2005), which also takes into account the extended exposure plateau of cytotoxicity:

(32)

where E is cell inhibition (1 minus surviving fraction), x is DNA-bound drug-time product (area under the curve, or AUC), A and m are phenomenologically fit parameters, and N is a function of nutrient n ranging from 0 to 1 used to mimic the effect of hypoxia and hypoglycemia (Sinek et al., in prep).

While the model retains unavoidable phenomenological features, the exponential component allows for further development as cell-cycle kinetics and the governing signalling pathways are incorporated. It is anticipated that a much more powerful pharmacodynamics model will result from a detailed stochastic formulation of relevant intracellular events and mechanisms. This description is critical for predicting lesion response to nanovectored drug, where exposure times are longer and concentrations more uniform than for intravenously administered free drug. While it is hoped that complete cancer regression will result in vivo under such conditions, the possibility remains that some cells will dynamically adjust their phenotype to increase drug resistance during exposure.

In silico investigations by Sinek et al. (in prep.) have provided insight into the behavior of the well-known anticancer agents doxorubicin (see Figure 10) and cisplatin. Results indicate that despite well-known penetration difficulties (Tannock et al., 2002, Zheng et al., 2001, Durand et al., 1990), doxorubicin clinically performs better than might be expected. This appears to be due to drug retention in tissue remote from vasculature, resulting in a more homogeneous AUC (area under curve) distribution, and causing significant cell inhibition. While this retention phenomenon and its consequences had been established earlier (Durand et al., 1990), its observation using only basic tissue/cell/drug parameters (e.g., membrane/drug permeability and drug interstitial diffusivity) as input to an in silico simulation underscores the powerful potential of multiscale modeling of tumor growth and angiogenesis fully integrated with experimentally observed biological phenomena.

Figure 10. Top: Doxorubicin is systemically administered via a 2-hour bolus. Drug distribution through the tumor at 0, 6, 12, and 18 hours after the bolus is shown from left to right, respectively, decreasing accordingly and as indicated by darker color in the left frame and lighter color in the right frame. Bottom: The tumor is shown at days 0, 8, and 16 from left to right, and undergoing regression due to the cytotoxic effects of doxorubicin. Black solid line represents tumor boundary, radial lines represent tumor vasculature, and dark regions represent necrotic tissue.

2.9 Genotype

In order to study molecular mechanisms of drug resistance, the multiscale model herein needs to have not only multispecies but also multiphase characteristics. This enables the description of multiple, coexistent clones in a tumor microenvironment responsible for drug resistance, which may be linked to molecular mechanisms that enhance cell survival, such as decreased drug uptake, increased drug efflux, diminished apoptosis, DNA damage repair, and alterations in drug target, drug metabolism and cell cycle checkpoint mediators (Dalton & Salmon, 1992, Gottesman et al., 2002). Resistance is favored by changes in protein expression (Poland et al., 2002), localization (Oloumi et al., 2000), and altered gene expression, such as amplification or mutation of genes encoding efflux pumps (Cadman, 1989, Desoize & Jardillier, 2000, Knowles & Phillips, 2001, Laderoute et al., 1992, Olive & Durand, 1994, Oloumi et al., 2002, Sutherland, 1988, Theodorescu et al., 1993, Timmins et al., 2004, Wartenberg et al., 1998). Resistance beyond that intrinsic to and acquired by particular cells is introduced by the tumor mass itself as a three-dimensional physical object (Kobayashi et al., 1993) affecting micro-scale conditions such as cellular density and arrangement, DNA conformation (Olive & Durand, 1985), and cell cycle stage (Olive & Durand, 1994). The extracellular matrix may provide structurally based resistance (Kobayashi et al., 1993), as could inter-cellular interactions, perhaps via adhesion molecules or tissue /cellular architecture (Graham et al., 1994, Olive & Durand, 1994, Rainaldi et al., 1999) supporting apoptosis suppression through increased cell-cell adhesion (Bates et al., 1994). Drugs may be sequestered in outlying tumor regions (Durand, 1990, Erlansson et al., 1992), although active efflux by cells in inner layers may minimize tissue penetration barriers (Wartenberg et al., 1998).

Resistance introduced by tumor tissue in 3-D may stem from substrate gradients in the cellular microenvironment producing hypoglycaemia, hypoxia, and acidosis (Fernandez et al., 2000), which may select for more resistant phenotypes (Raghunand et al., 2003), decrease drug sensitivity (Tannock, 2001), alter cell cycle kinetics (Durand & Vanderbyl, 1989), reduce cell proliferation (Sutherland, 1988), and foment a toxic cellular environment favoring upregulation of efflux mechanisms (Wartenberg et al., 1998). Diffusion gradients may enforce drug resistance at the cellular-scale through various molecular mechanisms. Hypoglycaemia, hypoxia, and acidosis can lead to accumulation of unfolded or malfolded proteins in the endoplasmic reticulum (ER), triggering transcription of several genes as part of the unfolded protein response: GRP78 (Lee, 2005, Park et al., 2004, whose expression has been linked to drug resistance (Dong et al., 2005), GRP94, whose elevated basal levels have been detected in cells with reduced proliferative capacity (Gazit et al., 1999), and HSP27, which has been shown to cause drug resistance (Fuqua et al., 1994). Glucose starvation causes the glucose-regulated stress response (Lee, 1987), which has been correlated with resistance to topoisomerase II-directed chemotherapeutic drugs such as doxorubicin (Tomida & Tsuruo, 1999, Fernandez et al., 2000, Li & Lee, 2006) through a decreased expression of topoisomerase II (Shen et al., 1989, Yun et al., 1995). Glucose deprivation further induces cellular oxidative stress (Lee et al., 1998, Spitz et al., 2000). Persistent oxidative stress at sublethal levels, caused by endogenous mechanisms, exposure to drug, and glucose deprivation, promotes cancer cell viability and resistance to apoptosis (Brown & Bicknell, 2001). Drug resistance is also abetted by gradients of growth factors (e.g. serum), reducing cell proliferation (Lee et al., 1997, Wosikowski et al., 1997, Hug et al., 1986) and leading to proliferation gradients (Zhang et al., 2005) with outlying tumor cells having highest mitotic activity and cells near necrotic regions having reduced proliferative rates or being mostly quiescent (Wartenberg & Acker, 1995, Carlsson, 1977). Cell cycle arrest can minimize cytotoxicity because some drugs (e.g., doxorubicin) induce antineoplastic and toxic effects through DNA intercalation and inhibition of DNA and RNA polymerases (Zunino et al., 1975), DNA alkylation (Taatjes et al., 1996), as well as interaction with topoisomerase II (Osheroff et al., 1994). Drug gradients might induce further resistance (Durand, 1981) due to physical barriers associated with tissue morphology that drug molecules must overcome to reach their nuclear target (Durand, 1990, Sutherland, 1988). Acidosis due to pH gradients can positively charge a drug such as doxorubicin, hindering its passive translocation over the cellular membrane (Mahoney et al., 2003, Kozin et al., 2001), and by ion gradients on drug distribution and ion trapping, where weakly basic drugs will concentrate in more acidic cellular compartments instead of reaching their intracellular target (Raghunand & Gillies, 2000). In addition, the effect of basic drugs on topoisomerase II activity is optimal at alkaline pH (Gieseler et al., 1996).

These effects at different scales seamlessly influence each other in tumor tissue, with genetic abnormalities triggering tumorigenicity and local environmental conditions inducing genetic change. For example, the unrestrained growth of certain tumors generates an increased nutrient requirement and leads to hypoxic tumor regions and necrosis, which may in turn trigger angiogenesis (Principi et al., 2003) and enhanced invasive capability (Young et al., 1988, Young & Hill, 1990, Cairns et al., 2001, Postovit et al., 2002, Rofstad et al., 2002). A multiscale model of nanovectored chemotherapy needs the capability to include a tumor cell population containing varying genotypes, in order to model the effects on tumor growth and morphology, such as those, for example, caused by expression of oncogenes and tumor suppressor genes affecting tumor growth and survival. Since proliferation of cells expressing a more aggressive phenotype may also be unwittingly facilitated by therapeutic strategies employing drugs that are more lethal to less malignant cell species, heterogeneity in the microenvironment through non-uniform cell death and selection for resistant clones also needs to be included.

A genotype module can be implemented that specifies transformation through random mutation from a less aggressive Species 1 to a more aggressive Species 2 (Frieboes et al., in review). For the latter, a set of oncogenes that boosts cell proliferation is assumed to double nutrient consumption and mitosis rate, and a set of mutated tumor suppressor genes is assumed to down-regulate the apoptosis rate. A cell population undergoing genetic mutations is thus assumed to contain at least two transformed species, each characterized by a genotype consisting of a set of tumor suppressor genes and a set of oncogenes. We have integrated this simplified genotype model with the multi-dimensional diffuse-interface component described earlier (Wise et al., in review), thus enabling quantitative description of evolution and progression of multiple clones in the multiscale model. A simple genotype is described by the matrix M=[o_g,t_g], where each matrix element describes possible states for a set of oncogenes and tumor-suppressor genes, while 0 and 1 indicate silent and active gene sets, respectively. This approach enables the specification of constitutive links between genotype and tumor growth.

Simulations were created to investigate how perturbations due to genetic mutations may lead to changes in morphology. Mutations were randomly introduced with constraints consistent with a time progression from relatively normal cells with o_g=0 and t_g=1 to more abnormal clones, the most malignant having o_g=1 and t_g=0, thus proliferating without restraints (see Table 1). Results revealed that as more malignant species appear, they perturb the nutrient environment by increased uptake, thus leading to tumor morphological instability and infiltrative fingering of glioma branches into normal brain tissue. Results matched glioblastoma characteristics, where areas of necrosis are an indication of aggressive cell proliferation as observed in vivo through MRI.

GENOTYPE M	MITOSIS RATE lM	APOPTOSIS RATE lA	NUTRIENT UPTAKE hn
[0,1]	Unchanged	Unchanged	Unchanged
[0,0]	Unchanged	Downregulated	Unchanged
[1,1]	Upregulated	Unchanged	Upregulated
[1,0]	Upregulated	Downregulated	Upregulated

Table 1. Definition of a genotype matrix M to specify differences between Species 1 and Species 2 as a function of glioma genotype, and the corresponding implications for model parameters.

This model also predicted (Frieboes et al., in review) that application of therapy, even in a uniform manner as might be achieved with nanovectors, could destabilize tumor morphology and lead to increased invasive potential. Therapy results would be short-lived because by ten months, the tumor would resume growth and infiltration of surrounding tissue, driven by aggressive proliferation of both species. One of the consequences of therapy was an increased invasiveness of the two species, shown as extensive fingering into healthy tissue and perhaps due to the therapy-induced exacerbation of hypoxia (Figure 11). This result corresponds to experimental observations in vitro and in vivo (Cairns et al., 2001, Rofstad et al., 2002, Young et al., 1988), in which hypoxia was determined to trigger aggressive tumor invasion.

Figure 11. Modeling of glioblastoma growth and response to treatment with a tumor containing two different genetic species. Left: Nutrient field shows necrosis (darker areas, lower left) where the more aggressive species proliferates faster. Right: After simulated treatment, tumor mass resumes growth after initial regression and exhibits infiltrative fingering of host tissue. Note increase in necrosis (darker areas) as a result of therapy, coinciding with this infiltration. Adapted with permission from Frieboes et al. (In review).

3. Model Integration and Parameter Measurement

While most of the mathematical models of tumor growth and angiogenesis developed to date provide worthwhile insight into cancer-related processes at a particular time and length scale, they fail to tackle the issue of integration across these scales. Phenomena occurring at various time and length scales must be appropriately coupled to capture the dynamics involved. Multiscale systems modeling complex biological processes such as cancer (Alarcon et al., 2003 and 2005, Ribba et al., 2006, Jiang et al., 2005) and various phenomena related to developmental biology (Sharp et al., 1993, Chaturvedi et al., 2005) have been modeled. Jiang et al. (2005) and Alarcon et al. (2003, 2005) present frameworks for building multiscale cancer progression models capable of integrating a hierarchy of processes at varying time and length scales. Most cancer models and multiscale systems (Araujo & McElwain, 2004, Alarcon et al., 2003 and 2005, Ribba et al., 2006, Jiang et al., 2005, Thomlinson & Gray, 1955, Burton, 1966, Greenspan, 1972, Byrne & Chaplain, 1995 and 1996, Maggelakis & Adam, 1990) mainly produce one- and two-dimensional simulations, which their ability to represent the complex three-dimensional in vivo tumor microenvironment.

The components reviewed in Section 2 are integrated here to form the multiscale model of nanovectored therapeutics delivery to tumors as shown schematically in Figure 12. In the following discussion, parameters and field variables are distinct, parameters having been measured before and remaining more or less fixed during a simulation. Variables (e.g., nutrient and pressure) result from operation of the simulator with its given parameters. This distinction is sometimes fuzzy, since we expect that some “parameters” on the coarser scales will be determined from activities on the finer scales. Still, this distinction serves as a guide. Tumor mass evolution starts with Tumor Growth specifying a pressure p and boundaries for viable and necrotic tumor tissues. Pressure becomes an input to Nutrient Transport, and the tumor boundaries become an input to Angiogenesis and Pharmacokinetics. In turn, Angiogenesis determines vascular boundary d for Nutrient Transport, Pharmacokinetics and Nanovector, and nutrient n and vascular transfer function n_n for Nutrient Transport, as function of flux Q in the vasculature. Pharmacokinetics DNA-bound drug s₃ is input to Pharmacodynamics and Genotype. Pharmacodynamics computes and submits apoptosis rate l_A to Tumor Growth. Genotype determines the mutation (transition) of one cell species to another. To each species corresponds a set of survival (mitotic, apoptotic), nutrient uptake, motility, and drug sensitivity parameters, which specify constitutive links from the molecular to the tumor scale. Nanovector supplies nanovector vascular transfer function n_s to Pharmacokinetics. Finally, Nutrient Transport calculates nutrient level n for Tumor Growth.

Figure 12. The multiscale model. Components are shown in rectangular boxes. Oval boxes indicate experimental methods to obtain values for model parameters, which are then input into the appropriate components as shown by thick arrows. Thin arrows represent input and output of variables and functions between model components

Accurate rendering of in vivo tumor growth depends upon correct parameter calibration, for which a wide variety of experimental techniques exists. It should be noted, however, that a naïve reliance on “what you see is what you get” is neither sufficient nor practicable. Rarely can parameters be determined directly from experimental data; rather a model is posited and parameters are fit to provide the best representation of the data. This approach is true even in the simplest cases, such as determining membrane permeability to a specific substrate. Permeability is not itself witnessed, but rather the concentration-time relation of substrate on one or both sides of the membrane. Often, it is difficult to isolate parameters for measurement for at least two reasons. First, in any complex system, performance is greater than the sum of component behaviors. Equivalently, parameter values often change when measured out of context, as would be expected when attempting to determine mitotic proclivity l_M of in vivo cells from in vitro monolayers. Second, the determinants that are subsumed in a particular model parameter are usually multiple and incompletely understood. As an example, nuclear drug affinity k₃₂ in the pharmacokinetics system (Eq. 21) is influenced not only by direct chemical binding between drug and DNA molecules, but also by protein vaults, which entrap and remove drug from the nucleus. While each of these may conceivably be modeled and independently measured, the question becomes where to stop. At some point it is necessary to measure the parameter in its full context, accepting the fact that a well-defined number may be representing the cumulative effects of poorly understood mechanisms. In extreme cases, parameters are entirely phenomenological, not reflecting any recognized substance or mechanism, but necessary nonetheless to produce a “good fit.” This is amply illustrated by the parameters P₁, P₂, n, and m from Eq. 22 and shown in the Hill-type equation in Figure 9. The challenge, therefore, is not only to develop experimental techniques for acquiring data, but also analytical methods for inferring parameter values from it. Such methods should, at a minimum, reliably provide key statistics such as means and deviations, and should look beyond the common (and often unsubstantiated) assumption of normality. Ordinary least squares, generalized method of moments, and maximum likelihood estimation will figure prominently in this endeavor.

Deriving parameters is only a first phase of a process that must necessarily end with model validation. We envision a rigorous process whereby parameters are first derived by applying the methods discussed above to training data, and goodness-of-fit is gauged relative to this data. The process is completely analogous to the well-known statistical practice of performing linear regression. At this stage, a poor fit of model performance to training data will require further model development, just as a low correlation coefficient indicates that a linear regression model is inappropriate. Once a good fit is achieved, the process of validation, or testing, may be entered. This is analogous to using the results of a linear regression to predict values outside the range of the fitting data. The tumor model must be tested in new circumstances, and its performance compared to the in vitro or in vivo system that it is simulating. Success at this stage builds confidence in the model, while failure again indicates design improvement is required.

Important parameters needing estimation are listed in Table 2 and discussed below along with proposed experiments to provide training data. Common in vitro experiments revolve around the use of cell monolayers and suspensions. Monolayers have been widely utilized for measuring mitosis l_M, apoptosis l_A, and nutrient metabolism h_n under conditions of varying nutrient and oxygen. Some information can also be gleaned regarding the nutrient threshold n_N at which cells necrose, as well as the rate l_N at which this happens. Both monolayers and suspensions can be used to fix cellular drug pharmacokinetic parameters (k_i’s and k_ij’s) (Dordal et al., 1995; Sadowitz et al., 2002). In particular, intracellular platinum quantification typically relies upon atomic absorbance spectroscopy, while doxorubicin quantification relies upon its fluorescent properties. Monolayers have also long been used for fitting Hill-type pharmacodynamic parameters (Troger et al., 1992; Levasseur et al., 1998).

In vitro spheroids provide a more realistic environment in which to measure the above parameters, since they more closely approximate in vivo conditions (Sutherland, 1988). They additionally provide a means for measuring parameters pertaining to tissue-level determinants, such as the lysis rate l_L (which reflects fluid efflux through the spheroid surface), nutrient and drug interstitial diffusivities, and some information regarding cell mobility m and adhesion g. Measurements of surface tension (cell adhesion) g as reflected in tumor morphology cohesiveness can also be obtained by varying glucose concentration and growth medium concentration. In vitro scaffolds may be another experimental platform, with greater power to measure mobility and adhesion, since the shape of the cell mass may be controlled and tensions may be manipulated. The methods used for extracting parameter values from spheroid experiments rely more heavily on model assumptions and construction than do the methods associated with monolayers and suspensions. For example, the simple steady-state equation governing nutrient n (or oxygen) diffusion and uptake h through a perfect spheroid of radius r is , where D is the diffusion constant. This is, strictly speaking, a model assumption, albeit widely accepted. The data available would typically be oxygen tension from surface into the center (Mueller-Klieser & Sutherland, 1982). This data is not synonymous with the parameters values desired, and only through fitting it to the solution of the differential equation can some information about the parameters be obtained, namely, the ratio h/D. In a similar fashion, drug, proliferation, and pH gradients, as well as the depth of the necrotic interface yield information from which associated parameters can be obtained. Robust statistical methods as discussed earlier become increasingly important for deriving parameters from three-dimensional cell cultures, and are essential when using in vivo experiments, as discussed next.

Since angiogenesis is difficult to reproduce in most in vitro settings, in vivo clinical data, animal (e.g. mouse) xenograft models with subcutaneous tumor implantation, and ex-vivo models (e.g. chick CAM assays) could be used. The in vivo mouse model includes a window that enables, for example, direct observation and measurement of vessel formation along with TAF and fibronectin distributions (R. Gatenby, in press). These data can be used to infer values of endothelial cell diffusivity D_e, fibronectin production n_f and uptake h_f, TAF uptake h_c, and chemotaxis a_c and haptotaxis a_f strengths Both the in vivo mouse and ex vivo chick models can provide data on vascular hemodynamics and intravascular nutrient and nanovector transport and extravasation. Combining this with directly measureable pO₂ and pH distributions and spatially variable tumor growth provides a means of determining vascular transfer coefficients, n_n and n_s, of nutrient and nanovectors. The clinical application of dynamic contrast-enhanced magnetic resonance imaging (DC MRI) provides another means for obtaining data regarding tumoral vasculature, including intravascular nanovector transport and extravasation (Su et al., 1995, Su et al., 1998). Bolus injections of contrast agents allow the tracking of tumor saturation in time, yielding curves to which permeability and vascular density functions can be fit.

Finally, parameters such as DNA repair rate and cell permeability could be manipulated at the genetic level by introducing agents that modify expression of cell proteins. Effects of these agents on cellular pathways could be modeled to study different nanovector designs to accomplish it. Genomic profiling provides an additional technique to elucidate molecular mechanisms of drug resistance (e.g. Francia et al., 2004). These experiments would be done first on established cell lines to perfect experimental methods and later on patient tissue biopsies.

Parameter	Meaning	Section
Tumor Growth
m	Cell Mobility	2.1
lM	Mitosis Rate	2.1
lA	Apoptosis Rate	2.1
lN	Necrosis Rate	2.1
lL	Lysis Rate	2.1
nN	Necrotic Thershold	2.1
g	Surface Tension (Cell Adhesion)	2.1
Angiogenesis
De	Endothelial Cell Diffusivity	2.2
nf	Rate of Fibronectin Production	2.2
hf	Rate of Fibronectin Uptake	2.2
hc	Rate of TAF Uptake	2.2
ac	Strength of Chemotaxis	2.2
af	Strength of Haptotaxis	2.2
Nutrient Transport
nn	Vascular Transfer Function for Nutrient	2.4
Dn	Interstitial Nutrient Diffusivity	2.4
hn	Rate of Nutrient Uptake	2.4
Nanoparticle Targeting, Drug Release, and Pharmaco-kinetics/-dynamics
k^+	Rate of Bond Formation	2.6
k^-	Rate of Bond Dissolution	2.6
ns	Vascular Transfer Function for Nanovectors	2.8
Ds	Drug Diffusivity	2.8
ki and kij	Compartment Transfer & Release Rates of Drug	2.8
P1, P2, m1, m2	Pharmacodynamic Parameters	2.8

Table 2. Important parameters and functions of the multiscale model needing estimation from experiments.

4. Summary and Future Work

We have presented a multiscale model of nanovectored therapeutics delivery to cancerous lesions and tumor response to treatment. This model covers a wide range of scales, from that of introduction of nanovectors into whole blood and subsequent removal via the reticuloendothelial system, through the scale of particle navigation in tumoral vasculature and cellular pharmacokinetics, down to the molecular scale of drug molecule release, binding to DNA, and alteration of gene expression. This is accomplished through a modular design in which each model component is dedicated to one process, passing requisite information to the others as needed. The primary component, around which all others are built, is the growth module, which capitalizes on technically demanding numerics to accurately render lesion response.

Here we list a few of the main areas in need of further development. The Tumor Growth module should include an interstitial fluid phase and distinguish between interstitial hydrostatic pressure and mechanical pressure due to cell-cell interactions within tumor tissue (Jain, 1990, Jain, 2001, Padera et al., 2004). Furthering this idea, the explicit representation of the extracellular matrix needs to be included. A true representation of cell-cell adhesion throughout the domain should be developed to supersede the current model of surface tension at the tumor/host interface. Surface tension would then naturally arise due to adhesion. Tumor Growth also assumes that cell mass density is uniform in the tumor (Chaplain, 1996), and that regions become necrotic where nutrient and oxygen concentration fall below some specified minimum. In reality, density is not uniform, different species can coexist at various locations within tumoral tissue, and necrosis is not instantaneous. Furthermore, only one concentration is described, without making distinction between oxygen and cell nutrients. In addition, the effects of pH in the microenvironment should be included. A more detailed model of cellular function at the genetic and protein pathway level needs to be developed, and its functioning constitutively linked to tissue-level behavior. This model could serve as a basis for a much more mechanistically based pharmacodynamics component. An expanded quantification of the genotype component would also aid in the description of carcinogenesis and the role of mitochondria, and provide insight into the link between genotype and an invasive phenotype. Future subcomponents could include DNA repair, gene expression pathways, and cellular apoptosis mechanisms. Finally, a more sophisticated genotype matrix could incorporate the effects of mutations on multiple model parameters.

The current simplified implementation of the multiscale model has already been successfully employed to investigate lesion response to both free and nanovectored therapeutics, the results of which underscore the impact of nutrient and drug heterogeneity (Sinek et al., 2004). Further simulations have shown that these heterogeneities may be a powerful driving force behind “morphological instability,” enhancing tumoral invasiveness (Cristini et al., 2005; Frieboes et al., 2006). More importantly, since the role of nutrient heterogeneity is at least as important as that of drug, the implication is that nanovectored therapy may not, by itself, be powerful enough to mitigate adverse morphological response to treatment. Indeed, there is ample opportunity and need to model delivery of nanovectored thereapeutics to cancerous tumors, in turn testing and suggesting hypothetical strategies. This “hypothesis testing” capability of the multiscale model gives it the power to assist not only in designing carrier-based delivery systems and application, but also drug development and clinical trials design.

Another strongly anticipated role of the simulator is that of clinical “therapy optimization.” As development proceeds, results will acquire greater quantitative authority to the extent that reliable, real-time individualized clinical prediction will be possible. In order to achieve this level of performance, it will be necessary to develop methods to accurately measure parameters of an individual patient’s tumor. The design of these methods, in itself, will require careful modeling and validation, and includes genetic analysis, in vitro assays, and in vivo imaging, all of whose output can be supplied to the model to yield a recommendation of optimal therapy (Figure 13). The predictive capabilities of multiscale computational models such as the one presented here, combined with advances in nanoscale delivery systems, drug design, and genetic analysis, will hopefully accelerate the day when cancer treatment can be optimized to maximize therapeutic benefit on an individual patient basis.

Figure 13. Information flow enabling the multiscale model to be used clinically for formulation of individualized optimal therapy protocols. Images and tissue from a patient’s tumor, subject to a battery of analyses, would yield the necessary information for input to the model, which would then provide the prognosis under differing therapeutic strategies.

REFERENCES

B. Alberts, A. Johnson, J. Lewis, M. Raff, K. Roberts, and P. Walter, Molecular Biology of the Cell (4th Ed., Garland Science, New York, 2002).

H. Acker, J. Carlsson, R. Durand, R.M. Sutherland, Eds, Spheroids in Cancer Research: Methods and Perspectives (Springer-Verlag, New York, 1984).

H. Acker, J. Carlsson, W. Mueller-Klieser, R.M. Sutherland, Comparative pO₂ measurements in cell spheroids cultured with different techniques. Br. J. Cancer 56, 325-327 (1987).

T. Alarcon, H.M. Byrne, P.K. Maini, A cellular automaton model for tumour growth in inhomogeneous environment. J Theor Biol. 225, 257-274 (2003).

T. Alarcon, H.M. Byrne, P.K. Maini, A multiple scale model for tumor growth. Multiscale Model. Simul. 3, 440-475 (2005).

V. Ananthakrishnan, W.N. Gill, A.J. Barduhn, Laminar dispersion in capillaries: Part I – Mathematical analysis. Am. Inst. Chem. Engrs. J. 11, 1063-1072 (1965).

A. Anderson, M. Chaplain, Continuous and discrete mathematical models of tumor-induced angiogenesis. Bull. Math. Biol. 60, 857-900 (1998).

A. Anderson, X. Zheng, V. Cristini, Adaptive unstructured volume remeshing-I: The method. J. Comput. Phys. 208, 616-625 (2005).

R.P. Araujo, D.L.S. McElwain, A history of the study of solid tumour growth: the contribution of mathematical modeling. Bull. Math. Biol. 66, 1039-1091 (2004).

R.P. Araujo, D.L.S. McElwain, A mixture theory for the genesis of residual stresses in growing tissues I: a general formulation. SIAM J. Appl. Math. 65, 1261-1284 (2005).

R. Aris, On the dispersion of a solute in a fluid flowing through a tube. Proc. R. Soc. Lond. Series A, Math. Phys. Sci. 235, 67-77 (1954).

D. Arnold, F. Brezzi, M. Fortin, A stable finite-element method for the stokes equations. Calcolo 21, 337 (1984).

J.W. Baish, Y. Gazit, D.A. Berk, M. Nozue, L.T. Baxter and R.K. Jain, Role of tumor vascular architecture in nutrient and drug delivery: an invasion percolation-based network model. Microvasc. Res. 51, 327 (1996).

J.W. Baish, P.A.Netti, R.K.Jain, Transmural coupling of fluid flow in microcirculatory network and interstitium in tumors. Microvasc. Res. 53, 128-41(1997).

A.N. Barclay, Concluding remarks and the challenge from the immune system, Faraday Discuss. 111, 345-50, (1998).

R.C. Bates, A. Buret, D.F. van Helden, M.A. Horton, G.F. Burns, Apoptosis induced by inhibition of intercellular contact. J Cell Biol 125, 403-415 (1994).

D. Bazile, C. Prud Homme, M. Bassoullet, M. Marlard, G. Spenlehauer and M. Veillard, Stealth Me-PEG–PLA nanoparticles avoid uptake by the mononuclear phagocyte system. J. Pharm. Sci. 84, 493 (1995).

G.I. Bell, Models for the specific adhesion of cells to cells. Science 200, 618-627 (1978).

G.I. Bell, M. Dembo, P. Bongrand, Cell adhesion. Competition between nonspecific repulsion and specific bonding; Biophys. J. 45, 1051-64, (1984).

L. Bello, V. Lucini, F. Costa, M. Pluderi, C. Giussani, F. Acerbi, G. Carrabba, M. Pannacci, D. Caronzolo, S. Grosso, S. Shinkaruk, F. Colleoni, X. Canron, G. Tomei, G. Deleris, A. Bikfalvi, Combinatorial Administration of Molecules That Simultaneously Inhibit Angiogenesis and Invasion Leads to Increased Therapeutic Efficacy in Mouse Models of Malignant Glioma. Clin. Can. Res. 10, 4527–4537 (2004 ).

N. Bellomo, L. Preziosi, Modelling and mathematical problems related to tumor evolution and its interaction with the immune system. Math. Comp. Model. 32, 413–542 (2000).

A. Bredel-Geissler, U. Karbach, S. Walenta, L. Vollrath, W. Mueller-Klieser, Proliferation-associated oxygen consumption and morphology of tumor cells in monolayer and spheroid culture. J. Cell. Physiol. 153, 44-52 (1992).

DD. Breimer, Future challenges of drug delivery. J. Control. Release 62, 3-6 (1999).

C. Breward, H. Byrne, C. Lewis, The role of cell-cell interactions in a two-phase model for avascular tumour growth. J. Math. Biol. 45, 125–152 (2002).

C. Breward, H. Byrne, C. Lewis, A multiphase model describing vascular tumour growth. Bull. Math. Biol. 65, 609–640 (2003).

N.S. Brown, R. Bicknell R, Hypoxia and oxidative stress in breast cancer. Oxidative stress: its effects on the growth, metastatic potential and response to therapy of breast cancer. Breast Can Res 3, 323-327 (2001).

A.C. Burton, Rate of growth of solid tumours as a problem of diffusion. Growth.30, 157-176 (1966).

H. Byrne, M. Chaplain, Growth of nonnecrotic tumors in the presence and absence of inhibitors. Math. Biosci. 130, 151–181 (1995).

H. Byrne, M. Chaplain, Growth of necrotic tumors in the presence and absence of inhibitors. Math. Biosci. 135, 187–216 (1996).

H. Byrne, M. Chaplain, Free boundary value problems associated with the growth and development of multicellular spheroids. Eur. .J Appl. Math. 8, 639 (1997).

H. Byrne, J. King, A two-phase model of solid tumour growth. Appl. Math. Letters 16, 567–573.

Araujo, R., McElwain, D., 2004. A history of the study of solid tumour growth: The contribution of mathematical modelling. Bull. Math. Biol. 66, 1039–1091 (2003).

E.C. Cadman, The selective transfer of drug-resistant genes in malignant cells. In: D. Kessel, ed., Resistance to Antineoplastic Agents (Boca Raton, Fla: CRC Press, 1989), pp. 168-183.

R.A. Cairns, T. Kalliomaki, R.P. Hill, Acute (cyclic) hypoxia enhances spontaneous metastasis of KHT murine tumors. Cancer Res. 61, 8903-8908 (2001).

J. Carlsson, A proliferation gradient in three-dimensional colonies of cultured human glioma cells. Int. J. Cancer 20, 129-136 (1977).

J. Carlsson, H. Acker, Relations between pH, oxygen partial pressure and growth in cultured cell spheroids. Int. J. Cancer 42, 715-720 (1988).

J.J. Casciari, S.V. Sotirchos, R.M. Sutherland RM, Glucose diffusivity in multicellular tumor spheroids. Cancer Res. 48, 3905-3909 (1988).

M.A.J. Chaplain, Avascular growth, angiogenesis and vascular growth in solid tumours: the mathematical modelling of the stages of tumour development. Math. Comp. Model. 23, 47-87 (1996).

R. Chaturvedi, C. Huang, B. Kazmierczak, T. Schneider, J.A. Izaguirre, T. Glimm, H.G.E. Hentschel, J.A. Glazier, S.A. Newman, M.S. Alber, On multiscale approaches to three-dimensional modeling of morphogenesis. J. R. Soc. Interface 2, 237-253 (2005).

L.A. Clark, G.T. Ye, R.Q. Snurr, Molecular traffic control in a nanoscale system. Phys. Rev. Lett. 84, 2893-2896 (2000).

C. Cosentino, F. Amato, R. Walczak, A. Boiarski, M. Ferrari, Dynamic model of biomolecular diffusion through two-dimensional nanochannels. J. Phys. Chem. B. 109, 7358-7364 (2005).

C. Cozens-Roberts, J.A.Quinn, D.A.Lauffenburger, Receptor-mediated cell attachment and detachment kinetics. I. Probabilistic model and analysis. Biophys. J. 58, 841-856 (1990).

V. Cristini, J. Blawzdziewicz, M. Loewenberg, An adaptive mesh algorithm for evolving surfaces: simulations of drop breakup and coalescence. J. Comp. Phys. 168:445-463 (2001).

V. Cristini, H.B. Frieboes, R. Gatenby, S. Caserta, M. Ferrari, J. Sinek, Morphologic instability and cancer invasion. Clin. Can. Res. 11, 6772-6779 (2005).

D.J.A Crommelin, H. Schreier, Liposomes. In: Colloidal Drug Delivery Systems, J. Kreuter, ed. (Marcel Dekker, Inc., New York/Basel/Hong Kong, 1994).

W.S. Dalton, S.E. Salmon, Drug resistance in myeloma: mechanisms and approaches to circumvention. Hematol Oncol Clin North Am 6, 383-393 (1992).

P. Decuzzi, Lee S, Decuzzi M, Ferrari M, Adhesion of microfabricated particles on vascular endothelium: a parametric analysis, Ann. Biomed. Eng. 32, 793-802 (2004).

P. Decuzzi, Lee S, Bhushan B, Ferrari M, theoretical model for the margination of particles within blood vessels. .Ann. Biomed. Eng. 33, 179-90 (2005).

P. Decuzzi, F.Causa, M.Ferrari, P.A. Netti, The Effective Dispersion of NanoVectors within the MicroVasculature. Ann. Biomed. Eng. 34, 633-41(2006).

P. Decuzzi, M. Ferrari, The adhesive strength of multivalent biomimetic surfaces. In preparation.

M.W. DeGregorio, G.M. Lui, B.A. Macher, J.R. Wilbur, Uptake, metabolism, and cytotoxicity of doxorubicin in human Ewing's sarcoma and rhabdomyosarcoma cells. Cancer Chemother. Pharmacol. 12, 59-63 (1984).

E.J.F. Demant, E. Friche, Kinetics of anthracycline accumulation in multidrug-resistant tumor cells: relationship to drug lipophilicity and serum albumin binding. Biochem. Pharm. 56, 1209-1217 (1998).

M. Demoy, J.P. Andreux, C. Weingarten, B. Gouritin, V. Guilloux and P. Couvreur, Spleen capture of nanoparticles: influence of animal species and surface characteristics. Pharm. Res. 16, 37 (1999).

B. Desoize, J. Jardillier, Multicellular resistance: a paradigm for clinical resistance? Crit. Rev. Oncol. Hematol. 36, 193-207 (2000).

D. Dong, B. Ko, P. Baumeister, S. Swenson, F. Costa, F. Markland, C. Stiles, J.B. Patterson, S.E. Bates, A.S. Lee AS, Vascular targeting and antiangiogenesis agents induce drug resistance effector GRP78 within the tumor microenvironment. Cancer Res 65, 5785-5791 (2005).

M.S. Dordal, J.N. Winter, A.J. Atkinson Jr., Kinetic analysis of P-glycoprotein-mediated doxorubicin efflux. J. Pharm. Exp. Ther. 263, 762-766 (1992).

M.S. Dordal, A.C. Ho, M. Jackson-Stone, Y.F. Fu, C.L. Goolsby, J.N. Winter, Flow cytometric assessment of the cellular pharmacokinetics of fluorescent drugs. Cytometry 20, 307-14 (1995).

R. Duncan, The dawning era of polymer therapeutics. Nature Rev. Drug Discovery 2, 347-360, (2003).

R.E. Durand, Flow cytometry studies of intracellular adriamycin in multicell spheroids in vitro. Cancer Res 41, 3495-3498 (1981).

R.E. Durand, Slow penetration of anthracyclines into spheroids and tumors: a therapeutic advantage? Cancer Chemother Pharmacol 26, 198-204 (1990).

R.E. Durand, S.L. Vanderbyl SL, Tumor resistance to therapy: a genetic or kinetic problem? Cancer Commun 1, 277-283 (1989).

A.W. El-Kareh, T.W. Secomb, A mathematical model for cisplatin cellular pharmacodynamics. Neoplasia 5, 161-169 (2003).

A.W. El-Kareh, T.W. Secomb, Two mechanism peak concentration model for cellular pharmacodynamics of doxorubicin. Neoplasia 7, 705-713 (2005).

M. Erlansson, E. Daniel-Szolgay, J. Carlsson, Relations between the penetration, binding and average concentration of cytostatic drugs in human tumour spheroids. Cancer Chemother. Pharmacol. 29, 343-353 (1992).

E.A. Evans, Detailed mechanics of membrane-membrane adhesion and separation. I. Continuum of molecular cross-bridges. Biophys. J. 48, 175-83, (1985).

E.A. Evans, A. Leung, D. Hammer, and S. Simon, Chemically distinct transition states govern rapid dissociation of single L-selectin bonds under force. PNAS USA 98, 3784, (2001).

S.S. Feng, S. Chien, Chemotherapeutic engineering: application and further development of chemical engineering principles for chemotherapy of cancer and other diseases. Chem. Eng. Sci. 58, 4087-4114 (2003).

P.M. Fernandez, S.O. Tabbara, L.K. Jacobs, F.C.R. Manning, T.N. Tsangaris, A.M. Schwartz, K.A. Kennedy, S.R. Patierno, Overexpression of the glucose-regulated stress gene GRP78 in malignant but not benign human breast lesions. Breast Cancer Res Treatment 59, 15-26 (2000).

M. Ferrari, Cancer Nanotechnology: Opportunities and Challenges. Nature Rev. Cancer 5, 161-171 (2005).

G. Francia, S. Man, B. Teicher, L. Grasso, R.S. Kerbel, Gene expression analysis of tumor spheroids reveals a role for suppressed DNA mismatch repair in multicellular resistance to alkylating agents. Mol. Cell. Biol. 24, 6837-6849 (2004).

A.J. Franko, R.M. Sutherland, Oxygen diffusion distance and development of necrosis in multicell spheroids. Radiat. Res. 79, 439-453 (1979).

J.P. Freyer, R.M. Sutherland, Selective dissociation and characterization of cells from different regions of multicell tumor spheroids. Cancer Res. 40, 3956-3965 (1980).

H.B. Frieboes, J.P. Sinek, O. Nalcioglu, J.P. Fruehauf, V. Cristini, Nanotechnology in cancer drug therapy: a biocomputational approach. In: A.P. Lee, J. Lee, M. Ferrari, eds., BioMEMS and Biomedical Nanotechnology, Vol. 1: Biological and Biomedical Nanotechnology (Springer-Verlag, New York, 2006) , p. 441-466.

H.B. Frieboes, X. Zheng, C.-H. Sun, B. Tromberg, R. Gatenby, V. Cristini, An integrated computational/experimental model of tumor invasion. Cancer Res. 66, 1597-1604 (2006).

H.B. Frieboes, S.M. Wise, V. Cristini, Three-dimensional diffuse-interface simulation of multispecies tumor growth—II: Investigation of Tumor Invasion. Bull. Math. Biol. (In review).

H.B. Frieboes, J.P. Fruehauf, R. Gatenby, V. Cristini, An integrated computational/experimental model of tumor drug response. (In preparation).

A. Friedman, F. Reitich, Analysis of a mathematical model for the growth of tumors. J. Math..Biol. 38, 262 (1999).

S.A.W. Fuqua, S. Oesterreich, S.G. Hilsenbeck, D.D. Von Hoff, J. Eckardt, C.K. Osborne, Heat shock proteins and drug resistance. Breast Cancer Res Treat 32, 67-71 (1994).

A. Gabizon, H. Shmeeda, A.T. Horowitz and S. Zalipsky, Tumor cell targeting of liposome-entrapped drugs with phospholipid-anchored folic acid-PEG conjugates. Adv. Drug Deliv. Rev. 56, 1177 (2004).

W.F. Ganong, Review of Medical Physiology (Lange Medical Books/McGraw-Hill, Medical Publishing Division, New York, 2003).

S.N. Gardner, Modeling multi-drug chemotherapy: tailoring treatment to individuals. J. Theor. Biol. 214, 181-207 (2002).

R.A. Gatenby, E.T. Gawlinski, A.F. Gmitro, B. Kaylor, R.J. Gillies, Acid-mediated tumor invasion: a multidisciplinary study. Cancer Res. 66, 5216-5223 (2006).

G. Gazit, J. Lu, A.S. Lee, De-regulation of GRP stress protein expression in human breast cancer cell lines. Breast Cancer Res Treatment 54, 135-146 (1999).

F. Gieseler, A. Glasmacher, D. Kampfe, H. Wandt, V. Nuessler, S. Valsamas, J. Kunze, K. Wilms, Topoisomerase II activities in AML and their correlation with cellular sensitivity to anthracyclines and epipodophyllotoxines. Leukemia 10, 1177-1180 (1996).

W. Gill, A note on the solution of transient dispersion problems. Proc. R. Soc. Lond. Series A, Math. Phys. Sci. 298, 335-339 (1967).

W. Gill, R. Sankarasubramanian, Exact analysis of unsteady convection. Proc. Roy. Soc. London. 316, 341-350 (1970).

M.M. Gottesman, T. Fojo, S.E. Bates, Multidrug resistance in cancer: role of ATP-dependent transporters. Nat Rev Cancer 2, 48-58 (2002).

T. Graeber, C. Osmanian, T. Jacks, D.E. Housman, C.J. Koch, S.W. Lowe, A.J. Giaccia, Hypoxia-mediated selection of cells with diminished apoptotic potential in solid tumors. Nature 379, 88–91 (1996).

C.H. Graham, H. Kobayashi, K.S. Stankiewicz, S. Man, S.J. Kapitain, R.S. Kerbel, Rapid acquisition of multicellular drug resistance after a single exposure of mammary tumor cells to antitumor alkylating agents. J Natl Cancer Inst 86, 975-982 (1994).

H.P. Greenspan, Models for the growth of a solid tumor by diffusion. Stud. Appl. Math. 52, 317-340 (1972).

H.P. Greenspan, On the growth and stability of cell cultures and solid tumors. J. Theor. Biol. 56, 229–242 (1976).

V. Gupta, S.S. Nivarthi, D. Keffer, A.V. McCormick, H.T. Davis, Evidence of single-file diffusion in zeolites. Science 274, 164 (1996).

K. Hahn, J. Karger, V. Kukla, Single-file diffusion observation. Phys. Rev. Lett. 76, 2762-2765 (1996).

G. Hamilton, Multicellular spheroids as an in vitro tumor model. Cancer Letters 131:29-34 (1998).

D.A. Hammer, Apte SM., Simulation of cell rolling and adhesion on surfaces in shear flow: general results and analysis of selectin-mediated neutrophil adhesion. Biophys. J. 63, 35-57, (1992).

Z. Haroon, K.G. Peters, C.S. Greenberg, M.W. Dewhirst, Angiogenesis and blood flow in the solid tumors. In: B. Teicher, ed., Antiangiogenic Agents in Cancer Therapy (Humana Press, Totowa, New Jersey, 1999), pp. 3-21.

T. Higuchi, Rate of release of medicaments from ointment bases containing drugs in suspension. J. Pharm. Sci. 50, 874-875 (1961).

S.K. Hobbs, W.L. Monsky, F. Yuan, W.G. Roberts, L. Griffith, V.P. Torchilin, R.K. Jain, Regulation of transport pathways in tumor vessels: Role of tumor type and microenvironment. PNAS 95, 4607-4612 (1998).

M. Holz, A. Fahr, Compartment modeling. Adv. Drug Deliv. Rev. 48(2-3), 249-264 (2001).

V. Hug, D. Johnston, M. Finders, G. Hortobagyi, Use of growth-stimulatory hormones to improve the in vitro therapeutic index of doxorubicin for human breast tumors. Cancer Res 46, 147-152 (1986).

S.J. Hurwitz, M. Terashima, N. Mizunuma, C.A. Slapak, Vesicular anthracycline accumulation in doxorubicin-selected U-937 cells: Participation of lysosomes. Blood 89, 3745-3754 (1997).

L. Illum and S.S. Davis, Effect of the nonionic surfactant poloxamer 338 on the fate and deposition of polystyrene microspheres following intravenous administration. J. Pharm. Sci. 72, 1086 (1983).

L. Illum and S.S. Davis, The organ uptake of intravenously administered colloidal particles can be altered using a non-ionic surfactant (poloxamer 338). FEBS Letts. 167, 79 (1984).

L. Illum, M.A. Khan, E. Mak, S.S. Davis, Evaluation and Carrier Capacity and Release Characteristics for Poly (butyl 2-cyanoacrylate) Nanoparticles. Int. J. Pharm. 30, 17-28 (1986).

J. Israelachvili, Intermolecular and Surface Forces, 2nd ed. (Academic Press, New York, 1992).

T.L. Jackson, Intracellular accumulation and mechanism of action of doxorubicin in a spatio-temporal tumor model. J. Theoret. Biol. 220, 201-13 (2003).

R.K. Jain, Determinants of tumor blood flow: a review. Cancer Res. 49, 2641-2658 (1988).

R.K. Jain, Physiological barriers to delivery of monoclonal antibodies and other macromolecules in tumors. Cancer Res. (Suppl.) 50, 814s- 819s (1990).

R.K. Jain, Normalizing tumor vasculature with anti-angiogenic therapy: A new paradigm for combination therapy. Nature Med. 7, 987-989 (2001).

R.K. Jain, Delivery of molecular medicine to solid tumors: lessons from in vivo imaging of gene expression and function. J. Control. Rel. 74, 7-25 (2001).

A.P. Jekunen, D.R Shalinsky, D.K Hom, K.D.Albright, D. Heath, and S.B. Howell, Modulation of cisplatin cytotoxicity by permeabilization of the plasma membrane by digitonin in vitro. Biochem. Pharm. 45, 2079-2085 (1993).

Y. Jiang, J. Pjesivac-Grbovic, C. Cantrell, J.P. Freyer, A multiscale model of avascular tumor growth. Biophys. J. 89, 3884-3894 (2005).

N.P. Johnson, P. Lapetoule, H. Razaka, G. Villani, Biological and biochemical effects of DNA damage caused by platinum compounds. In D.C.H. McBrien and T.F. Slater, Eds., Biochemical Mechanisms of Platinum Antitumour Compounds (IRL Press, Oxford, 1986), pp. 1-28.

J. Kim, A continuous surface tension force formulation for diffuse-interface models. J. Comput. Phys. 204, 784-804 (2005).

H.J. Knowles, R.M. Phillip, Identification of differentially expressed genes in experimental models of the tumor microenvironment using differential display. Anticancer Res. 21, 2305-2311 (2001).

H. Kobayashi, S. Man, C.H. Graham, S.J. Kapitain, B.A. Teicher, R.S. Kerbel, Acquired multicellular-mediated resistance to alkylating agents in cancer. PNAS 90, 3294-3298 (1993).

K. Kosmidis, P. Argyrakis, P. Macheras, A reappraisal of drug release laws using Monte Carlo simulations: the prevalence of the Weibull function. Pharmaceutical Research 20, 988-995 (2003).

S.V. Kozin, P. Shkarin, L.E. Gerweck, The cell transmembrane pH gradient in tumors enhances cytotoxicity of specific weak acid chemotherapeutics. Cancer Res 61, 4740-4743 (2001).

J. Kreuter, Evaluation of nanoparticles as drug-delivery systems. I. Preparation methods. Pharm. Acta Helv. 58, 217 (1983).

J. Kreuter, in Drug Targeting, P. Buri and A. Gumma, Eds. (Elsevier, Amsterdam, 1985), p. 51.

J. Kreuter, Nanoparticles. In J. Kreuter, ed., Colloidal Drug Delivery Systems (Marcel Dekker, Inc., New York/Basel/Hong Kong, 1994).

V. Kukla, J. Kornatowski, D. Demuth, I. Girnus, H. Pfeifer, L.V.C. Rees, S. Schunk, K.K. Unger, J. Karger, NMR studies of single-file diffusion in unidimensional channel zeolites. Science 272, 702-704 (1996).

P. Kunkel, U. Ulbricht, P. Bohlen, M.A. Brockmann, R. Fillbrandt, D. Stavrou, M. Westphal, K. Lamszus, Inhibition of Glioma Angiogenesis and Growth in Vivo by Systemic Treatment with a Monoclonal Antibody against Vascular Endothelial Growth Factor Receptor-2. Cancer Res. 61, 6624-8 (2001).

L.A. Kunz-Schughart, Multicellular tumor spheroids: intermediates between monolayer culture and in vivo tumor. Cell Biol. Int. 23, 157-161 (1999).

D.A. La Van, T. McGuire, R. Langer, Small-scale systems for in vivo drug delivery. Nature Biotech. 21, 1184-1191 (2003).

K.R. Laderoute, B.J. Murphy, S.M. Short, T.D. Grant, A.M. Knapp, R.M. Sutherland, Enhancement of transforming growth factor-alpha synthesis in multicellular tumour spheroids of A431 squamous carcinoma cells. Br. J. Cancer 65, 157-162 (1992).

K. Lamszus, P. Kunkel, M. Westphal, Invasion as limitation to anti-angiogenic glioma therapy. Acta Neurochir. Suppl. 88, 169-77 (2003).

R. Langer, N.A. Peppas, Chemical and physical structure of polymers as carriers for controlled release of bioactive agents: A review. Rev. Macromol. Chem. Phy. C23, 61-126 (1983).

J. Lankelma, R.F. Luque, H. Dekker, H.M. Pinedo, Simulation model of doxorubicin activity in islets of human breast cancer cells. Biochim. Biophys. Acta 1622, 169-178 (2003).

A.S. Lee, Coordinated regulation of a set of genes by glucose and calcium ionophores in mammalian cells. Trends Biochem Sci 12, 20-23 (1987).

A.S. Lee, The ER chaperone and signaling regulator GRP78/BiP as a monitor of endoplasmic reticulum stress. Methods 35, 373-381 (2005).

J.S. Lee, S. Scala, Y. Matsumoto, B. Dickstein, R. Robey, Z. Zhan, G. Altenberg, S.E. Bates, Reduced drug accumulation and multidrug resistance in human breast cancer cells without associated P-Glycoprotein or MRP overexpression. J Cell Biochem 65, 513-526 (1997).

Y.J. Lee, S.S. Galoforo, C.M. Berns, J.C. Chen, B.H. Davis, J.E. Sim, P.M. Corry, D.R. Spitz, Glucose deprivation-induced cytotoxicity and alterations in mitogen-actived protein kinase activation are mediated by oxidative stress in multidrug-resistant human breast carcinoma cells. J Biol Chem 273, 5294-5299 (1998).

J.R. Less, M.C. Posner, T.C. Skalak, N. Wolmark, R.K. Jain, Geometric Resistance to Blood Flow and Vascular Network Architecture in Human Colon Carcinoma. Abstract 356, Houston, TX (1991).

L.M. Levasseur, H.K. Slocum, Y.M. Rustum, W.R. Greco, Modeling of the time-dependency of in vitro drug cytotoxicity and resistance. Cancer Res. 58, 5749-5761 (1998).

J. Li, A.S. Lee, Stress induction of GRP78/BiP and its role in cancer. Curr Molecular Medicine 6, 45-54 (2006).

L. Liotta,, E. Kohn, The microenvironment of the tumour-host interface. Nature 411, 375–379 (2001).

S.A. Maggelakis, J.A. Adam, Mathematical model of prevascular growth of a spherical carcinoma. Math Comput. Modelling 13, 23-38 (1990).

B.P. Mahoney, N. Raghunand, B. Baggett, R.J. Gillies, Tumor acidity, ion trapping and chemotherapeutics. I. Acid pH affects the distribution of chemotherapeutic agents in vitro. Biochem Pharmacol 66, 1207-1218 (2003).

N.V. Mantzaris, S. Webb, H.G. Othmer HG, Mathematical modeling of tumor-induced angiogenesis. J. Math. Biol. 49, 111-187 (2004).

F. Martin, R. Walczak, A. Boiarski, M. Cohen, T. West, C. Cosentino, M. Ferrari, Tailoring width of microfabricated nanochannels to solute size can be used to control diffusion kinetics. J. Control. Release 102, 123-133 (2005).

Y. Matsumura, M. Gotoh, K. Muro, Y. Yamada, K. Shirao, Y. Shimada, M. Okuwa, S. Matsumoto, Y. Miyata, H. Ohkura, K. Chin, S. Baba, T. Yamao, A. Kannami, Y. Takamatsu, K. Ito and K. Takahashi, Phase I and pharmacokinetic study of MCC-465, a doxorubicin (DXR) encapsulated in PEG immunoliposome, in patients with metastatic stomach cancer. Ann. Oncol. 15, 517 (2004).

S.R. McDougall, A.R. Anderson, M.A.J. Chaplain, J.A. Sherratt, Mathematical modelling of flow through vascular networks: Implications for tumour-induced angiogenesis and chemotherapy strategies. Bull. Math. Biol. 64, 673-702 (2002).

S.R. McDougall, A.R. Anderson, M.A. Chaplain MA, Mathematical modeling of dynamic adaptive tumour-induced angiogenesis: clinical implications and therapeutic targeting strategies. J. Theor. Biol. (In press).

D.A. McQuarrie, Kinetics of small systems. I, J Chem. Phys. 38, 433-435, (1963).

T. Meersmann, J.W. Logan, R. Simonutti, S. Caldarelli, A. Comotti, P. Sozzani, L.G. Kaiser, A. Pines, Exploring single-file diffusion in one-dimensional nanochannels by laser-polarized 129-Xe NMR spectroscopy. J. Phys. Chem. 104, 11665-11670 (2000).

S.M. Moghimi and T.A. Gray, Single dose of IV-injected poloxamine coated long-circulating particle triggers macrophage clearance of subsequent doses in rats. Clin. Sci. 93, 371 (1997).

R. Montesana, K. Matsumoto, T. Nakamura, L. Orci, Identification of a fibroblast-derived epithelial morphogen as hepatocyte growth factor. Cell 67, 901-908 (1991).

J. Moreira, A. Deutsch, Cellular automaton models of tumor development: a critical review. Advances in Complex Systems 5(2-3), 247-267 (2002).

W.F. Mueller-Klieser, R.M. Sutherland, Oxygen tensions in multicell spheroids of two cell lines, Br. J. Cancer 45, 256-264 (1982).

W.F. Mueller-Klieser, Multicellular spheroids. A review on cellular aggregates in cancer research. J. Cancer Res. Clin. Oncol. 113, 101-122 (1987).

P.L. Olive, R.E. Durand, Effect of intercellular contact on DNA conformation, radiation-induced DNA damage, and mutation in Chinese hamster V79 cells. Radiat Res 101, 94-101 (1985).

P.L. Olive, R.E. Durand, Drug and radiation resistance in spheroids: cell contact and kinetics. Cancer and Metastasis Rev. 13, 121-138 (1994).

A. Oloumi, S.H. MacPhail, P.J. Johnston, J.P. Banath, P.L. Olive, Changes in subcellular distribution of toposiomerase IIalpha correlate with etoposide resistance in multicell spheroids and xenograft tumors. Cancer Res. 60, 5747-5753 (2000).

A. Oloumi, W. Lam, J.P. Banath, P.L. Olive, Identification of genes differentially expressed in V79 cells grown as multicell spheroids. Int. J. Radiat. Biol. 78, 483-492 (2002).

N. Osheroff, A.H. Corbett, M.J. Robinson, Mechanism of action of topoisomerase II-targeted antineoplastic drugs. Adv Pharmacol 29, 105-126 (1994).

T.P. Padera, B.R. Stoll, J.B. Tooredman, D. Capen, E. di Tomaso, R.K. Jain, Cancer cells compress intratumour vessels. Nature 427, 695 (2004).

J.C. Panetta, A mathematical model of breast and ovarian cancer treated with paclitaxel. Math. Biosci. 146, 89-113 (1997).

H.-R. Park, A. Tomida, S. Sato, Y. Tsukumo, J. Yun, T. Yamori, Y. Hayakawa, T. Tsuruo, K. Shin-ya, Effect on tumor cells of blocking survival response to glucose deprivation. J Nat Cancer Inst 96, 1300-1310 (2004).

C. Paul, C. Peterson, G. Gahrton, D. Lockner, Uptake of free and DNA-bound daunorubicin and doxorubicin into human leukemic cells. Cancer Chemo. Pharmacol. 2, 49-52 (1979).

S. Pennacchietti, P. Michieli, M. Galluzzo, M. Mazzone, S. Giordano, P. Comoglio, Hypoxia promotes invasive growth by transcriptional activation of the met protooncogene. Cancer Cell 3, 347-361 (2003).

N.A. Peppas, Analysis of Fickian and non-Fickian drug release from polymers. Pharm. Acta. Helv. 60, 110-111 (1985).

G. Peskir, On the diffusion coefficient: the Einstein relation and beyond. Stoch. Model. 19, 383-405 (2003).

J.W. Piper, R.A. Swerlick, C. Zhu, Determining force dependence of two-dimensional receptor-ligand binding affinity by centrifugation. Biophys. J. 74, 492-513 (1998).

M.J. Plank, B.D. Sleeman, A reinforced random walk model of tumour angiogenesis and anti-angiogenic strategies. Math. Med. Biol. 20, 135-181 (2003).

M.J. Plank, B.D. Sleeman, Lattice and non-lattice models of tumour angiogenesis. Bull. Math. Biol. 66, 1785-1819 (2004).

C.P. Please, G.J. Pettet, D.L. McElwain, A new approach to modelling the formation of necrotic regions in tumours. Appl. Math. Lett. 11, 89-94 (1998).

C.P. Please, G.J. Pettet, D.L. McElwain, Avascular tumour dynamics and necrosis. Math. Models Methods Appl. Sci. 9, 569-579 (1999).

J. Poland, P. Sinha, A. Siegert, M. Schnolzer, U. Korf, S. Hauptmann, Comparison of protein expression profiles between monolayer and spheroid cell culture of ht-29 cells revealed fragmentation of ck18 in three-dimensional cell culture. Electrophoresis 23, 1174-1184 (2002).

L.M. Postovit, M.A. Adams, G.E. Lash, J.P. Heaton, C.H. Graham, Oxygen-mediated regulation of tumor cell invasiveness. Involvement of a nitric oxide signaling pathway. J. Biol. Chem. 277, 35730–35737 (2002).

A.R. Pries, T.W. Secomb, P. Gaehtgens, Structural adaptation and stability of microvascular networks: theory and simulations. Am. J. Physiol. 275, 349-360 (1998).

A.R. Pries, T.W. Secombl, Microcirculatory network structures and models. Annals of Biomedical Engineering 28, 916-922 (2000).

A.R. Pries, B. Beglin, T.W. Secomb, Structural adaptation of vascular networks: role of the pressure response. Hypertension 38, 1476-1479 (2001).

A.R. Pries, T.W. Secomb, Control of blood vessel structure: insights from theoretical models. Am J Physiol Heart Circ Physiol. 288, 1010-1015 (2005).

M. Principi, M. Italiani, A. Guiducci, L. Aprile, M. Muti, G. Giulianelli, P. Ottaviano, Perfusion MRI in the evaluation of the relationship between tumour growth, necrosis and angiogenesis in glioblastomas and grade 1 meningiomas. Neuroradiology 45, 205-211 (2003).

N. Raghunand, R.J. Gillies, pH and drug resistance in tumors. Drug Resist Updat 3, 39-47 (2000).

N. Raghunand, R.A. Gatenby, R.J. Gillies, Microenvironmental and cellular consequences of altered blood flow in tumours. Br J Radiol 76, S11-S22 (2003).

G. Rainaldi, A. Calcabrini, G. Arancia, M.T. Santini, Differential expression of adhesion molecules (CD44, ICAM-I and LFA-3) in cancer cells grown in monolayer or as multicellular spheroids. Anticancer Res 19, 1769-1778 (1999).

B. Ribba, T. Colin, S. Schnell, A multiscale mathematical model of cancer, and its use in analyzing irradiation therapies. Theor. Biol. Med. Model. 3(7) (2006).

A. Rieber, H.J. Brambs, A. Gabelmann, V. Heilmann, R. Kreienberg, T. Kuhn, Breast MRI for monitoring response of primary breast cancer to neo-adjuvant chemotherapy. Eur. Radiol. 12, 1711-1719 (2002).

J.J. Roberts, R.J. Knox, F. Friedlos, D.A. Lydall, DNA as the target for the cytotoxic and antitumour action of platinum co-ordination complexes: comparative in vitro and in vivo studies of cisplatin and carboplatin. In: D.C.H. McBrien and T.F. Slater, eds., Biochemical Mechanisms of Platinum Antitumour Compounds (IRL Press, Oxford, 1986), pp. 28-64.

E. Rofstad, H. Rasmussen, K. Galappathi, B. Mathiesen, K. Nilsen, B.A. Graff, Hypoxia promotes lymph node metastasis in human melanoma xenografts by up-regulating the urokinase-type plasminogen activator receptor. Cancer Res. 62, 1847-1853 (2002).

J.L. Rubenstein , J. Kim, T. Ozawa, M. Zhang, M. Westphal, D.F. Deen, M.A. Shuman, Anti-VEGF antibody treatment of glioblastoma prolongs survival but results in increased vascular cooption. Neoplasia 2, 306-14 (2000)

S. Rudt and R.H. Muller, In vitro phagocytosis assay of nano- and microparticles by chemiluminescence. III. Uptake of differently sized surface-modified particles, and its correlation to particle properties and in vivo distribution. Eur. J. Pharm. Sci. 1, 31 (1993).

E. Saatdjian, Transport Phenomena Equations and Numerical Solutions (Wiley & Sons, New York, 2000), Chap. 3.

P.D. Sadowitz, B.A. Hubbard, J.C. Dabrowiak, J. Goodisman, K.A. Tacka, M.K. Aktas, M.J. Cunningham, R.L. Dubowy, and A.-K. Souid, Kinetics of cisplatin binding to cellular DNA and modulations by thiol-blocking agents and thiol drugs. Drug Metabolism and Disposition 30, 183-190 (2002).

S. Sanga, J.P. Sinek, H.B. Frieboes, M. Ferrari, J.P. Fruehauf, V. Cristini, Mathematical Modeling of Cancer Progression and Response to Chemotherapy: Towards the Development of a Multiscale Computer Simulator. Expert Rev. Anticancer Ther. (In press).

P. Sapra and T.M. Allen, Ligand-targeted liposomal anticancer drugs. Prog. Lipid. Res. 42, 439 (2003).

K.S. Soppimath, T.M. Aminabhavi, A.R. Kulkarni and W.E. Rudzinksi, Biodegradable polymeric nanoparticles as drug delivery devices. J. Control. Release 70, 1 (2001).

T.W. Secomb, R. Hsu, Motion of red blood cells in capillaries with variable cross-sections. J. Biomech. Eng. 118, 538-544 (1996).

D.H. Sharp, J. Reinitz, E. Mjolsness, Multiscale modeling of developmental processes. Open Systems Information Dynamics 2, 67-76 (1993).

J.W. Shen, J.R. Subjeck, R.B. Lock, W.E. Ross, Depletion of topoisomerase II in isolated nuclei during a glucose-regulated stress response. Mol Cell Biol 9, 3284-3291 (1989).

J. Siepmann, A. Goepferich, Mathematical modeling of bioerodible, polymeric drug delivery systems. Advanced Drug Delivery Reviews 48, 229-247 (2001).

J.P. Sinek, H.B. Frieboes, X. Zheng, V. Cristini, Two-dimensional chemotherapy simulations demonstrate fundamental transport limitations involving nanoparticles. Biomed. Microdev. 6, 297-309 (2004).

J.P. Sinek, H.B. Frieboes, B. Sivaraman, S. Sanga, V. Cristini, Mathematical and Computational Modeling: Towards the development and application of nanodevices for drug delivery. In: Nanotechnology for Life Sciences, Vol. 4: Nanodevices for Life Sciences (Wiley-VCH, New York, In press).

J.P. Sinek, S. Sanga, X. Zheng, V. Cristini, Predicting drug pharmacokinetics and effect in vascularized tumors using computer simulation. (In preparation).

S.A. Skinner, P.J.M. Tutton, P.E. O’Brien, Microvascular architecture of experimental colon tumors in the rat. Cancer Res. 50, 2411–2417 (1990).

B. Sleeman, I.P. Wallis, Tumour induced angiogenesis as a reinforced random walk: modeling capillary network formation without endothelial cell proliferation. Math. Comp. Model. 36, 339-358 (2002).

C.W. Song, J.C. Lyons, Y. Luo, Intra- and extracellular pH in solid tumors: influence on therapeutic response. In: B.A. Teicher, ed., Drug resistance in Oncology (Marcel Dekker, New York, 1993), pp. 25-51.

K.S. Soppimath, T.M. Aminabhavi, A.R. Kulkarni, W.E. Rudzinski, Biodegradable polymeric nanoparticles as drug delivery devices. J. Control. Release 70, 1-20 (2001).

D.R. Spitz, J.E. Sim, L.A. Ridnour, S.S. Galoforo, Y.J. Lee, Glucose deprivation-induced oxidative stress in human tumor cells. Annals NY Acad Sci 899, 349-362 (2000).

A. Stephanou, S.R. McDougall, A.R. Anderson, M.A. Chaplain, Mathematical modeling of flow in 2D and 3D vascular networks: applications to anti-angiogenic and chemotherapeutic drug strategies. Mathematical and Computer Modelling 41, 1137-1156 (2005).

A. Stephanou, S.R. McDougall, A.R. Anderson, M.A. Chaplain, Mathematical modeling of the influence of blood rheological properties upon adaptive tumour-induced angiogenesis. Mathematical and Computer Modelling 44, 96-123 (2006).

M.-Y. Su, A.A. Najafi, O. Nalcioglu, Regional comparison of tumor vascularity and permeability parameters measured by albumin-Gd-DTPA and Gd-DTPA. Magn. Reson. Med. 34, 402-411 (1995).

M.-Y. Su, A. Muhler, X. Lao, O. Nalcioglu, Tumor characterization with dynamic contrast-enhanced MRI using MR contrast agents of various molecular weights. Magn. Reson. Med. 39, 259-69 (1998).

M.-Y. Su, H. Yu, J.-Y. Chiou, J. Wang, J.P. Fruehauf, O. Nalcioglu, R.S. Mehta, C.H. Baick, Measurements of Volumetric Changes and Vascular Changes with Dynamic Contrast Enhanced MRI for Cancer Therapy Monitoring. Tech. Cancer Res. Treat. 1, 479-488 (2002).

R.M. Sutherland, R.E. Durand, Radiation response of multicell spheroids – an in vitro tumor model. Curr. Top. Radiat. Res. Q 11, 87-139 (1976).

R.M. Sutherland, Cell and environment interactions in tumor microregions: the multicell spheroid model. Science 240, 177-184 (1988).

D.J. Taatjes, G. Gaudiano, K. Resing, T.H. Koch, Alkylation of DNA by the anthracycline, antitumor drugs adriamycin and daunomycin. J Med Chem 39, 4135-4138 (1996).

I.F. Tannock, Principles of cell proliferation: cell kinetics, 1^st ed.  In: Current Cancer Therapeutics (Princeton Academic Press, Princeton, 1994), pp. 3-13.

I.F. Tannock, Tumor physiology and drug resistance. Cancer Metas. Rev. 20, 123-132 (2001).

I.F. Tannock, C.M. Lee, J.K. Tunggal, D.S.M. Cowan, M.J. Egorin, Limited penetration of anticancer drugs through tumor tissue: a potential cause of resistance of solid tumors to chemotherapy. Clin. Cancer Res. 8, 878-884 (2002).

G. Taylor, Dispersion of soluble matter in solvent flowing slowly through a tube. Proc. Roy. Soc. London. 219, 186-203 (1953).

H.F. Teutsch, A. Goellner, W. Mueller-Klieser, Glucose levels and succinate and lactate dehydrogenase activity in EMT6/Ro tumor spheroids. Eur. J. Cell. Biol. 66, 302-307 (1995).

D. Theodorescu, C. Sheehan, R.S. Kerbel, TGF-beta gene expression depends on tissue architecture. In Vitro Cell. Dev. Biol. 29, 105-108 (1993).

R.H. Thomlinson, L.H. Gray, The histological structure of some human lung cancers and the possible implications of radiotherapy. Br. J. Cancer 9, 539-549 (1955).

N.E. Timmins, T.L. Maguire, S.M. Grimmond, L.K. Nielsen, Identification of three gene candidates for multicellular resistance in colon carcinoma. Cytotechnology 46, 9-18 (2004).

Tomida A, Tsuruo T. Drug resistance mediated by cellular stress response to the microenvironment of solid tumors. Anti-Cancer Drug Design 14, 169-177 (1999).

V.P. Torchilin, V.G. Omelyanenko, M.I. Papisov, A.A. Bogdanov, Jr., V.S. Trubetskoy, J.N. Herron and C.A. Gentry, Poly(ethylene glycol) on the liposome surface: on the mechanism of polymer-coated liposome longevity. Biochimica et Biophysica Acta 1195, 11 (1994).

V. Troger, J.L. Fischel, P. Formento, J. Gioanni, G. Milano, Effects of prolonged exposure to cisplatin on cytotoxicity and intracellular drug concentration. Eur. J. Cancer 28, 82-86 (1992).

S.D. Tröster, U. Muller and J. Kreuter, Modification of the body distribution of poly(methyl methacrylate) nanoparticles in rats by coating with surfactants. Int. J. Pharm. 61, 85 (1990).

S.D. Tröster and J. Kreuter, Influence of the surface properties of low contact angle surfactants on the body distribution of 14-C-poly(methyl methacrylate) nanoparticles. J..Microencapsul. 9, 19 (1992).

U. Trottenberg, C. Oosterlee, A. Schueller, Multigrid. (Academic Press, New York, 2001).

S. Wang, R.J. Lee, G. Cauchon, D.G. Gorenstein and P.S. Low, Delivery of antisense oligodeoxyribonucleotides against the human epidermal growth factor receptor into cultured KB cells with liposomes conjugated to folate via polyethylene glycol. Proc. Nat. Acad. Sci. 92, 3318 (1995).

J.P. Ward, J.R. King, Mathematical modelling of avascular-tumour growth. Math. Med. Biol. 14, 39-69 (1997).

J.P. Ward, J.R. King, Mathematical modelling of avascular-tumour growth II: modelling growth saturation. Math. Med. Biol. 16, 171-211 (1999).

M. Wartenberg, H. Acker, Quantitative recording of vitality patterns in living multicellular spheroids by confocal microscopy. Micron 26, 395-404 (1995).

M. Wartenberg, C. Frey, H. Diedershagen, J. Ritgen, J. Hescheler, H. Sauer, Development of an intrinsic P-glycoprotein-mediated doxorubicin resistance in quiescent cell layers of large, multicellular tumor spheroids. Int. J. Cancer 75, 855-863 (1998).

Q. Wei, C. Bechinger, P. Leiderer, Single-file diffusion of colloids in one-dimensional channels. Science 287, 625-627 (2000).

W. Weibull, A statistical distribution of wide applicability. J. Appl. Mechan. 18, 293-297 (1951).

D.J. Wilkins and P.A. Myers, Studies on the relationship between the electrophoretic properties of colloids and their blood clearance and organ distribution in the rat. Brit. J. Exp. Pathol. 47, 568 (1966).

D.J. Wilkins, The biological recognition of foreign from native particle as a problem in surface chemistry. J. Colloid. Interf. Sci. 25, 84 (1967).

D.J. Wilkins, Interaction of charged colloids with the RES. The Reticuloendothelial System and Arteriosclerosis. N.R. Di Luzio and R. Paoletti, eds. (Plenum Press, New York, 1967), p. 25.

S.M. Wise, H.B. Frieboes, V. Cristini, Three-dimensional diffuse-interface simulation of multispecies tumor growth—I: Numerical method. Bull. Math. Biol. (In review).

K. Wosikowski, D. Schuurhuis, G.J.P.L. Kops, M. Saceda, S.E. Bates, Altered gene expression in drug-resistant human breast cancer cells. Clin Cancer Res 3, 2405-2414 (1997).

H.S. Yoo and T.G. Park, Folate-receptor-targeted delivery of doxorubicin nano-aggregates stabilized by doxorubicin-PEG-folate conjugate. J. Control. Rel. 100, 247 (2004).

S.D. Young, R.S. Marshall, R.P. Hill, Hypoxia induces DNA overreplication and enhances metastatic potential of murine tumor cells. PNAS 85, 9533-9537 (1988).

S.D. Young, R.P. Hill, Effects of reoxygenation of cells from hypoxic regions of solid tumors: anticancer drug sensitivity and metastatic potential. J. Natl. Cancer Inst. 82, 338–9 (1990).

F. Yuan, M. Dellian, D. Fukumura, M. Leunig, D.A. Berk, V.P. Torchilin, R.K. Jain, Vascular permeability in a human tumor xenograft: molecular size dependence and cutoff size. Cancer Res. 55, 3752-3756 (1995).

J. Yun, A. Tomida, K. Nagata, T. Tsuruo, Glucose-regulated stresses confer resistance to VP-16 in human cancer cells through a decreased expression of DNA topoisomerase II. Oncology Res 7, 583-590 (1995).

A. Zakrzewicz, T.W. Secomb, A.R. Pries, Angioadaptation: keeping the vascular system in shape. News Physiol. Sci. 17, 197-201 (2002).

M. Zhang, Z. Yan, L.L. Chow, C.H. Wan, Simulation of drug release from biodegradable polymeric microspheres with bulk and surface erosions. J. Pharm. Sci. 92, 2040-2056 (2003).

X. Zhang, W. Wang, W. Yu, Y. Xie, X. Zhang, Y. Zhang, X. Ma, Development of an in vitro multicellular tumor spheroid model using microencanpsulation and its application in anticancer drug screening and testing. Biotechnol Prog 21, 1289-1296 (2005).

J.H. Zheng, C.-T. Chen, J.L.-S. Au, M.G. Wientjes, Time- and concentration-dependent penetration of doxorubicin in prostate tumors. AAPS Pharmsci. 3, 1-9 (2001).

X. Zheng, S.M. Wise, V. Cristini, Nonlinear simulation of tumor necrosis, neo-vascularization and tissue invasion via an adaptive finite-element/level-set method. Bull. Math. Biol. 67, 211-259 (2005).

G. Zunino, R. Gambetta, A. Di Marco, The inhibition in vitro of DNA polymerase and RNA polymerase by daunomicyn and adriamycin. Biochem Pharmacol 24, 309-311 (1975).