Validation and parameter optimization of a hybrid embedded/homogenized solid tumor perfusion model

The goal of this paper is to investigate the validity of a hybrid embedded/homogenized in-silico approach for modeling perfusion through solid tumors. The rationale behind this novel idea is that only the larger blood vessels have to be explicitly resolved while the smaller scales of the vasculature are homogenized. As opposed to typical discrete or fully-resolved 1D-3D models, the required data can be obtained with in-vivo imaging techniques since the morphology of the smaller vessels is not necessary. By contrast, the larger vessels, whose topology and structure is attainable non-invasively, are resolved and embedded as one-dimensional inclusions into the three-dimensional tissue domain which is modeled as a porous medium. A sound mortar-type formulation is employed to couple the two representations of the vasculature. We validate the hybrid model and optimize its parameters by comparing its results to a corresponding fully-resolved model based on several well-defined metrics. These tests are performed on a complex data set of three different tumor types with heterogeneous vascular architectures. The correspondence of the hybrid model in terms of mean representative elementary volume blood and interstitial fluid pressures is excellent with relative errors of less than 4%. Larger, but less important and explicable errors are present in terms of blood flow in the smaller, homogenized vessels. We finally discuss and demonstrate how the hybrid model can be further improved to apply it for studies on tumor perfusion and the efficacy of drug delivery.


INTRODUCTION
Mathematical modeling of blood flow and mass transport is of increasing importance to study a number of highly relevant biomedical questions in health and disease. Computational tools offer the possibility to gain new insight into physiologically relevant processes such as the transport of nutrients, oxygen or drugs across the vascular system and inside the tissue microenvironment. These methods can ultimately lead to a new rationale for developing and non-invasive testing of novel therapies. 1 Concurrently, such in-silico models will make the design of drugs both cheaper and faster.
In this paper, we are concerned with the simulation of blood flow and tissue perfusion at the scale of the microcirculation with a special focus on solid tumors where transport processes can be decisive for disease progression and treatment efficacy. This includes, first, the vasculature, which is embedded in the surrounding tissue, second, passage across the blood vessel walls into the surrounding extravascular space and, third, flow of the fluid filling this space, namely the interstitial fluid (IF). Subsequently, we will distinguish between three different modeling strategies for these transport processes, namely discrete, continuum and hybrid approaches. All strategies have been developed and implemented in the context of our vascular multiphase tumor growth model. 2,3 For the discrete or fully-resolving variant, we follow a common modeling approach, where the vasculature is represented by a network of one-dimensional blood vessel segments embedded in the encompassing three-dimensional tissue domain which is modeled as a porous medium. A 1D partial differential equation (PDE) is employed to model mass transport in the vasculature while a corresponding 3D PDE governs the surrounding IF. Both domains are coupled via source terms which account for the exchange across the blood vessel wall. Such models are well-studied and first contributions include the so-called embedded multiscale method developed by D' Angelo and Quarteroni 4,5,6 and the Green's function method of Secomb et al. 7,8 More recent approaches with such a philosophy include drug delivery 9,10 , hyperthermia treatment 11,12 and a combination of a numerical framework with optical imaging to predict fluid and species mass transport through whole tumors with heterogenous blood vessel architecture. 13,14 These approaches are commonly termed discrete models in distinction from continuum models which involve a homogenization procedure. Thereby, the vasculature is approximated as a homogeneous porous medium resulting in two distinct pore spaces which are the aforementioned interstitium and the homogenized vasculature. Flow in both domains is modeled via the Darcy equation and suitable exchange terms are defined. 2,15,16,17,18,19,20 These two distinct approaches have different use cases: Discrete models can and should be applied when the entire structure of the vasculature including the smallest scales, i.e., the capillaries, is known and its resolution is needed for the question at hand. This is usually restricted to small domain sizes of an order of several mm 3 . By contrast, continuum models are used to simulate mass transport at larger scales, e.g., through whole organs. Both approaches have advantages and disadvantages: On the one hand, the computational cost of continuum models is usually smaller than for discrete ones which makes the application to larger domains possible in the first place. On the other hand, the information about the exact morphology of the vascular network is lost such that blood flow can only be described in an averaged sense. Discrete models, however, are computationally more expensive. Furthermore, they require the full structure of the part of the vasculature under consideration. This is usually realized via a graph whose edges are assigned the radius of the blood vessel segments between nodes. Such high-resolution data including blood vessel radii, connectivity and positions can at present only be acquired through ex-vivo imaging. 21 In addition, the acquisition of high-quality data is still challenging and error-prone especially on the finest scales. 22 By contrast, in-vivo imaging is currently only possible for larger vessels and flow therein. 21,23 Therefore, discrete models rely on data which is not available via noninvasive imaging. An additional difficulty is the assignment of blood pressure or flow boundary conditions which can only be estimated for large networks. 14,24 In any case, validation of these models is usually only performed on macroscopic quantities such as tissue perfusion 13 since measuring flow or pressures inside single micro-vessels is not possible. 14 This has motivated the development of hybrid methods which are especially suited for cases where the full vascular morphology is unknown or too large to be modeled with a discrete approach. The idea behind them is to explicitly resolve the larger vessels through a discrete model and to use a homogenized approach for the capillary bed. Next to our own work, 3 such hybrid approaches have also been developed by Vidotto et al. 25 , Shipley et al. 21 and Kojic et al. 26 Compared to pure homogenized formulations, their advantage is that the structure of the larger vessels is retained and, therefore, the heterogeneity of blood flow and pressure in the major vessel branches is better represented. Moreover, compared to discrete models, less anatomic data is needed since the morphology of the smallest vessels is not required. This could also have the additional advantage of a smaller computational cost and make them applicable to larger domains. Also quantities typically needed for validation such as tissue perfusion, blood flow or pressures at the resolution of current imaging techniques can equally be acquired from hybrid models. A related approach, where no homogenization of the capillaries is needed, is to generate a discrete surrogate network of the smaller scales based on the oxygen demand of the tissue. 22 We have previously incorporated a hybrid method for coupling discrete, one-dimensional blood vessels with a homogenized representation of the vasculature 3 into our vascular multiphase tumor growth model. 2 Therein, we couple a discrete representation of the pre-existing vasculature with a homogenized representation of the neovasculature which is formed during angiogenesis. For that, we employ constraint enforcement strategies which are well-known from solid mechanics. In our previous paper, the main focus was on modeling vascular tumor growth. In this contribution, we validate the applicability of the hybrid embedded/homogenized approach for the study of perfusion through solid tumors. Here, accurate models of fluid mass transport are of high relevance since efficient drug delivery to cancerous tissue relies on the fact that the drug reaches a large fraction of the tumor cells. Therefore, physiological characteristics such as microvascular flow, the structure of the extracellular matrix or the IF pressure profile may influence the transport of drugs through tumor tissue and, hence, ultimately the success of treatment. 1 For instance, increased interstitial pressure due to highly permeable vessels and inefficient lymphatic drainage has  been identified as an obstacle for successful drug delivery. 27,28,29 Novel nanoparticle-based therapies aim to exploit these properties of the tumor vasculature for more specific targeting of tumor sites. 30 Appropriate models of these transport phenomena can provide additional insight into and guidelines for drug design. In the context of cancer, this paradigm shift is described by the concept of transport oncophysics with the objective to engineer drugs with optimized transport properties. 31,32 Sound computational models are required to achieve this goal. We therefore took great care in the development of our hybrid model, i.e., both in the theoretical basis and its implementation. In this paper the main focus is on the validation of our hybrid embedded/homogenized scheme with three complex tumor-specific vascular networks based on large tissue samples containing more than 100 000 blood vessels. 13,14 We put a special emphasis on the extraction of the larger vessels from the fully-resolved network data such that it qualitatively matches the topology and distribution of larger vascular structures inside tumors available via in-vivo imaging. Thus, we make sure that the hybrid approach is investigated for cases which closely resemble real-life scenarios where the structure of the considered part of the microcirculation is not entirely known. Here, the complete topology of the vasculature in the given tissue domain is available which allows us to generate reference solutions with a fully-resolved model and to quantify the error introduced by the homogenization in the hybrid model. We evaluate the error by means of several well-defined metrics involving the agreement of pressures and flow between the two models. Concurrently, the parameters of the hybrid model are identified such that the correspondence of the models is maximized. Evaluating the error of the hybrid model in comparison to a fully-resolved one is a first and indispensable step towards realistic hybrid models of tumor perfusion relying only on non-invasively available physiological data. For a full validation and parameter optimization, similar methods as applied herein need to be combined with advanced in-vivo imaging techniques. The comparison of two purely numerical approaches and the inverse identification of the optimal parameters allows us to investigate the hybrid model in a controlled environment unaffected by any further influences such as uncertainties in experimental or clinical data.
The remainder of this work is structured as follows: We introduce both the hybrid and the fully-resolved model in Section 2. The employed tumor vasculature data sets as well as the setup of the models including the assignment of boundary conditions and the extraction of the hybrid model from the fully-resolved one are described in Section 3. Numerical experiments to compare the accuracy of the hybrid model w.r.t. the full model and to evaluate its main errors are conducted in Section 4. We illustrate some possible improvements of the hybrid model in Section 5 before summarizing our findings in Section 6.

