Mixture theory modeling for characterizing solute transport in breast tumor tissues

Background Tumor numerical models have been used to quantify solute transport with a single capillary embedded in an infinite tumor expanse, but measurements from different mammalian tumors suggest that a tissue containing a single capillary with an infinite intercapillary distance assumption is not physiological. The present study aims to investigate the limits of the intercapillary distance within which nanoparticle transport resembles solute extravasation in a breast tumor model as a function of the solute size, the intercapillary separation, and the flow direction in microvessels. Methods Solute transport is modeled in a breast tumor for different vascular configurations using mixture theory. A comparison of a single capillary configuration (SBC) with two parallel cylindrical blood vessels (2 BC) and a lymph vessel parallel to a blood vessel (BC_LC) embedded in the tissue cylinder is performed for five solute molecular weights between 0.1 kDa and 70 kDa. The effects of counter flow (CN) versus co-current flow (CO) on the solute accumulation were also investigated and the scaling of solute accumulation-decay time and concentration was explored. Results We found that the presence of a second capillary reduces the extravascular concentration compared to a single capillary and this reduction is enhanced by the presence of a lymph vessel. Varying the intercapillary distance with respect to vessel diameter shows a deviation of 10–30% concentration for 2 BC and 45–60% concentration for BC_LC configuration compared to the reference SBC configuration. Finally, we introduce a non-dimensional time scale that captures the concentration as a function of the transport and geometric parameters. We find that the peak solute concentration in the tissue space occurs at a non-dimensional time, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {T}_{peak}^{\ast } $$\end{document}Tpeak∗ = 0.027 ± 0.018, irrespective of the solute size, tissue architecture, and microvessel flow direction. Conclusions This work suggests that the knowledge of such a unique non-dimensional time would allow estimation of the time window at which solute concentration in tissue peaks. Hence this can aid in the design of future therapeutic efficacy studies as an example for triggering drug release or laser excitation in the case of photothermal therapies.


Background
The total cost of cancer care in the United States is projected to increase by 39% from 2010 to 2020 [1]. Primary areas of cancer research involve improving the efficacy of chemotherapeutic agents at the tumor sites and minimizing their toxic side effects in the non-target sites [2][3][4]. Conventional chemotherapeutic agents [5] are non-specifically distributed in the body which limits the effectiveness of the drug dose and increases toxicity in normal cells. Drug carriers with hydrodynamic diameter 3-200 nm accumulate preferentially in tumors owing to the enhanced permeability and retention (EPR) effect [6] exploiting the wider pores in tumor vessels and the impaired lymphatic drainage in diseased tissues. The transport mechanism of these nanoparticles in tumors is a function of the hemodynamics, nanoparticle transport parameters (solute permeability, solute diffusivity, reflection coefficient) as well as the extravascular matrix properties (porosity, hydraulic conductivity). Before the binding/uptake by the cancer cell these particles overcome three major transport barriers: transport through microvasculature, translocation across the endothelial wall, and diffusion within the extracellular tissue matrix. Using a numerical model to investigate the nanoparticle transport mechanics could enable determination of the exact time interval between nanoparticle introduction and drug release to achieve desired therapeutic efficacy based on patient specific tumor measurements.
The majority of the existing multiscale models use the Darcy's law, Starling's law, and Poiseuille's law to analyze extravascular, trans-capillary, and intravascular transport respectively [7][8][9][10][11][12]. Poiseuille's law cannot account for variations in capillary diameter and the inhomogeneous nature of blood. Deviations from Starling's law are expected when osmotic terms would include other endogenous solutes in addition to proteins. Darcy's law does not include the dependence of interstitial flow on local fluid chemical potential. Schuff et al. [13,14] used mixture theory equations in an axisymmetric tissue geometry containing a concentric blood vessel and showed the dependence of extravascular fluid transport on chemical gradients in addition to hydrostatic pressure which was previously suggested [15,16] and observed [17] but not commonly accounted for in previous transport models.
Recently developed tumor numerical models [11,18,19] focus on complex capillary distribution in the tissue, multi-stage drug delivery systems and dynamically changing tumor geometry but for the nutrient transport they rely on the advection-diffusion equation only instead of considering each type of particle (solute, fluid, cells etc) as a mixture component. In the present work, we focus on nanoparticle transport only, neglecting, drug pharmacokinetics, dynamic tumor geometry variation and complex capillary network to delineate a characteristic metric applicable to all kinds of solute and fluid transport in tumors. The mixture theory model is implemented in dual-tissue geometries to predict nanoparticle distribution in cancerous breast tissues over a wide range of particle sizes (0.5-15 nm) and molecular weights (0.1-70 kDa).
We hypothesize that nanoparticle distribution in breast tumors is a function of solute size, intercapillary separation, and flow direction and there exists a characteristic nondimensional time, T Ã peak , for which solute concentration in the tissue space is maximum. We test this hypothesis by investigating the transport mechanisms of five solute types (0.1, 3, 10, 40 and 70 kDa) in tumor systems containing a single vessel (SBC) and compare with tumors possessing dual-vessel (blood capillaries only (2 BC), blood capillary and a lymph capillary (BC_LC)) tissue systems with varying intercapillary separation. The mixture theory equations are used for the first time and their predictive capability validated with measurements of dextran transport in an in vitro tumor platform containing multiple blood vessels.

