Validation and parameter optimization of a hybrid embedded/homogenized solid tumor perfusion model
Johannes Kremheller, Sebastian Brandstaeter, Bernhard A. Schrefler, Wolfgang A. Wall
NNoname manuscript No. (will be inserted by the editor)
Validation and parameter optimization of a hybridembedded/homogenized solid tumor perfusion model
Johannes Kremheller · Sebastian Brandstaeter · Bernhard A. Schrefler · Wolfgang A. Wall
Received: date / Accepted: date
Abstract
The goal of this paper is to investigate thevalidity of a hybrid embedded/homogenized in-silico approach for modeling perfusion through solid tumors.The rationale behind this novel idea is that only thelarger blood vessels have to be explicitly resolved whilethe smaller scales of the vasculature are homogenized.As opposed to typical discrete or fully-resolved 1D-3Dmodels, the required data can be obtained with in-vivo imaging techniques since the morphology of the smallervessels is not necessary. By contrast, the larger vessels,whose topology and structure is attainable non-inva-sively, are resolved and embedded as one-dimension-al inclusions into the three-dimensional tissue domainwhich is modeled as a porous medium. A sound mortar-type formulation is employed to couple the two repre-sentations of the vasculature. We validate the hybridmodel and optimize its parameters by comparing its re-sults to a corresponding fully-resolved model based onseveral well-defined metrics. These tests are performedon a complex data set of three different tumor types
J. Kremheller · S. Brandstaeter · W. A. WallInstitute for Computational Mechanics, Technical Universityof Munich,Germany, D-85748E-mail: [email protected],E-mail: [email protected],E-mail: [email protected]. A. SchreflerInstitute for Advanced Study, Technical University of Mu-nich,Germany, D-85748B. A. SchreflerDepartment of Civil, Environmental and Architectural Engi-neering, University of Padova,ItalyE-mail: bernhard.schrefl[email protected] with heterogeneous vascular architectures. The corre-spondence of the hybrid model in terms of mean repre-sentative elementary volume blood and interstitial fluidpressures is excellent with relative errors of less than4 %. Larger, but less important and explicable errorsare present in terms of blood flow in the smaller, ho-mogenized vessels. We finally discuss and demonstratehow the hybrid model can be further improved to ap-ply it for studies on tumor perfusion and the efficacy ofdrug delivery.
Keywords microcirculation · tissue perfusion · homogenization · hybrid models · Mathematical modeling of blood flow and mass trans-port is of increasing importance to study a number ofhighly relevant biomedical questions in health and dis-ease. Computational tools offer the possibility to gainnew insight into physiologically relevant processes suchas the transport of nutrients, oxygen or drugs acrossthe vascular system and inside the tissue micro-environ-ment. These methods can ultimately lead to a new ra-tionale for developing and non-invasive testing of noveltherapies (Dewhirst and Secomb 2017). Concurrently,such in-silico models will make the design of drugs bothcheaper and faster.In this paper, we are concerned with the simulationof blood flow and tissue perfusion at the scale of themicrocirculation with a special focus on solid tumorswhere transport processes can be decisive for diseaseprogression and treatment efficacy. This includes, first,the vasculature, which is embedded in the surroundingtissue, second, passage across the blood vessel walls intothe surrounding extravascular space and, third, flow a r X i v : . [ c s . C E ] F e b Johannes Kremheller et al. of the fluid filling this space, namely the interstitialfluid (IF). Subsequently, we will distinguish betweenthree different modeling strategies for these transportprocesses, namely discrete, continuum and hybrid ap-proaches. All strategies have been developed and imple-mented in the context of our vascular multiphase tumorgrowth model (Kremheller et al. 2018; Kremheller et al.2019). For the discrete or fully-resolving variant, we fol-low a common modeling approach, where the vascu-lature is represented by a network of one-dimensionalblood vessel segments embedded in the encompassingthree-dimensional tissue domain which is modeled asa porous medium. A 1D partial differential equation(PDE) is employed to model mass transport in the vas-culature while a corresponding 3D PDE governs thesurrounding IF. Both domains are coupled via sourceterms which account for the exchange across the bloodvessel wall. Such models are well-studied and first con-tributions include the so-called embedded multiscalemethod developed by D’Angelo (2007), D’Angelo andQuarteroni (2008) and D’Angelo (2012) and the Green’sfunction method of Secomb and Hsu (1994) and Sec-omb et al. (2004). More recent approaches with such aphilosophy include drug delivery (Cattaneo and Zunino2015, 2014), hyperthermia treatment (Nabil et al. 2015;Nabil and Zunino 2016) and a combination of a numeri-cal framework with optical imaging to predict fluid andspecies mass transport through whole tumors with het-erogenous blood vessel architecture (D’Esposito et al.2018; Sweeney et al. 2019). These approaches are com-monly termed discrete models in distinction from con-tinuum models which involve a homogenization proce-dure of the vascular network. Thereby, the vasculatureis approximated as a homogeneous porous medium re-sulting in two distinct pore spaces which are the afore-mentioned interstitium and the homogenized vascula-ture. Flow in both domains is modeled via the Darcyequation and suitable exchange terms are defined (Chap-man et al. 2008; Shipley and Chapman 2010; Tully andVentikos 2011; Kremheller et al. 2018).These two distinct approaches have different usecases: Discrete models can and should be applied whenthe entire structure of the vasculature including thesmallest scales, i.e., the capillaries, is known and itsresolution is needed for the question at hand. This isusually restricted to small domain sizes of an orderof several mm . By contrast, continuum models areused to simulate mass transport at larger scales, e.g.,through whole organs. Both approaches have advan-tages and disadvantages: On the one hand, the com-putational cost of continuum models is usually smallerthan for discrete ones which makes the application tolarger domains possible in the first place. On the other hand, the information about the exact morphology ofthe vascular network is lost such that blood flow canonly be described in an averaged sense. Discrete mod-els, however, are computationally more expensive. Fur-thermore, they require the full structure of the partof the vasculature under consideration. This is usuallyrealized via a graph whose edges are assigned the ra-dius of the blood vessel segments between nodes. Suchhigh-resolution data including blood vessel radii, con-nectivity and positions can at present only be acquiredthrough ex-vivo imaging (Shipley et al. 2019). In addi-tion, the acquisition of high-quality data is still chal-lenging and error-prone especially on the finest scales(K¨oppl et al. 2020). By contrast, in-vivo imaging is cur-rently only possible for larger vessels and flow therein(Shipley et al. 2019; Li et al. 2020). Therefore, dis-crete models rely on data which is not available vianon-invasive imaging. An additional difficulty is the as-signment of blood pressure or flow boundary condi-tions which can only be estimated for large networks(Sweeney et al. 2019; Fry et al. 2012). In any case, val-idation of these models is usually only performed onmacroscopic quantities such as tissue perfusion (D’Es-posito et al. 2018) since measuring flow or pressures in-side single micro-vessels is not possible (Sweeney et al.2019).This has motivated the development of hybrid meth-ods which are especially suited for cases where the fullvascular morphology is unknown or too large to be mod-eled with a discrete approach. The idea behind them isto explicitly resolve the larger vessels through a dis-crete model and to use a homogenized approach for thecapillary bed. Next to our own work (Kremheller et al.2019), such hybrid approaches have also been devel-oped by Kojic et al. (2017), Vidotto et al. (2019) andShipley et al. (2019). Compared to pure homogenizedformulations, their advantage is that the structure ofthe larger vessels is retained and, therefore, the hetero-geneity of blood flow and pressure in the major vesselbranches is better represented. Moreover, compared todiscrete models, less anatomic data is needed since themorphology of the smallest vessels is not required. Thiscould also have the additional advantage of a smallercomputational cost and make them applicable to largerdomains. Also quantities typically needed for valida-tion such as tissue perfusion, blood flow or pressures atthe resolution of current imaging techniques can equallybe acquired from hybrid models. A related approach,where no homogenization of the capillaries is needed, isto generate a discrete surrogate network of the smallerscales based on the oxygen demand of the tissue (K¨opplet al. 2020). alidation and parameter optimization of a hybrid embedded/homogenized solid tumor perfusion model 3 We have previously incorporated a hybrid methodfor coupling discrete, one-dimensional blood vessels witha homogenized representation of the vasculature (Krem-heller et al. 2019) into our vascular multiphase tumorgrowth model (Kremheller et al. 2018). Therein, we cou-ple a discrete representation of the pre-existing vascula-ture with a homogenized representation of the neovas-culature which is formed during angiogenesis. For that,we employ constraint enforcement strategies which arewell-known from solid mechanics. In our previous pa-per, the main focus was on modeling vascular tumorgrowth. In this contribution, we validate the applica-bility of the hybrid embedded/homogenized approachfor the study of perfusion through solid tumors. Here,accurate models of fluid mass transport are of high rel-evance since efficient drug delivery to cancerous tissuerelies on the fact that the drug reaches a large fractionof the tumor cells. Therefore, physiological character-istics such as microvascular flow, the structure of theextracellular matrix or the IF pressure profile may influ-ence the transport of drugs through tumor tissue and,hence, ultimately the success of treatment (Dewhirstand Secomb 2017). For instance, increased interstitialpressure due to highly permeable vessels and inefficientlymphatic drainage has been identified as an obsta-cle for successful drug delivery (Baxter and Jain 1989;Jain 1994; Heldin et al. 2004). Novel nanoparticle-basedtherapies aim to exploit these properties of the tumorvasculature for more specific targeting of tumor sites(Matsumura and Maeda 1986). Appropriate models ofthese transport phenomena can provide additional in-sight into and guidelines for drug design. In the contextof cancer, this paradigm shift is described by the con-cept of transport oncophysics with the objective to engi-neer drugs with optimized transport properties (Ferrari2010; Nizzero et al. 2018).Sound computational models are required to achievethis goal. We therefore took great care in the develop-ment of our hybrid model, i.e., both in the theoret-ical basis and its implementation. In this paper themain focus is on the validation of our hybrid embed-ded/homogenized scheme with three complex tumor-specific vascular networks based on large tissue samplescontaining more than 100 000 blood vessels (D’Espositoet al. 2018; Sweeney et al. 2019). We put a special em-phasis on the extraction of the larger vessels from thefully-resolved network data such that it qualitativelymatches the topology and distribution of larger vascularstructures inside tumors available via in-vivo imaging.Thus, we make sure that the hybrid approach is inves-tigated for cases which closely resemble real-life scenar-ios where the structure of the considered part of themicrocirculation is not entirely known. Here, the com- plete topology of the vasculature in the given tissue do-main is available which allows us to generate referencesolutions with a fully-resolved model and to quantifythe error introduced by the homogenization in the hy-brid model. We evaluate the error by means of severalwell-defined metrics involving the agreement of pres-sures and flow between the two models. Concurrently,the parameters of the hybrid model are identified suchthat the correspondence of the models is maximized.Evaluating the error of the hybrid model in comparisonto a fully-resolved one is a first and indispensable steptowards realistic hybrid models of tumor perfusion rely-ing only on non-invasively available physiological data.For a full validation and parameter optimization, simi-lar methods as applied herein need to be combined withadvanced in-vivo imaging techniques. The comparisonof two purely numerical approaches allows us to investi-gate the hybrid model in a controlled environment un-affected by any further influences such as uncertaintiesin experimental or clinical data.The remainder of this work is structured as follows:We introduce both the hybrid and the fully-resolvedmodel in Section 2. The employed tumor vasculaturedata sets as well as the setup of the models includingthe assignment of boundary conditions and the extrac-tion of the hybrid model from the fully-resolved one aredescribed in Section 3. Numerical experiments to com-pare the accuracy of the hybrid model w.r.t. the fullmodel and to evaluate its main errors are conducted inSection 4. We illustrate some possible improvements ofthe hybrid model in Section 5 before summarizing ourfindings in Section 6.
In this section, we describe the mathematical modelswe employ to solve the interaction between microcir-culation and interstitial tissue perfusion including theirmain simplifications. We outline both our fully-resolvedand our hybrid approach and their discretization bymeans of the finite element method (FEM).2.1 Problem settingAs in other publications (Secomb and Hsu 1994; Sec-omb et al. 2004; D’Angelo 2007; D’Angelo and Quar-teroni 2008; D’Angelo 2012; Cattaneo and Zunino 2014,2015; Nabil et al. 2015; Nabil and Zunino 2016; D’Espo-sito et al. 2018; Sweeney et al. 2019; Kremheller et al.2019), topology and structure of the microcirculationis described by a graph with straight edges, i.e., bloodvessel segments. The segments connect the nodes of
Johannes Kremheller et al.
ΩΛ ∂ΩΛ S Λ L (a) Fully-resolved ΩΛ L Ω v ∂Ω v ∂Ω (b) Hybrid Fig. 1: Notation for domains and boundariesthe network. A radius R k is assigned to each segment Λ k . Available experimental data including the one em-ployed here is also commonly provided in this format.Therefore, the vascular domain is given as the unionof straight cylinders which are embedded in the three-dimensional tissue domain Ω , see also Figure 1a. Basedon the flow physics and the quantities of interest, it isobvious to employ a one-dimensional blood flow modelfor this part of the vascular system. In the following,we will denote the 1D embedded blood vessel networkwith Λ . Similar to Vidotto et al. (2019), we further di-vide it into two subsets Λ L and Λ S which correspond tothe larger and smaller vessels in the network such that Λ L := (cid:91) k ∈ I L Λ k , Λ S := (cid:91) k ∈ I S Λ k and Λ = Λ L ∪ Λ S (1)with the index sets of large and small blood vessel seg-ments I L and I S , respectively. We will show in detailhow this partition is realized in Section 3.3. Whereasthe larger vessels are kept in the hybrid model, thesmaller scale vessels are replaced by an appropriate ho-mogenized representation as a porous network occupy-ing the domain Ω v ⊂ Ω , cf. Figure 1b.2.2 Fully-resolved 1D-3D modelGravity, inertial effects and the pulsatility of blood floware neglected which are valid assumptions since we dealwith microcirculatory flow. The balance of mass in the1D vasculature domain Λ is then given by the followingequation − ∂∂s (cid:18) πR µ ˆ v ∂p ˆ v ∂s (cid:19) = − ˆ v → l M leak ρ ˆ v on Λ, (2) where we have applied the Hagen-Poiseuille law for flowin cylindrical pipes and assumed constant blood den-sity ρ ˆ v . Here and in the following, quantities definedon the 1D vasculature domain are denoted by super-script ˆ v . In the previous equation, R is the blood vesselradius, p ˆ v the pressure inside the vasculature, µ ˆ v theblood viscosity and s the arc length coordinate alongthe 1D blood vessel segment. To account for the non-Newtonian behaviour of blood, we employ the algebraicrelationship developed by Pries and Secomb (2005) for in-vivo blood viscosity depending on blood vessel diam-eter and hematocrit. As in Sweeney et al. (2019), we fixthe hematocrit to 0 .
45, thus, the blood viscosity µ ˆ v ineach individual blood vessel segment depends only onits diameter. Finally, the right-hand-side term ˆ v → l M leak = ρ l · πR · L p, ˆ v · (cid:0) p ˆ v − p l − σ (cid:0) π b − π l (cid:1)(cid:1) (3)is employed to model leakage of fluid across the bloodvessel wall into the interstitium. For that, we use Star-ling’s law with hydraulic conductance L p, ˆ v , density ofblood plasma ρ l , oncotic reflection coefficient σ and theoncotic pressures of blood π b and the interstitial fluid(IF) π l . In summary, the transvascular flux from thevascular network into the interstitial fluid is propor-tional to the pressure difference between vasculatureand IF whose pressure in (3) is denoted as p l . It haslong been known that blood vessels inside tumors areleakier than normal ones which, in combination with anon-functional lymphatic system, leads to increased in-terstitial pressure inside solid tumors and, concurrently,resistance to efficient drug delivery (Baxter and Jain1989; Jain 1994; Heldin et al. 2004). Note that our datasets are whole-tumor blood vessel networks where also alidation and parameter optimization of a hybrid embedded/homogenized solid tumor perfusion model 5 larger vessels are leaky (D’Esposito et al. 2018; Sweeneyet al. 2019) which is why we apply the transvascular ex-change term (3) also on the subset of larger vessels Λ L .As in related works, the tissue domain Ω is modeledas a porous medium. Therefore, flow in the interstitialfluid is accounted for by the following Darcy equation − ∇ · (cid:32) k l µ l ∇ p l (cid:33) = δ Λ · ˆ v → l M leak ρ l in Ω (4)with (isotropic) permeability k l = k l · I and IF viscosity µ l . Hence, the primary variable for fluid flow throughthe tissue is the IF pressure p l . The right hand siderepresents the counterpart of the leakage of fluid fromthe vasculature into the IF from (2). As proposed by(D’Angelo 2007; D’Angelo and Quarteroni 2008; D’An-gelo 2012) this mass transfer term is concentrated asa Dirac measure δ Λ along the centerline Λ of the vas-culature resulting in a 1D-3D coupled problem. Themathematical properties including reduced convergencerates due to the introduced singular line source in the3D pressure field are extensively studied by D’Angeloand Quarteroni (2008), D’Angelo (2012) and K¨oppl andWohlmuth (2014). Alternative 2D-3D coupled approa-ches, where the mass exchange is evaluated at the lat-eral surfaces of the cylindrical blood vessel segments,have been proposed to increase the regularity of the so-lution (K¨oppl et al. 2018, 2020). However, in our datasets the diameter D is smaller than the element size h in the 3D domain, see also Table 1. Therefore, we mod-ify the approach of (D’Angelo and Quarteroni 2008;D’Angelo 2012), which would involve taking the aver-age value p l of the pressure in the IF at the outer sur-face of the cylindrical vessels in the exchange term (3),and instead take the IF pressure value at the centerline Λ , which is a reasonable approach for the case h > D (D’Angelo 2007; Kremheller et al. 2019). This has re-cently also been investigated in the analogous solid me-chanics problem of embedding thin 1D structures, i.e.,beams, into 3D solid volumes (Steinbrecher et al. 2020).The weak form of the 1D-3D coupled problem may bewritten as (cid:18) ∂δp ˆ v ∂s , πR µ ˆ v ∂p ˆ v ∂s (cid:19) Λ + δp ˆ v , ˆ v → l M leak ρ ˆ v Λ = 0 (5a) (cid:32) ∇ δp l , k l µ l ∇ p l (cid:33) Ω − δp l , ˆ v → l M leak ρ l Λ = 0 (5b)with test functions δp ˆ v defined on the 1D domain and δp l defined on the 3D domain. Our approach allows fornon-matching 1D discretizations Λ h and 3D discretiza-tions Ω h such that the two domains can be meshed independently of each other. This requires the numeri-cal integration of products of 1D shape functions with3D shape functions and products of 3D shape functionswith 3D shape functions along the one-dimensional dis-cretization Λ h which we realize via a segment-basedline integration approach (Kremheller et al. 2019; Stein-brecher et al. 2020) to avoid integration over kinks ofshape functions. After space discretization, the nodalprimary variables of both domains are p ˆ v ∈ R n nodes ,Λ and p l ∈ R n nodes ,Ω , (6)that is, the nodal blood pressure in the discretized 1Ddomain and the nodal IF pressure in the discretized 3Ddomain, which consist of n nodes ,Λ and n nodes ,Ω , respec-tively. Details on the employed boundary conditions aregiven in Section 3.2.1.Finally, we arrive at the global system of equations,which may be written as a 2 × (cid:20) K ˆ v ˆ v G ˆ vl H l ˆ v K ll (cid:21) (cid:20) p ˆ v p l (cid:21) = (cid:20) F ˆ v F l (cid:21) . (7)Herein, the main diagonal blocks K ii comprise contri-butions from the diffusive term and the exchange termin (5a) and (5b) while the off-diagonal submatrices G ˆ vl and H l ˆ v contain the ”mixed” contributions from theexchange term. The right hand side terms F i representthe constant contribution to the exchange term stem-ming from the oncotic pressures in (3). To solve thecoupled linear system (7) we employ a GMRES itera-tive solver in combination with the AMG(BGS) blockpreconditioner presented in Verdugo and Wall (2016).2.3 Hybrid 1D-3D modelThe main idea behind our hybrid 1D-3D model, basedon our previous work introduced in Kremheller et al.(2019) in the context of our vascular multiphase tumorgrowth model (Kremheller et al. 2018), is the follow-ing: The full resolution of the larger vessels Λ L is kept,i.e., these are still modeled as a 1D embedded vascu-lature. Consequently, the hierarchy, topology and vas-cular properties such as individual blood vessel radiiand viscosities of each segment are retained, see alsoFigure 1b. The smaller vessels Λ S , for which this high-resolution data might either not be available throughnon-invasive imaging techniques or susceptible to er-rors, are instead represented as an additional porousnetwork. This results in a double-porosity formulationwhere the first porous network is, as before, the intersti-tial space and the second one the smaller vessels occu-pying the domain Ω v . In the following, we will present Johannes Kremheller et al. the governing equations and the space discretization ofthis formulation.As stated above, the model for the larger vesselsdoes not change. Therefore, the mass balance equationinside the large vessels is given by − ∂∂s (cid:18) πR µ ˆ v ∂p ˆ v ∂s (cid:19) = − ˆ v → l M leak ρ ˆ v on Λ L (8)with the only difference to (2) being that it holds onlyon the subset Λ L ⊂ Λ of bigger vessels. The mass bal-ance equation for the smaller vessels Λ S is replaced bya homogenized Darcy equation in the vascular domain Ω v , which we formulate as − ∇ · (cid:18) k v µ v ∇ p v (cid:19) = − v → l M leak ρ v in Ω v . (9)The unknown in this equation is the blood pressure p v in the homogenized part of the vasculature whichis now defined in the entire 3D domain Ω v , therebyreplacing the blood pressure of the smaller vessels inthe 1D domain Λ S as the variable governing flow insidethe smaller vessels. For simplicity, in a first step weconsider an isotropic permeability tensor k v = k v · I for the additional porous network. This permeabilityand the averaged blood viscosity µ v are the two modelparameters governing this equation together with theright hand side term v → l M leak = (cid:40) ρ l · L p,v ( S/V ) Λ S · (cid:0) p v − p l − σ (cid:0) π b − π l (cid:1)(cid:1) in Ω v Ω \ Ω v . (10)This term replaces the outflow of fluid from the smallervessels into the IF by a homogenized representationof the Starling equation (3) involving the surface-to-volume ratio of the smaller blood vessels ( S/V ) Λ S asan additional parameter. The mass balance equation ofthe IF for the fully-resolved model (4) is adapted as − ∇ · (cid:32) k l µ l ∇ p l (cid:33) = 1 ρ l (cid:18) δ Λ L · ˆ v → l M leak + v → l M leak (cid:19) in Ω (11)in the homogenized formulation. Comparing the twoequations, it is obvious that leakage from the large ves-sels is still treated equivalently, i.e., the large vesselsare still embedded as 1D inclusions in the tissue with aDirac measure (now defined only on Λ L ). By contrast,leakage from the smaller blood vessels is replaced by thehomogenized mass transfer term (10) from the vasculardomain Ω v into the interstitium, i.e., from (9) into (11). So far, this procedure is analogous to other hybridapproaches (Shipley et al. 2019; Vidotto et al. 2019).The main difference to our methodology lies in the cou-pling between the larger vessels Λ L and the homoge-nized vasculature Ω v . In the aforementioned publica-tions, this was realized at the free ends of the largervessels, i.e., as an outflow at the tips of the 1D dis-cretization into the homogenized 3D vasculature do-main. This was possible since the employed data setshad a clear vascular hierarchy with larger arterioles andvenules connected to smaller capillaries. Our vascularnetworks, which we will describe in detail in Section 3.1,have been segmented from solid tumors and, therefore,have a much more complex, disorganized structure in-cluding variable vessel lengths and diameters as well asdead ends. All this is typical for tumor-specific vascu-lature (Carmeliet and Jain 2000; Baluk et al. 2005). Asshown in detail in Section 3.3 for our data and the em-ployed methodology to distinguish between large, flow-carrying vessels and smaller ones, another approach ismore sensible: We enforce the coupling between largervessels and the homogenized vasculature along the en-tire 1D representation of larger vessels Λ L with a line-based coupling instead of a point-based coupling atthe tips of the larger vessels flowing into the capillarybed as described before. Compared to these hybrid ap-proaches, our proposed method has the advantage thatno additional parameter – apart from the penalty pa-rameter – is involved for the coupling of the two repre-sentations.For that, we formulate a constraint of equal pres-sures in Λ L and Ω v as g = p ˆ v − p v = 0 on Λ L , (12)which enforces a coupling between pressures p ˆ v in theone-dimensional, large vessel domain Λ L and homoge-nized pressures p v in the 3D domain Ω v . We aim toreproduce the fact that the pressure in a smaller ves-sel branching from a larger vessel at a specific node isequal to the pressure at the same node. If this smallervessel is homogenized and, thus, removed from the 1Drepresentation, we want to enforce these equal pres-sures between the resolved part and the homogenizedpart of the vasculature along the 1D vessel domain Λ L .In Section 3.3, we justify formulating this constraintalong the entire 1D domain Λ L considering the connec-tivity between larger and smaller vessels in our cases.We have previously employed a similar strategy in ourhybrid treatment of the vasculature in a multiphase tu-mor growth model (Kremheller et al. 2019) and therelated solid mechanics problem of beam-to-solid meshtying (Steinbrecher et al. 2020). We follow the same ap-proach as in the two aforementioned publications and alidation and parameter optimization of a hybrid embedded/homogenized solid tumor perfusion model 7 incorporate the constraint with an additional Lagrangemultiplier (LM) field into the weak form of our hybridmodel, which reads as (cid:18) ∂δp ˆ v ∂s , πR µ ˆ v ∂p ˆ v ∂s (cid:19) Λ L + δp ˆ v , ˆ v → l M leak ρ ˆ v Λ L + (cid:0) δp ˆ v , λ (cid:1) Λ L = 0 (13a) (cid:18) ∇ δp v , k v µ v ∇ p v (cid:19) Ω v + δp v , v → l M leak ρ v Ω v − ( δp v , λ ) Λ L = 0 (13b) (cid:32) ∇ δp l , k l µ l ∇ p l (cid:33) Ω − δp l , ˆ v → l M leak ρ l Λ L − δp l , v → l M leak ρ l Ω v = 0 (13c) (cid:0) δλ, (cid:0) p ˆ v − p v (cid:1)(cid:1) Λ L = 0 (13d) Therein, the first line is the weak form of flow in thelarger vessels (8) which is coupled to the weak form offlow in the homogenized vasculature domain, i.e., thesecond line (13b) with a continuous LM field λ definedalong the blood vessel center line. The third line is theweak form of flow in the IF. Compared to the fully-resolved model, conf. eqn. (5b), the additional masstransfer term arises due to leakage from the homog-enized part of the vasculature into the IF. The fourthline represents the variational form of the coupling con-straint (12). Conveniently, the LM field employed to en-force this constraint can then be interpreted as a masstransfer term from the 1D resolved bigger vessels intothe 3D homogenized vasculature, i.e., λ = ˆ v → v M . Alter-natively, a Gauss-point-to-segment scheme could alsobe employed but suffers from over-constraining of thesystem for large penalty parameters (Kremheller et al.2019; Steinbrecher et al. 2020). Spatial discretizationof the weak form (13a)-(13d) leads to a saddle-pointproblem with nodal primary variables p ˆ v ∈ R n nodes ,Λ L , λ ∈ R n nodes ,Λ L , p l ∈ R n nodes ,Ω , p v ∈ R n nodes ,Ωv , (14)that is, nodal pressures and nodal LMs in Λ L ,h , nodalIF pressures in Ω h and nodal blood pressures of thehomogenized vasculature in Ω v,h . In the following, wewill specifically focus on the discretization of the termsarising due to the LM method. Approximating thosecontributions with a finite element interpolation yieldsa mortar-type formulation where the nodal LMs areadditional degrees of freedom, condensed out with adual approach (Wohlmuth 2000; Popp et al. 2010) ora penalty regularization of the mortar method is em-ployed to remove the additional degrees of freedom andthe saddle-point structure (Yang et al. 2005). Here, wefollow the latter approach just as in our previous workon 1D-3D type couplings (Kremheller et al. 2019; Stein-brecher et al. 2020). The contributions to the weak formof the mass balance equations, i.e., the two last terms in (13a) and (13b) can be written as δΠ LM ,h = n nodes ,Λ L (cid:88) j =1 n nodes ,Λ L (cid:88) k =1 λ j D jk δp ˆ vk − n nodes ,Λ L (cid:88) j =1 n nodes ,Ωv (cid:88) l =1 λ j M jl δp vl (15)with the so-called mortar matrices D [ j, k ] = D jk = (cid:90) Λ L ,h ˆ Φ j ˆ N k ds (16)and M [ j, l ] = M jl = (cid:90) Λ L ,h ˆ Φ j N l ds. (17)The entries of these matrices involve integrals of prod-ucts of LM shape functions ˆ Φ j defined on the discretized1D domain Λ L ,h with 1D shape functions ˆ N k and with3D shape functions N l defined in the 3D domain Ω v .Hence, these terms are again evaluated using a segment-based approach. We choose linear shape functions forboth primary variables and the LM interpolation, i.e.,ˆ Φ j = ˆ N j . The weak form of the constraint (13d) maythen be written in discretized form as δΠ λ,h = δ λ T (cid:0) Dp ˆ v − Mp v (cid:1) =: δ λ T g (cid:0) p ˆ v , p v (cid:1) , (18)where we have defined a weighted pressure gap g ateach node in Λ L ,h . This gap is then further used for thepenalty regularization of the mortar method to explic-itly define the nodal LMs in terms of 1D and 3D nodalblood pressures as λ = (cid:15) κ − g (cid:0) p ˆ v , p v (cid:1) . (19)Hence, the LMs are no longer independent variables inthe system but depend on the primary variables p ˆ v and p v . This overcomes the two major drawbacks of theLM method, namely, the increased system size and itssaddle-point structure. Depending on the penalty pa-rameter (cid:15) >
0, the constraint g = is relaxed and theexact solution is only recovered for (cid:15) → ∞ . Addition-ally, the nodal LM in (19) has been scaled with theinverse of the diagonal matrix κ [ j, j ] = (cid:90) Λ L ,h ˆ Φ j ds. (20)As proposed by Yang et al. (2005) this removes thedependency of the nodal LM on its ”gap”, i.e., in ourcase it makes its entries independent of the elementlengths associated with its corresponding node. Thiscan now be used to replace the LM vector such that Johannes Kremheller et al. the matrix-vector form of our hybrid model emerges as K ˆ v ˆ v + (cid:15) D T κ − D G ˆ vl − (cid:15) D T κ − MH l ˆ v K ll J lv − (cid:15) M T κ − D L vl K vv + (cid:15) M T κ − M p ˆ v p l p v = F ˆ v F l F v . (21)As in the fully-resolved model (7), main diagonal blocksare denoted as K ii and the coupling blocks G ˆ vl and H l ˆ v stem again from the transvascular 1D-3D exchangeterm. Additionally, the coupling blocks J lv and L vl ac-count for exchange between homogenized vasculatureand IF. The terms involving the mortar matrices D , M and κ couple blood flow in the larger vessels with thehomogenized vasculature using our mortar penalty ap-proach. Obviously, the LMs are no longer part of thesystem which is, consequently, not of saddle-point typeanymore. The drawback, however, is that the choiceof the penalty parameter influences the accuracy withwhich the constraint is fulfilled. Large penalty parame-ters yield better accuracy in terms of constraint fulfill-ment but can lead to an ill-conditioning of the systemmatrix. We will comment on the choice of the penaltyparameter in Remark 4. Remark 1
The concrete implementation of the hybridmodel is slightly different than described here for illus-trative purposes. The equations for IF flow and bloodflow are evaluated simultaneously on the 3D domainand not assembled into two separate block matrices aswritten in (21). This means that the degrees of freedomare actually re-ordered in a node-wise manner com-pared to (21) such that one row corresponding to thenodal IF pressure at a node j is followed by a row cor-responding to the homogenized blood pressure at thisnode j . Therefore, we actually solve a system which isblocked with 2 × This section describes the setup of our fully-resolvedand of our hybrid model. We first analyse the real-worldtumor data sets which we will employ for all our numer-ical tests. Subsequently, the assignment of boundaryconditions in both models is described. Then, we illus-trate how we create the hybrid model with homogenized vasculature starting from the full topology of the vas-cular networks. Finally, the definition of representativeelementary volumes for homogenization is introduced.3.1 Analysis of real-world tumor data setsWe have obtained three different vasculature data setsfrom REANIMATE (D’Esposito et al. 2018; Sweeneyet al. 2019), which is a framework combining mathe-matical modeling with high-resolution imaging data topredict transport through tumors. We only briefly de-scribe the experimental procedure here. More detailsare given in the two aforementioned papers. Two dif-ferent colorectal cell lines, namely SW1222, LS174T,and one glioma cell line, GL261, were grown subcuta-neously for 10 to 14 days in mice, resected and opticallycleared, and finally imaged using optical projection to-mography. The data was then segmented to obtain thecomplete blood vessel networks inside the tumors in thegraph format as discussed before.The topologies and blood vessel radii of the threedistinct cases are illustrated in Figure 2 together withrepresentative results of blood vessel and IF pressureof the fully-resolved model. Further network data hasbeen collected in Table 1: All three networks containmore than 100 000 blood vessel segments and nodes.The SW1222 case is the biggest tumor both in net-work size and tissue volume. The latter has been calcu-lated by approximating the hull of the tumor using the alphaShape function of Matlab R2018b (MathWorksInc. (2018)). The hull is then smoothed, remeshed usingGmsh (version 4.4.1, Geuzaine and Remacle 2009) andslightly enlarged to encompass all vasculature nodes.Its enclosed region is integrated to give the tumor vol-ume, see also Figure 3 for the SW1222 tumor. Notethat this tumor domain corresponds to the domain Ω v on which the additional porous network of smaller ves-sels is present in the hybrid model and where its ad-ditional governing equation (9) is defined and solved.Furthermore, all topologies are rotated such that theirprincipal axes align with the coordinate axes. The threedifferent cases show distinct vascular architectures, forinstance, the SW1222 network is much denser with ahigher blood vessel volume fraction and blood vesselsurface-to-volume ratio than the two other types. Inaddition, its blood vessel diameters are generally largerand have a much higher variability. Finally, we have an-alyzed the boundary nodes, i.e., the tips which are onlyconnected to one other node. All topologies have a com-parable number of boundary nodes lying on the afore-mentioned enclosing alpha shape whereas the GL261and SW1222 tumors have a much higher number of tipsinside the domain than the LS174T tumor. alidation and parameter optimization of a hybrid embedded/homogenized solid tumor perfusion model 9(a) SW1222(b) LS174T(c) GL261 Fig. 2: Full topology and structure of the vascular networks (left, colour-coded by the respective radii), represen-tative results for simulated blood pressures (middle) and IF pressures (right) in the fully-resolved model (Samespatial scale is used for all three cases)3.2 Assignment of boundary conditionsThe assignment of physiologically reasonable boundaryconditions on large vascular networks is quite challeng-ing since flows or pressures cannot be measured on thelevel of individual micro-vessels. Sweeney et al. (2019)developed an algorithm (Sweeney et al. 2018) to ap-ply boundary conditions which match in-vivo measure- ments of perfusion for the present data set. We will re-use this framework here to generate the boundary con-ditions for the fully-resolved case and briefly describeit in Section 3.2.1. Boundary conditions for the hybridmodel are detailed in Section 3.2.2.
LS174T GL261 SW1222 UnitNo. of segments (elements) of 1D network 186 092 120 340 419 198 − No. of nodes of 1D network 178 592 110 062 385 218 − Tumor volume 190 . . . Tumor dimensions 4 . × . × .
88 2 . × . × .
32 6 . × . × .
08 mm × mm × mmBlood vessel volume fraction 1 .
13 4 .
01 14 .
90 %Blood vessel surface-to-volume ratio 1 . × − . × − . × − µ m − Mean blood vessel diameter ± std. dev. 22 . ± . . ± . . ± . µ mMean blood vessel segment length ± std. dev. 27 . ± . . ± . . ± . µ mNo. of boundary nodes of 1D network on tumor hull 1559 2419 1933 − No. of boundary nodes of 1D network inside domain 1855 6599 13 772 − No. of elements of 3D domain 15 955 142 15 141 173 13 231 813 − No. of nodes of 3D domain 2 660 273 2 524 666 2 207 655 − Mean element size in Ω v . . . µ mEdge length of REV 1250 750 1500 µ m Table 1: Details on tumor vasculature data sets and discretization
Ω ∂ΩΩ v ∂Ω v Fig. 3: Mesh of three-dimensional domain for SW1222 tumor. Tumor domain (equivalent to the domain Ω v onwhich the additional porous network of smaller vessels is present in the hybrid model) is depicted in red and hasbeen obtained as the alpha shape of nodes of the vascular network. For the fully-resolved model, boundary conditions forthe blood pressure p ˆ v in the 1D network and the IFpressure p l need to be assigned. For the blood ves-sel pressure, we re-use the approach of Sweeney et al.,which has been made publicly available (Sweeney et al.2018) and is based on earlier work by Fry et al. (2012).Thereby, boundary conditions are assigned on the tipsof the network, i.e., on the boundary nodes of the 1Drepresentation of the vasculature both on the tumor hulland inside the tumor as given in Table 1. The followingalgorithm is applied: First, a high or low pressure of5999 . . all end points ofthe 1D network have been assigned such a high/lowpressure. Additionally, the method prevents that points which are in close proximity to each other are assignedthe far apart pressure values which would produce un-physiologically large flows. Second, a no-flux boundarycondition is randomly assigned to the interior bound-ary nodes until 33 % of all boundary nodes have thiscondition. This value is consistent with the fraction ofdead ends in tumor vasculature estimated from exper-imental studies (Morikawa et al. 2002). Third, the re-maining 62 % of boundary nodes are given as unknownsto the flow optimization scheme of Fry et al. (2012).This scheme aims at solving a constrained optimizationproblem for incomplete pressure boundary data by min-imizing the error of pressures and wall shear stress w.r.t.pre-defined target values. D’Esposito et al. (2018) andSweeney et al. (2019) have shown that this procedurefor assignment of boundary conditions ensures that tu-mor perfusion is in good agreement with experimentaldata. Note that the entire algorithm is not determinis- alidation and parameter optimization of a hybrid embedded/homogenized solid tumor perfusion model 11 tic due to the random selection of nodes for high/lowboundary conditions on the external surface of the tu-mor and of nodes for no-flux boundary condition in itsinterior. Therefore, the analyses in the following sec-tions will be performed on five different sets of pressureboundary conditions on the 1D network per tumor case.Concerning the IF pressure, we want to prescribethe far-field pressure for the IF as p l ∞ = 0 Pa followingSweeney et al. (2019). In order to achieve this withinour finite element approach, we enlarge the domain Ω radially to a sphere of radius 80 mm as shown in Fig-ure 3 for the SW1222 case. This allows us to set aDirichlet boundary condition of p l = 0 Pa on its bound-ary ∂Ω and, thereby, to mock the far-field pressure. Wevalidated this approach in the following way for all threevascular networks: We solved the fully-resolved modeland compared the IF pressure solution (for one specificset of pressure boundary conditions on the 1D network)with a case where the domain was only enlarged to asphere with radius 40 mm (with corresponding zero IFpressure Dirichlet boundary condition on its outer sur-face). No visible differences in the IF pressure distri-bution in our domain of interest inside and around thetumor domain could be detected. This indicates thatthe enlargement is big enough insofar as the solution inthe domain of interest is not influenced by the size of theenlargement any more. We can also gradually coarsenthe mesh when moving away further from the vasculardomain as depicted in Figure 3 since the IF pressuregradient flattens and tends to zero further away fromthe center of the domain. This enables the use of a suf-ficiently fine mesh for the region surrounding the em-bedded vascular network while the computational costfor extending the domain is not too high. In addition to boundary conditions for the IF pressure p l and the blood pressure p ˆ v , the hybrid model requiresboundary conditions for the pressure in the homoge-nized vasculature p v . The IF pressure is treated as inthe fully-resolved model and we set it to zero at theboundary of the domain ∂Ω . In the following, we will al-ways compare the accuracy of the hybrid variant w.r.t.the fully-resolved one for one specific set of pressureboundary conditions on the 1D network obtained withthe procedure described in the previous sectoin. Thus,to perform this comparison the pressure boundary con-ditions on the 1D network are transferred from thefully-resolved model to the hybrid model in the follow-ing manner: The boundary conditions of blood pressure p ˆ v on the larger vessels Λ L can directly be taken fromthe boundary conditions of the fully-resolved model. If a node with a Dirichlet boundary condition in thefully-resolved vasculature Λ is part of the larger ves-sels Λ L we simply keep this boundary condition on the1D discretization also in the hybrid model. Dirichletboundary conditions on the smaller vessels Λ S cannotbe assigned on the 1D discretization since smaller ves-sels are homogenized. However, we can employ them toassign boundary conditions for p v on the boundary ofthe domain of homogenized vessels ∂Ω v as depicted inFigures 1b and 3. Similar to Vidotto et al. (2019), wesmooth these values to account for the homogenizationof the smaller vessels: Each condition belonging to anode of the smaller vessels Λ S at the tumor surface isassigned to all 3D nodes lying on the surface ∂Ω v withina distance of less than 400 µ m for the SW1222 and theLS174T tumor and less than 200 µ m for the GL261 tu-mor. Nodes of the 3D mesh which lie within this dis-tance of multiple boundary nodes on Λ S are assignedthe mean pressure value of all these boundary nodes.On the rest of the surface ∂Ω v we set a no-flux bound-ary condition. We also do not set a boundary conditionfor p v on nodes of the 3D mesh in close proximity to endnodes of the 1D network since this would mean settingdifferent boundary conditions on nodes whose pressuresshould be coupled due to the constraint on pressures p ˆ v and p v and, thus, would lead to an overconstrainmentof the system. The resulting distribution of boundaryconditions p v over ∂Ω v is illustrated in Figure 4 forthree exemplary cases.3.3 Distinction between fully-resolved and hybridmodelAs previously stated, we envision that our hybrid modelcould be applied in cases where the full structure ofthe vascular network is unknown such that only thetopology of the larger vessels can be acquired via non-invasive imaging. However, in our data sets we actuallyhave the full structure available. In line with the maingoals of this paper, namely, to validate the hybrid ap-proach, to quantify the error with respect to the fully-resolved case and to determine its optimal parametersfor perfusion through solid tumors, we artificially cre-ate the hybrid model from the fully-resolved one. In thehybrid approaches of Shipley et al. (2019) and Vidottoet al. (2019) this was realized by a radius-based crite-rion. Their employed data sets had a clear hierarchytypical for the microcirculation with larger arteriolesbranching into smaller capillaries which in turn con-nect and form larger venules. Thus, it was possible toexploit the hierarchical structure of the vasculature bykeeping only the larger vessels in the set I L . Fig. 4: Exemplary distributions of boundary conditions for homogenized pressure p v on boundary of vasculardomain ∂Ω v in the hybrid model variant for all three casesFor our tumor vasculature data sets this is not asstraightforward. While there are some thicker vascularbranches, especially in the SW1222 case, no clear hier-archical vascular architecture can be extracted from thetopologies in Figure 2 with a radius-based criterion. Toillustrate this fact, we compare the full architecture ofthe SW1222 network with a network where only the top10 % of vessels with the largest radius are kept in Fig-ures 5a and 5c. Many small unconnected clusters of sev-eral blood vessel segments appear due to the heteroge-neous, extremely variable distribution of the radius andlack of vascular hierarchy. Branches connecting thesecluster which have a smaller radius are removed. Ap-plying our or any hybrid model on this topology wouldnot be possible as hybrid approaches also rely on a ”sen-sible” topology for Λ L which preserves the structure ofthe entire network via one or several connected sub-graphs of larger vessels which feed respectively drainthe smaller, homogenized vessels. Only then, the 1Dblood flow model and corresponding boundary condi-tions can reasonably be applied on Λ L together withsuitable exchange terms into the smaller vessels. Thus,we instead distinguish between smaller and larger ves-sels based on the flow within the vessels. This yieldsa better preservation of the network architecture forthe hybrid case, see Figures 5a and 5b. Now, connectedsubgraphs of larger vessels Λ L emerge which connectin- and outlets of the main flow-carrying vessels withthe smaller vessels.Hence, our strategy to obtain Λ L is as follows: Wefirst solve the fully-resolved model (using the bound-ary conditions described in Section 3.2.1). Then, all ele-ments except the ones with the highest flow are deletedfrom the vascular graph, e.g., the top 10 % with the highest flow are kept. However, there are still some verysmall clusters consisting of only a few segments presentin the graph. We additionally delete connected compo-nents from the graph whose overall length is smallerthan 250 µ m, i.e., sub-components which are smallerthan ten segments with the average segment lengthsgiven in Table 1. By that, we delete an additional 0 . − . I L of largervascular branches which are kept in the hybrid ap-proach as exemplarily featured in Figure 5b. Here, weshow only the SW1222 case but equivalent results holdfor the other two network topologies. In the following,we will denote cases where the top X % of elements withhighest flow are kept and the small connected compo-nents are removed according to the procedure describedabove as ”case X %”.Recall that the assignment of pressure boundaryconditions on the fully-resolved vascular network is notdeterministic. Moreover, different boundary conditionswill produce distinct flow patterns in the vasculatureand, hence, also different sets of large and small vesselsin our procedure and a different topology for Λ L . There-fore, the following analysis will always be performed forfive sets of pressure boundary conditions on the 1D net-work with corresponding distinct sets of large and smallvessels I L and I S .In Table 2 the mean diameters D Λ L and D Λ S oflarger and smaller vessels are compared. It is obviousthat the diameters in the set of small vessels I S whichare removed from the hybrid model are considerablysmaller than the diameters of the large vessels. This be-haviour is most pronounced for the SW1222 topologywhere for the case 5 % the mean diameters in Λ L are 2 . alidation and parameter optimization of a hybrid embedded/homogenized solid tumor perfusion model 13(a) full topology (b) 10 % with highest flow (c) 10 % with biggest radius Fig. 5: Extraction of large vessels Λ L from the entire network Λ – Comparison between flow-based criterion (andsorting out of small connected components) and radius-based criteriontimes bigger than in Λ S . Naturally, this ratio drops forall topologies when a higher percentage of segments iskept in the large vessel set. For the LS174T and GL261data sets the difference in blood vessel diameters is notas large but this can be attributed to the fact that thediameters are less dispersed than in the SW1222 topol-ogy, see also the mean and standard deviation of thediameters in Table 1. Also in these cases, the diametersin Λ L are larger by approximately one standard devi-ation of the diameter of the entire vasculature (as inthe SW1222 case). In summary, our approach incorpo-rates mainly the vessels with larger radii in the set I L whereas also some segments with smaller radii are keptto preserve the main topology of the networks in thehybrid model. Therefore, there is also a significant con-gruence of the sets of large vessels I L between differentpressure boundary condition cases. For instance, in thecase where 10 % of the blood vessels are kept in the hy-brid model, the average percentage of identical retainedsegments between two different pressure boundary con-dition cases is 45 % for the LS174T tumor, 51 % forthe GL261 tumor and 78 % for the SW1222 tumor. InRemark 3 we further comment on how the obtainedtopologies of larger vessels Λ L relate to real in-vivo tu-mor imaging data.Next, we justify our line-based coupling approachbetween the large vessels Λ L and the homogenized vas-culature. For that, we have analyzed the connectiv-ity between larger and smaller vessels for the fully-resolved topologies in Table 3. Here, we denote by ϕ = n nodes ,Λ L ∩ Λ S /n nodes ,Λ L the proportion of nodes of thelarger vessels Λ L which have a direct connection to anode of the smaller vessels Λ S . The presented data il-lustrates that for the GL261 and the SW1222 tumor almost every third to every fourth node of the mainbranches Λ L is directly connected to a node of thesmaller blood vessel segments Λ S , i.e., at every thirdto fourth node a smaller vessel branches away from Λ L .For the LS174T network, the connectivity is slightlysmaller. Here, only 13 −
18 % of nodes in larger ves-sels are connected to smaller vessels. In all cases, thesenumbers obviously again drop when keeping a largerportion of the entire network in the set I L .In the hybrid approach, information about thesesmaller branching vessels is lost since they are removedfrom the 1D representation of the vasculature. As statedabove, we want to enforce equal pressures between largerand smaller vessels as this equality also holds in thefully-resolved model at branching points. The high con-nectivity between the two network parts supports ourline-based mortar penalty coupling between the resolvedand homogenized part of the vasculature in which weactually couple the entire network of big vessels Λ L with the homogenized vasculature. Of course, we knowthe connecting nodes between larger and smaller ves-sels here as we know the full topology of all networks,so we could also enforce the coupling between resolvedand homogenized part in a point-based manner at theselocations. However, in the more realistic case when onlythe architecture of larger vessels is known without theexact locations where smaller vessels branch away, thisis not the case. Therefore, we adopt our line-based cou-pling within the hybrid model hereafter to compare theresults with the fully-resolved reference solution. Notethat the network tips of Λ L (both in the interior of thedomain and on the tumor hull) are actually also cou-pled with the homogenized vasculature since the dis-crete constraint of a vanishing weighted pressure gap is µ m] D Λ L D Λ S D Λ L D Λ S D Λ L D Λ S D Λ L D Λ S SW1222 104 . . . . . . . . . . . . . . . . . . . . . . . . Table 2: Comparison of mean blood vessel radius in larger vessels Λ L and smaller vessels Λ S (all values indicate themean taken over five different sets of pressure boundary conditions on the 1D network produced by the methodologydescribed in Section 3.2.1, ”case X %” denotes the case where X % of the 1D blood vessels are retained in thehybrid approach) case 5 % case 10 % case 15 % case 20 % ϕ CV D CV | Q | ϕ CV D CV | Q | ϕ CV D CV | Q | ϕ CV D CV | Q | SW1222 0 .
29 0 .
39 2 .
40 0 .
28 0 .
45 2 .
14 0 .
26 0 .
50 1 .
77 0 .
24 0 .
55 1 . .
18 0 .
23 0 .
88 0 .
16 0 .
24 0 .
85 0 .
15 0 .
24 0 .
83 0 .
13 0 .
24 0 . .
30 0 .
35 1 .
21 0 .
29 0 .
37 1 .
16 0 .
28 0 .
38 1 .
13 0 .
27 0 .
40 1 . Table 3: Analysis of connectivity between fully-resolved and homogenized part of vasculature: ϕ is the fraction ofnodes of larger vessels with a direct connection to smaller vessels, CV D and CV | Q | are measures of the variability ofthe diameter and flow, respectively, in the segments connecting larger and smaller vessels (data includes the meantaken over five different sets of pressure boundary conditions on the 1D network produced by the methodologydescribed in Section 3.2.1, ”case X %” denotes the case where X % of the 1D blood vessels are retained in thehybrid approach)enforced along the entire 1D discretization and, thus,also at the end nodes.Finally, we analyze also the elements connectinglarger and smaller vessels, i.e., those 1D elements ofthe smaller vessels where one node is part of Λ L andthe other part of Λ S . We gather all these elements andcompute their mean diameter and mean absolute flowvalue. Then, we calculate the coefficient of variation ofthese quantities, CV D and CV | Q | as the ratio of stan-dard deviation of the diameter, resp. flow to its mean inthese connectivity elements. The results are collected inTable 3. Obviously, the SW1222 case shows the highestvariability in both flow and diameter followed by theGL261 and the LS174T case. For all cases, the variabil-ity of the flow is larger than for the diameter since thevolumetric flow in a segment depends on the fourthpower of the diameter due to the employed Hagen-Poiseuille relationship. These results are consistent withthe topology of the entire network where the variabil-ity of the blood vessel diameter is also larger for theSW1222 tumor than the GL261 and the LS174T tu-mor, see Table 1. In Section 4.3 we will show that thishigher variability makes it harder to match the flowfrom large into small vessels between the two models. Remark 2
We believe that our hybrid approach is alsoapplicable to more organized, hierarchical networks as,for example, the topology used by Vidotto et al. (2019).In this publication the network was partitioned by aradius-based threshold, see Figure 1 therein. The larger vascular structures contain very short branches goingaway from the main vessels. At the tips of these shortbranches, the node-based coupling is performed. If oneinstead removed these very short branches and left onlythe major, flow-carrying vessels in Λ L , a line-based cou-pling along these vessels could again be implemented.3.4 Determination of representative elementaryvolume sizeThe existence of a representative elementary volume(REV) is an important concept for different homoge-nization procedures (Davit et al. 2013). In general, sucha volume should be big enough to smooth out fluctua-tions of spatial heterogeneities yet small enough to re-solve the physical effects of interest. In this section, weinvestigate the choice of REVs in the context of ourmodel and the employed data sets. Naturally, we willinvestigate the properties of the smaller vessels Λ S inthe following since this is the part of the vasculaturewhich is homogenized and treated as a porous contin-uum in the hybrid approach. Furthermore, five differentsets of pressure boundary conditions on the 1D networkare studied. This is again due to the fact that differentpressure boundary conditions on the 1D network willlead to different flow patterns in the vascular networkand, therefore, also different sets I L and I S of large andsmall vessels (potentially with a different distributionthroughout the domain) with the employed flow-basedcriterion. alidation and parameter optimization of a hybrid embedded/homogenized solid tumor perfusion model 150 1 ,
000 2 ,
000 3 ,
000 4 , . . . l REV l [ µ m] ε v Λ S [ − ] ,
000 2 ,
000 3 ,
000 4 , . . · − l REV l [ µ m] ( S / V ) Λ S [ µ m − ] (a) SW12220 1 ,
000 2 ,
000 3 ,
000 4 , · − l REV l [ µ m] ε v Λ S [ − ] ,
000 2 ,
000 3 ,
000 4 , · − l REV l [ µ m] ( S / V ) Λ S [ µ m − ] (b) LS174T0 500 1 ,
000 1 ,
500 2 , · − . l REV l [ µ m] ε v Λ S [ − ] ,
000 1 ,
500 2 , . . · − l REV l [ µ m] ( S / V ) Λ S [ µ m − ] (c) GL261 Fig. 6: Determination of representative elementary volume (REV) size – evolution of blood vessel volume fraction ε vΛ S and surface-to-volume ratio ( S/V ) Λ S of smaller blood vessels Λ S is shown for increasing possible REV sizesFor this purpose, we have devised the following pro-cedure:1. For each network topology we create five differentcases with a different set of pressure boundary con-ditions on the 1D network for the fully-resolved modelas described in Section 3.2.1.2. We partition all cases into the distinct sets of largeand small vessels as described in Section 3.3. Wehere investigate the case 10 % for all different topolo-gies but equivalent results have been obtained forleaving the top 5 %, 15 % or 20 % of vessels with thelargest flow in the system.3. We select random positions in the vasculature do-main Ω v in the range [ x min +0 . · l x , x max − . · l x ],[ y min + 0 . · l y , y max − . · l y ] and [ z min + 0 . · l z , z max − . · l z ], where l i denotes the domain lengths in the respective coordinate directions and( · ) min and ( · ) max the minimum and maximum co-ordinate value in each direction in Ω v . In this way,the random positions are chosen such that they donot lie too close to the boundaries of the domain.4. For each of the random positions within the domainwe define a cube with edge length l edge = 1 / · max( l x , max( l y , l z )). The random position is chosenas the center of that cube.5. The size of the cubes is successively increased in allcoordinate directions by l edge while keeping theircenters fixed. The blood vessel volume fraction ε vΛ S and the surface-to-volume ratio ( S/V ) Λ S of smallerblood vessels Λ S is computed for each cube at eachsize. If a cube protrudes from the domain duringthis enlargement, these quantities are calculated onthe intersection of the cube with the domain Ω v . Per case with different boundary conditions this is per-formed for ten randomly generated cube centers. Theresults are shown in Figure 6 for only three REVs perboundary condition case, that is, in total 15 cases tonot clutter the plots. The evolution of the blood vesselvolume fraction ε vΛ S and of the surface-to-volume ratio( S/V ) Λ S of the smaller blood vessels Λ S for increasingthe edge length of the cubes is illustrated. Therein, wedenote the length scale as l = (cid:112) V cube ∩ Ω v to accountfor cases when a larger cube protrudes from the do-main Ω v . All three topologies exhibit similar features: ε vΛ S and ( S/V ) Λ S fluctuate strongly for smaller lengths.Then, most curves stabilize and remain almost station-ary while increasing the size of the averaging volume.Finally, for even larger volumes the curves slowly con-verge to the values of these quantities across the entiredomain. This behaviour can be expected in porous me-dia (Davit et al. 2013) and we consequently define thelength of the REVs l REV at the point where the initialoscillations of too small control volumes fade out andthe values stabilize.Splitting the domains into these REVs of equal sizeis not an easy task due to their irregular, elliptic shape.We first defined a regular grid of REV centers andperformed an initial Voronoi tesselation based on thisgrid. Due to the shape of the domain this resultedin too small or too large REVs. Therefore, we per-formed an optimization of the Voronoi tesselation wherethe objective function had the goal to define REVsof equal volume and equal dimensions. The resultingREVs are visualized in Figure 7. The mean deviationof the REVs from the previously determined volumeand lengths from Figure 6 in all three coordinate direc-tions is less than 5 % for the domains of all three tumortypes.Finally, we employ these REVs to study the distri-bution of blood vessels inside the domain. For that, wedefine the non-dimensionalized radial distance of eachREV ˜ r REV as the distance of the center of the REVto the center of the domain divided by the distance ofthe center of the domain to the tumor hull in directionof the center of the REV. Again, this analysis is per-formed for all three tumor types for five different sets ofpressure boundary conditions on the 1D network sincethose influence the flow in the 1D vasculature and, con-sequently, also the composition of Λ L and Λ S as previ-ously mentioned. The results for the volume fraction ofbig vessels ε vΛ L , small vessels ε vΛ S and the entire vascu-lature ε vΛ are shown in Figure 8. The clearest structureis evident for the the SW1222 case: towards the tu-mor hull, ε vΛ S and ε vΛ L and, thus, also the sum of theformer two, ε vΛ , gradually increase. Close to the cen-ter of the domain, there is still a significant amount of smaller blood vessels while almost no larger blood ves-sels are present. This is consistent with experimentaldata showing higher blood vessel density and perfusionin the tumor periphery (D’Esposito et al. 2018; Forsteret al. 2017) with only a few major vessels penetratinginto the center of the tumor (Holash et al. 1999). Thesetrends are also present in the LS174T tumor, albeit, farless pronounced than for the SW1222 tumor. By con-trast, the GL261 vascular network shows a completelydifferent behaviour. While the vascular density of largevessels remains almost constant over the tumor radius,the one of the smaller blood vessels Λ S drops and, thus,also the overall volume fraction ε vΛ . Remark 3
The validity of the obtained topologies anddistributions for Λ L and Λ S and the applicability ofthe proposed hybrid approach is supported by state-of-the art optoacoustic in-vivo imaging techniques (Liet al. 2020). The currently attainable spatial resolutionis less than 50 µ m throughout the tumor domain whichis in the range of the diameter of larger vessels fromTable 2. Furthermore, the larger vessels which are re-tained in the hybrid model are more concentrated at thetumor periphery (at least for the SW1222 and LS174Tcase) and are, thus, more accessible to imaging. Qualita-tively, the topology of the larger vessels from Figure 5bis in good agreement with corresponding imaging datafrom tumors (Li et al. 2020) where larger feeding vesselsare visible at the tumor rim. From these experiments,one can extract a similar topology of larger vessels Λ L to apply our hybrid model. Hence, we conclude thatthe employed methodology of splitting into larger andsmaller vessels yields a valid scenario resembling realexperimental data and can, therefore, be used to inves-tigate our hybrid embedded/homogenized approach forsolid tumor perfusion. In this section we perform several numerical experi-ments to evaluate the performance of the hybrid modelin comparison to the fully-resolved one. We first definea comparison metric and optimize the parameters ofthe hybrid model such that the best possible correspon-dence between the models is achieved according to thismetric. Subsequently, we study several other quantitiesto compare the two models and present a further im-provement of the hybrid model via a vascular volumefraction dependent permeability for the homogenizedvessels.The tumor hull is smoothed and triangulated us-ing Gmsh (version 4.4.1, Geuzaine and Remacle 2009)and its enclosed volume is meshed with linear tetrahe- alidation and parameter optimization of a hybrid embedded/homogenized solid tumor perfusion model 17(a) SW1222, n REV = 78 (b) LS174T, n REV = 114 (c) GL261, n REV = 55
Fig. 7: Representative elementary volumes of all three tumor domains . . . . · − . . . .
25 ˜ r REV [ − ] ε v [ − ] ΛΛ L Λ S (a) SW1222 0 . . . . . . . · − ˜ r REV [ − ] ε v [ − ] ΛΛ L Λ S (b) LS174T 0 . . . . · − ˜ r REV [ − ] ε v [ − ] ΛΛ L Λ S (c) GL261 Fig. 8: Dependency of volume fraction of big vessels Λ L , small vessels Λ S and entire vasculature Λ over non-dimensionalized radial distance from center of domain (data is taken from five different sets of 1D blood pressureboundary conditions if 10 % of 1D blood vessels are retained in hybrid model for each network structure, dashedlines indicate linear least squares fits)dral elements using Trelis (2020). An exemplary meshfor the SW1222 topology is shown in Figure 3 and pa-rameters of the 3D mesh are given in Table 1. Notethat the 3D mesh is completely independent of the dis-cretization of the 1D networks, that is, the nodes ofthe two meshes do not match which is an advantageousfeature provided by our recently introduced hybrid ap-proach (Kremheller et al. 2019). Both the fully-resolvedand the hybrid FEM model have been implemented inthe in-house parallel multi-physics research code BACI (2021). Parameters for both models are listed in Ta-ble 4. Remark 4
In preliminary simulations, we determinedthe proper range for the penalty parameter (cid:15) . As acompromise between accuracy and a well-conditionedsystem matrix, we defined the following criterion: δ = 1 n nodes ,Λ L n nodes ,Λ L (cid:88) i =1 (cid:12)(cid:12) κ − [ i, i ] g [ i ] (cid:12)(cid:12) p ˆ v [ i ] < . (22)This rule states that the mean relative pressure error δ in terms of the length-independent nodal pressure dif- ference vector κ − g and the nodal pressure p ˆ v in the1D network is less than 1 %. The values for all cases aregiven in Table 4. Note that the penalty parameter hasunits of [length] / [time · pressure] such that the LM fieldrepresents a 1D-3D mass transfer term, or volumetricflow per length. This allows interpreting the penalty pa-rameter as very large permeability governing the masstransfer between resolved and homogenized vasculaturein the hybrid model. Remark 5
In our opinion, the main advantage of thehybrid model is not a reduction of computational costcompared to the full model, but the fact that it reliesonly on data available through non-invasive imaging.Nevertheless, we also did a first preliminary evaluationcomparing the computational costs of the two modelsand found that the hybrid model was not significantlyfaster than the fully-resolved one and in some cases evenslower. The effort for finding 1D-3D elements interact-ing with each other, building the integration segmentsand evaluating the coupling terms along the 1D vascu-lature is obviously smaller for the hybrid model sinceless 1D vessels are present. However, this is balancedor even outweighed by its increased effort in severalother aspects: The evaluation of the 3D elements ismore costly since two equations per node (in Ω v ) haveto be evaluated, the system size, which is dominatedby the number of 3D nodes, and, thus, the linear solvertime is increased and the condition of the system isworse compared to the fully-resolved case due to thepenalty approach, which in turn raises the linear solvertime. However, for all our studies we used the same 3Dmeshes for both hybrid and full model. The cost for thehybrid model could be greatly reduced by employing acoarser 3D mesh. Vidotto et al. (2019) showed that thisstill gave acceptable results in terms of REV pressuresfor their approach.4.1 Definition of metric for comparison of the twomodelsTo assess the performance of the hybrid model in pre-dicting microvascular flow and IF pressure inside solidtumors in comparison with the fully-resolved model, asuitable metric is warranted. Ideally, the hybrid modelshould match the fully-resolved one in terms of bloodand IF pressure as well as blood and IF flow to obtainan accurate representation of the perfusion through thetumor. Therefore, we define our metric as a combina-tion of these quantities. The first contribution is thecorrespondence of blood pressures in the large vessels Λ L which are present both in the fully-resolved and the hybrid model. We define the coefficient of determina-tion R in terms of nodal blood pressures in the largevessels between the two models as R = 1 − (cid:80) n nodes ,Λ L i =1 (cid:16) p ˆ v [ i ] (cid:12)(cid:12) full − p ˆ v [ i ] (cid:12)(cid:12) hyb (cid:17) (cid:80) n nodes ,Λ L i =1 (cid:0) p ˆ v [ i ] (cid:12)(cid:12) full − µ I L (cid:0) p ˆ v (cid:12)(cid:12) full (cid:1)(cid:1) , (23)where µ I L (cid:0) p ˆ v (cid:12)(cid:12) full (cid:1) is the mean blood pressure in thelarge vessels of the fully-resolved model. A value of R = 1 would mean a perfect correspondence of bothmodels while smaller values suggest larger deviations.A negative R indicates that the hybrid model performsworse than simply taking the mean value of the fully-resolved model. The second contribution to our metricis the correspondence of blood pressures in the smallvessels Λ S between the fully-resolved and the hybridmodel, which we calculate as R = 1 − (cid:80) n nodes ,Λ S i =1 (cid:16) p ˆ v [ i ] (cid:12)(cid:12) full − p v ( X [ i ]) (cid:12)(cid:12) hyb (cid:17) (cid:80) n nodes ,Λ S i =1 (cid:0) p ˆ v [ i ] (cid:12)(cid:12) full − µ I S (cid:0) p ˆ v (cid:12)(cid:12) full (cid:1)(cid:1) . (24)Since the smaller vessels Λ S are not retained in the hy-brid model, we compare nodal blood pressures in thesmaller vessels of the fully-resolved model with the ho-mogenized blood pressure field p v of the hybrid modelevaluated at the nodal positions X [ i ] of the smaller ves-sels. Again, this is formulated in terms of a coefficientof determination, now involving all nodes in the smallvessels and µ I S (cid:0) p ˆ v (cid:12)(cid:12) full (cid:1) is the mean blood pressure inthe small vessels of the fully-resolved model.Equivalently, the coefficient of determination of theIF pressure is given by R = 1 − (cid:80) n nodes ,Ω i =1 (cid:16) p l [ i ] (cid:12)(cid:12) full − p l [ i ] (cid:12)(cid:12) hyb (cid:17) (cid:80) n nodes ,Ω i =1 (cid:0) p l [ i ] (cid:12)(cid:12) full − µ (cid:0) p l (cid:12)(cid:12) full (cid:1)(cid:1) (25)with the mean IF pressure µ (cid:0) p l (cid:12)(cid:12) full (cid:1) of the full model inthe tissue domain Ω . Instead of the point-wise compar-ison of pressures in (24) and (25), one could also com-pare mean REV pressures of the two models. We willadditionally calculate and compare mean (blood andIF) pressures inside the REVs in Section 4.3. With theprevious three equations, the metrics for blood and IFpressure have been defined. Also the flow in the largervessels Λ L is covered since larger vessels present in bothmodels have the same diameter, length and blood vis-cosity. Therefore, if the nodal pressures match, also theflow between the nodes, i.e., inside the elements is iden-tical. The same applies for flow in the IF if the same3D mesh and hydraulic conductivity k l /µ l is employedin both models which we will assume hereafter. What alidation and parameter optimization of a hybrid embedded/homogenized solid tumor perfusion model 19 Quantity Symbol Value Unit Reference Eqns.Density of blood ρ ˆ v , ρ v − Formaggia et al. (2009) (2),(8),(9)Viscosity of blood µ ˆ v [a] Pa s Pries and Secomb (2005) (2),(8)Density of interstitial fluid and blood plasma ρ l − known (3),(4),(11),(10)Hydraulic conductivity of interstitial fluid k l /µ l . × − µ m Pa − s − Boucher et al. (1998) (4),(11)Hydraulic conductivity for transvascular flow L p, ˆ v , L p,v . × − µ m Pa − s − Baxter and Jain (1989) (3),(10)Oncotic reflection coefficient σ .
82 - Sweeney et al. (2019) (3),(10)Oncotic pressure of blood π b . π l . k v /µ v see Table 5 µ m Pa − s − - (10)Surface-to-volume ratio for transvascular flow ( S/V ) Λ S see Table 5 µ m − - (10)Penalty parameter (cid:15) µ m Pa − s − Remark 4 (19)SW1222: case 5 % 400SW1222: case 10 %, 15 %, 20 % 100LS174T: case 5 % 100LS174T: case 10 %, 15 %, 20 % 50GL261: case 5 % 100GL261: case 10 %, 15 %, 20 % 50 [a]
The value for blood viscosity is calculated separately in each 1D element using the algebraic relationship of Pries and Secomb (2005) withhematocrit value fixed to 0 . Table 4: Parameters and values is still missing, is a metric for comparison of blood flowinside the smaller blood vessels Λ S which are homoge-nized in the hybrid model. We define this measure asfollows R ,Λ S = 1 − (cid:80) n REV i =1 (cid:80) j =1 (cid:18) Q ˆ vj (cid:12)(cid:12) Λ S , full − Q vj (cid:12)(cid:12) Λ S , hyb (cid:19) (cid:80) n REV i =1 (cid:80) j =1 (cid:18) Q ˆ vj (cid:12)(cid:12) Λ S , full − µ (cid:18) Q ˆ v (cid:12)(cid:12) Λ S , full (cid:19)(cid:19) , (26)i.e., for each of the n REV
REVs we compare the vol-umetric flow Q j of the fully-resolved and the hybridmodel in all three coordinate directions and compare itwith each other via the coefficient of determination offlow in the smaller vessels R ,Λ S .Next, we will detail how we calculate the flows inthe REVs in both models. In the center of each REVwe define a square (cid:3) j with dimensions l REV × l REV suchthat its normal n j is aligned with coordinate direction j . The volumetric flow in the homogenized part of thevasculature in coordinate direction j is then given by Q vj (cid:12)(cid:12) Λ S , hyb = (cid:90) (cid:3) j − k v µ v n j · ∇ p v dA (27)as the surface integral of the flux through the square.For the fully-resolved model, we define it as Q ˆ vj (cid:12)(cid:12) Λ S , full = (cid:88) (cid:3) j ∩ Λ S − πR µ ˆ v ∂p ˆ v ∂s · sgn ( t · n j ) , (28)which is the sum of the volumetric flow of all segmentswhich are part of the smaller vessels and cut by thesquare (cid:3) j . Therein, t is the tangential vector of a seg-ment pointing from its first to its second node andsgn ( · ) denotes the sign function. Finally, we define the total coefficient of determina-tion between the two models as the sum of the contri-butions from blood pressure in large vessels (23), bloodpressure in small vessels (24), IF pressure (25) and flowin small vessels (26) as R = 14 (cid:0) R + R + R + R ,Λ S (cid:1) . (29)This metric, where all four contributions are weightedequally, will be employed to study the accuracy of thehybrid model w.r.t. the full model and to find the op-timal parameters of the hybrid model.4.2 Optimization of parameters of the hybrid modelCompared to the fully-resolved model, the hybrid onehas two additional parameters, which are the hydraulicconductivity of the homogenized vessels k v /µ v in (9)governing blood flow and the surface-to-volume ratio( S/V ) Λ S accounting for transvascular flow from the ho-mogenized vessels into the IF in (10). Our goal in thissection is to determine these parameters such that theagreement in terms of blood flow and blood and IF pres-sures of the hybrid model with the fully-resolved modelis maximized. For that we employ the total coefficientof variation (29) between the two models deduced in theprevious section and formulate the following optimiza-tion problem in terms of the parameters of the hybridmodel:argmax k v /µ v , ( S/V ) Λ S R , (30)that is, we aim to find the parameters of the hybridmodel, for which the correspondence of the two mod-els is optimized. With these optimal parameters we can then evaluate the accuracy of the hybrid model w.r.t.the fully-resolved one. For the optimization procedure,we parallelized the least-squares method of the SciPypackageSciPy (2020) (v1.5.2) and interfaced it to thesoftware framework QUEENS (Biehler et al. 2021). In-ternally, SciPy employs the Levenberg-Marquardt algo-rithm to solve the nonlinear least-squares problem (30).Derivatives of the metric (29) w.r.t. the parameters areapproximated using forward finite differences. This im-plies that the hybrid model has to be solved three timesper iteration step. In preliminary simulations, we con-firmed that different initial conditions (from a sensibleparameter range) converged to the same optimum.Since the full topology of the vasculature is avail-able, we could also obtain these parameters by a suit-able homogenization procedure for the permeability aspreviously done for other hybrid models (Shipley et al.2019; Vidotto et al. 2019). We did not follow this ap-proach here for the following reasons: First, the chaoticstructure of the blood vessel network implying also avery chaotic blood flow pattern typical for the solidtumors would make this very challenging. Second, wewant to create a best-case scenario by fitting the pa-rameters of the hybrid model such that possible errorsintroduced by a homogenization scheme are minimal.The general algorithm can be described as follows:1. Obtain a set of boundary conditions for the fullmodel as described in Section 3.2.1 and solve thefull model to generate a reference solution.2. Extract the topology of larger vessels for the hybridmodel from the full model, conf. Section 3.3, andapply boundary conditions on the hybrid model,conf. Section 3.2.2.3. Find the optimal parameters of the hybrid modelby maximizing the total coefficient of variation (30).During the optimization procedure repeated evalua-tions of the hybrid model with different parametersare required.Representative results of the optimization schemeare depicted in Figure 9 for all four contributions to thetotal coefficient of determination. Very good agreementbetween the two models in terms of nodal pressuresin the larger vessels p ˆ v can be observed in Figure 9a.This can be expected because the same boundary con-ditions on the large vessels are applied in both cases.Thus, large and small pressure values show very goodagreement, further away from these boundary condi-tions in the medium pressure range, deviations becomelarger. The clusters with the largest errors are separatebranches which are not directly connected to nodes ofthe 1D network carrying boundary conditions. The cor-respondence for the nodal IF pressures p l in Figure 9cis also very good. For low IF pressures this is again due to the zero pressure boundary condition assignedon ∂Ω for both cases, but also for higher IF pressuresinside the tumor, which is the actual domain of interest,the pressure differences are very small, in this case, themaximum absolute error is 237 . . p ˆ v (cid:12)(cid:12) full ≈ −
20 % of the largervessels are kept in the hybrid model. We found that fivesets of pressure boundary conditions on the 1D networkwere enough to study our hybrid model since randomlypicking only four out of the five boundary conditioncases changed the mean result by at most 8 %. More-over, taking the mean parameter of a case X % overall different boundary condition cases instead of theoptimal value for each specific case only changed thetotal coefficient of determination by less than 2 %. Fur-thermore, we compare the result of the optimizationprocedure for (
S/V ) Λ S with the calculated surface-to-volume ratio of the smaller vessels for each case. The alidation and parameter optimization of a hybrid embedded/homogenized solid tumor perfusion model 212 ,
000 3 ,
000 4 ,
000 5 ,
000 6 , , , , , , p ˆ v (cid:12)(cid:12) full [Pa] p ˆ v (cid:12)(cid:12) h y b [ P a ] (a) R = 0 .
983 2 ,
000 3 ,
000 4 ,
000 5 ,
000 6 , , , , , , p ˆ v (cid:12)(cid:12) full [Pa] p v ( X ) (cid:12)(cid:12) h y b [ P a ] (b) R = 0 . ,
000 2 ,
000 3 , , , , p l (cid:12)(cid:12) full [Pa] p l (cid:12)(cid:12) h y b [ P a ] (c) R = 0 . − − · − − · Q ˆ v (cid:12)(cid:12) Λ S , full [ µ m s − ] Q v (cid:12)(cid:12) Λ S , h y b [ µ m s − ] (d) R ,Λ S = 0 . Fig. 9: Exemplary comparison of hybrid model (with optimized parameters) with fully-resolved model for onespecific network topology (SW1222, 10 % of 1D blood vessels have been retained in hybrid model). In each sub-figure, solution of hybrid model is plotted over solution of fully-resolved model and the dashed line indicatesperfect agreement between the models with 1:1-correspondence. Comparison of blood pressure in large vessels isdepicted in subfigure a), blood pressure in small vessels in b), IF pressure in c) and flow in small vessels in d).Coefficient of determination for agreement between both variants is given for each quantity and overall coefficientof determination calculated according to (29) is R = 0 .
716 for this case.relative error E S/V is smaller than 5 % for all cases val-idating that the optimization procedure converges to aphysically reasonable result. The permeability is largestfor the SW1222 topology which can be expected con-sidering the much denser network of this case. For alltumors, it decreases if a larger proportion of the 1Dvessels is kept in the model, which is also sensible sincethe smaller the proportion of homogenized vessels, theless permeable these vessels.As already described above, all cases exhibit a verygood correspondence in terms of blood pressures inlarger vessels and IF pressures, proven by the valuesfor R and R in Table 5. If the fidelity of the hybridmodel is increased by resolving a larger proportion of the network structure, the agreement between the twomodels grows likewise. This is also the case for the coef-ficient of determination of blood pressure in smaller ves-sels R . Here, the SW1222 and the LS174T case exhibitcomparable accuracy whereas the GL261 case experi-ences a larger discrepancy. We can attribute this to thefact that this topology has the largest number of tips atthe tumor hull and also the largest number of dead endsconsidering that it is the smallest data set, see Table 1.Hence, the pressure error is largest due to non-matchingboundary conditions between fully-resolved and hybridmodel as mentioned above. However, we will show inSection 4.3 that in terms of REV blood pressures itsconformity with the hybrid model is as good as the Network Case k v /µ v [ µ m ] ( S/V ) Λ S [ µ m − ] E S/V [%] R R R R ,Λ S R SW1222 case 5 % 16 . ± .
840 (6 . ± . × − .
17 0.944 0.488 0.994 0.163 0.647case 10 % 3 . ± .
506 (5 . ± . × − .
93 0.988 0.654 0.998 0.176 0.704case 15 % 1 . ± .
247 (5 . ± . × − .
66 0.993 0.739 0.999 0.262 0.748case 20 % 0 . ± .
122 (4 . ± . × − .
84 0.989 0.792 0.999 0.193 0.743LS174T case 5 % 1 . ± .
105 (1 . ± . × − .
21 0.905 0.643 0.990 0.282 0.705case 10 % 1 . ± .
096 (1 . ± . × − .
18 0.916 0.683 0.991 0.260 0.713case 15 % 0 . ± .
032 (1 . ± . × − .
70 0.930 0.695 0.992 0.244 0.715case 20 % 0 . ± .
064 (1 . ± . × − .
30 0.944 0.718 0.993 0.207 0.715GL261 case 5 % 1 . ± .
288 (6 . ± . × − .
35 0.917 0.195 0.997 0.199 0.577case 10 % 0 . ± .
093 (5 . ± . × − .
63 0.927 0.233 0.996 0.107 0.566case 15 % 0 . ± .
073 (5 . ± . × − .
12 0.941 0.295 0.996 0.113 0.586case 20 % 0 . ± .
061 (4 . ± . × − .
21 0.950 0.346 0.996 0.134 0.607
Table 5: Results of the optimization procedure for hydraulic conductivity and surface-to-volume ratio of homog-enized vasculature in the hybrid model. Relative error w.r.t. calculated surface-to-volume ratio and R -values foragreement between both variants in terms of blood pressure in large vessels, blood pressure in small vessels, IFpressure and flow in small vessels is additionally provided. Overall coefficient of determination R between fully-resolved and hybrid model is calculated according to (29). (All data includes the mean taken over five differentsets of pressure boundary conditions on the 1D network produced by the methodology described in Section 3.2.1,”case X %” denotes the case where X % of the 1D blood vessels are retained in the hybrid approach)other cases. The difference in flow in the small vesselsis the largest source of error in all cases. Also takingmore 1D vessels into account for the hybrid model doesnot necessarily improve the behaviour. We believe thatthis is due to the chaotic flow patterns in the smallervessels and to the fact that we define the permeabilitytensor as isotropic and constant over the entire domain Ω v . Apparently, this is insufficient to resolve the flowin the homogenized vasculature in comparison to thefull model. We tried to increase the agreement by giv-ing a higher weight to the coefficient of determinationof flow in the smaller vessels R ,Λ S in the definitionof the total coefficient of determination (29) but couldnot achieve any significant improvements. However, theagreement in terms of flow in the entire (resolved andhomogenized) vasculature, which will be investigated inthe next section, is much better.4.3 Additional comparisons of results of both modelsWe further study the agreement of the hybrid modelwith the optimized parameters from the previous sec-tion in terms of several other quantities. For that, wedefine the mean REV IF pressure in the j -th REV as p l ( j ) = 1 V REV j (cid:90) REV j p l dV. (31)This is employed to study the absolute and relativemean IF pressure error between the two models in eachREV as E l abs ( j ) = abs (cid:16) p l ( j ) (cid:12)(cid:12) full − p l ( j ) (cid:12)(cid:12) hyb (cid:17) (32) and E l rel ( j ) = abs (cid:16) p l ( j ) (cid:12)(cid:12) full − p l ( j ) (cid:12)(cid:12) hyb (cid:17) p l ( j ) (cid:12)(cid:12) full . (33)Equivalently, we define the mean blood pressure in thehomogenized vasculature of the hybrid model in the j -th REV as p v ( j ) (cid:12)(cid:12) hyb = 1 V REV j (cid:90) REV j p v dV (34)and as p v ( j ) (cid:12)(cid:12) full = 1 n nodes ,Λ S ∩ REV j n nodes ,Λ S ∩ REV j (cid:88) i =1 p ˆ v [ i ] (35)for the smaller vessels of the fully-resolved model. Thelatter is the mean blood pressure of all n nodes ,Λ S ∩ REV j nodes of the smaller blood vessels which lie inside the j -th REV. This allows us to define the absolute and rel-ative mean blood pressure error (in the smaller vessels)between the two models in each REV as E v abs ( j ) = abs (cid:16) p v ( j ) (cid:12)(cid:12) full − p v ( j ) (cid:12)(cid:12) hyb (cid:17) (36)and E v rel ( j ) = abs (cid:16) p v ( j ) (cid:12)(cid:12) full − p v ( j ) (cid:12)(cid:12) hyb (cid:17) p v ( j ) (cid:12)(cid:12) full . (37)Furthermore, we denote by ( · ) the mean value of theseerror measures over all n REV
REVs. Note also that boththe mean REV blood and IF pressure vary consider-ably between different REVs. The pressure difference alidation and parameter optimization of a hybrid embedded/homogenized solid tumor perfusion model 23Network Case E v abs [Pa] E v rel [%] E l abs [Pa] E l rel [%] R ,Λ L → Λ S R ,Λ SW1222 case 5 % 90 . .
25 49 . .
66 0 .
192 0 . . .
43 32 . .
08 0 .
091 0 . . .
10 24 . .
83 0 .
142 1 . . .
05 21 . .
72 0 .
100 1 . . .
88 73 . .
83 0 .
640 0 . . .
53 68 . .
58 0 .
567 0 . . .
48 63 . .
30 0 .
513 0 . . .
33 59 . .
10 0 .
444 0 . . .
64 41 . .
16 0 .
472 0 . . .
52 44 . .
29 0 .
489 0 . . .
22 44 . .
33 0 .
435 0 . . .
02 45 . .
33 0 .
407 0 . Table 6: Additional error measures for the agreement of both models. Shown are the absolute and relative errorof the hybrid approach in terms of mean REV blood pressure in smaller vessels and mean REV interstitial fluidpressure and the R -values for agreement between both variants in terms of flow from large to small vessels andflow in the entire vasculature. (All data includes the mean taken over five different sets of pressure boundaryconditions on the 1D network produced by the methodology described in Section 3.2.1, ”case X %” denotes thecase where X % of the 1D blood vessels are retained in the hybrid approach)between single REVs varies in a range of 800 − − λ = ˆ v → v M interpreted as a mass transfer term, or volumetric flowper length, as detailed in Section 2.3. Note that thiscan represent both a flow from large 1D vessels intothe homogenized vasculature if locally the pressure inthe resolved vasculature is bigger than the homogenized pressure or, vice versa, a flow from the homogenizedvasculature into the larger vessels if the homogenizedpressure is bigger than the blood pressure in the 1D vas-culature. Consequently, for each REV the flow betweenthe two compartments is given by the integral of theLM field along the part of the larger vessels Λ L ∩ REV j inside the specific REV j , or ˆ v → v M ( j ) (cid:12)(cid:12) hyb = (cid:90) Λ L ∩ REV j λ ds. (38)In the full model we can directly evaluate the masstransfer between large and small vessels inside the con-necting elements of both sets, which are those elementsof the smaller vessels where one node is part of Λ L andthe other part of Λ S . Assuming that the first node ispart of the larger vessels and the second one part of thesmaller vessels, the flow between large and small vesselsin the j -th REV is given by ˆ v → v M ( j ) (cid:12)(cid:12) full = n ele ,Λ L → Λ S ∩ REV j (cid:88) i =1 − πR µ ˆ v ∂p ˆ v ∂s (39)as the sum of the volumetric flows in the elements con-necting large and small vessels which lie inside the spe-cific REV j . The number of these elements is denotedby n ele ,Λ L → Λ S ∩ REV j in the previous equation. To com-pare the mass transfer between large and small vesselsin both models, we again define a coefficient of deter-mination as R ,Λ L → Λ S = 1 − (cid:80) n REV j =1 (cid:18) ˆ v → v M ( j ) (cid:12)(cid:12) full − ˆ v → v M ( j ) (cid:12)(cid:12) hyb (cid:19) (cid:80) n REV j =1 (cid:18) ˆ v → v M ( j ) (cid:12)(cid:12) full − µ (cid:18) ˆ v → v M ( j ) (cid:12)(cid:12) full (cid:19)(cid:19) , (40) with the respective mass transfer terms for the hybridand the full model for each REV. Again, µ ( · ) denotesthe mean of the mass transfer between large and smallvessels of the full model over all n REV
REVs.The reference solution of the fully resolved modelfor this volumetric flow per REV varies considerablybetween the different REVs and both positive values,representing an overall outflow from the larger vesselsinto the smaller vessels in a specific REV, and nega-tive values, representing an overall inflow into the largervessels from the smaller vessels in a specific REV, arepresent. This indicates that in- or outflow from largerto smaller vessels is indeed a meaningful quantity de-scribing the spatially varying flow patterns inside thevascular network. To reproduce this behaviour in thehybrid model variant, a good agreement with the refer-ence solution is desirable. The results for the coefficientof determination R ,Λ L → Λ S are again assembled inTable 6. The LS174T case shows the best agreementwith the fully-resolved model while the SW1222 casedelivers the worst results. We believe that this can beattributed to the much higher dispersion of the diame-ter and, thus, also the flow in the connectivity elements,which we have already studied by the coefficient of vari-ability in Table 3. The LS174T case, which has the leastdispersed distribution of both values, performs best inmatching the flow between larger and smaller vessels inthe hybrid model. There is a small decline of the agree-ment for higher percentages of retained vessels in allcases. However, the flow between large and small ves-sels is not included in the parameter optimization pro-cedure. Hence, we assume that the better performancein terms of the other quantities is at the expense of thismetric.Finally, we study the correspondence between thetwo models in terms of the blood flow in the entirevasculature Λ . Previously, in Table 5 only flow in thesmaller vessels Λ S , respectively, the homogenized vas-culature was investigated. For the full model, the totalflow in Λ in each REV in coordinate direction j is cal-culated as in (28), but now both large and small vesselsare taken into account. For the hybrid model, the totalflow can be obtained as the sum of the flow in the ho-mogenized vessels as given by (27) and the flow in thelarger, resolved vessels, that is, eqn. (28) evaluated forthe larger vessels of the hybrid model. The two quan-tities are compared in Table 6 defining a coefficient ofdetermination for flow in the entire vasculature R ,Λ as in (26). Evidently, the agreement between the twomodels is very good and much better than the previ-ously reported agreement of flow in the smaller vessels R ,Λ S . This is due to the fact that, as expected, flowis dominated by the larger vessels, the values calculated for flow in the entire vasculature are one to two orders ofmagnitude larger than the in the small vessels depend-ing on the investigated case. As we are able to matchthe pressure in the large vessels very well and, thus, alsothe flow therein, very good accordance can be achievedfor the total flow in both small and big vessels. As flowin the big vessels is decisive for the overall perfusion ofthe domain and could also be more easily acquired withexperiments for further validation this is an encourag-ing result for the applicability of the hybrid approach.Nevertheless, we demonstrate how to enhance the cor-respondence of the hybrid model also in terms of flowin the smaller vessels in the following section. In this section, we discuss some possible improvementsfor the hybrid model and implement one of them. Themost straightforward one would be to define the per-meability of the homogenized vessels not as a constantover the entire domain Ω v but per REV. Instead of anisotropic permeability tensor, one could easily integrateanisotropic effects based on the blood vessel structureinside each REV. Both has been done in the hybridmodel of Vidotto et al. (2019), where a diagonal perme-ability tensor with different permeabilities in all threecoordinate directions was employed. This could poten-tially augment the agreement in terms of mass fluxes inthe homogenized vasculature, which is the main sourceof error in the hybrid model. However, we did not inte-grate this into our optimization procedure since we be-lieve that this would result in overfitting of the chaoticflow in the tumor such that we would meet every singleboundary condition case very well but with largely dif-ferent results for the permeability tensors between thecases with distinct flow patterns. With a single scalarpermeability the results for the permeability betweendifferent boundary condition cases did not fluctuategreatly. Moreover, in a real-world case where only thearchitecture of the larger vessels is known, it seems un-realistic to deduce the entire permeability field from thelimited amount of information.Instead we tried to enhance the model by takinginformation of volume fractions of the smaller vesselsinto account. Our rationale behind this approach is thatwhile the complete structure of the smaller vessels can-not be obtained non-invasively, regions with higher orsmaller microvascular density of small vessels could stillbe identified. This information could then be employedto enrich the hybrid model. The overall trend we ob-served in Table 5 is that the higher the volume frac-tion of the homogenized vessels, the larger their perme-ability. It also reasonable to assume that areas with a alidation and parameter optimization of a hybrid embedded/homogenized solid tumor perfusion model 255 · − . . · ε vΛ S [ − ] | Q ˆ v | (cid:12)(cid:12) Λ S , f u ll [ µ m s − ] (a) SW1222, R = 0 .
50 0 0 . . · − . · ε vΛ S [ − ] | Q ˆ v | (cid:12)(cid:12) Λ S , f u ll [ µ m s − ] (b) LS174T, R = 0 .
66 2 3 4 · − · ε vΛ S [ − ] | Q ˆ v | (cid:12)(cid:12) Λ S , f u ll [ µ m s − ] (c) GL261, R = 0 . Fig. 10: Dependency of absolute flow in small vessels Λ S of the full model on volume fraction of small vessels inREVs for one representative case per tumor topology where 10 % of 1D blood vessels have been retained in hybridmodel (dashed lines indicate linear least squares fits with corresponding R -values)higher vascular volume fraction are more permeable toblood flow. Therefore, we investigated the relationshipof the volume fraction of smaller vessels ε vΛ S in eachREV on the perfusion of the smaller blood vessels inthe full model. Results are shown in Figure 10. Here,the absolute volumetric flow in each coordinate direc-tion (calculated as in (28) but not taking the flow direc-tion into account) is plotted over the volume fractionof the smaller vessels ε vΛ S . The clearest picture emergesfor the LS174T topology with a good correlation of flowin smaller vessels with their volume fraction. A similar,yet less distinctive trend is present for the SW1222 casewhereas no relationship can be observed for the GL261tumor.Therefore, inside each REV j we defined the isotro-pic permeability tensor as k v µ v ( j ) · I = α · ε vΛ S ( j ) · I , (41)that is, a simple linear dependency of the permeabil-ity in the j -th REV on the volume fraction of smallervessels in the j -th REV with proportionality constant α . We also tested a nonlinear Kozeny-Karman law, butobtained slightly better results with the linear fit. Thus,we will exclusively study this linear dependency here-after. Next, the optimization of the nonlinear least-squares problem (30) is performed for the proportional-ity constant α . Results are shown in Table 7 for the case10 %, which can be compared with the cases with con-stant permeability over the entire domain from Table 5.We obtained a slightly better agreement in terms of flowin the smaller vessels R ,Λ S and, thus, also for the to-tal coefficient of determination R for the SW1222 andGL261 case. Compared to that, the correspondence offlow in the smaller vessels was markedly better than the constant permeability case for the LS174T topol-ogy. This is coherent with Figure 10 where the latternetwork showed the most evident correlation of bloodflow on volume fraction. Thus, one could expect thatno significant improvement was possible for the GL261case where volume fraction and flow seem to be de-coupled. However, also for the SW1222 topology, whichshowed at least a moderate dependency of blood flow onvolume fraction, the agreement could not be increasedsignificantly. Therefore, at least for one of our casesit was beneficial to include blood vessel volume frac-tion information into the hybrid model while it was notdetrimental for the other two.It also is conceivable that at least preferential di-rections of smaller vessels can be detected non-inva-sively even though their complete structure cannot beresolved. A further enhancement of the model couldbe achieved when taking this information about theanisotropy of smaller vessels into account. However, wewant to emphasize that our whole study is based on nu-merical results. Experimental findings indicate no de-pendency between blood vessel diameter and flow intumors (Dewhirst and Secomb 2017; Dewhirst et al.1989; Leunig et al. 1992) and a high vascular densitydoes not automatically imply efficient perfusion, nutri-ent supply and drug delivery for solid tumors (Jain2005). These properties could make it impossible todeduce permeabilities of blood vessels inside tumorsfrom macroscopic quantities such as blood vessel vol-ume fractions or preferential directions. By contrast,non-invasive measurements of perfusion (D’Esposito etal. 2018; Thomas et al. 2000) could prove helpful toenhance the hybrid model.Similarly, improvements are possible for flow fromthe larger into the smaller vessels. In this contribu- α [ µ m ] R R R R ,Λ S R SW1222 37 . ± . . ± . . ± . Table 7: Results of the optimization procedure for non-constant permeability depending on volume fraction ofsmaller vessels. 10 % of 1D blood vessels have been retained in hybrid model for each tumor topology. Shown arethe proportionality constant α relating permeability and blood vessel volume fraction of smaller vessels inside eachREV according to (41). R -values for agreement between both variants in terms of blood pressure in large vessels,blood pressure in small vessels, IF pressure and flow in small vessels is additionally provided. Overall coefficientof determination R between fully-resolved and hybrid model is calculated according to (29). (All data includesmean taken over five different sets of pressure boundary conditions on the 1D network per case)tion, we have assumed equal pressures in resolved andhomogenized vasculature and, thereby, infinite (or atleast a very large) permeability governing the flow be-tween the two compartments such that a constraint ofequal pressures holds along the resolved 1D vascula-ture. This has the major advantage that the couplingbetween resolved and homogenized vasculature is essen-tially parameter-free. Only the penalty parameter hasto be chosen large enough such that a sufficient accu-racy in the pressure constraint is achieved as describedin Remark 4. The GL261 and LS174T case had a lessdispersed distribution of the radius in connecting ele-ments and, thus, of the permeability between larger andsmaller vessels. For these topologies, our approach couldestimate the mass transfer between larger 1D vesselsand smaller homogenized vessels more accurately. TheSW1222 case had a much higher variability of radiusand flow between Λ L and Λ S . In this case, it could beadvantageous to employ finite permeabilities to modelthe mass transfer and assign higher permeabilities toREVs or regions along the larger vessels where manybranches go away from the main vessels. However, thiswould require additional parameterization of the modelas well as additional data on regions where a lot of flowfrom larger into smaller vessels can be expected.In addition, we have so far only employed very sim-ple algorithms to optimize the parameters of the hybridmodel. A much more powerful framework for coarse-graining physical models has been developed by Grigoand Koutsourelakis (2019) and tested for flow throughporous media. This could also be applied in our caseto infer the optimal parameters of the hybrid modelper REV. However, this would require much more mi-crostructural features, such as tortuosity, blood vesseldistances or radius data on the smaller homogenizedblood vessels to calibrate the hybrid model. Again, it isquestionable if this data can be acquired non-invasive-ly and if these parameters are determining blood flowthrough tumors. In this work, we have studied a hybrid embedded/homo-genized model for computational modeling of solid tu-mor perfusion. Its guiding principle is that the completemorphology of vascular networks including blood vesseldiameters and topology cannot be acquired with cur-rently available in-vivo imaging techniques. Thus, fully-resolved discrete models relying on this data cannot beapplied in ”real world” scenarios. If, however, the struc-ture of larger vessels constituting the main branches ofthe vasculature is available, our hybrid approach whereonly these larger branches are completely resolved isa sensible alternative. By contrast, the smaller scalevessels are homogenized such that their exact struc-ture is not required anymore. This results in a two-compartment or double porosity formulation where thelarger vessels are still represented as one-dimensionalembedded inclusions. The coupling between the resolvedand homogenized part of the vasculature is realized viaa line-based pressure constraint along the 1D larger ves-sels, which we enforce with a mortar-type formulationwith penalty regularization. This also has the advantagethat compared to previous hybrid models no additionalparameter – apart from the penalty parameter whichhas to be chosen large enough – is required to couplethe two distinct representations of the vasculature.The results of the hybrid model have been comparedwith reference solutions generated by a fully-resolved1D-3D model. For that, three different network topolo-gies extracted from three different tumor types grownin mice have been employed. These topologies consistof up to 420 000 vessel segments and have dimensions ofup to 6 mm × ×
11 mm. To date, this is the largestand most challenging test case for a hybrid model, es-pecially considering the abnormal and tortuous struc-ture of the networks typical for the vasculature insidetumors. We have further shown how we artificially gen-erate the hybrid from the fully-resolved model, definerepresentative elementary volumes and assign bound- alidation and parameter optimization of a hybrid embedded/homogenized solid tumor perfusion model 27 ary conditions. We are confident that the artificiallycreated topologies of larger vessels are representative ofreal in-vivo imaging data sets of larger vessels insidetumors such that they enable us to draw meaningfulconclusions for more realistic scenarios where the fulltopology is not available such that a hybrid approachis the only option.For comparison of the results of the two models,we have defined several rigorous metrics involving theblood pressure in both resolved and homogenized vas-culature, the pressure in the interstitial fluid and bloodflow in the homogenized vasculature. These metrics havethen been employed to obtain the optimal parametersfor the hybrid model and to study its accuracy w.r.t. thefully-resolved one. We have obtained very good agree-ment in terms of blood pressure in the larger vesselsand IF pressure. Larger deviations are present for bloodpressure and flow in the homogenized vasculature. How-ever, these limitations can be expected since the infor-mation on the smaller vessels is not retained in the hy-brid model. Overall, the best correspondence has beenachieved for the SW1222 case which also had the clear-est vascular structure and distinction between largerand smaller vessels. All topologies showed a very goodagreement in terms of REV IF pressure and REV bloodpressure in smaller vessels with mean deviations in arange of 20 −
70 Pa and 40 −
110 Pa resp. 0 . − . . − . −
10 % of allblood vessels segments by keeping them in the hybridmodel since there is only a marginal improvement of theagreement with the fully-resolved model in terms of allinvestigated metrics when retaining a higher percent-age (15 −
20 %) of blood vessels. Concerning the flowbetween smaller and larger vessels the error was mainlycaused by the large variability of diameter and flow inthe connectivity elements between large and small ves-sels for the SW1222 case. Possibly, this error could bereduced by allowing a varying permeability for couplingthe two compartments. By including information aboutthe blood vessel volume fraction of smaller blood vesselsinto the definition of their permeability tensor a bet-ter agreement with flow therein could be achieved forthe LS174T case. Nevertheless, the abnormal vascularstructure and blood flow patterns of tumor vasculaturecould impede this approach.Several other potential improvements have been dis-cussed and remain subject to future work. Furthermore,the inclusion of species transport to simulate drug deliv-ery or nutrient transport lies at hand. Species transportincluding the coupling between resolved and homoge-nized vasculature is possible within our hybrid multi-phase tumor growth model (Kremheller et al. 2019) andwe have already studied nanoparticle delivery to solid tumors employing the homogenized compartment only(Wirthl et al. 2020). These models could ultimately en-hance our understanding of the limitations of currentdrug delivery strategies and aid in devising more tar-geted therapies.The next step towards a more realistic scenario forhybrid computational models of tissue perfusion is todevise a strategy which combines data which is avail-able non-invasively (Li et al. 2020). For instance, thiscould be measurements on tissue perfusion, hypoxic ar-eas, REV pressure data, volume fractions of homoge-nized blood vessels or their preferential direction overthe entire domain. The methods and metrics developedhere could also be applied to either fully-resolved or hy-brid models in combination with this data. Inevitably,this implies that the large majority of the (homogenizedor resolved) pressure boundary conditions is unknown.In this case, the boundary conditions together with thepermeability of both porous networks would be the un-knowns of the optimization procedure to generate a re-alistic model of tissue perfusion through solid tumors.Subsequently, this could be employed for in-silico stud-ies of drug delivery or optimization of treatment strate-gies.
Acknowledgements
The authors gratefully acknowledge thesupport of the Technical University of Munich – Institute forAdvanced Study, funded by the German Excellence Initiativeand the T ¨UV S ¨UD Foundation. Research reported in thispublication was supported by the National Cancer Instituteof the National Institutes of Health under Award NumberU54CA210181. The content is solely the responsibility of theauthors and does not necessarily represent the official views ofthe National Institutes of Health. The software QUEENS wasprovided by the courtesy of AdCo Engineering GW GmbH,which is gratefully acknowledged. We would also like to thankthe authors of REANIMATE (D’Esposito et al. 2018; Sweeneyet al. 2019), especially Paul W. Sweeney for sharing theirdata and making their code publicly available. In addition,we gratefully thank Barbara Wirthl and Anh-Tu Vuong fordiscussions about various aspects of the two models.
References
BACI (2021). BACI: A Comprehensive Multi-Physics Simulation Framework. https://baci.pages.gitlab.lrz.de/website/ . Accessed: 2021-02-10.Baluk, P., Hashizume, H., and McDonald, D. M. (2005).Cellular abnormalities of blood vessels as targets incancer.
Curr. Opin. Genet. Dev. , 15(1):102–111.Baxter, L. T. and Jain, R. K. (1989). Transport of fluidand macromolecules in tumors. I. Role of interstitialpressure and convection.
Microvasc. Res. , 37(1):77–104.
Biehler, J., Nitzler, J., Brandstaeter, S., Wall, W. A.,and Gravemeier, V. (2021).
QUEENS – A SoftwarePlatform for Uncertainty Quantification, Physics-Informed Machine Learning, Bayesian Optimization,Inverse Problems and Simulation Analytics: UserGuide . AdCo Engineering GW .Boucher, Y., Brekken, C., Netti, P. A., Baxter, L. T.,and Jain, R. K. (1998). Intratumoral infusion of fluid:estimation of hydraulic conductivity and implicationsfor the delivery of therapeutic agents. Br. J. Cancer ,78(11):1442–1448.Carmeliet, P. and Jain, R. K. (2000). Angiogenesis incancer and other diseases.
Nature , 407(6801):249–257.Cattaneo, L. and Zunino, P. (2014). A computationalmodel of drug delivery through microcirculation tocompare different tumor treatments.
Int. J. Numer.Meth. Biomed. Engng. , 30(11):1347–1371.Cattaneo, L. and Zunino, P. (2015). Computationalmodels for fluid exchange between microcirculationand tissue interstitium.
Networks Heterog. Media ,9(1):135–159.Chapman, S. J., Shipley, R. J., and Jawad, R. (2008).Multiscale Modeling of Fluid Transport in Tumors.
Bull. Math. Biol. , 70(8):2334.D’Angelo, C. (2007).
Multiscale modelling ofmetabolism and transport phenomena in living tis-sues . PhD thesis, EPFL PP, Lausanne.D’Angelo, C. (2012). Finite Element Approximationof Elliptic Problems with Dirac Measure Terms inWeighted Spaces: Applications to One- and Three-dimensional Coupled Problems.
SIAM J. Numer.Anal. , 50(1):194–215.D’Angelo, C. and Quarteroni, A. (2008). On the cou-pling of 1D and 3D diffusion-reaction equations. Ap-plications to tissue perfusion problems.
Math. ModelsMethods Appl. Sci. , 18(08):1481–1504.Davit, Y., Bell, C. G., Byrne, H. M., Chapman, L.A. C., Kimpton, L. S., Lang, G. E., Leonard, K. H. L.,Oliver, J. M., Pearson, N. C., Shipley, R. J., Waters,S. L., Whiteley, J. P., Wood, B. D., and Quintard, M.(2013). Homogenization via formal multiscale asymp-totics and volume averaging: How do the two tech-niques compare?
Adv Water Resour. , 62:178–206.D’Esposi to, A., Sweeney, P. W., Ali, M., Saleh, M.,Ramasawmy, R., Roberts, T. A., Agliardi, G., Des-jardins, A., Lythgoe, M. F., Pedley, R. B., Shipley,R., and Walker-Samuel, S. (2018). Computationalfluid dynamics with imaging of cleared tissue and ofin vivo perfusion predicts drug uptake and treatmentresponses in tumours.
Nat. Biomed. Eng. , 2(10):773–787. Dewhirst, M. W. and Secomb, T. W. (2017). Transportof drugs from blood vessels to tumour tissue.
Nat RevCancer , 17:738.Dewhirst, M. W., Tso, C. Y., Oliver, R., Gustafson,C. S., Secomb, T. W., and Gross, J. F. (1989). Mor-phologic and hemodynamic comparison of tumor andhealing normal tissue microvasculature.
Int. J. Ra-diat. Oncol. , 17(1):91–99.Ferrari, M. (2010). Frontiers in cancer nanomedicine:Directing mass transport through biological barriers.
Trends Biotechnol. , 28(4):181–188.Formaggia, A., Quarteroni, A., and Veneziana, A.(2009).
Cardiovascular Mathematics: Modeling andsimulation of the circulatory system . Springer-Verlag.Forster, J. C., Harriss-Phillips, W. M., Douglass, M.J. J., and Bezak, E. (2017). A review of the devel-opment of tumor vasculature and its effects on thetumor microenvironment.
Hypoxia , 5:21–32.Fry, B. C., Lee, J., Smith, N. P., and Secomb, T. W.(2012). Estimation of Blood Flow Rates in Large Mi-crovascular Networks.
Microcirculation , 19(6):530–538.Geuzaine, C. and Remacle, J.-F. (2009). Gmsh: A 3-dfinite element mesh generator with built-in pre- andpost-processing facilities.
Int. J. Numer. MethodsEng. , 79(11):1309–1331.Grigo, C. and Koutsourelakis, P.-S. (2019). A physics-aware, probabilistic machine learning frameworkfor coarse-graining high-dimensional systems in theSmall Data regime.
J. Comput. Phys. , 397:108842.Heldin, C.-H., Rubin, K., Pietras, K., and Ostman, A.(2004). High interstitial fluid pressure [mdash] an ob-stacle in cancer therapy.
Nat Rev Cancer , 4(10):806–813.Holash, J., Wiegand, S. J., and Yancopoulos, G. D.(1999). New model of tumor angiogenesis: dynamicbalance between vessel regression and growth medi-ated by angiopoietins and VEGF.
Oncogene , 18:5356.Jain, R. K. (1994). Barriers to Drug Delivery in SolidTumors.
Sci. Am. , 271(1):58–65.Jain, R. K. (2005). Normalization of Tumor Vascula-ture: An Emerging Concept in Antiangiogenic Ther-apy.
Science , 307(5706):58 LP – 62.Kojic, M., Milosevic, M., Simic, V., Koay, E. J., Flem-ing, J. B., Nizzero, S., Kojic, N., Ziemys, A., andFerrari, M. (2017). A composite smeared finite ele-ment for mass transport in capillary systems and bi-ological tissue.
Comput. Methods Appl. Mech. Eng. ,324(Supplement C):413–437.K¨oppl, T., Vidotto, E., and Wohlmuth, B. (2020).A 3D-1D coupled blood flow and oxygen transportmodel to generate microvascular networks.
Int. J.Numer. Method. Biomed. Eng. , 36:e3386. alidation and parameter optimization of a hybrid embedded/homogenized solid tumor perfusion model 29
K¨oppl, T., Vidotto, E., Wohlmuth, B., and Zunino, P.(2018). Mathematical modeling, analysis and numer-ical approximation of second-order elliptic problemswith inclusions.
Math. Models Methods Appl. Sci. ,28(05):953–978.K¨oppl, T. and Wohlmuth, B. (2014). Optimal A Pri-ori Error Estimates for an Elliptic Problem withDirac Right-Hand Side.
SIAM J. Numer. Anal. ,52(4):1753–1769.Kremheller, J., Vuong, A.-T., Schrefler, B. A., andWall, W. A. (2019). An approach for vascular tumorgrowth based on a hybrid embedded/homogenizedtreatment of the vasculature within a multiphaseporous medium model.
Int. J. Numer. Method.Biomed. Eng. , 35:e3253.Kremheller, J., Vuong, A.-T., Yoshihara, L., Wall,W. A., and Schrefler, B. A. (2018). A monolithic mul-tiphase porous medium framework for (a-)vasculartumor growth.
Comput. Methods Appl. Mech. Eng. ,340:657–683.Leunig, M., Yuan, F., Menger, M. D., Boucher, Y.,Goetz, A. E., Messmer, K., and Jain, R. K. (1992).Angiogenesis, Microvascular Architecture, Microhe-modynamics, and Interstitial Fluid Pressure duringEarly Growth of Human Adenocarcinoma LS174T inSCID Mice.
Cancer Res. , 52(23):6553 LP – 6560.Li, J., Chekkoury, A., Prakash, J., Glasl, S., Vetschera,P., Koberstein-Schwarz, B., Olefir, I., Gujrati, V.,Omar, M., and Ntziachristos, V. (2020). Spatialheterogeneity of oxygenation and haemodynamicsin breast cancer resolved in vivo by conical multi-spectral optoacoustic mesoscopy.
Light Sci. Appl. ,9(1):57.MathWorks Inc. (2018). Documentation of matlab’s alphaShape function. https://de.mathworks.com/help/matlab/ref/alphashape.html . Accessed:2021-02-10.Matsumura, Y. and Maeda, H. (1986). A NewConcept for Macromolecular Therapeutics in Can-cer Chemotherapy: Mechanism of Tumoritropic Ac-cumulation of Proteins and the Antitumor AgentSmancs.
Cancer Res. , 46(12 Part 1):6387 LP – 6392.Morikawa, S., Baluk, P., Kaidoh, T., Haskell, A., Jain,R. K., and McDonald, D. M. (2002). Abnormalities inPericytes on Blood Vessels and Endothelial Sproutsin Tumors.
Am. J. Pathol. , 160(3):985–1000.Nabil, M., Decuzzi, P., and Zunino, P. (2015). Mod-elling mass and heat transfer in nano-based cancerhyperthermia.
R. Soc. Open Sci. , 2(10):150447.Nabil, M. and Zunino, P. (2016). A computationalstudy of cancer hyperthermia based on vascular mag-netic nanoconstructs.
R. Soc. Open Sci. , 3(9):160287. Nizzero, S., Ziemys, A., and Ferrari, M. (2018). Trans-port Barriers and Oncophysics in Cancer Treatment.
Trends in Cancer , 4(4):277–280.Popp, A., Gitterle, M., Gee, M. W., and Wall, W. A.(2010). A dual mortar approach for 3D finite defor-mation contact with consistent linearization.
Int. J.Numer. Methods Eng. , 83(11):1428–1465.Pries, A. R. and Secomb, T. W. (2005). Microvascu-lar blood viscosity in vivo and the endothelial sur-face layer.
Am. J. Physiol. - Hear. Circ. Physiol. ,289(6):H2657–H2664.SciPy (2020). Documentation ofscipy.optimize.least squares (version 1.5.2). https://docs.scipy.org/doc/scipy-1.5.2/reference/generated/scipy.optimize.least_squares.html .Accessed: 2021-02-10.Secomb, T. W. and Hsu, R. (1994). Simulation of O2transport in skeletal muscle: diffusive exchange be-tween arterioles and capillaries.
Am. J. Physiol. -Hear. Circ. Physiol. , 267(3):H1214–H1221.Secomb, T. W., Hsu, R., Park, E. Y. H., and Dewhirst,M. W. (2004). Green’s Function Methods for Anal-ysis of Oxygen Delivery to Tissue by MicrovascularNetworks.
Ann. Biomed. Eng. , 32(11):1519–1529.Shipley, R. J. and Chapman, S. J. (2010). MultiscaleModelling of Fluid and Drug Transport in VascularTumours.
Bull. Math. Biol. , 72(6):1464–1491.Shipley, R. J., Smith, A. F., Sweeney, P. W., Pries,A. R., and Secomb, T. W. (2019). A hybrid discrete-continuum approach for modelling microcirculatoryblood flow.
Math. Med. Biol.
Steinbrecher, I., Mayr, M., Grill, M. J., Kremheller, J.,Meier, C., and Popp, A. (2020). A mortar-type finiteelement approach for embedding 1D beams into 3Dsolid volumes.
Comput. Mech.
Sweeney, P. W., D’Esposito, A., Walker-Samuel, S., andShipley, R. J. (2019). Modelling the transport offluid through heterogeneous, whole tumours in silico.
PLoS Comput. Biol. , 15(6):e1006751.Sweeney, P. W., Walker-Samuel, S., and Shipley, R. J.(2018). Vascular and interstitial flow solver for dis-crete microvascular networks.Thomas, D. L., Lythgoe, M. F., Pell, G. S., Calamante,F., and Ordidge, R. J. (2000). The measurementof diffusion and perfusion in biological systems us-ing magnetic resonance imaging.
Phys. Med. Biol. ,45(8):R97–R138.Trelis (2020). Trelis (Version 17.0) [Computer software].American Fork, UT: csimsoft. Retrieved from http://csimsoft.com .Tully, B. and Ventikos, Y. (2011). Cerebral water trans-port using multiple-network poroelastic theory: ap-plication to normal pressure hydrocephalus.
J. Fluid
Mech. , 667:188–215.Verdugo, F. and Wall, W. A. (2016). Unified computa-tional framework for the efficient solution of n-fieldcoupled problems with monolithic schemes.
Comput.Methods Appl. Mech. Eng. , 310:335–366.Vidotto, E., Koch, T., K¨oppl, T., Helmig, R., andWohlmuth, B. (2019). Hybrid models for simulat-ing blood flow in microvascular networks.
MultiscaleModel. Simul. , 17(3):1076–1102.Wirthl, B., Kremheller, J., Schrefler, B. A., and Wall,W. A. (2020). Extension of a multiphase tumourgrowth model to study nanoparticle delivery to solidtumours.
PLOS ONE , 15(2):1–22.Wohlmuth, B. I. (2000). A Mortar Finite ElementMethod Using Dual Spaces for the Lagrange Mul-tiplier.
SIAM J. Numer. Anal. , 38(3):989–1012.Yang, B., Laursen, T. A., and Meng, X. (2005). Two di-mensional mortar contact methods for large deforma-tion frictional sliding.