MATHEMATICAL MODELS AND NUMERICAL METHODS
In this section, we describe the mathematical models we employ to solve the interaction between microcirculation and interstitial tissue perfusion including their main simplifications. We outline both our fully-resolved and our hybrid approach and their discretization by means of the finite element method (FEM).

Problem setting
As in other publications 3,4,5,6,7,8,9,10,11,12,13,14 , topology and structure of the microcirculation is described by a graph with straight edges, i.e., blood vessel segments. The segments connect the nodes of the network. A radius is assigned to each segment Λ . Available experimental data including the one employed here is also commonly provided in this format. Therefore, the vascular domain is given as the union of straight cylinders which are embedded in the three-dimensional tissue domain Ω, see also Figure 1a. Based on the flow physics and the quantities of interest, it is obvious to employ a one-dimensional blood flow model for this part of the vascular system. In the following, we will denote the 1D embedded blood vessel network with Λ. Similar to Vidotto et al. 25 , we further divide it into two subsets Λ L and Λ S which correspond to the larger and smaller vessels in the network such that with the index sets of large and small blood vessel segments L and S , respectively. We will show in detail how this partition is realized in Section 3.3. Whereas the larger vessels are kept in the hybrid model, the smaller scale vessels are replaced by an appropriate homogenized representation as a porous network occupying the domain Ω ⊂ Ω, cf. Figure 1b.

Fully-resolved 1D-3D model
Gravity, inertial effects and the pulsatility of blood flow are neglected which are valid assumptions since we deal with microcirculatory flow. The balance of mass in the 1D vasculature domain Λ is then given by the following equation where we have applied the Hagen-Poiseuille law for flow in cylindrical pipes and assumed constant blood density ̂ . Here and in the following, quantities defined on the 1D vasculature domain are denoted by superscript̂ . In the previous equation, is the blood vessel radius, ̂ the pressure inside the vasculature, ̂ the blood viscosity and the arc length coordinate along the 1D blood vessel segment. To account for the non-Newtonian behaviour of blood, we employ the algebraic relationship developed by Pries and Secomb 33 for in-vivo blood viscosity depending on blood vessel diameter and hematocrit. As in 14 , we fix the hematocrit to 0.45, thus, the blood viscosity ̂ in each individual blood vessel segment depends only on its diameter. Finally, the right-hand-side term̂ is employed to model leakage of fluid across the blood vessel wall into the interstitium. For that, we use Starling's law with hydraulic conductance ,̂ , density of blood plasma , oncotic reflection coefficient and the oncotic pressures of blood and the interstitial fluid (IF) . In summary, the transvascular flux from the vascular network into the interstitial fluid is proportional to the pressure difference between vasculature and IF whose pressure in (3) is denoted as . It has long been known that blood vessels inside tumors are leakier than normal ones which, in combination with a non-functional lymphatic system, leads to increased interstitial pressure inside solid tumors and, concurrently, resistance to efficient drug delivery. 27,28,29 Note that our data sets are whole-tumor blood vessel networks where also larger vessels are leaky 13,14 which is why we apply the transvascular exchange term (3) also on the subset of larger vessels Λ L . As in related works, the tissue domain Ω is modeled as a porous medium. Therefore, flow in the interstitial fluid is accounted for by the following Darcy equation with (isotropic) permeability = ⋅ and IF viscosity . Hence, the primary variable for fluid flow through the tissue is the IF pressure . The right hand side represents the counterpart of the leakage of fluid from the vasculature into the IF from (2). As proposed by D'Angelo and Quarteroni 4,5,6 this mass transfer term is concentrated as a Dirac measure Λ along the centerline Λ of the vasculature resulting in a 1D-3D coupled problem. The mathematical properties including reduced convergence rates due to the introduced singular line source in the 3D pressure field are extensively studied in 5,6,34 . Alternative 2D-3D coupled approaches, where the mass exchange is evaluated at the lateral surfaces of the cylindrical blood vessel segments, have been proposed to increase the regularity of the solution. 22,35 However, in our data sets the diameter is smaller than the element size ℎ in the 3D domain, see also Table 1. Therefore, we modify the approach of 5,6 , which would involve taking the average value of the pressure in the IF at the outer surface 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 ℎ > . 4,3 This has recently also been investigated in the analogous solid mechanics problem of embedding thin 1D structures, i.e., beams, into 3D solid volumes. 36 The weak form of the 1D-3D coupled problem may be written as with test functions ̂ defined on the 1D domain and defined on the 3D domain. Our approach allows for non-matching 1D discretizations Λ ℎ and 3D discretizations Ω ℎ such that the two domains can be meshed independently of each other. This requires the numerical integration of products of 1D shape functions with 3D shape functions and products of 3D shape functions with 3D shape functions along the one-dimensional discretization Λ ℎ which we realize via a segment-based line integration approach 3,36 to avoid integration over kinks of shape functions, see Appendix A. After space discretization, the nodal primary variables of both domains are ̂ ∈ ℝ nodes,Λ and ∈ ℝ nodes,Ω , that is, the nodal blood pressure in the discretized 1D domain and the nodal IF pressure in the discretized 3D domain, which consist of nodes,Λ and nodes,Ω , respectively. Details on the employed boundary conditions are given in Section 3.2.1. Finally, we arrive at the global system of equations, which may be written as a 2 × 2 block matrix Herein, the main diagonal blocks comprise contributions from the diffusive term and the exchange term in (5a) and (5b) while the off-diagonal submatrices ̂ and ̂ contain the "mixed" contributions from the exchange term. The right hand side terms represent the constant contribution to the exchange term stemming from the oncotic pressures in (3). To solve the coupled linear system (7) we employ a GMRES iterative solver in combination with the AMG(BGS) block preconditioner presented in Verdugo and Wall. 37