Results
Numerical simulations with baseline values of breast tumor transport (Table 1) show a strong influence of vascular configurations on the solute accumulationdecay temporal history in the tissue space. A schematic of the vascularized breast tumor configuration along with the transport pathways is shown in Fig. 1. The blood vessel in Fig. 1 allows both intravasation and extravasation depicted by blue and red arrows respectively. The lymph vessel allows intravasation (blue arrow) only and drains the lymphatic fluid out of the tissue. All the equations used in the simulations are explained in the Appendix 1.
Test 1: effect of solute size in three configurations for fixed intercapillary separation The effect of solute size in the double blood capillary (2 BC) and the blood capillary -lymph vessel (BC_LC) configurations compared to the single capillary (SBC) is shown in Fig. 2. The tissue volume surrounding each capillary is equal.
The concentration-time history essentially shows an initial accumulation period until the solute reaches its maximum concentration in the tissue space, then it is followed by a concentration decay.
It is observed in the accumulation phase of the SBC configuration ( Fig. 2a inset)

Test 2: effect of flow direction in microvessels
The microvessel flows considered in test1 are in the same axial direction and are called co-current (CO) flows. They are compared with oppositely directed axial flows in the microvessels which are called counter current (CN) flows. Figure 3 compares the CN flow with the CO flow for 2 BC configuration and BC_LC configuration respectively. In both configurations there is no difference between the two flow types during the solute accumulation phase in the tissue. As the tissue concentration decays, the extravascular solute concentration is less in CN flow at later times compared to CO flow. The percentage reduction in concentration is more pronounced for larger (10 kDa: 0.5%; 40 kDa: 1.7%; 70 kDa: 2.5%) solutes in 2 BC configuration (Fig. 3a, b). A similar observation is made for heavier (10 kDa: 0.12%; 40 kDa: 0.4%; 70 kDa: 0.5%) solutes but the difference is less in BC_LC configuration (Fig. 3c, d) compared to the 2 BC configuration. The faster drainage in presence of additional capillaries delays the solute accumulation. Thus, it takes a longer time for the solute to reach its peak concentration in a tissue location and compared to SBC configuration the magnitude of the peak accumulation is lower.
Test 3: effect of intercapillary separation on transport of 3 kDa and 10 kDa solutes The intercapillary separation (L) was varied with respect to the vessel diameter (d) in the next set of tests for the 2 BC and BC_LC configurations. We calculated the surface averaged extravascular concentration at a radial distance 0.2 L from the blood vessel wall in SBC configuration (C SBC ) as well as in the 2 BC configuration (C 2BC ) where L is the intercapillary separation (Fig. 4a).
The dotted lines represent L/d values of 5 (yellow) and 21.5 (purple) which are the lower and upper limits of normalized intercapillary separation in breast tumors. L/d values above 21.5 are typically found in normal (non-diseased) tissues.

