Correction: Interstitial Fluid Flow and Drug Delivery in Vascularized Tumors: A Computational Model

an error where incorrect symbols were used in the equations throughout the paper. Please download the PDF again to view the corrected article. The originally published, uncorrected article and the republished, corrected article are provided here for reference. Copyright: ß 2014 The PLOS ONE Staff. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.


Introduction
Cancer is a complex disease which involves phenomena across different scales from the molecular genetic level to the tissue as a whole. Cancerous cells of solid tumors have undergone mutations all of which combined lead to cancer [1]. These involve a dysfunctional control of proliferation, the ability to survive under low nutrient conditions and the stimulation of increased vascularization through angiogenesis [2]. This leads to an advantage in the competition over space and nutrients whereby cancer cells are also able to evade the immune systems which would otherwise kill malfunctioning cells. Solid tumors grow as compact masses. In order to grow larger than a few millimeters they must acquire additional nutrient supply through a blood vessel network. In response to inadequate supply cells produce signaling substances called growth factors which diffuse through the tissue and stimulate sprouting of new blood vessels from preexisting host vessels (angiogenesis). In tumors this angiogenic activity is located within a few hundred micrometers from the tumor rim. Fueling further growth, the resulting neovasculature is progressively coopted together with the original blood vessels by the expanding tumor mass while also pushing the neovascularization zone further into normal tissue. Chemical signaling by the tumor is however abnormal, leading to chaotic non-hierarchical vascular organiza-tion. Behind the invasive edge, angiogenic activity ceases. Further proliferation of endothelial cells instead leads to circumferential growth and dilated tumor vessels. Also vessel walls degenerate via detachment of structural support cells like pericytes and smooth muscle cells. In conjunction with decreased blood flow rates vessels become prone to collapse leading to large unvascularized regions. Some surviving vessels thread the tumor, distal to which the tumor tissue becomes necrotic due to the lack of nutrients. As a whole such a typical tumor vasculature is characterized by tortuous vessels, chaotic connectivity and heterogeneous distribution as well as a compartmentalization into a zone with high micro-vascular density (MVD) near the invasive edge and a rapid density drop towards the center.
The interstitial fluid (IF), which is a solution that bathes and surrounds the human cells and originates from blood plasma extravasating from capillaries through pores and intercellular clefts in the vessel wall, plays an important role in the development and treatment of tumors. Due to degenerate walls many tumor vessels are leaky leading to a stronger coupling of the interstitial fluid pressure (IFP) with the blood pressure. This leads to an interstitial hypertension which can be elevated up to the blood pressure. The resulting decreased pressure difference across vessel walls is believed pose a barrier to drug delivery due to decreased convective trans-vascular transport or even back-flow [3][4][5][6][7][8].
Moreover functional lymphatics which would normally drain the superfluous fluid are absent in most tumors aggravating the IFP increase [2]. As a result the IFP profile assumes a plateau in the center of the tumor and drops off rapidly across the boundary. This gradient drives a strong outward directed convective flow at around 0.1 mm=s. Signaling chemicals or tumor cells can therewith be transported into normal tissue or into lymphatics, promoting invasive behavior and metastatization [8,9]. Indeed high IFP is associated with a negative prognosis. High IFP also has negative implications for chemotherapeutic treatment. Through the outward convection, drug may be removed from the peripheral regions.
Mathematical modeling of interstitial fluid flow and delivery was first approached in a radially symmetric geometry with homogeneously distributed source and sink terms using a porous media model for the interstitial flow velocity [3]. This predicted in agreement with experimental results an elevated IFP profile and the corresponding velocity profile. Convective (drug) extravasation was virtually limited to a small peripheral region. Later, extensions were developed using explicit representations of the blood vessel network in two dimensions. In [10] individual vessels were arranged as rectangular grid and their blood pressure was coupled to the IFP at their lattice sites, including the effect of fluid loss through the vessel walls. A similar approach was taken in [11] but with a tumor vasculature which was generated from in-growth from two parent vessels. Except for [3] these studies did not consider drug transport. Recently, simulations of IF flow and drug transport were conducted based on imaging data from real tissues [12]. An analysis of biophysical parameters governing the distribution of the local drug concentration was performed in [13] primarily focusing on the effects of varying tissue permeabilities for diffusing drugs. The modeling incorporated a tumor vasculature, realistic tumor lesions and cellular uptake and binding. However convective transport was neglected. In [14] a model was introduced in which interstitial fluid interacts with a growing tumor, also incorporating a vascular network that evolves dynamically from an initial capillary grid. IFF and hence convective transport of macro-molecules depend crucially on the spatial distribution and strength of IF sources and sinks within the tumor, which are determined by the spatial arrangement of blood vessels together and their local blood pressure. Even when lymphatics are absent within the tumor, leaky vessels with low blood pressure represent also sinks for IFF inducing non-trivial flow patterns inside the tumor with unexpected effects for the convective transport of macro-molecules. It is clear that the predictive power of a computational model for IFF and drug delivery depends critically on the physiological relevance of the underlying model for the tumor vasculature. In the present paper we present for the first time a IFF and drug delivery study with a realistic, hierarchically organized arteriovenous initial vasculature, circumferential growth of tumor vessels and IF back flow into tumor vessels.
In earlier work we developed a mathematical model for vascularized tumor growth which involved an initial vasculature [15,16] arranged as a grid and updating rules representing angiogenesis, dilation and collapse. More recently, it was extended with an arteriovenous initial vasculature where quite realistically few arteries and veins branch out in a tree-like manner down to the lowest level where they are connected by capillaries [17][18][19]. Using this framework, biophysical aspects of tumor blood flow and the spatial distribution of tumor blood vessel were analyzed but it did not involve the IF explicitly nor the presence of drugs.
In this paper we want to compute the interstitial fluid pressure and flow within a tumor and use this information to predict drug delivery within the tumor and its various dependencies on physiological parameters. These parameters include the blood vessel network morphology as opposed to simplified vasculature models, blood flow characteristics, blood pressure, permeabilities of the vessel walls, the interstitium and lymphatic walls, the mass of drug particles, i.e. the ratio of convection vs. diffusion. The tumor phenotype that we consider is a vascularized solid tumor for example a melanoma or glioma which grow in their natural environment in a human host with characteristic features as described above. Samples of such tumors were studied experimentally in [20,21,32].
One particular question is how far an elevated IFP is an obstacle to drug delivery. The general consensus is that elevated IFP reduces the convective flux through the vessel walls, due to the lowered pressure difference. However, our results indicate that this does not necessarily need to be true and that the relation between IFP and IFF is more complex, also involving vessel wall and tissue conductivity. Moreover our model predicts Peclet numbers (ratio of diffusion to convection) of the order of 1 in the tumor periphery. There the IF flow is largest and almost perpendicular to the boundary into normal tissue. Hence neither diffusion nor convection can be neglected. Finally our model also allows to study extreme cases like the delivery of very heavy drug particles, which are transported purely by convection. Their distribution is harder to predict than for a highly diffusive drug since it is dependent on long range transport along chaotic IF flow patterns which are eventually governed by vascular morphology.
This paper is organized as follows: Our mathematical model is defined in the next section. We first define the representation of the vessel network followed by remodeling rules during tumor growth and the procedure to construct the initial network. Then we define continuous parts, including the representation of tissue, modeling of tumor growth, IF flow, chemical concentration fields involved in tumor growth and finally drug transport. Finally a brief overview of our numerical implementation and a derivation of parameters are given. The subsequent results section comprises a discussion of a typical base case including a brief presentation of results obtained from tumor growth simulations and an in-depth analysis of IF flow and drug delivery. After that various other cases are analyzed before the paper is finally concluded.

Model
The model has been developed for simulations in three dimensions in Cartesian space. It is a hybrid model with discrete (vessels) and continuous parts (everything except vessels), see Figure 1 for an illustration. Continuous distributions are defined in the spatial domain V that we choose to be a cubic box. Discrete vessels are defined on a face centered cubic (FCC) lattice L, which has 60 0 branching angles between parent and child branches. Varying branching angles would require modeling the vessel network in continuous space, which is computationally much more demanding but would not change the large scale morphology of the resulting network. V overlaps with L and both are centered at the origin. The lateral size of L is 8 mm and that of V is 4.5 mm. The size of L is chosen to be larger to reduce boundary effects. V is initially filled with normal tissue and contains a small tumor nucleus located in its center.
We study IF flow and drug delivery for static tumor configurations. This means the tumor growth is simulated up to a specific time without explicit involvement of IF or modeling effects of drugs. Then IFP and IF flow are computed, and finally the spatiotemporal distribution of the drug concentration with the tumor frozen in time. A coupling would be interesting in the context of studying various therapeutic protocols, which we defer to forthcoming publications.