Hybrid 1D-3D model
The main idea behind our hybrid 1D-3D model, based on our previous work introduced in 3 in the context of our vascular multiphase tumor growth model 2 , is the following: The full resolution of the larger vessels Λ L is kept, i.e., these are still modeled as a 1D embedded vasculature. Consequently, the hierarchy, topology and vascular properties such as individual blood vessel radii and viscosities of each segment are retained, see also Figure 1b. The smaller vessels Λ S , for which this high-resolution data might either not be available through non-invasive imaging techniques or susceptible to errors, are instead represented as an additional porous network. This results in a double-porosity formulation where the first porous network is, as before, the interstitial space and the second one the smaller vessels occupying the domain Ω . In the following, we will present the governing equations and the space discretization of this formulation.
As stated above, the model for the larger vessels does not change. Therefore, the mass balance equation inside the large vessels is given by with the only difference to (2) being that it holds only on the subset Λ L ⊂ Λ of bigger vessels. The mass balance equation for the smaller vessels Λ S is replaced by a homogenized Darcy equation in the vascular domain Ω , which we formulate as The unknown in this equation is the blood pressure in the homogenized part of the vasculature which is now defined in the entire 3D domain Ω , thereby replacing the blood pressure of the smaller vessels in the 1D domain Λ S as the variable governing flow inside the smaller vessels. For simplicity, in a first step we consider an isotropic permeability tensor = ⋅ for the additional porous network. This permeability and the averaged blood viscosity are the two model parameters governing this equation together with the right hand side term This term replaces the outflow of fluid from the smaller vessels into the IF by a homogenized representation of the Starling equation (3) involving the surface-to-volume ratio of the smaller blood vessels ( ∕ ) Λ S as an additional parameter. The mass balance equation of the IF for the fully-resolved model (4) is adapted as in the homogenized formulation. Comparing the two equations, it is obvious that leakage from the large vessels is still treated equivalently, i.e., the large vessels are still embedded as 1D inclusions in the tissue with a Dirac measure (now defined only on Λ L ). By contrast, leakage from the smaller blood vessels is replaced by the homogenized mass transfer term (10) from the vascular domain Ω into the interstitium, i.e., from (9) into (11). So far, this procedure is analogous to other hybrid approaches. 21,25 The main difference to our methodology lies in the coupling between the larger vessels Λ L and the homogenized vasculature Ω . In the aforementioned publications, this was realized at the free ends of the larger vessels, i.e., as an outflow at the tips of the 1D discretization into the homogenized 3D vasculature domain. This was possible since the employed data sets had a clear vascular hierarchy with larger arterioles and venules connected to smaller capillaries. Our vascular networks, 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 including variable vessel lengths and diameters as well as dead ends. All this is typical for tumor-specific vasculature. 38,39 As shown in detail in Section 3.3 for our data and the employed methodology to distinguish between large, flow-carrying vessels and smaller ones, another approach is more sensible: We enforce the coupling between larger vessels and the homogenized vasculature along the entire 1D representation of larger vessels Λ L with a line-based coupling instead of a point-based coupling at the tips of the larger vessels flowing into the capillary bed as described before. Compared to these hybrid approaches, our proposed method has the advantage that no additional parameter -apart from the penalty parameter -is involved for the coupling of the two representations.
For that, we formulate a constraint of equal pressures in Λ L and Ω as which enforces a coupling between pressures ̂ in the one-dimensional, large vessel domain Λ L and homogenized pressures in the 3D domain Ω . We aim to reproduce the fact that the pressure in a smaller vessel branching from a larger vessel at a specific node is equal to the pressure at the same node. If this smaller vessel is homogenized and, thus, removed from the 1D representation, we want to enforce these equal pressures between the resolved part and the homogenized part of the vasculature along the 1D vessel domain Λ L . In Section 3.3, we justify formulating this constraint along the entire 1D domain Λ L considering the connectivity between larger and smaller vessels in our cases. We have previously employed a similar strategy in our hybrid treatment of the vasculature in a multiphase tumor growth model 3 and the related solid mechanics problem of beam-to-solid mesh tying. 36 We follow the same approach as in the two aforementioned publications and incorporate the constraint with an additional Lagrange multiplier (LM) field into the weak form of our hybrid model, which reads as , Therein, the first line is the weak form of flow in the larger vessels (8) which is coupled to the weak form of flow in the homogenized vasculature domain, i.e., the second line (13b) with a continuous LM field defined along the blood vessel center line. The third line is the weak form of flow in the IF. Compared to the fully-resolved model, conf. eqn. (5b), the additional mass transfer term arises due to leakage from the homogenized part of the vasculature into the IF. The fourth line represents the variational form of the coupling constraint (12). Conveniently, the LM field employed to enforce this constraint can then be interpreted as a mass transfer term from the 1D resolved bigger vessels into the 3D homogenized vasculature, i.e., =̂ → . Alternatively, a Gauss-point-to-segment scheme could also be employed but suffers from over-constraining of the system for large penalty parameters. 3,36 Spatial discretization of the weak form (13a)-(13d) leads to a saddle-point problem with nodal primary variables ̂ ∈ ℝ nodes,Λ L , ∈ ℝ nodes,Λ L , ∈ ℝ nodes,Ω and ∈ ℝ nodes,Ω , that is, nodal pressures and nodal LMs in Λ L,ℎ , nodal IF pressures in Ω ℎ and nodal blood pressures of the homogenized vasculature in Ω ,ℎ . In the following, we will specifically focus on the discretization of the terms arising due to the LM method.
Approximating those contributions with a finite element interpolation yields a mortar-type formulation where the nodal LMs are additional degrees of freedom, condensed out with a dual approach 40,41 or a penalty regularization of the mortar method is employed to remove the additional degrees of freedom and the saddle-point structure. 42 Here, we follow the latter approach just as in our previous work on 1D-3D type couplings. 3,36 The contributions to the weak form of the mass balance equations, i.e., the two last terms in (13a) and (13b) can be written as with the so-called mortar matrices and The entries of these matrices involve integrals of products of LM shape functionsΦ defined on the discretized 1D domain Λ L,ℎ with 1D shape functionŝ and with 3D shape functions defined in the 3D domain Ω . Hence, these terms are again evaluated using a segment-based approach, see Appendix A. We choose linear shape functions for both primary variables and the LM interpolation, i.e.,Φ =̂ . The weak form of the constraint (13d) may then be written in discretized form as where we have defined a weighted pressure gap at each node in Λ L,ℎ . This gap is then further used for the penalty regularization of the mortar method to explicitly define the nodal LMs in terms of 1D and 3D nodal blood pressures as Hence, the LMs are no longer independent variables in the system but depend on the primary variables ̂ and . This overcomes the two major drawbacks of the LM method, namely, the increased system size and its saddle-point structure. Depending on the penalty parameter > 0, the constraint = is relaxed and the exact solution is only recovered for → ∞. Additionally, the nodal LM in (19) has been scaled with the inverse of the diagonal matrix As proposed by Yang et al. 42 this removes the dependency of the nodal LM on its "gap", i.e., in our case it makes its entries independent of the element lengths associated with its corresponding node. This can now be used to replace the LM vector such that the matrix-vector form of our hybrid model emerges as  As in the fully-resolved model (7), main diagonal blocks are denoted as and the coupling blocks ̂ and ̂ stem again from the transvascular 1D-3D exchange term. Additionally, the coupling blocks and account for exchange between homogenized vasculature and IF. The terms involving the mortar matrices , and couple blood flow in the larger vessels with the homogenized vasculature using our mortar penalty approach. Obviously, the LMs are no longer part of the system which is, consequently, not of saddle-point type anymore. The drawback, however, is that the choice of the penalty parameter influences the accuracy with which the constraint is fulfilled. Large penalty parameters yield better accuracy in terms of constraint fulfillment but can lead to an ill-conditioning of the system matrix. We will comment on the choice of the penalty parameter in Remark 4.

Remark 1.
The concrete implementation of the hybrid model is slightly different than described here for illustrative purposes. The equations for IF flow and blood flow are evaluated simultaneously on the 3D domain and not assembled into two separate block matrices as written in (21). This means that the degrees of freedom are actually re-ordered in a node-wise manner compared to (21) such that one row corresponding to the nodal IF pressure at a node is followed by a row corresponding to the homogenized blood pressure at this node . Therefore, we actually solve a system which is blocked with 2×2 submatrices, where the upper part corresponds to the resolved part of the vasculature and the lower part to the IF and the homogenized vasculature. For this system, we again employ the AMG(BGS) preconditioner 37 with the GMRES iterative solver.

SETUP OF COMPUTATIONAL MODELS
This section describes the setup of our fully-resolved and of our hybrid model. We first analyse the real-world tumor data sets which we will employ for all our numerical tests. Subsequently, the assignment of boundary conditions in both models is described. Then, we illustrate how we create the hybrid model with homogenized vasculature starting from the full topology of the vascular networks. Finally, the definition of representative elementary volumes for homogenization is introduced.

Analysis of real-world tumor data sets
We have obtained three different vasculature data sets from REANIMATE 13,14 , which is a framework combining mathematical modeling with high-resolution imaging data to predict transport through tumors. We only briefly describe the experimental procedure here. More details are given in the two aforementioned papers. Two different colorectal cell lines, namely SW1222, LS174T, and one glioma cell line, GL261, were grown subcutaneously for 10 to 14 days in mice, resected and optically cleared, and finally imaged using optical projection tomography. The data was then segmented to obtain the complete blood vessel networks inside the tumors in the graph format as discussed before.
The topologies and blood vessel radii of the three distinct cases are illustrated in Figure 2 together with representative results of blood vessel and IF pressure of the fully-resolved model. Further network data has been collected in Table 1: All three networks contain more than 100 000 blood vessel segments and nodes. The SW1222 case is the biggest tumor both in network

FIGURE 2
Full topology and structure of the vascular networks (left, colour-coded by the respective radii), representative results for simulated blood pressures (middle) and IF pressures (right) in the fully-resolved model (Same spatial scale is used for all three cases) size and tissue volume. The latter has been calculated by approximating the hull of the tumor using the alphaShape function of Matlab (MathWorks Inc., Natick, MA) 43 . The hull is then smoothed, remeshed using Gmsh (version 4.4.1) 44 and slightly enlarged to encompass all vasculature nodes. Its enclosed region is integrated to give the tumor volume, see also Figure 3 for the SW1222 tumor. Note that this tumor domain corresponds to the domain Ω on which the additional porous network of smaller vessels is present in the hybrid model and where its additional governing equation (9) is defined and solved. Furthermore, all topologies are rotated such that their principal axes align with the coordinate axes. The three different cases show distinct vascular architectures, for instance, the SW1222 network is much denser with a higher blood vessel volume fraction and blood vessel surface-to-volume ratio than the two other types. In addition, its blood vessel diameters are generally larger and have a much higher variability. Finally, we have analyzed the boundary nodes, i.e., the tips which are only connected to one other node. All topologies have a comparable number of boundary nodes lying on the aforementioned enclosing alpha shape whereas the GL261 and SW1222 tumors have a much higher number of tips inside the domain than the LS174T tumor.

Assignment of boundary conditions
The assignment of physiologically reasonable boundary conditions on large vascular networks is quite challenging since flows or pressures cannot be measured on the level of individual micro-vessels. Sweeney et al. 14 developed an algorithm 45 to apply boundary conditions which match in-vivo measurements of perfusion for the present data set. We will re-use this framework here to generate the boundary conditions for the fully-resolved case and briefly describe it in Section 3.2.1. Boundary conditions for the hybrid model are detailed in Section 3.2.2.