Double blood capillary (2 BC) embedded in tissue cylinder
The 2 BC configuration shows minimum deviation (3 kDa: 1.5%; 10 kDa: 2.5%) from the SBC configuration for L/d =  (Fig. 4b). This is because the spacing between two blood capillary walls is so small that they effectively function as a single capillary with twice the original capillary diameter. So, the solute accumulation almost the resembles that in a SBC configuration. The solute concentration deviation is substantial (9-33%) for L/d = 5-25 which is the range of interest as reported in Table 2. The isolated capillary assumption will not hold true for extravascular solute accumulation in this regime. The deviation (3 kDa: 4.7%; 10 kDa: 8.5%) reduces for L/d = 125 because the large spacing between capillary walls minimizes the cumulative effect of the two blood vessels on the peak tissue concentration.

Blood capillary and lymph capillary (BC_LC) embedded in tissue cylinder
The BC_LC configuration (Fig. 4b) shows minimum deviation (3 kDa: 7.4%; 10 kDa: 7.3%) from the SBC configuration for L/d = 125 due to the same reason as the 2 BC configuration. But with decreasing L/d the sink action of the lymph vessel become increasingly dominant resulting in 44-55% deviation of the maximum solute concentration in tissue volume from that in the corresponding SBC configuration in the L/d regime of 1 to 25.

Non-dimensional time vs peak non-dimensional concentration analysis
The results discussed in the previous cases have shown that the variation of solute size, microvessel arrangement, number of microvessels and intercapillary separation all contribute differently to the solute accumulation time vs solute decay time in the tissue space.
Hence the non-dimensional extravascular solute concentration (C*) and non-dimensional time (T*) were calculated for all test cases to account for variations of these four parameters. The non-dimensional profiles for fixed L/d = 1, 5, 10, 25, 125 were plotted in Fig. 5a. All concentration peaks lie within T* = 0.1 shown by the dotted black line. So, in Fig. 5b the scaled concentration time-histories from T* = 0 to T* = 0.1 were analyzed. The peaks were extracted and plotted in Fig. 5c. The solid line red curves that correspond to the largest intercapillary distance (L/d = 125) have T Ã peak values which are one order of magnitude less than the average T Ã peak . This is because for a large L, the second vessel does not contribute to the solute accumulation in the measurement location which is at a distance 0.2 L from the first blood vessel. The concentration gradients across each capillary wall dynamically change the extravascular flux across the wall and for a large L the solute may get trapped within a certain distance of the capillary. The non-dimensional equations do not account for these and hence the deviation of T Ã peak for L/d = 125. It was concluded that the average non-dimensional time at which the peak concentration occurs in all configurations for all solutes is T Ã peak =0.0 27± 0.018 (Fig. 5c).