Blood Vessel Network
Let V~ffvg, fngg be the graph which formally describes vessels (edges) fvg and junctions (nodes) fng. Vessels coincide with bonds from L, starting at one site of L and ending at another. They can span multiple bonds but must be straight. Each vessel carries biophysical properties like radius or blood flow. We will introduce them as needed.
Blood flow is an essential part of our model. For vessels we compute the flow rate q v (volume through its cross-section per time) and shear stress on the vessel wall f v and the blood pressure p v at the endpoints. The indices are dropped in the following. We assume ideal pipe flow within the vessels, obeying Hagen-Poisseuille's law where Dp is the pressure difference between the vessel ends, l is the vessel length, g the blood viscosity and r the vessel radius. g is composed of the blood plasma viscosity 4 : 10 {6 kPa s times the relative viscosity g r (r,H) which is a function of the hematocrit H and the radius. For g r we use a formula based on in vivo experimental data [22]. For simplicity, we assume that H~0:45, the average in the human body. Mass conservation at each node requires that the flow rates of attached vessels sum to zero: S v q v~0 (Kirchhof's law). Together with appropriate boundary conditions a system of linear equations for the nodal pressures arises which is solved numerically. As boundary condition the pressure is fixed at the arterial and venous roots of the vascular trees. These boundary pressures are determined as function of the vessel radius, also guided by experimental data [23] (see Supplement S1 (1)). Note that we do not incorporate the extravasated fluid into the mass balance, which is justified since, as we will demonstrate in the results section, the amount of extravasated liquid is orders of magnitudes smaller than the total vascular blood flow. In the rest of the paper q and f will denote the absolute value of the flow and shear force within a vessel -above they carried a sign. Blood vessel network remodeling, the process in which the hierarchically structured initial network is reorganized by the growing tumor is defined by a set of stochastic and continuous processes which model angiogenesis, dilation, degeneration and collapse. They are implemented as updating rules which are applied consecutively in each time step. As a result, vessels are created, deleted or they change their properties. These rules are adopted straight forward from the 2d case [18] and presented here again for completeness.
Our time stepping scheme advances the vessel network in fixed steps of width Dt. Assuming that the frequency of a stochastic event is determined by a rate parameter k we approximate the probability for its occurrence in one time step as p~kDt. We chose Dt sufficiently small such that pv1.The time evolution of continuous processes described by differential equations of first order in time is handled by Euler's method with time step Dt. In the following we describe the individual vascular remodeling processes that are incorporated into our model, for an illustration of theses processes (see Supplement S1).
Sprout initiation models the event when endothelial cells (ECs) leave the parent vessel in order to grow a new sprout. It is a stochastic process that adds new vessel segments to the existing network. Lattice sites occupied by the existing network are visited in random order and at each of these sites a new segment is attached with probability Dt=t (sprout) EC provided that the following conditions are met: the growth factor concentration is non-zero, the distance to the next branching point is larger than l (spr) and the time spent within the tumor is less than t (switch)

EC
. The new segments are created along neighboring lattice edges where the growth factor gradient is maximal and where no other vessels are already present. A vessel is tagged as ''within the tumor'' if at least one of the endpoints is within the tumor, which is true where the level set function h(x)v0 (see below). Vessels also have a property which can tag them as sprouts and tell for how long they have been sprouts. We denote this ''life-time as sprout'' as t.
Sprout migration is the process in which initial sprout vessels continue to grow. The probability is Dt=t (sprout) EC for vessels which are tagged as sprout. A growth event is realized by appending a vessel segment along a single lattice edge in the same direction as the existing sprout. Sprout vessels are untagged and become normal vessels if the tip fuses with another vessel such that blood can flow or if their twt (migr) EC , where t (migr) EC is a parameter which defines the maximal sprout growth time. If the tip fuses with another sprout without creating a conducting branch then it remains tagged as sprout. Sprout initiation can also start from sprouts which emulates tip splitting as observed in-vivo and invitro. Sprouts are excluded from the collapse, degeneration and circumferential growth mechanisms.
Wall degeneration models the detachment and disintegration of cell layers and membranes around the vessel lumen. Therefore we implement the property w which reflects the vessel wall thickness for normal vessels, and continuously decreases for tumor vessels with the rate Dw until zero. For values smaller than e.g. the size of an EC, w becomes an abstracted representation of the stability and tightness which the remaining EC layer provides. w is initialized (sprouts and initial vessels) with the wall thickness of real healthy vessels in dependence on their radius [24] (see Supplement S14).
Vessel collapse models pinch off of blood flow and complete disintegration of the vessel. It is a stochastic process where a vessel can be removed with probability Dt=t (coll) EC under the condition that its wall stability variable wvw (coll) and its wall shear stress f vf (coll) . Thereby r (max) , t (switch) EC , w (coll) and f (coll) are model parameters. Recently the Ang-Tie system was modeled in a similar context [25]. This is straight forward to include in future work. Here we model the effects of it, rather than the system directly.
Vessel dilation models the switch to circumferential growth within the tumor [21]. During circumferential growth the vessel radius increases continuously with the rate Dr. The requirement is hereby that rvr (max) , the average growth factor concentration over the segment is non-zero and the time spent within the tumor is larger than t (switch)

Initial Blood Vessel Network Construction
To our knowledge there are no data sets from real networks available that cover a few millimeters of tissue and represent the complete vasculature including micro vessels in a form which is convertible to a ''network of pipes'' as required for our modeling purposes. Therefore we decided to generate it algorithmically. Our aim is to maximize the lattice occupation with a network which exhibits a hierarchical topology and homogeneously distributed capillaries.
A well known method is constraint constructive optimization (CCO) [26] in which a tree is grown by successively adding branch segments at locations given by some optimality criterion e.g. minimal total surface area. The constraints are that there is no geometrical overlap of the branch segments and that new segments must reach certain previously unperfused tissue blocks. The radii at branching points r b ,r c behave according to Murray's law [27] r a a~r a b zr a c , where r a is the radius of the parent branch and a is an exponent. a has been found to range between 2:55 and 3. The latter is a common choice which is also taken here.
In [28] a method was presented in which vascular trees on a lattice are stochastically grown and remodeled. This has the advantage of being relatively simple and being capable of building connected networks comprising arteries, capillaries and veins, whereas in CCO one typically has ''dangling'' terminal branches where capillaries should connect to. In [18] we adopted this method in order to obtain initial networks for 2d simulations. Later we applied it on a cubic lattice [19] and here we apply it on a FCC lattice. An illustration of the steps can be found in Supplement S1.
The initial network construction is based on a relatively coarse lattice with a lattice constant that corresponds to the mean intercapillary distanceh h L~8 0mm. After the construction of the network on the coarser network it has to be mapped on the finer lattice withh h L~1 0mm for use in the subsequent simulation. FCC lattices can be subdivided not unlike cubic lattices, meaning that sites and edges of the coarse lattice coincide with site and edges on the fine lattice. The tedious details are omitted here.
The construction is initialized by placing nodes which serve as roots for the trees onto boundary sites of the lattice. The type of these nodes is either arterial or venous, placed in alternating order.
The subsequent construction is then carried out in two stages. In stage one, trees are grown by a stochastic process in which ''structural elements'' are successively appended to one of the current tree leafs. As structural element we take either single vessels or a Y-shaped aggregate of three vessels. The element, its orientation and the leaf are determined by randomly (see Supplement S1). Eventually the lattice is filled but arterial and venous side are not interdigitating sufficiently to yield a homogeneous capillary distribution. This is corrected in a second remodeling stage. Capillaries are temporarily inserted in-between neighboring arterial and venous terminals. We set the capillary radii to 4mm. Radii of terminal branches are set to 4:5mm for arteries and 5mm for veins. Radii of higher level vessels are determined by Murray's law. As a result an intermediate functional vascular tree is obtain for which blood flow is computed. Shear-stress dependent growth and shrinkage is carried out by stochastic removal and attachment of vessels from or to terminal branches. High shear-stress means higher probability to grow and vice versa. The idea originates from the observation that high shear stress indeed promotes vessel survival and stability [29]. We repeat this stage until the number of capillaries reaches a plateau. Trees can potentially grow from each of the root nodes. A few of these trees establish themselves while most of them regress and disappear.
For this paper we extended the algorithm from [18,19] with an ''outer'' loop producing increasingly fine resolved networks in a hierarchical fashion. This effectively reduces the tortuousity of major vessels. The first level (coarsest network) is constructed as described above. Then the lattice is refined, halving the lattice spacing and doubling the number of sites in each direction. Arteriovenous trees are kept in place and capillaries are discarded. Each vessel segment now occupies two lattice bonds. The lattice spacing is then reset to its former value. Hence, the spatial extend and segment lengths are effectively doubled. Now the random growth and remodeling steps are executed as above, where the previous terminal nodes now serve as new roots. This up-scaling and growth procedure is repeated a preset number of times. The results shown here were generated from a 25|31|31 base lattice and 2 up-scaling steps.