Fully-resolved model
For the fully-resolved model, boundary conditions for the blood pressure ̂ in the 1D network and the IF pressure need to be assigned. For the blood vessel pressure, we re-use the approach of Sweeney et al., which has been made publicly available 45 and is based on earlier work by Fry et al. 24 Thereby, boundary conditions are assigned on the tips of the network, i.e., on the boundary nodes of the 1D representation of the vasculature both on the tumor hull and inside the tumor as given in Table 1.
The following algorithm is applied: First, a high or low pressure of 5999.4 Pa or 1999.8 Pa (corresponding to 45 mmHg or 15 mmHg) is randomly applied to the boundary points on the tumor surface until 5 % of all end points of the 1D network have been assigned such a high/low pressure. Additionally, the method prevents that points which are in close proximity to each other are assigned the far apart pressure values which would produce unphysiologically large flows. Second, a no-flux boundary condition is randomly assigned to the interior boundary nodes until 33 % of all boundary nodes have this condition. This value is consistent with the fraction of dead ends in tumor vasculature estimated from experimental studies. 46 Third, the remaining 62 % of boundary nodes are given as unknowns to the flow optimization scheme of Fry et al. 24 This scheme aims at solving a constrained optimization problem for incomplete pressure boundary data by minimizing the error of pressures and wall shear stress w.r.t. pre-defined target values. D'Esposito et al. 13 and Sweeney et al. 14 have shown that this procedure for assignment of boundary conditions ensures that tumor perfusion is in good agreement with experimental data. Note that the entire algorithm is not deterministic due to the random selection of nodes for high/low boundary conditions on the external surface of the tumor and of nodes for no-flux boundary condition in its interior. Therefore, the analyses in the following sections will be performed on five different sets of pressure boundary conditions on the 1D network per tumor case.
Concerning the IF pressure, we want to prescribe the far-field pressure for the IF as ∞ = 0 Pa following Sweeney et al. 14 In order to achieve this within our finite element approach, we enlarge the domain Ω radially to a sphere of radius 80 mm as shown in Figure 3 for the SW1222 case. This allows us to set a Dirichlet boundary condition of = 0 Pa on its boundary Ω and, thereby, to mock the far-field pressure. We validated this approach in the following way for all three vascular networks: We solved the fully-resolved model and compared the IF pressure solution (for one specific set of pressure boundary conditions on the 1D network) with a case where the domain was only enlarged to a sphere with radius 40 mm (with corresponding zero IF pressure Dirichlet boundary condition on its outer surface). No visible differences in the IF pressure distribution in our domain of interest inside and around the tumor domain could be detected. This indicates that the enlargement is big enough insofar as the solution in the domain of interest is not influenced by the size of the enlargement any more. We can also gradually coarsen the mesh when moving away further from the vascular domain as depicted in Figure 3 since the IF pressure gradient flattens and tends to zero further away from the center of the domain. This enables the use of a sufficiently fine mesh for the region surrounding the embedded vascular network while the computational cost for extending the domain is not too high.

Hybrid model
In addition to boundary conditions for the IF pressure and the blood pressure ̂ , the hybrid model requires boundary conditions for the pressure in the homogenized vasculature . The IF pressure is treated as in the fully-resolved model and we set it to zero at the boundary of the domain Ω. In the following, we will always compare the accuracy of the hybrid variant w.r.t. the fully-resolved one for one specific set of pressure boundary conditions on the 1D network obtained with the procedure described in the previous sectoin. Thus, to perform this comparison the pressure boundary conditions on the 1D network are transferred from the fully-resolved model to the hybrid model in the following manner: The boundary conditions of blood pressure ̂ on the larger vessels Λ L can directly be taken from the boundary conditions of the fully-resolved model. If a node with a Dirichlet boundary condition in the fully-resolved vasculature Λ is part of the larger vessels Λ L we simply keep this boundary condition on the 1D discretization also in the hybrid model. Dirichlet boundary conditions on the smaller vessels Λ S cannot be assigned on the 1D discretization since smaller vessels are homogenized. However, we can employ them to assign boundary conditions for on the boundary of the domain of homogenized vessels Ω as depicted in Figures 1b and 3. Similar to Vidotto et al. 25 , we smooth these values to account for the homogenization of the smaller vessels: Each condition belonging to a node of the smaller vessels Λ S at the tumor surface is assigned to all 3D nodes lying on the surface Ω within a distance of less than 400 µm for the SW1222 and the LS174T tumor and less than 200 µm for the GL261 tumor. Nodes of the 3D mesh which lie within this distance of multiple boundary nodes on Λ S are assigned the mean pressure value of all these boundary nodes. On the rest of the surface Ω we set a no-flux boundary condition. We also do not set a boundary condition for on nodes of the 3D mesh in

TABLE 2
Comparison of mean blood vessel radius in larger vessels Λ L and smaller vessels Λ S (all values indicate the mean taken over five different sets 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) close proximity to end nodes of the 1D network since this would mean setting different boundary conditions on nodes whose pressures should be coupled due to the constraint on pressures ̂ and and, thus, would lead to an overconstrainment of the system. The resulting distribution of boundary conditions over Ω is illustrated in Figure 4 for three exemplary cases.

Distinction between fully-resolved and hybrid model
As previously stated, we envision that our hybrid model could be applied in cases where the full structure of the vascular network is unknown such that only the topology of the larger vessels can be acquired via non-invasive imaging. However, in our data sets we actually have the full structure available. In line with the main goals of this paper, namely, to validate the hybrid approach, to quantify the error with respect to the fully-resolved case and to determine its optimal parameters for perfusion through solid tumors, we artificially create the hybrid model from the fully-resolved one. In the hybrid approaches of 21,25 this was realized by a radius-based criterion. Their employed data sets had a clear hierarchy typical for the microcirculation with larger arterioles branching into smaller capillaries which in turn connect and form larger venules. Thus, it was possible to exploit the hierarchical structure of the vasculature by keeping only the larger vessels in the set L . For our tumor vasculature data sets this is not as straightforward. While there are some thicker vascular branches, especially in the SW1222 case, no clear hierarchical vascular architecture can be extracted from the topologies in Figure 2 with a radiusbased criterion. To illustrate this fact, we compare the full architecture of the SW1222 network with a network where only the top 10 % of vessels with the largest radius are kept in Figures 5a and 5c. Many small unconnected clusters of several blood vessel segments appear due to the heterogeneous, extremely variable distribution of the radius and lack of vascular hierarchy. Branches connecting these clusters which have a smaller radius are removed. Applying our or any hybrid model on this topology would not be possible as hybrid approaches also rely on a "sensible" topology for Λ L which preserves the structure of the entire network via one or several connected subgraphs of larger vessels which feed respectively drain the smaller, homogenized vessels. Only then, the 1D blood flow model and corresponding boundary conditions can reasonably be applied on Λ L together with suitable exchange terms into the smaller vessels. Thus, we instead distinguish between smaller and larger vessels based on the flow within the vessels. This yields a better preservation of the network architecture for the hybrid case, see Figures 5a and 5b. Now, connected subgraphs of larger vessels Λ L emerge which connect in-and outlets of the main flow-carrying vessels with the smaller vessels.
Hence, our strategy to obtain Λ L is as follows: We first solve the fully-resolved model (using the boundary conditions described in Section 3.2.1). Then, all elements except the ones with the highest flow are deleted from the vascular graph, e.g., the top 10 % with the highest flow are kept. However, there are still some very small clusters consisting of only a few segments present in the graph. We additionally delete connected components from the graph whose overall length is smaller than 250 µm, i.e., sub-components which are smaller than ten segments with the average segment lengths given in Table 1. By that, we delete an additional 0.1 − 0.8 % of segments which are part of these smaller sub-components. This methodology gives the set L of larger vascular branches which are kept in the hybrid approach as exemplarily shown in Figure 5b. Here, we show only the SW1222 case but equivalent results hold for the other two network topologies. In the following, we will denote cases where the top X % of elements with highest flow are kept and the small connected components are removed according to the procedure described above as "case X %".
Recall that the assignment of pressure boundary conditions on the fully-resolved vascular network is not deterministic. Moreover, different boundary conditions will produce distinct flow patterns in the vasculature and, hence, also different sets of large and small vessels in our procedure and a different topology for Λ L . Therefore, the following analysis will always be performed for five sets of pressure boundary conditions on the 1D network with corresponding distinct sets of large and small vessels L and S .
In Table 2 the mean diameters Λ L and Λ S of larger and smaller vessels are compared. It is obvious that the diameters in the set of small vessels S which are removed from the hybrid model are considerably smaller than the diameters of the large vessels. This behaviour is most pronounced for the SW1222 topology where for the case 5 % the mean diameters in Λ L are 2.5 times bigger than in Λ S . Naturally, this ratio drops for all topologies when a higher percentage of segments is kept in the large vessel set. For the LS174T and GL261 data sets the difference in blood vessel diameters is not as large but this can be attributed to the fact that the diameters are less dispersed than in the SW1222 topology, see also the mean and standard deviation of the diameters in Table 1. Also in these cases, the diameters in Λ L are larger by approximately one standard deviation of the diameter of the entire vasculature (as in the SW1222 case). In summary, our approach incorporates mainly the vessels with larger radii in the set L whereas also some segments with smaller radii are kept to preserve the main topology of the networks in the hybrid model. Therefore, there is also a significant congruence of the sets of large vessels L between different pressure boundary condition cases. For instance, in the case where 10 % of the blood vessels are kept in the hybrid model, the average percentage of identical retained segments between two different pressure boundary condition cases is 45 % for the LS174T tumor, 51 % for the GL261 tumor and 78 % for the SW1222 tumor. In Remark 3 we further comment on how the obtained topologies of larger vessels Λ L relate to real in-vivo tumor imaging data.
Next, we justify our line-based coupling approach between the large vessels Λ L and the homogenized vasculature. For that, we have analyzed the connectivity between larger and smaller vessels for the fully-resolved topologies in Table 3. Here, we denote by = nodes,Λ L ∩Λ S ∕ nodes,Λ L the proportion of nodes of the larger vessels Λ L which have a direct connection to a node of the smaller vessels Λ S . The presented data illustrates that for the GL261 and the SW1222 tumor almost every third to every fourth node of the main branches Λ L is directly connected to a node of the smaller blood vessel segments Λ S , i.e., at every third to fourth node a smaller vessel branches away from Λ L . For the LS174T network, the connectivity is slightly smaller. Here, only 13 − 18 % of nodes in larger vessels are connected to smaller vessels. In all cases, these numbers obviously again drop when keeping a larger portion of the entire network in the set L .
In the hybrid approach, information about these smaller branching vessels is lost since they are removed from the 1D representation of the vasculature. As stated above, we want to enforce equal pressures between larger and smaller vessels as this equality also holds in the fully-resolved model at branching points. The high connectivity between the two network parts supports our line-based mortar penalty coupling between the resolved and homogenized part of the vasculature in which we actually couple the entire network of big vessels Λ L with the homogenized vasculature. Of course, we know the connecting nodes between larger and smaller vessels here as we know the full topology of all networks, so we could also enforce the coupling between resolved and homogenized part in a point-based manner at these locations. However, in the more realistic case when only the architecture of larger vessels is known without the exact locations where smaller vessels branch away, this is not the case. Therefore, we  TABLE 3 Analysis of connectivity between fully-resolved and homogenized part of vasculature: is the fraction of nodes of larger vessels with a direct connection to smaller vessels, and | | are measures of the variability of the diameter and flow, respectively, in the segments connecting larger and smaller vessels (data includes the mean taken over five different sets 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) adopt our line-based coupling within the hybrid model hereafter to compare the results with the fully-resolved reference solution. Note that the network tips of Λ L (both in the interior of the domain and on the tumor hull) are actually also coupled with the homogenized vasculature since the discrete constraint of a vanishing weighted pressure gap is enforced along the entire 1D discretization and, thus, also at the end nodes.
Finally, we analyze also the elements connecting larger and smaller vessels, i.e., those 1D elements of the smaller vessels where one node is part of Λ L and the other part of Λ S . We gather all these elements and compute their mean diameter and mean absolute flow value. Then, we calculate the coefficient of variation of these quantities, and | | as the ratio of standard deviation of the diameter, resp. flow to its mean in these connectivity elements. The results are collected in Table 3. Obviously, the SW1222 case shows the highest variability in both flow and diameter followed by the GL261 and the LS174T case. For all cases, the variability of the flow is larger than for the diameter since the volumetric flow in a segment depends on the fourth power of the diameter due to the employed Hagen-Poiseuille relationship. These results are consistent with the topology of the entire network where the variability of the blood vessel diameter is also larger for the SW1222 tumor than the GL261 and the LS174T tumor, see Table 1. In Section 4.3 we will show that this higher variability makes it harder to match the flow from large into small vessels between the two models.