Discussion
We present a numerical model for breast tumor that can predict passive transport of nanoparticles across a multilayer barrier when the tissue architecture and nanoparticle properties are specified. The solute size (molecular weight and hydrodynamic diameter) influences its extravascular concentration dynamically across time. Smallest solutes are cleared faster from the tissue but they are also susceptible to getting trapped in the recirculation zone set up by countercurrent blood flow in adjacent vessels [56,57]. High clearance rates measured for solutes < 10 kDa during in vivo measurements in tumors validate this observation. The therapeutic outcome of breast cancer drugs (Doxorubicin: 0.54 kDa; Cisplatin: 0.3 kDa) having similar molecular weight as the smallest solute investigated here can now be predicted for patient-specific tumor biopsies. Heavier nanoparticles (50-200 kDa) are preferred vehicles for the tumor location specific targeting and drug delivery [57,58]. According to the results, heavier solutes with hydrodynamic diameter (5-14 nm) take longer to attain maximum accumulation at a specific tissue location and are also removed slowly [59,60]. Presence of adjacent blood vessels with counter-current flow accelerate the clearance process owing to drainage from both ends of the tissue. Presence of a lymph vessel reduces the magnitude of their peak concentration considerably owing to high permeability cross the lymph capillary wall. In various in vitro/in vivo studies the lymph wall is shown to allow unidirectional flow only [61][62][63][64] that facilitates better drainage. Easy removal of all sizes of drugs through the lymphatic system contributes to role of lymphatics in worse response [65].
Blood vessels have been mostly modelled as nonporous structures permeable to solutes embedded in the flow [66,67]. To our knowledge, this is the first numerical model that account for endothelial porosity directly measured from dextran transport in an in vitro breast tumor model. Additionally, this work shows inclusion of multiple vessels in a tumor numerical model is necessary to accurately measure transport phenomena. The SBC assumption works only for tissues where capillaries are so close (L/d = 1) that they act as a single vessel, e.g. when nearby lymph vessels collapse [61,62,64,65] during metastasis of some cancers or if they are so far Human colorectal tumor center 177 17.7 [55] apart (L/d = 125) that the solute flux from one does not reach the other. The dextran transport investigated using a fabricated 3D microfluidic platform measured tissue porosity, endothelial porosity, nanoparticle permeability and nanoparticle diffusivity. Simulations driven by these parameters showed a close correspondence of numerical and experimental concentration-time histories. These parameters, when reported in literature, span over several orders of magnitude. The wide range can be attributed to the complex in vivo measurements whose intrusive nature would perturb the tissue microenvironment [17,68]. Thus, the ability to measure these parameters in vitro can be used to design future non-invasive transport investigation studies.
Condensing all the effects of tissue architecture, solute and fluid transport properties, we suggest there exists a unique time T Ã peak at which nanoparticle concentration in the tissue is maximum. Previous analytical solutions of a convective-dispersive solute transport equation with time-dependent inlet boundary condition [69,70] have shown the dependence of time constant on the input timescale, convection timescale, decay timescale and the diffusion timescales but did not account for varying intercapillary separation. Later Chapman et al. and other researchers [8,11,12,71] modeled transport in tumors characterized by intercapillary separation but did not account for the dynamic change of hydraulic permeability as a function of concentration gradients which in turn modulates the extravascular solute flux. The presented work, for the first time, analyzes the solute concentration in the tissue in the light of mixture theory equations for varying solute types, two parallel microvessels, differing flow directions in parallel vessels and tissue architecture and proposes a non-dimensional time at which solute concentration is maximum in the tissue.
Since this approach non-dimensionalizes the intercapillary separation (L) with the vessel diameter (d), T Ã peak can be predicted for tissues ranging over several scales and also for different disease stages (cancer vs normal). This prediction would aid in efficient endothelial targeting, triggering drug release and laser excitation for photothermal therapies [72][73][74]. The T Ã peak estimation can hugely impact the clinical landscape as it would customize treatment based on tumor specificity. Future studies with varying nanoparticle design, changing dosage, presence of magnetic targeting, receptor binding can all be implemented first to the mixture theory model whose predictions would increase the efficacy of the targeted drug delivery in patient specific tumors.
The major limitation associated with the study are as follows: The complex vascular network have been simplified. The tortuosity and diameter variation of the microvessels were neglected. The extracellular matrix was considered stationary and not allowed to deform. A zero flux boundary condition was prescribed at each microvessel outlet which deviates from the physiological condition where a constant solute flux is drained to other organs like the liver from the microvessel outlets.