The Continuum Model for Tissue
Our tissue model is based on the framework developed in [30] which describes the tissue as a mixture of various tissue constituents. A mathematical model is formulated in terms of smoothed fields of quantities such as density, velocity, stress, etc. Several constituents can coexist at one material point due to smoothing. Assuming incompressibility, one can describe the composition in terms of volume fractions with the constraint that the fractions sum up to one at every point in space. For brevity, we just give the final set of equations. A derivation can be found in [30], see also [31]. The result is a system of partial differential equations of the diffusion convection reaction type.
First, let us denote tissue constituents and their volume fractions:  T zw D , inV T , w N zw D , inV\V T This is analogous to immiscible liquids, where cell-cell adhesion forces correspond to the atomic forces in the liquids. We assume however that the adhesion forces are very weak, which allows us to neglect the surface tension term which would normally appear in the momentum balance equation. It will be included in future work.
All w constituents move with the same velocity field v w which is driven by the gradient of a solid pressure (the isotropic component of the stress tensor) p w . It is based on the assumption that inertia is dominated by friction against a rigid ECM through which cells flow like a liquid through a porous medium. Therefore we have where (1) is the condensed momentum balance, K w is a mobility constant and Q w and Q D are source terms to be defined below.
Note that this set of equations is applicable to tumors that have a clearly delineated rim as for instance rat C6 gliomas and human glioblastomas [21], human malignant melanoma [20], leiomyomata [32], etc. It is not valid for non-solid cancers like Leukemia and highly invasive tumors which do not have such a clearly delineated rim. The motion of LV T is formally defined by (4). In practice we use the level set method [33] to represent V T and LV T and introduce an auxiliary field h(x) which gives the closest distance to LV T . It is signed so that h(x)v0 for x[V T . Over time it evolves according to the advection equation We can now define w T~( 1{H h (h))(w{w D ) and w N~Hh (h)(w{w D ), where h is the lattice spacing of the numerical grid (see below) H denotes a smoothed Heaviside step function with width (see Supplement S1).
For the pressure we take p w (w)~max½E(w{w w 0 ),0: For simplicity and the lack of better knowledge we use a linear elastic law with elastic modulus E. Also, we assume that cells do not exert pressure upon each other when w is less than the volume fraction in a fully relaxed statew w 0 .
The source terms are composed of contributions from T,N and D as follows where Q + a stands for proliferation and apoptosis of phase a and Q aD stands for necrosis.
We assume proliferation depends on packing density [34], i.e. volume fraction w, and on available nutrients c o . Cells do not proliferate in regions with high density where apoptosis reduces the density towards the so called homeostatic (equilibrium) density w 0,a . At w~w 0,a , and for sufficiently high c o , apoptosis and proliferation rates cancel so that the net production rate Q z a zQ { a vanishes. Moreover Q z a varies linearly with w{w 0,a . Under low nutrient conditions proliferative activity stops, i.e. Q z a~0 for c o vf z , where f z is a threshold parameter. Consequently, apoptosis and possibly necrosis reduce the cell density. The difference between these events is that apoptosis leaves no debris as cells are deconstructed in an orderly fashion, i.e. the respective cellular material vanishes. Necrosis occurs under very low nutrient conditions if c o vf x vf z , where f x is also a threshold parameter. The fraction of cells undergoing necrosis is transferred to w D via the rates Q x a . In total this behavior is summarized in the following formulas: for a~T,N, where c z a , c { a and c x a are constant rate coefficients (proliferation, apoptosis and necrosis), s w determines the sensitivity to density variations and H is the Heaviside function. Note, the use of ''min'' in conjunction with the Heaviside function. It limits the proliferation rate by either c z a or 0 (no proliferation) depending on nutrients.

Interstitial Fluid Flow
IF is commonly modeled as a liquid within a porous medium, e.g. [11,12,14,35]. We follow this approach and assume that cells and ECM collectively constitute the porous medium. Here we consider only stationary states, with a static tumor and a rigid medium, thus Ll=Lt~0. Mass balance for the IF fraction l requires that with its velocity v l and source distribution Q l . Neglecting inertia terms one obtains the momentum balance equation where T l is the stress tensor of the liquid, (+) ij~L T ij =Lx j and f l an interaction force with the other constituents. We use the results in [31] and [30] where constitutive relations for T l and f l were obtained for the case of a solid-fluid mixture. Assuming that the interstitial fluid is an ideal inviscid liquid, its stress tensor consists only of the contribution from the interstitial fluid pressure pressure p i or in the following just p.
The interaction force is defined in such a way that we later obtain a variant of Darcy's law The first term represents friction with cells and ECM fibers, where K l is a tissue dependent permeability coefficient. Substitution of equations (8) and (9)  which leads to an elliptic equation for the pressure Note that l 2 K l is the classical conductivity of the porous medium. We define K l~Kl (h) so that it smoothly interpolates between parameter values for tumor K l,T and normal tissue K l,N . K l,T and K l,N are chosen so that the conductivity in the bulk assumes experimentally determined values. Note that l is almost constant distal to the tumor boundary and varies over a small value range since w [ ½0:5,0:6. The source term is composed of contributions from the vessel network Q lv and lymphatic sinks Q ll so that Q l~Qlv zQ ll . Both are determined by the flux across the channel walls. For vessels, this flux is driven by the pressure difference p v {p i and an osmotic contribution s½p v {p i (Starlings equation) [36]. For lymphatics we assume an analogous relation but neglect osmosis due to the lack of data.
where p L is the lymphatic pressure, L (v) l and L (L) l are permeabilities, S=V and S (L) =V are the channel surface area densities per volume and p i and p v denote the so called oncotic pressures. s, the reflection coefficient, is a tissue specific value.
The standard approach for modeling exchange with vessels on a small scale would use boundary conditions at the vessel walls, while tessellating the surrounding space with a fine grained mesh. However this would make the large length scale which we are interested in inaccessible due to the size (we have of the order of 10 6 vessel segments). Instead we integrate the flux approximately over the vessel surfaces within each numerical grid cells and add it as source term. An approximation inherent to this method is that the space covered by the vasculature is not excluded from the interstitial space.
Hence (12) is not applied in this exact form. The source flux is implemented as superposition of smoothed delta functions d e (see Supplement S1 ). Their locations y are generated from a stochastic uniform sampling of the surfaces of the cylindrical pipes which make up the vessel network. We write this formally as y [ v, where v symbolizes a vessel. For a numerical grid cell with index i and center x i , Q lv then becomes where h is the grid spacing, L v l the wall permeability, S v,y is the area corresponding to a sample on v, and p v (y) is the blood pressure in v at position y, linearly interpolated between the nodes. Different degrees of vessel leakiness are incorporated based on the maturity state w. We assume that w reflects the vessel wall thickness for sufficiently large vessels and that the wall's resistance (L v l ) {1 increases proportionally to the wall's thickness. This eventually leads us to where l l,T and l l,N are experimentally determined permeabilities for capillaries in tumor and normal tissue, respectively, and w (init) (r) is a formula based on experimental data [24] from which we obtain the physiologically normal thickness of the vessel wall depending on the radius (see Supplement S1 (2)). For small w the identification with the wall thickness breaks down and it becomes a mere abstract quantity inversely related to the amount of leakiness. In order to obtain realistic permeabilities for tumor vessels as well, we are therefore free to bound L v l from above by l l,T . Lymphatics on the other hand are modeled as continuous sink distribution, where their surface area S (L) depends on the tissue type via h analogous to K l and moreover L (L) l and p L are assumed to be constant. Hence we can use (13) directly in the numerical implementation.

Chemical Concentration Fields
The basis for the description of dissolved chemicals is the following diffusion convection reaction equation which determines the evolution of the concentration c (a) in constituent a [30].
where D (a) c are effective diffusion constants (assumed to be scalar) and Q (a) c a source term. For nutrients and growth factors, we approximate the concentration as the equal in all phases c (a) :c under the assumption that the exchange among constituents is very fast. Then, summation of (16) over all a gives where D c is the composite effective diffusion coefficient, and v c the composite velocity of the mixture. In the following we will use subscripts to c to denote specific chemical species: c o denote nutrients, c g are growth factors. For drug we distinguish concentrations in two different compartments s i :~c d,i for which i~1,2 denote the extra-and intra cellular space, respectively. Nutrients are represented by the most prominent one, namely oxygen with its concentration c o . The time scale on which c o relaxes after changes is negligible, of the order of seconds, and thus c o is assumed to be always in a quasi stationary state, instantaneously adapting to changes in the system. Convection is neglected due to the dominance of diffusion. Consequently we obtain where we already replaced Q with a particular form of the source term: The second term represents consumption with the tissue dependent rate c { o . The third term represents the diffusive flux across the vessel wall, which we treat analogous to the interstitial fluid source term (14), only L v l is replaced by L v o , and p v by the blood oxygen concentration c v o . Since we already assumed that the hematocrit is constant over the whole vasculature, we further assume for simplicity that the oxygen concentration c v o also constant over the perfused parts and zero in unperfused vessels.
Growth factors are collectively represented as a single diffusible species with its concentration c g . A prominent representative is VEGF which is over expressed in under-oxygenated tumor cells. We assume a constant production rate by tumor cells in locations where c o vf z and that it diffuses, binds and degrades everywhere equally. Instead of solving a diffusion equation we use a simpler and faster approximation based on a Green's function approach: every source site generates a linearly decaying contribution to c g with the cutoff or diffusion radius R g . Thus we have where we define G(x)!max(0,1{x=R g ), with a normalization constant so that Ð G(jxj)~1. Note that consequently, by definition of the angiogenesis rules, sprouting occurs within R g of oxygen deprived TCs and a c g gradient arises along which sprouts are oriented.
Transport and uptake of drug is modeled as diffusion advection process in the interstitial fluid and sequestration into the cell constituent. We distinguish between the concentrations s 1 in the IF and s 2 within cells as average over the solvent volume. The tissue volume average reads s~s 1 lzs 2 w with the volume fractions l and w as defined above. Following (16), we define specialized mass balance equations as l Ls 1 Lt z+ : (s 1 lv l )~+ : (lD s +s 1 )zQ s 1 zQ 12 ð20Þ w where Q 12 is the exchange rate between the two compartments, Q s 1 the source contributions from vessels and lymphatics and D s the diffusion coefficient in the IF. For a simple derivation of Q 12 , we assume the total flux of molecules across the cell-fluid interface area S within some volume V has the form J~S(k k 21 s 2 {k k 12 s 1 ) with the rate constantsk k ij which model the combined effect of diffusion through the cell membrane and intracellular binding and unbinding. We write S in terms of the single cell volume V c and surface area S c as S~lwVS c =V c , assuming that only the fraction l of the cell surface is in contact with the IF. Then we obtain with Furthermore the contributions from vascular and lymphatic exchange are given by where S has the original meaning of vessel surface area again. The diffusive permeability L v s is defined exactly like L v l in (14) and (15) with correspondingly exchanged subscripts including the permeabilities of tumor vessels l s,T and normal capillaries l s,N . Q z lv stands for extravasated fluid volume per mixture volume carrying the concentration s v which is the concentration within the vessels. We assume that s v is homogeneous over the whole network but time dependent where the dependency is given as closed formula e.g. an exponential decay after a hypothetical injection at t~0. Q { lv stands for fluid uptake by vessels. Analogously Q { ll for uptake through lymphatic. Since these terms represent flow out of the interstitial space, they are multiplied by s 1 in order to obtain the respective solute flux. We could define a Q z ll for symmetry but in practice fluid always flows into lymphatics, never in the opposite direction. We treat Q + lb analogously to Q lb , for b~l,v in (12) and (13) with the exception that only contributions are added where the blood or lymphatic pressure is larger (z) or lower ({) than the IF pressure. Indeed Q z lv zQ { lv zQ z ll zQ { ll~Q l .

Numerical Implementation
Solutions to partial differential equations are computed by finite difference schemes on a regular uniform staggered grid [37]. Numerical values for concentrations, volume fractions, etc. are defined on grid cells, while velocities and fluxes are defined on faces. The grid spacing h is 30 mm which corresponds approximately to two to three tissue cells. The diffusion terms are discretized by standard 9 point centered difference stencils. All system of linear equations are solved with (algebraic multigrid -if needed) preconditioned conjugate gradient method. Specifically, we use the solvers in [38]. Advection terms are treated by a central scheme for conservation laws [39]. In the time, the operator splitting technique [37] allows treatment of various sub-systems separately, i.e. sub-systems are advanced one by one, always using the newest available state. The cell volume fractions w and w D are updated simultaneously with the 2nd order improved Euler method. The level set function h is updated likewise. The vessel network is updated in 1 hour steps. In these periods for w,w D and h smaller sub-steps must be taken the length of which is dictated by the stability conditions of the time integration methods. In practice these steps are about 0:02h wide. Sometimes h must be ''redistanced'' in order to restore the distance function property j+hj~1. The WENO method presented in [40] works very well for our purposes. The computation of c g and c o as well as redistancing are not performed every step. We determine the time between updating these fields by the time it takes tissue cells to cross a numerical grid cell and also by the time scale of the source term, which gives the time min (h=v w ,(LQ w =Lt) {1 ), where we take the minimum over space and time since the last update. The numerical solution of the drug concentrations s 1 and s 2 is carried out using the central advection scheme [39] and the improved Euler method.

Parameters
A list of parameters for our base case system can be found in tables 1 and 2. The sprouting parameters t (sprout) EC and t (migr) EC are estimated from in-vitro endothelial cell (EC) migration experiments in [41]. It is known that angiogenesis is inhibited in ECs near existing branching points. For this l (spr)~2 0mm seem reasonable, which are about two nearby ECs. The vessel dilation rate Dr and maximal radius r (max) was extracted from [21] where the spatial compartmentalization of human melanoma was described. The wall thickness w is initialized at t~0 depending on the vessel radius (see Supplement S1 (2)) guided by experimental data [24]. The wall degradation rate was estimated from [21] based on tissue section at increasing stages of tumor progression. We simply observed how long it takes until the supporting wall structures of a vessel with a certain radius are removed. For the critical collapse shear stress f (coll) we assumed that it is a low percentage of the average shear stress within the system, also guided by comparison of predicted vascular density levels with data from [20].
The oxygen level in our model is dimensionless and normalized to 1 which is the level within vessels. We divide the diffusion equation (18) The precise number is arbitrary and non-crucial but reflects that tumor cells have a higher oxygen consumption rate leading to decreased tissue oxygen levels. We then tuned l o,N =D o so that the oxygen level in normal tissue is above ca. 1=2. For simplicity we assume that the permeabilities in tumor and normal tissue are equal l o,T~lo,N .
We assume that tumor and tissue cells have the same fastest possible proliferation rate (c z N and c z T ) of once per day. The time to live for normal cells is assumed to be 10 days after which they undergo apoptosis, yielding c { N . Tumor cells have acquired mutations which enable them to circumvent the apoptotic mechanisms. Therefore we set c { T~0 . Cells under severe hypoxia are assumed to die off relatively quickly within 48 h (c x T and c x N ) and become necrotic tissue.
The oxygen threshold below which cells become necrotic is f x~0 :03. Cells stop proliferating when the oxygen level is below f z~0 :3 which is ca. 60% of the lowest level in normal tissue. Since only tumor cells are ever exposed to low oxygen levels, we also do not distinguish between tumor and normal tissue here.
Our cell volume fractions reflect a high-density prototypical tissue. We assume that tumor cells have become less sensitive to solid pressure from nearby cells and so their homeostatic level is 0.6 (w 0,T ) while it is 0.4 (w 0,N ) for normal cells. Here we follow [34], where the idea for this pressure regulated proliferation originated. A similar model was also employed in [42] but with different parameters.
Parameters for interstitial fluid flow and drug transport are summarized in table 2. The permeability coefficients K l,a and l l,a as well as osmosis parameters p i , p v and s are obtained from [35] and the references therein. Where the actual permeability is provided, e.g. l 2 K l , we compute the coefficient by dividing with the typical l&0:2 in the respective tissue. For lymphatics we assume a wall permeability (L (L) l ) which is of the same order of magnitude as for capillaries (l l,N ). The lymphatic surface area per volume (S (L) N =V ) is estimated by assuming that there is a grid-like network with a channel every 100mm and a radius of 10mm. We leave the tumor without lymphatics (S (L) T~0 ), since these are absent or dysfunctional in tumors (see [5] and the references therein).
The drug distribution model is guided by experimental data on the pharmacokinetics of Doxorubicin, which has been used for a long time to treat various cancers. For the diffusion coefficient D s and exchange rates k ij we follow [13] who presented a similar model with additional cellular compartments. The diffusion constant stems from [43] where it was estimated as ca. 1000mm 2 =min, which we use as well. Given an isolated system with two compartments and transition rates k 12 and k 21 , the steady state concentration ratio equals the ratio of the rates. In experiments with cell cultures [44], intra-cellular to medium ratios of ca. 100 were observed. It seems reasonable to associate k 12 , the cell uptake rate with diffusion through the vessel wall. In [13] this rate was estimated as 5:4min {1 . Thus we simply use Note that k 21 is close to the estimated lysosomal release rate which is also the slowest rate in the model proposed in [13], so this may be identified as bottleneck for release. It remains to determine the vascular permeabilities. We estimate these based on the diffusivity in plasma. As an approximation we write the permeability of a planar sheet of thickness L as D=L, where D is the diffusion constant. We take 180mm 2 =s for D in plasma and we identify L~10mm with the thickness of the capillary wall. We assume that drug only diffuses through the gaps between endothelial cells (ECs). Assuming that in very leaky tumor vessels, there would be a circular gap with 1:5mm radius per EC [35], we arrive at the fraction of gaps over the vessel surface G~0:017, assuming 10 2 mm 2 per EC, thus l s,T~D G=L.
Assuming the difference between tumor and normal capillaries is due to leakiness, we determine l s,N by requiring that the ratios are equal: l s,N =l s,T~ll,N =l l,T .

Tumor Growth
Snapshots from a simulation are displayed in Figure 2. We performed 15 simulation runs, producing 15 final states which differ in their initial blood vessel networks (and the seeds for the random number generator used for the stochastic events during the simulation). For a video visualizing the spatiotemporal evolution of the model see Supplements S14 and S14.
Initially, the tumor is prepared as a small sphere in which the tissue consists of tumor cells instead of normal cells. We define the distance function h at t~0 as the signed distance from the sphere boundary. The tumor is located in the center of the simulation box and has a radius of 0.5 mm. Increased oxygen consumption leads to decreased oxygen levels within the tumor which leads to expression of growth factors which again stimulates angiogenesis within R g . Eventually blood-perfused neovasculature raises the oxygen level in the tumor periphery and enables further tumor growth. The first snapshot in Figure 2 shows the system after 100 h. At this point the system is in a state displaying the typical compartmentalization into high micro vascular density (MVD) rim, decreasing MVD toward the tumor center, isolated vessels threading the tumor, necrotic regions associated with unvascularized regions and tumor proliferation confined to its rim. The tumor continues to grow by vascularizing and pushing into the surrounding tissue, leaving a torturous chaotic tumor network behind. The final snapshot is taken at t = 700 h where the tumor has reached the edge of the continuum domain V. Its final radius is ca. 2.5 mm. By design of the tumor-vessel interactions similar observations were reported in earlier work in [15], [16], [18], where much simpler tumor models were used. See also Figure 3 where important system variables as a function of the distance from the invasive edge h are shown.
To generate these plots we sorted data points based on their spatial position into bins, or shells surrounding the invasive edge according to their h value. The width of the bins is 30mm. Unless stated otherwise, we computed the averages of the binned values for each simulation run. The plotted data displays the means and standard deviations of the ensemble, not the spatial fluctuations. Spatial fluctuations can be seen in the map plots and are analyzed in more detail only for the distribution of drug. Note that we may define averages over (parts of the) vessel network formally as line integrals over the vessel center lines divided by the total length of respective parts. In practice we generate sampling points on center-lines and sort these into bins as we do with numerical grid values. This is applied e.g. in Figure 3B (vessel radius). The predictions described above are in good agreement with experimental data from [20] and [21] for human melanoma xenografts and gliomas respectively. Our new continuum model describes tissue more realistically by the incorporation of actual host tissue cells, cell motility and cell-cell adhesion. We will report a more detailed analysis of the resulting morphological aspects elsewhere.

Interstitial Fluid Flow
Pressure, velocity and source terms were computed numerically for the final tumor configurations at t = 700 h. For the IF flow studies we assumed that the tumors are static from this time on. Generally the motion of the IF is coupled to the motion of the other tissue constituents since ''empty'' spaces are filled with the IF. However in our case the velocity of the fluid is orders of magnitude larger than the velocity of cells, for which reason we can neglect these interactions. Figure 4 shows slices through simulation data of one sample. Figure 4A displays the vessel volume fraction h v . The data are generated by superposition of smoothed delta functions which are distributed stochastically within the cylindrical volumes comprising the network edges. In Figure 4C we plotted the source term Q l which is the IF volume flowing in or out of the interstitial space per volume and time. By definition, lymphatics are absent within the tumor, thus therein the only sources and sinks are blood vessels, which appear as lengthy blobs with positive (extravasation) or negative (uptake) contributions. Uptake is possible since the blood pressure can also be lower than the local IFP. At the tumor rim we see a significant amount of fluid being taken up, since there is a strong outward flux from the tumor which is absorbed into lymphatics and potentially also into parts of the neovascular plexus.
The IFP profile is elevated within the tumor and decays rapidly over its boundary. See Figure 4D and Figure 5A. The peak pressure in the tumor center is ca 6 kPa (45 mmHg). Outside it is 0.5 kPa (3.75 mmHg). We set the lymphatic pressure to 20.5 kPa, and the average blood pressure is ca 6.25 kPa (47 mmHg). In   models using pure capillary networks a pressure range from 15 to 25 mmHg is commonly either directly set or imposed via boundary conditions, e.g. in [10,12,14,15]. In the presence of higher level arteries (as in our model) a much higher IFP can be observed, because of the elevated pressure in arteries and connected vessels. Mean in vivo tumor IFPs are reported in [35] from 0 to 5 kP (38 mmHg) for various tumor types. Most tumors exhibit a high degree of heterogeneity, with deviations from the mean of 100%. The elevated IFP has a finite ''penetration depth'' into normal tissue (see below). Beyond that the IFP appears relatively homogeneous. Fine grained fluctuations are introduced by fluctuations of the capillary surface area per grid cell.
The order of magnitude of the interstitial fluid velocity is considered to be 0.1 to 1mm=s [4,[45][46][47]. In particular, it is found to be highest close to the tumor boundary where the pressure gradient is steepest, leading to a strong outward flux. Our results, shown in Figure 4A and Figure 5B are in agreement with this experimental observation. The velocity peak in the region is at 0:2mm=s.
Furthermore the velocity patterns in our data show a significant amount of fluid being transported in between tumor vessels. In Figure 4D this can be directly observed. For the whole ensemble we measured the outward projection of the velocity vector v jj :v l : +h=j+hj together with the velocity magnitude jvj:jv l j, which are plotted in Figure 5. Spatial distributions from one simulation can be seen in Figure 4E and F. Within the tumor center, v jj vanishes whereas jvj decreases to ca. 0.05 mm=s. This means that the flow becomes increasingly isotropic toward the center, where the IF apparently flows in between isolated tumor vessels instead of to the tumor boundary. This flow is more than an order of magnitude faster than IFF in normal tissue.
In the tumor periphery fluid uptake by lymphatics is significantly increased compared to normal tissue further away ( Figure 5C, Q ll ) because since lymphatic vessel are absent in the tumor, the extravasated fluid must cross the tumor boundary to be absorbed in normal tissue. The determining equation for the IF pressure p, and the equations for oxygen and growth factors, have a structure like + 2 u{cu~r, for some constant c, dependent variable u(x) and arbitrary distribution r(x). For this equation any local change in r causes an exponentially decaying disturbance in u. The length scale of the decay is 1= ffiffi ffi c p . Since c~L (L) l S (L) =V (K l l 2 ) {1 for the IFP, we see a ''penetration depth'' of 98mm.
Furthermore, in addition to the fluid that originates from the tumor interior, we find that the neovascular plexus extravasates a huge amount of fluid ( Figure 5C, Q z lv ). The collective surface area of these vessels is large due to the amount of vessels, they are very permeable and the p v {p i difference is relatively large, so this is not surprising but it has implications for the escape of tumor cells from the rim into the lymphatic system, and subsequently metastasis.
At this point we should discuss the validity of our approximation that neglects the loss of blood plasma from the vasculature due to extravasation. Apart from the study of IF flow this approximation is standard but in particular for tumor blood flow the coupling through leaky thread-like vessels it is a more severe simplification, where one would ideally solve for the IFP and blood pressure simultaneously and fully coupled. To clarify this we first determined the fraction of extravasated fluid relative to the total vascular blood flow into the tumor. To be precise we computed where V\LV T symbolizes the set of all vessels which intersect LV T with blood flow directed into the tumor. We obtain Q z(rel) l~4 +1:6 : 10 {4 over our base case states, which suggest that a fully coupled solution would not differ significantly from our solution. For completeness, the absolute values are presented in Table 3.
Moreover, we estimate the length scale over which IFP coupling would cause a decrease of the blood pressure along a single isolated vessel. For this purpose let q(x) be the flow rate as volume per time through a blood vessel with the axial coordinate x. It is determined by q~{Ap 0 , i.e. Poiseuille's law, where 0 denotes the derivative with respect to x, A a conductivity constant and p~p(x) the blood pressure. The fluid loss through the gaps of the vascular walls is the derivative of q, i.e. q 0~B (p i {p), whereby we have incorporated the osmosis contribution s(p v {p i ) into an effective blood pressure p. B is defined as B~2prL v l , where L v l is the wall permeability constant. p i denotes the interstitial pressure.
If we now assume that p i :const, we can easily derive a characteristic length scale l over which p approaches p i . By combining the above equations, we solve for p(x) and obtain p(x)~exp (x=l)C 1 z exp ({x=l)C 2 zC 3 xzC 4 , where C 1 ,C 2 ,C 3 and C 4 are constants which must be determined by boundary conditions. For l one obtains l~ffi ffiffiffiffiffiffiffiffi ffi A=B p i.e. the root of the ratio of flow to wall conductivity. Note that for l?? one can recover the standard Poiseuille's law. Based on our model parameters, the order of magnitude of l is estimated to lie within 10 3 to 10 6 mm, depending on the vessel. Actual measured values are shown in Figure 6D. Within the tumor l is actually longer than the system size. For capillaries l may be around 1 mm, which is still much longer than the length of typical capillaries. The value of l for major normal vessels is similar to tumor vessels. Other relevant variables namely B, q 0 and q are shown in Figure 6A, B, and C, respectively. In spite of the simplifications used we think this justifies the uncoupled evaluation of the interstitial fluid flow.
In the following we further comment on blood flow and extravasation. In silico data for B~2prL v l in Figure 6A shows a 20 fold increase from normal tissue to tumor tissue. Remarkably we can distinguish between an increase of r and L v l . The clustering with short ramp-up is associated with initially thin vessels which dilate to the maximal radius r (max) . Their permeability L v l increases simultaneously and eventually hits is upper bound l l,T . A plateau of maximal B forms where r and L v l reached their bounds which is here about 250mm behind the tumor boundary. The lower ramp corresponds to thicker vessels, too thick to dilate (because rwr (max) ), but for which L v l increases. The flow rate q, shown in Figure 6C, displays a clear distinction between tumor and normal tissue as do most other variables. In the tumor center it is more uniform and orders of magnitude higher than in normal tissue. q varies with the radius like r 4 , thus dilated vessels provide very well conducting pathways acting as arteriovenous shunts. Most of the data points in normal tissue stem from capillaries which have respectively slow flow rates. Going up the vascular hierarchy, we find increasing flow rates beyond those of tumor vessels and of course less vessels. In the neovascular plexus close to the tumor boundary we observe lower-than-normal flow rates, which is plausible since the blood volume is distributed over more vessels, which implies slower flow velocities in order to to satisfy mass conservation. See also our results and discussion in [18].
The flow rate through the vessel wall q 0 , shown in Figure 6B, correlates well with the wall permeability B weighted by the vessel circumference. Its magnitude within the tumor is about an order of magnitude larger than in normal tissue. This contradicts the common hypothesis that increased interstitial pressure hinders fluid extravasation. However, the permeability increase dominates the decreased IF pressure difference, which is only halved within the tumor compared to normal tissue.

Drug Transport
Here we analyze the spatiotemporal evolution of the concentration distribution s of some substance over the time frame of 96  hours. s was computed numerically according to (20) based on the data from previous simulations of tumor growth and interstitial fluid flow. Thereby the tumor is considered static and the IF flow is in a stationary state. Normally a tumor would grow further during this time frame. However since we do not model pharmacodynamics (i.e. cell killing) it seems reasonable to make this simplification. Parameters for our base case are derived from data for Doxorubicin, a commonly used chemotherapeutic drug. assume that initially the tissue is ''clean'' i.e. without drug. A bolus injection into the hosts body is modeled by the time varying blood plasma concentration s v (t). For an injection we take the exponential function s v (t)~exp({t=t sv ) with the time scale t sv1 h. Our results are presented with unit-less normalized concentrations for which s v (0)~1.
The distribution of Doxorubicin is usually observed in vivo with the aid of fluorescence imaging, e.g. in [43,44,48]. Typically one observes exponentially decaying concentration profiles around tumor vessels, at least during the first few hours. Eventually eventually drug is distributed relatively evenly. When the blood is cleared of drug, molecules diffuse back into the blood stream and so the concentration near vessels decreases again. The overall drug level decreases over the course of a few days until all drug is cleared from the tissue.
Our simulations results agree well with this observation. To illustrate this, Figure 7 shows the spatiotemporal evolution of the concentration s in a sequence of snapshots. The data are taken from one of the 15 systems we considered as the base case. See Supplement S16 for a video.
Upon closer inspection of these figures and comparison with Figure 4 A it becomes apparent that some vessels are not releasing drug. The explanation is that close to the tumor rim we have two classes of vessels with comparable radii but different permeabilities: (i) Mature vessels stemming from arterioles or venules which release little drug and (ii) dilated capillaries. Toward the tumor center these differences vanish due to wall degeneration.
In Figure 8A drug concentrations are plotted as average SsT h against h in the same way as Figure 3, i.e. as average over shells of constant h. Initially these drug profiles show a strong similarity with the profiles of the vessel volume fraction which is also plotted. Both have a peak at the tumor rim and decay into the tumor to significantly lower levels than in normal tissue. Of course, we expected such a correlation, because the ''bulk'' transvascular flux is proportional to the vessel surface area. This proportionality also implies a faster drug uptake by diffusion once blood is cleared from drug. Therefore one might naively expect that the tumor rim is cleared fastest of drug, afterwards normal tissue, and finally the tumor center. The actual result at t = 96 h is different, namely that the profile is relatively flat and monotonously increasing into normal tissue. Following the discussion below, it is clear that convection fulfills a significant contribution to transport, moving drug outwards and flattening the profile. This is studied in more detail considering a diffusion dominated system in the the case (iii) (see below).
The ratio of convection to diffusion is quantifiable by the Peclet number which is defined as Pe L~L V =D, where L is a characteristic length, V the velocity and D the diffusion constant. Pe L w1 means that transport is convection dominated and Pe L v1 means diffusion domination. It depends on L, so we analyze our system by determining L while requiring that Pe L~1 whereby we denote respective L as L dc~D =V . The mean over the tumor SL dc T T is 150+4mm. This practically the same length as the diffusion range of oxygen, i.e. the distance from drug sources (vessels) up to which viable tumor cells are present. Hence diffusion and convection are predicted to be equally relevant. In contrast L dc is two orders of magnitude larger within normal tissue which means that transport to spaces in-between capillaries is strongly diffusion dominated (see Supplement S13 Figure 1A for a spatial L dc map).
In Figure 8B the mean drug concentration SsT r is shown against the distance from tumor vessels r. The profiles were generated exactly like the plots against h by defining r as the (signed) distance from the region where the vessel volume fraction w v w0:01 which captures all vessels. Another interpretation of r is a penetration depth into the cuffs surrounding vessels and further into necrotic regions. It shows that over the first few hours Ss(r)T is very well fitted by an exponential decay function ! exp (r=r 0 (t)). Where the length scale r 0 is approximately 30mm, which is in good agreement with [43,44,48]. During the 3 to 24 hours period the exponential behavior vanishes. Thereby the level at the ''tail'' increases while the level in vascularized regions For the further analysis and quantification of drug delivery we have to introduce an appropriate metric that represents the timeindependent spatial distribution of drug doses. Quantities which commonly enter pharmocodynamical models are the maximal concentration over time and the area under curve (AUC) which is the time integral of the concentration. Even in the case of Doxorubicin, which has been known for a long time, models using either one have emerged, see e.g. the discussion in [49]. Mathematically those quantities correspond to the jj:jj ? and jj:jj 1 norms (over time), respectively. jjf jj 1 can be bounded by jjf jj ? ƒjjf jj 1 ƒtjjf jj ? where f is a proper function and t is the observation duration. This bound is of course not very strict so it seems justified to consider both. We use intracellular concentrations since we assume that drug has to enter the cell to e.g. bind to DNA in order to efficacious. Our denotations are ICMAX for jjs 2 jj ? and ICAUC for jjs 2 jj 1 .
As expected we find the highest values near tumor vessels from where the it decreases into the unvascularized (necrotic) regions further away, see Figure 9A and B. ICMAX and ICAUC are qualitatively similar and also display strong similarities with early concentration distributions at ca. t = 3 h. ICMAX is comparably sharp and in fact approximately s=w at t = 3 h. The extracellularly dissolved drug is negligible since its concentration is two orders of magnitude lower. Hence, s&ws 2 . Obviously ICMAX is very sensitive to the high initial drug levels, whereas ICAUC gives more weight to later smoother distributions and thus appears more blurry. Their penetration depth is about 200mm so the whole viable tumor mass obtains significant (meaning non-zero) contributions. The ratio of local maximum to minimum for the AUC is 23+ whereas it is 122+2 for ICMAX. In Figure 9C and D, ensemble averages of ICMAX and ICAUC are plotted as radial distributions together with the vessel volume fraction. The correlations are obvious and can be explained by the same reasoning as for the time-dependent concentration distributions. The displayed compartmentalization in decreasing exposure towards the center and peripheral peak is qualitatively retained for most of the parameter variations except for extreme cases with drastically increased interstitial transport rates such as case (v) (see below).

Variations
In this section we discuss a number of physiologically relevant variations of the base case scenario. We consider the following cases separately: For all cases we produced the corresponding data shown in Figures 3, 4, 5, 6, 7, 8, and 9. They are compiled in Supplements S2, S3, S4, S5, S6, S7, S8, S9, S10, S11, and S12. Videos visualizing the spatiotemporal evolution of the concentration distributions in cases (i), (ii) and (iii) are provided as Supplements S17, S18 and S19.
Below we discuss the physiological implications of these results. In Figure 10 we present a comparison of the mean S:T j and standard deviation (STD) std j (:) of the ICAUC and ICMAX distributions. j stands for a region over which mean and STD are taken. V denotes the viable tumor where w t w0:5, TB the boundary where {250mmƒhƒ0 and TI the interior where hv{250mm. This serves as an assessment of drug delivery efficiency. Higher mean, and smaller STD (i.e. less spatial fluctuations) means better delivery. Histograms representing actual probability density functions for ICAUC and ICMAX are shown in comparison in Supplement S13.
(i) Heavier drug particles. In this case we consider a drug with a molecular weight of 10 5 g=mol. This corresponds to the application of viruses or nano particles as delivery system, which we realize by the adjustment of diffusion related parameters, namely l s,N ,l s,T ,D s and k ij . For a particle performing Brownian motion, the diffusion constant scales with 1= ffiffiffiffi m p of the particle mass m. Hence we scale those parameters by ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 540=10 5 p , where 540g=mol is the molecular weight of Doxorubicin from the base case. This leads to a strongly convection dominated transport where L dc is of the order of 1mm. We see that for heavy particles the IF flow clearly dominates the drug distribution as seen in Figure 11 (see also Supplement S2, and the video provided in Supplement S17). As a result the tumor contains many small, isolated regions where a small amount of drug is delivered. This is reflected by the PDFs (p p(ICMAX), p p(ICAUC)) which are broader and lower values are more common. Interestingly, the p auc distribution completely changes from bell shaped to almost box shaped, see Supplement S13 Figure 2. The mean values SICMAXT V and SICAUCT V decrease by ca. 50% which can be expected due to the decreased transvascular diffusivity.
(ii) Prolonged infusions. Chemotherapy often follows a complicated schedule with several prolonged infusions, which is important to avoid toxicity to normal tissue. The basic idea is that drug given in low concentrations accumulates in the tumor, whereas the concentration in normal tissue remains tolerable. As a simple case we consider the administration of a prolonged infusion. In (ii)a during a 24 h period and in (ii)b continuously during 96 h. The ICMAX and ICAUC distributions, are surprisingly similar to the base case (see Supplements S3, S4 and the video Supplement S18). The difference lies in the scale. The average concentration SsT T increases approximately linearly in time for as long as the infusion is active.
That means even by the end of our simulation at t = 96 h the tumor is far from saturated. We expect this to occur when influx equals removal of drug. Assuming that s 2 =s 1~k12 =k 21 in a quasi steady state the time scale for a purely diffusing particle to move 1 mm can be estimated via the diffusion law SxT 2~D t as 1700 h, which is a reasonable estimate for the saturation time scale. Note that our results likely overestimate concentrations and underestimate the speed of the transport due to the lack of saturable cellular components.
(iii) Neglected Convection. This case serves to gain insight in the role of convective transport of drug through the interstitium.
For this purpose the convective term in (20) was neglect. Figure 12 shows that initially (t~1=2 and 3h) this case is indistinguishable from the base case. But later, SsT h profiles remain peaked within the tumor, leading to significantly increased peripheral drug concentrations. In fact the average concentration in the center remains nearly constant up to a ca. 250mm wide peripheral shell (see also Supplement S5, and the video Supplement S19). Thus, convection has the effect of ''flattening'' the profile apparently by driving additional drug drug into the neovascularized rim where it is reabsorbed once the blood stream is cleared of drug.
In the following cases (vi), (v), (iv) and (vii) we analyze the effect of varying permeabilities. In [3] it was shown that the IFP profile depends only on the ratio of vascular to interstitial hydraulic conductivity. Here we have spatially varying coefficients but the same scaling law is expected. Nonetheless we vary both parameters since their effect on the velocity and thus drug transport is different. Beside the IFP it seems appropriate to consider the actual flow through the tumor. We quantify this by the mean values S:T T of the following components of the source term in (12) : SQ l T T measures the amount of fluid that must leave the tumor through its boundary in interstitial space. SQ z l T T is the amount of extravasated fluid since vessels are the only sources. SQ { l T T is the amount of fluid taken up by vessels or lymphatics. Note that SQ { l T T is typically an order of magnitude less than SQ z l T T , which is clear since there should be very little back flow into vessels. Their response to parameter variations is shown in Figure 13.
(iv) Vascular permeability. We vary l l,T between 1/1000 and 10 times the b.c. value for leaky tumor vessels. Note that this does not simply scale all permeabilities (i.e. L v l ) equally, rather l l,T is the cutoff for L v l which increases up to this value for w?0 (see (15)). Also note the conductivity of normal capillaries is l l,N~1 =100l l,T (b:c:).   Increasing l l,T has little effect on the IFP and IFV profile since the IFP already approaches rapidly the blood pressure, see Figure 14 A. The plots of the source terms (Figure 13 A show that for increasing permeabilities most of the little additional flow is reabsorbed into vessels (SQ { l T T ). The flow that crosses the boundary is insignificantly altered. For lower permeabilities the uptake SQ { l T T decreases much more rapidly than the influx SQ z lv T T . The latter varies weakly within one order of magnitude over the whole l l,T range, implying that the hydraulic resistance of other components, i.e. interstitium and lymphatics limits the flow. Interestingly, setting l l,T to the level of normal vessels is not enough to lower the IFP to a normal level. It is only reduced by about 50%. Further reduction to 1=1000l l,T (b:c:) yields near zero IFP (& 0.1 kPa) but an outward gradient persists.
Predictions for drug delivery were made for l l,T~1 0,1=10 and 1=100l l,T of the base case. See Figure 10 (iv)-A to C (see also Supplements S6, S7 and S8, respectively). In conjunction with l l,T we also varied the diffusive permeability l s,T by the same factor. The mean ICMAX and ICAUC level in V TB are invariant since this regions exhibits rather normal vessels since leakiness increases gradually toward the center. In the interior, SICMAXT TI and SICAUCT TI increase with the permeability as expected. Here too, for (iv)-A we find ICMAX levels in the tumor which are comparable to normal tissue.
(v) Hydraulic conductivity of interstitium. The base case considers a dense tissue which implies a low conductivity due to the dependence on the available interstitial volume. Sparser tissues were estimated to be orders of magnitude more permeable [50]. To analyze such a situation we scale the tissue permeability coefficient K l,a up to 100 times simultaneously for a~N and T. For brevity, a is omitted in the following. Increasing K l alone produces unrealistic results where the IFP within normal tissue rises well above 0. Hence the lymphatic permeability L (L) l , was also increased by the same factor.
With increasing K l we observe a decreasing IFP which drops to ca. 1=7 of the base case value for a 100 times increase of K l , as shown in Figure 14B. The IFV decreases sub-proportionally, up to a factor of 20. Interestingly, the region where most fluid is extravasated shifts from close to the boundary to the tumor center. Concomitantly one finds a shift of the peak in the v(h) toward the tumor center and less back flow through vessels as compared to the base case as indicated by lowered SQ { l T T . For the total flow SQ l T T we can identify two regimes: Up to a 10 times increase we see an approximately linear variation, and for larger K l a logarithmic behavior.
Drug delivery is analyzed for an increase of K l by a factor of 10, where we additionally upscale the diffusion constant D s by the same amount (see Supplement S9 for additional figures). This is rather ad-hoc but based on the assumption that the free volume fraction and the amount of ECM components varies, on which the effective diffusion constant and IF permeability linearly depend on as a first order approximation. As a result we obtain a significantly more homogeneous ICMAX and ICAUC distributions. See the comparison Figure 10B,C-(v) vs. (b.c.). The effect is most drastic for ICAUC but also the ICMAX fluctuations over V TB are significantly reduced. Remarkably, the tumor interior now shows higher drug concentrations ( Figure 10A-(v)) than the exterior. The fact that the drug delivery to the interior is comparable to normal tissue is surprising since permeabilities were adjusted also for normal tissue not just the tumor.
(vi) Amount of normal lymphatics. The lymphatic system is not well documented due to the lack of specific markers for its  channel walls. In our parameter determination we assumed that it has a similar capacity as the capillary bed. Moreover tumors can induce lymph angiogenesis similar to normal angiogenesis [9] thereby increasing the amount of lymphatics nearby. Hence we consider variations of the source term coefficient L (L) l S (L) =V from 100 to 1/100 times the base case (b.c.) value.
As discussed above, the IFP has a ''penetration depth'' across the tumor boundary which is proportional to ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi K l,N l 2 =L (L) l S (L) q and indeed we see a variation of this magnitude in the radial profile. Moreover with increasing L (L) l S (L) =V the IFP drops asymptotically to the lymphatic pressure p L , which can be explained by the analogy of an electrical resistor. A L (L) l S (L) =V decrease on the other hand drives the IFP in unrealistically high. The central tumor IFP is thereby relatively invariant to these parameter changes, see Figure 14C. The IF velocity varies according to the gradient across the boundary and can be increased by up to 50% in the extreme case. The global flow SQ l T T increases by about the same magnitude. The average back flow SQ lv {T T decreases insignificantly within the error bars. On the other hand for low L (L) l S (L) =V more back flow is observed, which is reasonable since the uptake capacity of lymphatics is decreased.
For drug transport we only consider the case with a factor of 10. Qualitatively the ICMAX and ICAUC distributions appear similar to the base case (figures are provided in Supplement S10). The increased flow apparently leads to higher drug levels all over the tumor, but the effect is strongest in the TB region. Spatial fluctuations are also insignificantly changed. Thus we can conclude that the delivery slightly improved. Given that the fluid transport rate SQ l T increased by 50% this is surprising since case (i) demonstrated that a larger convection-diffusion ratio can have a negative effect.
(vii) Tumor lymphatics. In this hypothetical case we assume that functional lymphatics exist within the tumor. We model this by a non-zero lymphatic sink coefficient L (L) l S (L) =V as a fraction of the coefficient in normal tissue. This differs from taking a higher tissue conductivity by the assumption of the underlying channel organization. The latter case implies a grid like structure whereas the former implies a hierarchical structure where exchange takes place over the thinnest channels analogously to blood vessel network. Those lymphatic capillaries can be expected to exhibit approximately equal blood pressure over the whole tissue. Hence the use of homogeneously distributed sinks with lymphatic pressure p L .
Resulting IFPs and IFV plots are shown in Figure 14D. As can be seen an increasing presence of lymphatics ''pulls'' down the IFP to p L . The outward velocity component v jj decreases proportionally as can be expected based on the IFP gradient. But the actual velocity magnitude jvj increases with the amount of lymphatics. This can be expected as well since the velocity field is increasingly superimposed by flow away from vessels to nearby locations where fluid is absorbed. This is also reflected by increasing transvascular flow SQ l zT T , decreased flow across the boundary, SQ l T T and the back flow SQ { lv T T as shown in Figure 13C. We computed the drug transport for systems with (a) 0.1 and (b) 1 times the normal lymphatics (see Supplements S11 and S12, respectively). What we would expect is that the changed flow patterns (directed away from vessels) facilitate the delivery to more distant regions. As Figure 10A and B-(vii) show this is indeed the case. Even in (a) a significant improvement is achieved. The ICMAX and ICAUC levels in the rather unrealistic case (b) are increased dramatically. This must also be attributed to leaky tumor vessels. Thereby spatial fluctuations are either unchanged (std : (ICMAX)) or significantly reduced std : (ICAUC).

Discussion
Based upon an extension of our model for remodeling of tumor vasculature introduced in [15,16,18,19] we studied interstitial fluid flow (IFF) and drug transport.
Since the spatial details of IFF depends crucially the spatial arrangement of blood vessels and their blood pressure we considered, for the first time, a model for IFF and drug transport that involves a realistic arteriovenous initial vascular network evolving dynamically with the tumor. As a result tumor vessels are connected to a wide variety of original host vessels -from capillaries to ca. 50 mm arterioles and veins. Blood flow is determined for the whole system including the whole remaining host vasculature which comprises not only a single parent vessels but many vessels covering the entire simulation domain. It should be emphasized that an arteriovenous initial vasculature produces blood flow patterns and blood pressure fields that differs substantially from grid-like arrangements [15,18]: the latter usually have a fixed blood pressure gradient direction imposed upon them which can also (unrealistically) impose a preferred direction on tumor vessel growth and skew the IFP distribution. For a tumor blood vessel network emerging from an arteriovenous initial vaculature we obtain much higher blood ow rates through individual vessels due to arteriovenous shunts and circumferential growth, and more irregular (spatial) blood pressure distributions. Due to transvascular coupling these circumstances are also relevant for the IFP, IFF and drug transport.
A first remarkable result is that in spite of an expected IFP plateau within the tumor IFF does not cease and still allows for substantial convective transport, which is opposite to the currently prevailing view [3][4][5][6][7][8][9] that states that increased IFP poses a barrier to successful drug delivery within tumors. The physical explanation is simply that it is misleading to consider the pressure drop along the vessel wall alone as the driving force for IFF -in principle the complete network of hydraulic resistors has to be taken into account to obtain reliable predictions. Qualitatively it is sufficient to bear in mind that for IFF the tumor is series of resistors with a fixed potential outside the tumor (the lymphatics of healthy tissue) -decreasing one resistor, for instance by increasing the permeability of the vessel walls, generally increases the flow, in spite of the lower potential drop along the decreased resistor.
Another interesting result is that heavy macro-molecules are still distributed more or less evenly into viable areas in the tumor perimeter in spite of a pronounced outward IFP gradient there, which one could naively expect to remove drug before an efficacious dose is achieved. The reason is IFF between tumor internal vessels, that transports macro-molecules convectively from leaky high pressure vessels through the tumor tissue into neighboring low pressure vessels. In the following we discuss our results quantitatively in detail.
It is already established that leaky tumor vessels and lack of tumor lymphatics lead to a drastically increased hydraulic pressure of the interstitial fluid (IFP) in the tumor and that the resulting gradient drives fluid out of the tumor with a velocity of 0.1 to 1mm=s [4,51] and our results agree with these experimental observations very well. Our model predictions of the central IFP is ca. 6 kPa, which is at the upper limit of experimentally observed values. This value is also higher than in normal capillaries due to the coupling with higher level arterioles in which the blood pressure is naturally higher. The vessel walls of arteries and veins within the tumor become leaky [21] and thus increase their conductivity, which means tighter coupling between IFP and blood pressure which as a result causes the IFP to approach the level of nearby arteries or veins. Indeed, the predicted mean pressure difference between blood and interstitium vanishes towards the tumor center. Local fluctuations produce flow into, out of and in-between vessels which is not necessarily directed outwards. The outward component of the velocity vector dominates in the tumor periphery. Around the tumor perimeter we find a thin layer where where the fluid is absorbed into lymphatics. Absorption into blood vessels frequently occurs deep inside the tumor. It should be noted that the fluctuations among samples are very large, about 100%. Consequently, depending on their location in the microenvironment, seemingly by chance some tumors could be more likely to metastasize through the blood stream than others of the same kind. Although the transvascular flow of tumor vessels is elevated by an order of magnitude compared to normal capillaries, we estimated that only a small fraction of the order of 0.01 percent of the blood flow which enters the tumor is lost into the interstitium. The biophysical factors upon this ratio depends are the vascular morphology and the permeability of vessel walls, interstitium, and lymphatics. Here, dilated vessels and arteriovenous shunts lead generally to elevated flow rates (q) in tumor vessels. Note that q depends strongly on the vessel radius r, i.e. q!r 4 .
In addition to the base case we considered scenarios in which the conductivities of vessel walls, lymphatics and the interstitium are varied individually. We observe universally that the influx through vessels, and therewith the IFF through the tumor, increases with the conductivity. The sensitivity to these changes is rather low, even for variations of several orders of magnitude. Lower IFP does thereby not correlate with increased flow. For example the reduction of IFP can be achieved by increasing the tissue conductivity or decreasing the vessel wall conductivity. In the first case the IFF is increased while in the second it is decreased. This relationship between IFF and conductivity can be easily understood on the basis of analogy of our flow equations with an electrical network of ohmic resistors: As an extreme idealization let us consider a linear chain of three resistors, representing the walls of the tumor vasculature, the interstitium and lymphatic walls, whereby the potentials are fixed at the ends. Then of course the current is determined by the total resistance and the voltage drop over one resistor grows with its resistance. An increase of the first resistor (decrease of vessel wall conductance) as well as decrease of the second resistor (increase of tissue conductivity) both imply a lower voltage drop at the second resistor (lower IFP). But the current (the IFF) decreases in the 1st case and increases in the 2nd.
We quantify the local exposure to drug with the help of both the time integrated intracellular concentration (denoted as ICAUC for intracellular area under curve) as well as the maximal concentration (denoted as ICMAX for intracellular maximal concentration), which we consider separately. Our model predicts that on a large scale, drug delivery is compartmentalized similar to the vasculature. Let us subdivide the system into concentric shells with increasing distance from the invasive edge and consider averaged quantities over such shells. Close to the invasive edge we typically find an exposure peak, the location of which is not exactly aligned with the MVD peak due to outward convection. Towards the tumor center the decreasing vessel density leads to a sharp drop down to a plateau at a certain lower bound which is for our base case ca. 50% of the level at the tumor periphery. The height of this plateau depends on various factors including extravasation rate, transport rate through the interstitium, but predominantly vessel density.
The tumor center is threaded by many isolated vessels, and our model predicts exponentially decaying concentration gradients around them in agreement with experimental results e.g. [43,48]. Over longer time scales of many hours to days, the initial distributions smear out akin to the behavior of freely diffusing particles, and also decrease globally due to reabsorption. The magnitude of these time scales depend on the transport rates, cellular uptake and binding dynamics.
Although we did not study it in detail, we want to stress that cellular uptake and retention dynamics also govern the total amount of extravasated drug to a certain extend because the faster the uptake rate, the lower the interstitial concentration, the larger transvascular gradients which drive diffusive fluxes and counteract re-absorption by convection.
It has been suspected for over a decade that IFF could wash out drug from the tumor. Our model predicts this effect but as drastic in magnitude as expected. Simulation runs in which convective transport was neglected show that a ca. 2.5 times higher drug concentration is retained in the periphery whereas the radial profile of the time-independent exposure measures (ICAUC, ICMAX) merely exhibit a smoother decay from periphery to center. In the opposite extreme case of non-diffusing particles under normal convection the the interstitial flow causes sufficient flux to achieve significant drug delivery which is on a coarse scale comparable to the base case. Here we do however observe islands in the viable tumor region where no drug is delivered, implying that such tumor fragments could remain viable after treatment. In cases with diffusion, at least small amounts of drug are delivered.
Different tissue types as well as the effect of therapies that improve drug delivery can be described within the model by changing various permeability constants. We varied diffusion constants simultaneously with conductive constants by the same factors as first approximation to hypothetical changes e.g. of the intercellular channel geometry. As a result, global exposure levels and the amount of extravasated fluid correlate well with these permeabilities as well as with each other. Unfortunately relative local fluctuations, i.e. the STD over space of ICAUC and ICMAX are not correlated with these variations, making them unsuited as basis for achieving are uniform efficacy or dose.
Within the biologically relevant parameter space of which we only considered a small part, analysis of other cases will certainly lead to additional insight. For example in [14] tumor interstitium and capillary permeabilities (K l,T and l l,T in terms of our model) were varied simultaneously, leading to a flattened tumor IFP profile for increasing permeability. From such observations in an experimental setting one could draw conclusions about the nature of the tumor tissue.
With this in mind a tumor therapy that comprises a treatment that solely reduces the vessel leakiness appears not to be effective. A currently frequently discussed alternative is to improve the efficacy of drugs by turning the ill-formed tumor vasculature normal again [52] by pharmacological means during or before a conventional therapy. A step in this direction is obviously to reduce the leakiness of tumor vessels. However, our simulations predict an actual reduction of the global ICAUC and ICMAX levels while their relative fluctuations may even increase. On the other hand increased permeability leads to significantly improved delivery in the tumor center. Unfortunately this case is unlikely to be a therapeutic goal since it could increase the direct exposure of the blood stream to tumor cells increasing the chance of tumor cells entering it. A good treatment strategy would be to prioritize the maintenance, or fabrication of a dense tumor vasculature, rather than exclusively tightening the leaks in the already sparse vasculature. If a tumor is detected early enough the vasculature could thus be kept intact and deliver drugs effectively, simply due to the amount of functional capillaries, even though their walls are expected to conduct drug worse than leaky walls.
Vascular targeting therapies that take the opposite direction namely that aim to destroy the remaining tumor vasculature completely exists hand have proven to be promising, see e.g. [53]. A common problem in therapy is that tumor cell in vivo are more resistant to treatment than cell cultures. Many factors are involved but not all are purely of genetic nature. For example the lack of oxygen plays an important role for the development of radiation resistance. Or the drug delivery through the sparse tumor vasculature is insufficient, which is supported by our model prediction that the drug concentrations in the tumor are much lower than in normal tissue. If cells in the interior were killed indirectly by vascular targeting, a viable shell around the boundary would remain which could be effectively treated conventionally due to the better vascularization. It is plausible that the removal of the inner tumor vessels would lead to a reduction of the IFP and thereby lower IFF across the boundary. As our results show, the absence of this peripheral flow can improve the drug exposure of the boundary significantly potentially leading to better treatment in combination with vascular targeting.
Interestingly, the ICMAX and ICAUC levels increase with the interstitial permeability, lymphatic permeability, or when providing the tumor with a fraction of surviving lymphatics. In addition to that and perhaps more importantly, drug concentrations become more homogeneous as well. This effect is much more prominent for ICAUC than for ICMAX, where for the former the spatial STD decreases from over 0.4 to ca. 0.2. With a 10 times increased interstitial conductivity the ICAUC distribution even becomes completely smooth over the whole tumor at levels higher than in normal tissue. Unfortunately, in those cases the IFF is also increased which would likely increase the shedding of TCs into lymphatics, not to mention the danger of having TCs in direct contact with tumor lymphatics. Hence, despite significantly improved delivery the parameters above are certainly not a useful therapeutic target.
The final ICMAX and ICAUC distributions for continuous infusions are hardly distinguishable from the bolus case except for the scale. The reason is that the time scale to achieve a stationary state (estimated ca. 1700 h) is much longer than the infusion periods (24 h and 96 h) which we tried. For Doxorubicin cellular uptake is relatively fast (of the order of minutes) and release is slow (of the order of hours), hence the local concentrations are in a steady state where most of the drug is arrested in cells, lowering the effective transport rates. It can be expected that drugs with low uptake rates perform better in this respect.

Supporting Information
Supplement S1 Further details on the model. Supplement S14 Video of a growing tumor in silico. A quarter of the system is cut out to open the view into the interior. The presentation is analogous to Figure 2, i.e. the yellow mass depicts the viable tumor. Void spaces within the tumor are necrotic regions. The blood vessel network is color coded by blood pressure. Red is high (arteries), and blue is low (veins). (AVI) Supplement S15 Video of a growing tumor in silico. Here we visualized a ca. 200 mm thick slice through the center of the system. The presentation is analogous to Figure 2, i.e. the yellow mass depicts the viable tumor. Void spaces within the tumor are necrotic regions. The blood vessel network is color coded by blood pressure. Red is high (arteries), and blue is low (veins). (AVI) Supplement S16 Video of the drug concentration distribution for the base case.

(AVI)
Supplement S17 Video of the drug concentration distribution for case i.

(AVI)
Supplement S18 Video of the drug concentration distribution for case ii-b.

(AVI)
Supplement S19 Video of the drug concentration distribution for case iii. (AVI)