Remark 2.
We believe that our hybrid approach is also applicable to more organized, hierarchical networks as, for example, the topology used by Vidotto et al. 25 In this publication the network was partitioned by a radius-based threshold, see Figure  1 therein. The larger vascular structures contain very short branches going away from the main vessels. At the tips of these short branches, the node-based coupling is performed. If one instead removed these very short branches and left only the major, flow-carrying vessels in Λ L , a line-based coupling along these vessels could again be implemented.

Determination of representative elementary volume size
The existence of a resentative elementary volume (REV) is an important concept for different homogenization procedures. 47 In general, such a volume should be big enough to smooth out fluctuations of spatial heterogeneities yet small enough to resolve the physical effects of interest. In this section, we investigate the choice of REVs in the context of our model and the employed data sets. Naturally, we will investigate the properties of the smaller vessels Λ S in the following since this is the part of the vasculature which is homogenized and treated as a porous continuum in the hybrid approach. Furthermore, five different sets of pressure boundary conditions on the 1D network are studied. This is again due to the fact that different pressure boundary conditions on the 1D network will lead to different flow patterns in the vascular network and, therefore, also different sets L and S of large and small vessels (potentially with a different distribution throughout the domain) with the employed flow-based criterion.
For this purpose, we have devised the following procedure: 1. For each network topology we create five different cases with a different set of pressure boundary conditions on the 1D network for the fully-resolved model as described in Section 3.2.1.

2.
We partition all cases into the distinct sets of large and small vessels as described in Section 3.3. We here investigate the case 10 % for all different topologies but equivalent results have been obtained for leaving the top 5 %, 15 % or 20 % of vessels with the largest flow in the system.

The size of the cubes is successively increased in all coordinate directions by
while keeping their centers fixed. The blood vessel volume fraction Λ S and the surface-to-volume ratio ( ∕ ) Λ S of smaller blood vessels Λ S is computed for each cube at each size. If a cube protrudes from the domain during this enlargement, these quantities are calculated on the intersection of the cube with the domain Ω .
Per case with different boundary conditions this is performed for ten randomly generated cube centers. The results are shown in Figure 6 for only three REVs per boundary condition case, that is, in total 15 cases to not clutter the plots. The evolution of the blood vessel volume fraction Λ S and of the surface-to-volume ratio ( ∕ ) Λ S of the smaller blood vessels Λ S for increasing the edge length of the cubes is illustrated. Therein, we denote the length scale as = 3 √ ∩ Ω to account for cases when a larger cube protrudes from the domain Ω . All three topologies exhibit similar features: Λ S and ( ∕ ) Λ S fluctuate strongly for smaller Splitting the domains into these REVs of equal size is not an easy task due to their irregular, elliptic shape. We first defined a regular grid of REV centers and performed an initial Voronoi tesselation based on this grid. Due to the shape of the domain this resulted in too small or too large REVs. Therefore, we performed an optimization of the Voronoi tesselation where the objective function had the goal to define REVs of equal volume and equal dimensions. The resulting REVs are visualized in Figure 7. The mean deviation of the REVs from the previously determined volume and lengths from Figure 6 in all three coordinate directions is less than 5 % for the domains of all three tumor types.
Finally, we employ these REVs to study the distribution of blood vessels inside the domain. For that, we define the nondimensionalized radial distance of each REṼ REV as the distance of the center of the REV to the center of the domain divided by the distance of the center of the domain to the tumor hull in direction of the center of the REV. Again, this analysis is performed for all three tumor types for five different sets of pressure boundary conditions on the 1D network since those influence the flow in the 1D vasculature and, consequently, also the composition of Λ L and Λ S as previously mentioned. The results for the volume fraction of big vessels Λ L , small vessels Λ S and the entire vasculature Λ are shown in Figure 8. The clearest structure is evident for the the SW1222 case: towards the tumor hull, Λ S and Λ L and, thus, also the sum of the former two, Λ , gradually increase. Close to the center of the domain, there is still a significant amount of smaller blood vessels while almost no larger blood vessels are present. This is consistent with experimental data showing higher blood vessel density and perfusion in the tumor periphery 13,48 with only a few major vessels penetrating into the center of the tumor. 49 These trends are also present in the LS174T tumor, albeit, far less pronounced than for the SW1222 tumor. By contrast, the GL261 vascular network shows a completely different behaviour. While the vascular density of large vessels remains almost constant over the tumor radius, the one of the smaller blood vessels Λ S drops and, thus, also the overall volume fraction Λ . Remark 3. The validity of the obtained topologies and distributions for Λ L and Λ S and the applicability of the proposed hybrid approach is supported by state-of-the art optoacoustic in-vivo imaging techniques. 23 The currently attainable spatial resolution is less than 50 µm throughout the tumor domain which is in the range of the diameter of larger vessels from Table 2. Furthermore, the larger vessels which are retained in the hybrid model are more concentrated at the tumor periphery (at least for the SW1222 and LS174T case) and are, thus, more accessible to imaging. Qualitatively, the topology of the larger vessels from Figure 5b is in good agreement with corresponding optoacoustic imaging data from tumors 23 where larger feeding vessels are visible at the tumor rim. From these experiments, one can extract a similar topology of larger vessels Λ L to apply our hybrid model. A connected subgraph of larger peripheral feeder blood vessels as in Figure 5b rather than single clusters as in Figure 5c is attainable. Hence, we conclude that the employed methodology of splitting into larger and smaller vessels yields a valid scenario resembling real experimental data and can, therefore, be used to investigate our hybrid embedded/homogenized approach for solid tumor perfusion.

NUMERICAL EXPERIMENTS
In this section we perform several numerical experiments to evaluate the performance of the hybrid model in comparison to the fully-resolved one. We first define a comparison metric and optimize the parameters of the hybrid model such that the best possible correspondence between the models is achieved according to this metric. Subsequently, we study several other quantities to compare the two models and present a further improvement of the hybrid model via a vascular volume fraction dependent permeability for the homogenized vessels. The tumor hull is smoothed and triangulated using Gmsh (version 4.4.1) 44 and its enclosed volume is meshed with linear tetrahedral elements using Trelis 17.0. 50 An exemplary mesh for the SW1222 topology is shown in Figure 3 and parameters of the 3D mesh are given in Table 1. Note that the 3D mesh is completely independent of the discretization of the 1D networks, that is, the nodes of the two meshes do not match which is an advantageous feature provided by our recently introduced hybrid approach. 3 Both the fully-resolved and the hybrid FEM model have been implemented in the in-house parallel multi-physics research code BACI. 51 Parameters for both models are listed in Table 4.

Remark 4.
In preliminary simulations, we determined the proper range for the penalty parameter . As a compromise between accuracy and a well-conditioned system matrix, we defined the following criterion: This rule states that the mean relative pressure error in terms of the length-independent nodal pressure difference vector −1 and the nodal pressure ̂ in the 1D network is less than 1 %. The values for all cases are given in Table 4. Note that the penalty parameter has units of [length] 2 ∕[time ⋅ pressure] such that the LM field represents a 1D-3D mass transfer term, or volumetric flow per length. This allows interpreting the penalty parameter as very large permeability governing the mass transfer between resolved and homogenized vasculature in the hybrid model.