Conclusions
The study described in this paper focuses on quantification of solute transport across parallel blood vessels and initial lymph vessels in the light of mixture theory. Transport of nanoparticles to the targeted tumor volume is defined by the transport through microchannels, diffusion across endothelium and transport within the porous matrix, all of which were accounted for in the presented work. The results show that the solute size strongly influences its own rate of removal and rate of accumulation in the tissue. The flow physics in the extravascular space facilitate tissue drainage of nanoparticles depending on the solute size, the intercapillary separation and the microvessel arrangement in the tissue. A unique non-dimensional time T Ã peak was reported for the first time for peak solute accumulation in absence of pharmacokinetics. This is the time at which peak concentration of a nanoparticle occurs at any tissue location. The knowledge of the nanoparticle introduction time, tissue mechanical properties and solute dependent properties will allow, in future, to design in vitro tissue models testing various nanoparticle designs and concomitantly, predict for patient specific tumors, the appropriate time of drug release that can substantially improve drug efficacy.

Mixture theory model parameters
The mixture theory model requires fifteen input parameters that account for the mechanical properties of the porous matrix, type of the injected solute, and vascular geometry. The governing equations and boundary conditions for the mixture theory model have been derived in [13]. A sensitivity analysis of the input parameters, calibration and subsequent validation of the model was carried out by Schuff et al. in [14]. The mixture theory equations model the transport [13,14] of the fluid and solute in three distinct regions of a representative vascularized tumor geometry: a) in the intravascular space which consists of the plasma layer concentric with an inner core of red blood cells, b) across the capillary wall which is thin and semi-permeable and c) the extravascular space that comprises of the interstitial fluids and solutes flowing through a fibrous matrix. A finite element software package COMSOL 4.2 (COMSOL, Burlington, MA) was used to run the simulations. The simulation run time for each solute transport in a specific tissue configuration was 1.5 days.

Experimental validation of mixture theory model
The accuracy of the computational model was confirmed with experimental measurements performed in a 3D vascularized in vitro tumor microenvironment. Essential model parameters such as tissue porosity, vessel porosity, solute permeability, and solute diffusivity were measured using the in vitro platform and implemented in the model ( Table 3). The concentration-time histories were obtained from the mixture theory equations using the minimum, maximum, and mean values of the tissue parameters measured from in vitro platform. The simulation results were also compared with experimental measurements from dextran transport in the same vascularized in vitro platform (Fig. 6). Details explaining the fabrication and measurement of porosity, solute permeability and solute diffusivity from the fabricated platform can be found in the Appendix 2. To our knowledge, this is the first in vitro model that measured different porosity values in the extravascular and intravascular spaces.

Experimental and numerical comparison of concentrationtime histories
Numerical simulations for the model validation studies were separately processed with identical tissue properties and boundary conditions as obtained from the experiment. The intensity-time histories were spatially averaged at a radial location of 600 μ m. These were normalized by the maximum intensity inside the vessel at that time instant. For each of 3 kDa and 70 kDa solutes, transport was studied in N = 3 tissue samples with identical fabrication parameters. The normalized intensity profile of dextran particles from these experiments corresponds to the normalized concentration from the numerical simulations (Fig. 6). The error bars correspond to the experimental variability observed across 3 samples at each time instant. For 3 kDa, the 1e-12 --smaller solute, the experimental data till 1 h matches with the simulation curve from maximum values of input parameters. The deviation of experimental results from mean simulated values decreases with increasing time from 1 h to 2 h. The higher deviation of the 3 kDa numerical curves from the experimental curve can be attributed to the higher susceptibility of the small solutes (< 10 nm) to get trapped in local recirculation zones in a non-homogeneous tissue [56,64]. The experimental data for the 70 kDa solute almost coincide with the mean simulation curve and is closely enveloped by the maximum and minimum simulation curves.

Vascular tissue configurations and test matrix
In the current study three geometrical configurations of the vessels in the tissue are considered: the single blood capillary (SBC) configuration, the double blood capillary (2 BC) configuration, the blood capillary and the lymph vessel (BC_LC) configuration. The simulations were split into 3 tests listed in Table 4. The flow direction is same (CO) in parallel microvessels for tests 1 and 3. The intercapillary distance (L) in the dual-microchannel configurations is 100 μ m in tests 1 and 2 [29,48,52].