Remark 5.
In our opinion, the main advantage of the hybrid model is not a reduction of computational cost compared to the full model, but the fact that it relies only on data available through non-invasive imaging. Nevertheless, we also did a first preliminary evaluation comparing the computational costs of the two models and found that the hybrid model was not significantly faster than the fully-resolved one and in some cases even slower. The effort for finding 1D-3D elements interacting with each other, building the integration segments and evaluating the coupling terms along the 1D vasculature is obviously smaller for the hybrid model since less 1D vessels are present. However, this is balanced or even outweighed by its increased effort in several other aspects: The evaluation of the 3D elements is more costly since two equations per node (in Ω ) have to be evaluated, the system size, which is dominated by the number of 3D nodes, and, thus, the linear solver time is increased and the condition of the system is worse compared to the fully-resolved case due to the penalty approach, which in turn raises the linear solver time. However, for all our studies we used the same 3D meshes for both hybrid and full model. The cost for the hybrid model could be greatly reduced by employing a coarser 3D mesh. Vidotto et al. 25 showed that this still gave acceptable results in terms of REV pressures for their approach. (3),(10) Hydraulic conductivity of vasculature ∕ see Table 5 µm 2 Pa −1 s −1 -(10) Surface-to-volume ratio for transvascular flow ( ∕ ) Λ S see Table 5 µm −1 - [a] The value for blood viscosity is calculated separately in each 1D element using the algebraic relationship of Pries and Secomb 33 with hematocrit value fixed to 0.45.

Definition of metric for comparison of the two models
To assess the performance of the hybrid model in predicting microvascular flow and IF pressure inside solid tumors in comparison with the fully-resolved model, a suitable metric is warranted. Ideally, the hybrid model should match the fully-resolved one in terms of blood and IF pressure as well as blood and IF flow to obtain an accurate representation of the perfusion through the tumor. Therefore, we define our metric as a combination of these quantities. The first contribution is the correspondence 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 determination 2 in terms of nodal blood pressures in the large vessels between the two models as Next, we will detail how we calculate the flows in the REVs in both models. In the center of each REV we define a square □ with dimensions REV × REV such that its normal is aligned with coordinate direction . The volumetric flow in the homogenized part of the vasculature in coordinate direction is then given by as the surface integral of the flux through the square. For the fully-resolved model, we define it as which is the sum of the volumetric flow of all segments which are part of the smaller vessels and cut by the square □ . Therein, is the tangential vector of a segment pointing from its first to its second node and sgn (⋅) denotes the sign function. Finally, we define the total coefficient of determination between the two models as the sum of the contributions from blood pressure in large vessels (23), blood pressure in small vessels (24), IF pressure (25) and flow in small vessels (26) as This metric, where all four contributions are weighted equally, will be employed to study the accuracy of the hybrid model w.r.t. the full model and to find the optimal parameters of the hybrid model.

Optimization of parameters of the hybrid model
Compared to the fully-resolved model, the hybrid one has two additional parameters, which are the hydraulic conductivity of the homogenized vessels ∕ in (9) governing blood flow and the surface-to-volume ratio ( ∕ ) Λ S accounting for transvascular flow from the homogenized vessels into the IF in (10). Our goal in this section is to determine these parameters such that the agreement in terms of blood flow and blood and IF pressures of the hybrid model with the fully-resolved model is maximized. For that we employ the total coefficient of variation (29) between the two models deduced in the previous section and formulate the following optimization problem in terms of the parameters of the hybrid model: that is, we aim to find the parameters of the hybrid model, for which the correspondence of the two models 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 SciPy package (version 1.5.2) 54 and interfaced it to the software framework QUEENS. 55 Internally, SciPy employs the Levenberg-Marquardt algorithm to solve the nonlinear leastsquares problem (30). Derivatives of the metric (29) w.r.t. the parameters are approximated using forward finite differences. This implies that the hybrid model has to be solved three times per iteration step. In preliminary simulations, we confirmed that different initial conditions (from a sensible parameter range) converged to the same optimum.
Since the full topology of the vasculature is available, we could also obtain these parameters by a suitable homogenization procedure for the permeability as previously done for other hybrid or continuum models. 18,19,20,21,25 We did not follow this approach here for the following reasons: First, the chaotic structure of the blood vessel network implying also a very chaotic blood flow pattern typical for the solid tumors would make this very challenging. Second, we want to create a best-case scenario by fitting the parameters of the hybrid model such that possible errors introduced by a homogenization scheme are minimal.
The general algorithm can be described as follows: 1. Obtain a set of boundary conditions for the full model as described in Section 3.2.1 and solve the full model to generate a reference solution.
2. Extract the topology of larger vessels for the hybrid model from the full model, conf. Section 3.3, and apply boundary conditions on the hybrid model, conf. Section 3.2.2.
3. Find the optimal parameters of the hybrid model by maximizing the total coefficient of variation (30). During the optimization procedure repeated evaluations of the hybrid model with different parameters are required.
Representative results of the optimization scheme are depicted in Figure 9 for all four contributions to the total coefficient of determination. Very good agreement between the two models in terms of nodal pressures in the larger vessels ̂ can be observed in Figure 9a. This can be expected because the same boundary conditions on the large vessels are applied in both cases. Thus, large and small pressure values show very good agreement, further away from these boundary conditions in the medium pressure range, deviations become larger. The clusters with the largest errors are separate branches which are not directly connected to nodes of the 1D network carrying boundary conditions. The correspondence for the nodal IF pressures in Figure 9c is also very good. For low IF pressures this is again due to the zero pressure boundary condition assigned on Ω for both cases, but also for higher IF pressures inside the tumor, which is the actual domain of interest, the pressure differences are very small, in this case, the maximum absolute error is 237.1 Pa corresponding to a maximum relative error of 8.4 %. The pressure in the smaller vessels, resp. the homogenized vasculature in the hybrid model, exhibits larger errors, see Figure 9b. Overall, the agreement is still reasonable. We found that the error is largest for branches ending in tips with boundary conditions on the 1D vasculature either inside the domain or on the tumor hull. For instance, this is the case for the larger errors around ̂ | | f ull ≈ 3300 Pa. The boundary conditions on these tips inside the domain are not retained in the hybrid model and for the tips on the tumor hull, they are smeared over several 3D nodes as described in Section 3.2.2. Hence, while the error in the medium pressure range is distributed symmetrically, larger deviations at both ends of the pressure spectrum towards the smeared values are present. This error due to point-wise non-matching boundary conditions can also not be improved by the optimization of the parameters. However, in Section 4.3, we will show that averaged REV pressures of both models are in very good accordance. We believe that this is a more interpretable and fairer comparison metric as the hybrid model cannot be expected to exactly match the pressure distribution of the fully-resolved one (in particular on the boundary) since the information about the exact topology of the smaller vessels is not represented. Finally, the results for the flow in the smaller vessels are shown in Figure 9d. Here, the poorest agreement of the two models is present, especially, larger flows are not met properly.

FIGURE 9
Exemplary comparison of hybrid model (with optimized parameters) with fully-resolved model for one specific network topology (SW1222, 10 % of 1D blood vessels have been retained in hybrid model). In each subfigure, solution of hybrid model is plotted over solution of fully-resolved model and the dashed line indicates perfect agreement between the models with 1:1-correspondence. Comparison of blood pressure in large vessels is depicted 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 coefficient of determination calculated according to (29) is 2 tot = 0.716 for this case.
Further results for all cases have been collected in Table 5. For each tumor network, we generate five different sets of pressure boundary conditions on the 1D network from which different flow patterns and, therefore, also different sets of larger and smaller blood vessels emerge as discussed in Section 3.3. Then, we investigate different cases, where 5 − 20 % of the larger vessels are kept in the hybrid model. We found that five sets of pressure boundary conditions on the 1D network were enough to study our hybrid model since randomly picking only four out of the five boundary condition cases changed the mean result by at most 8 %. Moreover, taking the mean parameter of a case X % over all different boundary condition cases instead of the optimal value for each specific case only changed the total coefficient of determination by less than 2 %. Furthermore, we compare the result of the optimization procedure for ( ∕ ) Λ S with the calculated surface-to-volume ratio of the smaller vessels for each case. The relative error ∕ is smaller than 5 % for all cases validating that the optimization procedure converges to a physically reasonable result. The permeability is largest for the SW1222 topology which can be expected considering the much denser network of this case. For all tumors, it decreases if a larger proportion of the 1D vessels is kept in the model, which is also sensible since the smaller the proportion of homogenized vessels, the less permeable these vessels.
As already described above, all cases exhibit a very good correspondence in terms of blood pressures in larger vessels and IF pressures, proven by the values for 2 L and 2 IF in Table 5. If the fidelity of the hybrid model is increased by resolving a