Test 1: effect of solute size on transport
Five solutes (0.1 kDa, 3 kDa, 10 kDa, 40 kDa, 70 kDa) are injected at the inlet of the blood capillaries in each configuration to delineate the effect of solute size on transport. The geometry, hydraulic parameters and    Table 1. The transport properties of 0.1 kDa solute are estimated from the calibration model developed in [14].

Test 2: effect of flow direction in capillaries on transport
The second test investigates the influence of co-current flow and counter-current flow in both 2 BC and BC_LC tissue configurations and compares the solute transport mechanisms with the SBC configuration.

Test 3: effect of intercapillary separation on transport
In the third test, the capillary separation (L) is varied with respect to blood capillary diameter (d), (L/d = 1, 5, 10, 25, 125), to test the solute accumulation in the tissue space as compared to the SBC configuration. Due to long computational time for each solute, this test is run for two solutes, 3 kDa (representative of a therapeutic drug) and 10 kDa (representative of the size of drug carrying nanoparticle). The intercapillary separation, L, is defined as the shortest distance measured between all non-adjacent capillary pairs in a loop [46]. The tissue diameter (D) is twice the separation value (D = 2 L) to maintain the same volume of tissue around each microvessel with respect to the SBC configuration. The surface averaged extravascular concentration is calculated at a radial distance 0.2 L from the blood vessel wall in SBC configuration (C SBC ) as well as in the 2 BC configuration (C 2BC ). The percentage concentration deviation shown in Fig. 4a is compared across different values of intercapillary separation. The intercapillary distances for different tissue types in humans and small animals are recorded in Table 2. The blood capillary diameter, d, is 10 μm in accordance with the values reported in literature across humans, mice and rats [46,48,49]. The L/d for tissue types shown in Table 2 lies between 5 and 21.5. So, the analysis was performed for L/d = 1, 5, 10, 25 and 125 for two solute molecular weights of 3.0 kDa and 10 kDa. For all the tests, the average extravascular concentration in the tissue volume surrounding the blood capillary is measured for each configuration. They are normalized by the maximum intravascular concentration in the blood capillary volume and the concentration-time history for each test is analyzed in the Results section.

Non-Dimensionalization of the convection-diffusion-decay process
The final objective of this paper is to assimilate results of tests 1, 2 and 3 to provide a unique non-dimensional time at which tissue solute concentration is maximum. The solute concentration-time histories in specific radial locations of the tissue space are influenced by the drainage of the solutes in addition to the advection diffusion and decay processes as modelled by Eq. 1. S is the density of the micro vessels (unit is no of vessels/m) whose walls act as a sink in the tissue volume. The product of S and solute permeability (P s d ) is the decay constant k. We assume exponential decay and define τ as the time when the inlet concentration drops to 36.7% (1/e) of the maximum inlet concentration. The first form of Eq. 1 is non-dimensionalized to find the time scales of the other transport mechanisms involved. Defining the following where Uo is the average velocity of fluid in tissue. Substituting these in Eq. 1; we get the non-dimensional form as shown in Eq. 2.
The effect of τ on the solute accumulation in a tissue space is modified by the number of microvessels (n) present around it and the intercapillary separation (L/d) between them. The solute dependent timescales, obtained from Eq. 2 are the modified input time scale ( nτ L d ), the diffusion timescale ( L 2 D tissue ) and the decay timescale ( 1 k ). These values for each solute in SBC configuration are shown in Table 5. Since different timescales are dominant at different phases and radial locations of transport for differing solute types, a sum of all the solute dependent time scales is used to non-dimensionalize the time of solute accumulation and decay in tissue as shown in Eq. 3.
The extravascular concentration is rescaled to account for the variable solute molecular weight (Mw), solute density (ρ), tissue porosity ∅, inlet solute concentration (Co) and varying intercapillary separation (L/d) as defined by Eq. 4.