TABLE 5
Results of the optimization procedure for hydraulic conductivity and surface-to-volume ratio of homogenized vasculature in the hybrid model. Relative error w.r.t. calculated surface-to-volume ratio and 2 -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 coefficient of determination 2 tot between fully-resolved and hybrid model is calculated according to (29). (All data includes the mean taken over five different sets 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) larger proportion of the network structure, the agreement between the two models grows likewise. This is also the case for the coefficient of determination of blood pressure in smaller vessels 2 S . Here, the SW1222 and the LS174T case exhibit comparable accuracy whereas the GL261 case experiences a larger discrepancy. We can attribute this to the fact that this topology has the largest number of tips at the tumor hull and also the largest number of dead ends considering that it is the smallest data set, see Table 1. Hence, the pressure error is largest due to non-matching boundary conditions between fully-resolved and hybrid model as mentioned above. However, we will show in Section 4.3 that in terms of REV blood pressures its conformity with the hybrid model is as good as the other cases. The difference in flow in the small vessels is the largest source of error in all cases. Also taking more 1D vessels into account for the hybrid model does not necessarily improve the behaviour. We believe that this is due to the chaotic flow patterns in the smaller vessels and to the fact that we define the permeability tensor as isotropic and constant over the entire domain Ω . Apparently, this is insufficient to resolve the flow in the homogenized vasculature in comparison to the full model. We tried to increase the agreement by giving a higher weight to the coefficient of determination of flow in the smaller vessels 2 f low,Λ S Additional error measures for the agreement of both models. Shown are the absolute and relative error of the hybrid approach in terms of mean REV blood pressure in smaller vessels and mean REV interstitial fluid pressure and the 2 -values for agreement between both variants in terms of flow from large to small vessels and flow in the entire vasculature. (All data includes the mean taken over five different sets 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) Furthermore, we denote by (⋅) the mean value of these error measures over all REV REVs. Note also that both the mean REV blood and IF pressure vary considerably between different REVs. The pressure difference between single REVs varies in a range of 800 − 1200 Pa for the IF and a range of 800 − 2000 Pa for blood. The data of this analysis is collected in Table 6. Overall, a remarkable agreement of the mean REV pressures for both blood and IF can be observed in all cases. As in Table 5, the SW1222 tumor has the best agreement, but also the GL261 case which previously showed the biggest nodal blood pressure errors in the homogenized vessels is very accurate in terms of mean REV blood pressure. As described above, the error is located mainly on the tips of the smaller vessels, of which the GL261 has the most compared to its network size. Nevertheless, the average blood pressure in the REVs is still matched very well for this and all other cases even though locally the pressure error is larger. We can expect that these small-scale spatial fluctuations of blood pressures cannot be represented correctly in the homogenized vessels of the hybrid model while macroscopically the average REV pressures show good agreement. Anticipating a validation with experimental data, it is anyhow not possible to perform a point-wise comparison of (blood and IF) pressures such that the average REV pressure is the more relevant and meaningful metric. Additionally, we investigate the volumetric flow between large and small vessels and compare the results of both models.
In the hybrid model, the flow between large and small vessels is given by the LM field =̂ → interpreted as a mass transfer term, or volumetric flow per length, as detailed in Section 2.3. Note that this can represent both a flow from large 1D vessels into the homogenized vasculature if locally the pressure in the resolved vasculature is bigger than the homogenized pressure or, vice versa, a flow from the homogenized vasculature into the larger vessels if the homogenized pressure is bigger than the blood pressure in the 1D vasculature. Consequently, for each REV the flow between the two compartments is given by the integral of the LM field along the part of the larger vessels Λ L ∩ REV inside the specific REV , or In the full model we can directly evaluate the mass transfer between large and small vessels inside the connecting elements of both sets, which are those elements of the smaller vessels where one node is part of Λ L and the other part of Λ S . Assuming that the first node is part of the larger vessels and the second one part of the smaller vessels, the flow between large and small vessels in the -th REV is given bŷ as the sum of the volumetric flows in the elements connecting large and small vessels which lie inside the specific REV . The number of these elements is denoted by ele,Λ L →Λ S ∩REV in the previous equation. To compare the mass transfer between large and small vessels in both models, we again define a coefficient of determination as with the respective mass transfer terms for the hybrid and the full model for each REV. Again, (⋅) denotes the mean of the mass transfer between large and small vessels of the full model over all REV REVs. The reference solution of the fully resolved model for this volumetric flow per REV varies considerably between the different REVs and both positive values, representing an overall outflow from the larger vessels into the smaller vessels in a specific REV, and negative values, representing an overall inflow into the larger vessels from the smaller vessels in a specific REV, are present. This indicates that in-or outflow from larger to smaller vessels is indeed a meaningful quantity describing the spatially varying flow patterns inside the vascular network. To reproduce this behaviour in the hybrid model variant, a good agreement with the reference solution is desirable. The results for the coefficient of determination 2 f low,Λ L →Λ S are again assembled in Table 6. The LS174T case shows the best agreement with the fully-resolved model while the SW1222 case delivers the worst results. We believe that this can be attributed to the much higher dispersion of the diameter and, thus, also the flow in the connectivity elements, which we have already studied by the coefficient of variability in Table 3. The LS174T case, which has the least dispersed distribution of both values, performs best in matching the flow between larger and smaller vessels in the hybrid model. There is a small decline of the agreement for higher percentages of retained vessels in all cases. However, the flow between large and small vessels is not included in the parameter optimization procedure. Hence, we assume that the better performance in terms of the other quantities is at the expense of this metric.
Finally, we study the correspondence between the two models in terms of the blood flow in the entire vasculature Λ. Previously, in Table 5 only flow in the smaller vessels Λ S , respectively, the homogenized vasculature was investigated. For the full model, the total flow in Λ in each REV in coordinate direction is calculated as in (28), but now both large and small vessels are taken into account. For the hybrid model, the total flow can be obtained as the sum of the flow in the homogenized vessels as given by (27) and the flow in the larger, resolved vessels, that is, eqn. (28) evaluated for the larger vessels of the hybrid model. The two quantities are compared in Table 6 defining a coefficient of determination for flow in the entire vasculature 2 f low,Λ as in (26). Evidently, the agreement between the two models is very good and much better than the previously reported agreement of flow in the smaller vessels 2 f low,Λ S . This is due to the fact that, as expected, flow is dominated by the larger vessels, the values calculated for flow in the entire vasculature are one to two orders of magnitude larger than the in the small vessels depending on the investigated case. As we are able to match the pressure in the large vessels very well and, thus, also the flow therein, very good accordance can be achieved for the total flow in both small and big vessels. As flow in the big vessels is decisive for the overall perfusion of the domain and could also be more easily acquired with experiments for further validation this is an encouraging result for the applicability of the hybrid approach. Nevertheless, we demonstrate how to enhance the correspondence of the hybrid model also in terms of flow in the smaller vessels in the following section.

IMPROVEMENTS FOR THE HYBRID MODEL
In this section, we discuss some possible improvements for the hybrid model and implement one of them. The most straightforward one would be to define the permeability of the homogenized vessels not as a constant over the entire domain Ω but per REV. Instead of an isotropic permeability tensor, one could easily integrate anisotropic effects based on the blood vessel structure inside each REV. Both has been done in the hybrid model of Vidotto et al. 25 , where a diagonal permeability tensor with different permeabilities in all three coordinate directions was employed. This could potentially augment the agreement in terms of mass fluxes in the homogenized vasculature, which is the main source of error in the hybrid model. However, we did not integrate this into our optimization procedure since we believe that this would result in overfitting of the chaotic flow in the tumor such that we would meet every single boundary condition case very well but with largely different results for the permeability tensors between the cases with distinct flow patterns. With a single scalar permeability the results for the permeability between different boundary condition cases did not fluctuate greatly. Moreover, in a real-world case where only the architecture of the larger vessels is known, it seems unrealistic to deduce the entire permeability field from the limited amount of information.
Instead we tried to enhance the model by taking information of volume fractions of the smaller vessels into account. Our rationale behind this approach is that while the complete structure of the smaller vessels cannot be obtained non-invasively, regions with higher or smaller microvascular density of small vessels could still be identified. This information could then be employed to enrich the hybrid model. The overall trend we observed in Table 5 is that the higher the volume fraction of the homogenized vessels, the larger their permeability. It also reasonable to assume that areas with a higher vascular volume fraction are more permeable to blood flow. Therefore, we investigated the relationship of the volume fraction of smaller vessels Λ S in each REV on the perfusion of the smaller blood vessels in the full model. Results are shown in Figure 10. Here, the absolute volumetric flow in each coordinate direction (calculated as in (28) but not taking the flow direction into account) is plotted over the volume fraction of the smaller vessels Λ S . The clearest picture emerges for the LS174T topology with a good correlation of flow in smaller vessels with their volume fraction. A similar, yet less distinctive trend is present for the SW1222 case whereas no relationship can be observed for the GL261 tumor.
Therefore, inside each REV we defined the isotropic permeability tensor as that is, a simple linear dependency of the permeability in the -th REV on the volume fraction of smaller vessels in the -th REV with proportionality constant . We also tested a nonlinear Kozeny-Karman law, but obtained slightly better results with the linear fit. Thus, we will exclusively study this linear dependency hereafter. Next, the optimization of the nonlinear leastsquares problem (30) is performed for the proportionality constant . Results are shown in Table 7 for the case 10 %, which can be compared with the cases with constant permeability over the entire domain from Table 5. We obtained a slightly better agreement in terms of flow in the smaller vessels 2 f low,Λ S and, thus, also for the total coefficient of determination 2 tot for the  TABLE 7 Results of the optimization procedure for non-constant permeability depending on volume fraction of smaller vessels. 10 % of 1D blood vessels have been retained in hybrid model for each tumor topology. Shown are the proportionality constant relating permeability and blood vessel volume fraction of smaller vessels inside each REV according to (39). 2 -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 coefficient of determination 2 tot between fully-resolved and hybrid model is calculated according to (29). (All data includes mean taken over five different sets of pressure boundary conditions on the 1D network per case) SW1222 and GL261 case. Compared to that, the correspondence of flow in the smaller vessels was markedly better than the constant permeability case for the LS174T topology. This is coherent with Figure 10 where the latter network showed the most evident correlation of blood flow on volume fraction. Thus, one could expect that no significant improvement was possible for the GL261 case where volume fraction and flow seem to be decoupled. However, also for the SW1222 topology, which showed at least a moderate dependency of blood flow on volume fraction, the agreement could not be increased significantly. Therefore, at least for one of our cases it was beneficial to include blood vessel volume fraction information into the hybrid model while it was not detrimental for the other two.
It also is conceivable that at least preferential directions of smaller vessels or their tortuosity can be detected non-invasively even though their complete structure cannot be resolved. A further enhancement of the model could be achieved when taking this information about the anisotropy of smaller vessels or their tortuosity into account during the homogenization procedure. 18,19,20 However, we want to emphasize that our whole study is based on numerical results. Experimental findings indicate no dependency between blood vessel diameter and flow in tumors 1,56,57 and a high vascular density does not automatically imply efficient perfusion, nutrient supply and drug delivery for solid tumors. 58 These properties could make it impossible to deduce permeabilities of blood vessels inside tumors from macroscopic quantities such as blood vessel volume fractions or preferential directions. By contrast, non-invasive measurements of perfusion 13,59 could prove helpful to enhance the hybrid model.
Similarly, improvements are possible for flow from the larger into the smaller vessels. In this contribution, we have assumed equal pressures in resolved and homogenized vasculature and, thereby, infinite (or at least a very large) permeability governing the flow between the two compartments such that a constraint of equal pressures holds along the resolved 1D vasculature. This has the major advantage that the coupling between resolved and homogenized vasculature is essentially parameter-free. Only the penalty parameter has to be chosen large enough such that a sufficient accuracy in the pressure constraint is achieved as described in Remark 4. The GL261 and LS174T case had a less dispersed distribution of the radius in connecting elements and, thus, of the permeability between larger and smaller vessels. For these topologies, our approach could estimate the mass transfer between larger 1D vessels and smaller homogenized vessels more accurately. The SW1222 case had a much higher variability of radius and flow between Λ L and Λ S . In this case, it could be advantageous to employ finite permeabilities to model the mass transfer and assign higher permeabilities to REVs or regions along the larger vessels where many branches go away from the main vessels. However, this would require additional parameterization of the model as well as additional data on regions where a lot of flow from larger into smaller vessels can be expected.
In addition, we have so far only employed very simple algorithms to optimize the parameters of the hybrid model. A much more powerful framework for coarse-graining physical models has been developed by Grigo et al. 60 and tested for flow through porous media. This could also be applied in our case to infer the optimal parameters of the hybrid model per REV. However, this would require much more microstructural features, such as tortuosity, blood vessel distances or radius data on the smaller homogenized blood vessels to calibrate the hybrid model. Again, it is questionable if this data can be acquired non-invasively and if these parameters are determining blood flow through tumors.

CONCLUSION
In this work, we have studied a hybrid embedded/homogenized model for computational modeling of solid tumor perfusion. Its guiding principle is that the complete morphology of vascular networks including blood vessel diameters and topology cannot be acquired with currently available in-vivo imaging techniques. Thus, fully-resolved discrete models relying on this data cannot be applied in "real world" scenarios. If, however, the structure of larger vessels constituting the main branches of the vasculature is available, our hybrid approach where only these larger branches are completely resolved is a sensible alternative. By contrast, the smaller scale vessels are homogenized such that their exact structure is not required anymore. This results in a two-compartment or double porosity formulation where the larger vessels are still represented as one-dimensional embedded inclusions. The coupling between the resolved and homogenized part of the vasculature is realized via a line-based pressure constraint along the 1D larger vessels, which we enforce with a mortar-type formulation with penalty regularization. This also has the advantage that compared to previous hybrid models no additional parameter -apart from the penalty parameter which has to be chosen large enough -is required to couple the two distinct representations of the vasculature.
The results of the hybrid model have been compared with reference solutions generated by a fully-resolved 1D-3D model. For that, three different network topologies extracted from three different tumor types grown in mice have been employed. These topologies consist of up to 420 000 vessel segments and have dimensions of up to 6 mm × 8 mm × 11 mm. To date, this is the largest and most challenging test case for a hybrid model, especially considering the abnormal and tortuous structure of the networks typical for the vasculature inside tumors. We have further shown how we artificially generate the hybrid from the fullyresolved model, define representative elementary volumes and assign boundary conditions. We are confident that the artificially created topologies of larger vessels are representative of real in-vivo imaging data sets of larger vessels inside tumors such that they enable us to draw meaningful conclusions for more realistic scenarios where the full topology is not available such that a hybrid approach is the only option.
For comparison of the results of the two models, we have defined several rigorous metrics involving the blood pressure in both resolved and homogenized vasculature, the pressure in the interstitial fluid and blood flow in the homogenized vasculature. These metrics have then been employed to obtain the optimal parameters for the hybrid model and to study its accuracy w.r.t. the fully-resolved one. We have obtained very good agreement in terms of blood pressure in the larger vessels and IF pressure. Larger deviations are present for blood pressure and flow in the homogenized vasculature. However, these limitations can be expected since the information on the smaller vessels is not retained in the hybrid model. Overall, the best correspondence has been achieved for the SW1222 case which also had the clearest vascular structure and distinction between larger and smaller vessels. All topologies showed a very good agreement in terms of REV IF pressure and REV blood pressure in smaller vessels with mean deviations in a range of 20 − 70 Pa and 40 − 110 Pa resp. 0.7 − 3.8 % and 1.1 − 2.9 %. It is sufficient to resolve 5 − 10 % of all blood vessels segments by keeping them in the hybrid model since there is only a marginal improvement of the agreement with the fully-resolved model in terms of all investigated metrics when retaining a higher percentage (15 − 20 %) of blood vessels. Concerning the flow between smaller and larger vessels the error was mainly caused by the large variability of diameter and flow in the connectivity elements between large and small vessels for the SW1222 case. Possibly, this error could be reduced by allowing a varying permeability for coupling the two compartments. By including information about the blood vessel volume fraction of smaller blood vessels into the definition of their permeability tensor a better agreement with flow therein could be achieved for the LS174T case. Nevertheless, the abnormal vascular structure and blood flow patterns of tumor vasculature could impede this approach.
Several other potential improvements have been discussed and remain subject to future work. Furthermore, the inclusion of species transport to simulate drug delivery or nutrient transport lies at hand. Species transport including the coupling between resolved and homogenized vasculature is possible within our hybrid multiphase tumor growth model 3 and we have already studied nanoparticle delivery to solid tumors employing the homogenized compartment only. 61 These models could ultimately enhance our understanding of the limitations of current drug delivery strategies and aid in devising more targeted therapies.
The next step towards a more realistic or even clinical usage of hybrid computational models for tissue perfusion is to devise a strategy which combines data which is available non-invasively. 23 Faced with such a scenario, where a hybrid model is the only applicable option since the entire network topology is not known, the methods and metrics developed here could be applied in the following way: 1. Gather all physiological data, which can be accessed via in-vivo measurements for the specific case. For instance, this could be tissue perfusion, hypoxic areas, REV or point-wise measurements of pressure or flow, volume fractions of homogenized blood vessels or their preferential direction and the transport of tracer molecules through the domain. 3. Formulate an optimization problem similar to (30) to match the available information about transport, e.g., REV or pointwise flow and pressure data. However, not only the parameters of the homogenized vasculature would be part of the optimization as in this contribution but also the large majority of the (homogenized or resolved) pressure boundary conditions, which are additionally unknown. Note however, that far less boundary conditions compared to a fully-resolved setting need to be applied.
4. Employ the obtained flow state for in-silico studies of drug delivery or the optimization of treatment strategies.

APPENDIX A SEGMENT-BASED LINE INTEGRATION SCHEME FOR THE EVALUATION OF 1D-3D COUPLING TERMS
The numerical integration of the 1D-3D coupling terms is outlined in this appendix. These terms are the leakage terms from the 1D embedded vasculature into the IF occurring in the weak forms (5a) and (5b) respectively (13a) and (13c) and the mortar coupling matrices and from (16) and (17). After spatial discretization with appropriate shape functions, these terms comprise a line integral along the inclusion containing the product of shape functions either defined on the 1D or on the 3D domain. A one-dimensional segment-based integration for these types of terms has been proposed in our previous publications. 3,36 The goal of this procedure is to perform the integration in a general non-conforming scenario as sketched in Figure A1. Numerically, the 1D integrals are evaluated with Gauss quadrature. However, at first each 1D element is segmented by finding its 2D or 3D interaction partners, i.e., those elements of the discretization of the surrounding domain it intersects. This yields 1D pieces interacting with a single 2D/3D element on which also the shape functions of the respective 2D/3D element are continuous. Then, GPs are defined in the single segments and mapped from the element parameter space on the 1D centerline. In addition, the spatial coordinate of the respective GP ( ( )), is required to obtain the shape function values of the respective 2D/3D element at this position. The integral for a specific 1D element then emerges as the sum over the integrals of all its segments and the respective contributions are assembled into the DOFs of the 1D and the interacting 2D/3D elements.