Mixture theory model equations
The conservation of mass and linear momentum equations from which the mixture theory equations are derived are stated below. The equation derivation steps are described in [13,14] by Schuff et al. A representative tissue configuration is shown in Fig 7.

Balance of Mass
Assuming zero rate of mass production of constituent α with density ρ α and velocity u α in a mixture, the general expression for mass balance at each region is ∂ρ α ∂t þ ∇:ðρ α u α Þ ¼ 0. u α is related to the extravascular/ Darcy velocity of the same constituent as U α ¼ ϕ α u α where ϕ α is the volume fraction of constituent α. In each region the mixture was assumed to be saturated (i.e X α ϕ α ¼ 1). For the extravascular space, assuming incompressible flow and a stationary solid matrix, the mass conservation equation reduced to ∇:U f ¼ ∇:ðϕ f u f Þ ¼ 0. In our simulations due to assumption of dilute solution, ϕ f = 1, ϕ s = ϕ.

Balance of Linear Momentum
Assuming absence of body forces per unit volume of any tumor space, the momentum equation for the constitu- is the material derivative, T α ¼ −ϕ α PIþϕ αT α is the Cauchy stress tensor, P = pore pressure andT α is the stress tensor for a dense configuration of constituent α.

For fluids with viscosity
Assuming a dilute solution of solutes in the tumor, solute-solute interactions are negligible for constituent s, T s ¼ 0 . ∏ α accounts for the interaction terms between

Fluid transport Equations
Intravascular Space: The arterial pressure, Par, is specified at the inlet while the prescribed hydrostatic pressure difference, ∇P balances the viscous stresses to govern the fluid flow inside the capillaries according to Eq. 5.
Capillary Wall: At the capillary wall, jump boundary conditions are applied as ½½ϕ f ðu f Þ:n ¼ 0. The [[a]] denotes the difference of the quantity "a" across the walls where n is the unit normal outward to the wall. The fluid flux, qe ¼ ð ϕ f u f Þ:n, across the capillary wall (of radius Ro) as shown in the final form (Eq. 6) accounts for a) the hydrostatic pressure difference across the capillary wall due to the hydraulic conductivity,L p , b) a constant osmotic pressure gradient, σ p ðπ p i −π p e Þ , due to protein molecules; and c) a variable osmotic pressure that depends on the concentration difference of the injected solute (of molecular weight M s w and density ρ s T ) in the intravascular (i) and the extravascular (e) space of a fibrous matrix with porosity ∅.
Extravascular Space: The fluid transport in the extravascular space is influenced by a) the hydrostatic pressure difference in the tissue space, b) the hydraulic permeability of the tissue, k and c) the solute concentration gradients, ∇c s in the tissue space as depicted in Eq. 7.
The retardation factor, R s F , is the ratio of the solute velocity and the inline fluid velocity in the tissue space. A is the interaction coefficient of the solute "s". A constant pressure Po is applied to the tissue boundary.
Solute transport Equations Intravascular Space: The initial solute concentration, C o , is used to prescribe the concentration at the inflow according to Eq. 8. The concentration time history for the five solutes prescribed at the blood vessel inlet is shown in Fig. 8. The concentration has been normalized by the peak intravascular concentration of each solute. Capillary wall: The solute particles have the same jump boundary conditions as the fluid at the capillary wall. The solute particles carried by the fluid flux and those which permeate into the tissue space due to the concentration difference across the capillary wall constitute the solute flux (J s ) across it. Equation 10 is a modified version of Starling's law where P s d is the solute permeability coefficient and σ s is the reflection coefficient.
Js ¼ ϕ f c s u s :n ¼ P s d c s i jr ¼ Ro−c s e jr ¼ Ro Extravascular Space: In the porous matrix of the extravascular space the solute undergoes both advection and diffusion as shown in the governing transport Eq. 11.
The mixture theory equations are appropriate for this study since it shows the dependence of tissue mechanical properties like hydraulic conductivity on chemical gradients that is not captured by traditional transport models [13].

Difference of lymph vessel input parameters from that of blood vessel
Solute concentration is zero at the lymph capillary inlet. The osmotic pressure gradient due to protein molecules is absent in the lymph vessel, ðσ p ðπ p i −π p e Þ ¼ 0Þ. The solute permeability coefficient, P s d , across the lymph capillary wall is twice its value in blood capillary wall shown in Table 1 to account for the free permeability of lymph vessels to macromolecules. The fluid flux and the solute flux equations are modified to allow intravasation only.

Fabrication of microfluidic platform
Type I collagen was used as the extracellular matrix of the tumor with a single integrated endothelialized blood vessel. Excised tendons from rat tails were dissolved in a pH 2.0 HCl solution for 12 h at 23°C. The solution was centrifuged at 30000 g for 45 min and sterilized using 10% (v/v) chloroform for 24 h at 4°C. The mold for the in vitro tumor microfluidic platform with the embedded single vessel was fabricated as described in previous work [75,76]. PDMS housing was fabricated using common soft-lithography methods. Polydimethylsiloxane (PDMS) and curing agent was mixed with 10:1 ratio and baked at 75°C for 1 h. Hardened PDMS housing and the glass cover slip was plasma treated for 18 W for 30 s. Plasma treated surfaces were assembled to create a permanent bonding. The housing was treated with 1% (v/v) polyethyleneimine in dH 2 O for 10 min followed by 0.1% (v/v) glutaraldehyde in dH 2 O for 20 min and washed with dH 2 O twice. Collagen solution of 7 mg/ml was prepared by neutralizing stock solution with 1X DMEM, 10X DMEM, 1 N NaOH, and mixing 1 × 10 6 /ml MDA-MB-231 breast cancer cells uniformly in collagen which was then placed in the housing. A 22G (711 μm) needle was inserted into the mold and after polymerization and the needle removal, a cylindrical vascular channel was created within the collagen. 2 × 10 6 TIME cells were injected into the vascular channel and exposed to flow preconditioning protocols for 3 days to form a confluent, aligned endothelialized vessel. As a result, in vitro platform shown in Fig. 9 was fabricated.

Measuring mechanical properties from fabricated tissue platform
Porosity Vasculature porosity was measured using confocal microscopy and imaging mKate labeled TIME cells. Scanning electron microscopy (Zeiss, Super40) was used  to image 3 fibrous matrix samples under 15kx, 20kx, and 25kx magnifications at three vertical planes. Tissue porosity was measured by applying a Frangi filter on obtained images.
Solute Permeability Obtained intensity profiles of 3 kDa and 70 kDa Dextran particles were used to calculate solute permeability coefficient as shown previously [77]. 3 kDa and 70 kDa dextran particles were suspended in serum free endothelial basal medium at 10 μg/ml concentration and perfused in the vascular channel at 260 μL/min, which corresponds to 1 dyne/cm 2 physiological shear stress for tumor vasculatures at every 3 min for 2 h. The transport of these solutes was imaged using a confocal microscope (Leica SP8, 10X magnification). Normalized intensity profiles from this images as a function of time were used to compare with normalized concentration profiles from the equivalent numerical simulations.
Solute Diffusivity Fluorescence recovery after photobleaching (FRAP) technique was used as described previously by Voigt et al. to measure the diffusion coefficient for a range of dextran molecular weights (4 kDa-150 kDa) for varying pH values, collagen concentrations and temperatures [78,79]. We selected diffusivity values of 3 kDa and 70 kDa for pH 7.6, collagen concentration 7 mg/ ml at 37°C from the database of the mentioned study.
Hydraulic tissue permeability Hydraulic permeability for collagen at a concentration of 7 mg/ml as is used in fabrication of the platform were collected from the existing literature on vascularized in vitro experiments [42,80,81].