An approach for vascular tumor growth based on a hybrid embedded/homogenized treatment of the vasculature within a multiphase porous medium model

The aim of this work is to develop a novel computational approach to facilitate the modeling of angiogenesis during tumor growth. The preexisting vasculature is modeled as a 1D inclusion and embedded into the 3D tissue through a suitable coupling method, which allows for nonmatching meshes in 1D and 3D domain. The neovasculature, which is formed during angiogenesis, is represented in a homogenized way as a phase in our multiphase porous medium system. This splitting of models is motivated by the highly complex morphology, physiology, and flow patterns in the neovasculature, which are challenging and computationally expensive to resolve with a discrete, 1D angiogenesis and blood flow model. Moreover, it is questionable if a discrete representation generates any useful additional insight. By contrast, our model may be classified as a hybrid vascular multiphase tumor growth model in the sense that a discrete, 1D representation of the preexisting vasculature is coupled with a continuum model describing angiogenesis. It is based on an originally avascular model which has been derived via the thermodynamically constrained averaging theory. The new model enables us to study mass transport from the preexisting vasculature into the neovasculature and tumor tissue. We show by means of several illustrative examples that it is indeed capable of reproducing important aspects of vascular tumor growth phenomenologically.

of the most important aspects of cancer is that its governing mass transport phenomena may span multiple scales ranging across the organism to tissue, cellular, and molecular levels. Especially, transport of nutrients or growth signals across multiple scales and various biological barriers plays a crucial role in tumor growth. Concurrently, also drugs have to be transported across these scales and overcome these barriers which is still an inhibiting factor to efficient drug delivery and, hence, therapeutic success.
Tumor growth is commonly divided into three distinct stages, namely the avascular, the vascular, and the metastatic phase. 5 In the first avascular phase, tumor growth is confined to a size of approximately 1-2mm 36,7 since pure diffusion of nutrients is not enough to sustain the vast consumption by tumor cells (TCs). However, solid tumors can acquire the capability to trigger angiogenesis, which is the formation of new capillaries from the preexisting vasculature. A new blood vessel network evolves which supplies the tumor with nutrients and, thus, enables rapid growth. This so-called angiogenic switch 8,9 from avascular to vascular growth is crucial for the tumor to become harmful. Therefore, a sophisticated cancer model able to provide clinically relevant data has to take angiogenesis and vascular tumor growth into account.
The main objective of this work is the development of a novel coupling method to investigate angiogenesis and the multiscale transport of species from the preexisting vasculature through the tumor neovasculature (NV) to tumor tissue. We want to combine a discrete, 1D representation of the preexisting vasculature with the continuous or homogenized representation of the NV based on the theory of porous media, which we have recently introduced in our multiphase tumor growth model. 10 Most vascular tumor growth models incorporate a discrete angiogenesis model into a continuum representation of tumor tissue resulting in a hybrid vascular tumor growth model according to the common classification of models. 5,11,12 Note that we employ the term discrete also to describe angiogenesis models that aim for a full resolution of capillaries and blood flow therein. These models are mainly adapted from the work of Anderson and Chaplain, 13 see also the review papers [14][15][16] and the representative examples of such hybrid models. [17][18][19][20][21][22][23] More recent state-of-the-art fully resolved blood vessel network models 20,[24][25][26][27] are based on an integrative framework for vascular remodeling 16 including angiogenesis as well as vessel regression, dilation, and collapse during tumor progression. Thereby, the initial arterio-venous vasculature in host tissue develops into a tumor-specific vessel network.
The reason why we favor a combined approach over a full resolution of preexisting and NV is the abnormal structure and mechanics of tumor NV. It is well-established that tumor-induced angiogenesis results in a tortuous, dilated blood vessel network with variable vessel lengths and diameters and without the usual vascular hierarchy of arterioles, capillaries, and venules. 28,29 This causes highly heterogeneous blood flow. 30 While at first sight, discrete models seem to offer more insight into the formation of the specific network and its structure, this can and should also be challenged. It appears virtually impossible to resolve the complex morphology of the network and the full spatial and temporal heterogeneity of blood flow inside tumors where almost no relationship between vessel diameter and flow velocity is present. 31,32 Furthermore, it is more than questionable if a fully resolved blood vessel network offers more information on the actual quantities of interest especially when considering the inherently stochastic nature of angiogenesis and tumor vasculature remodeling which precludes predicting the specific in vivo network topology. Relevant quantities which could be obtained from in silico models and can actually be acquired through imaging are microvascular densities, hotspots of vascularization or very badly vascularized regions inside the tumor, and averaged transport properties to detect hypoxic and drug resistant areas. 33,34 To obtain this information from resolved models, averaging over several different simulations is necessary. 25 Therefore, we adopt a homogenized or smeared representation of the NV while resolving the preexisting larger vessels, for which information about structure and blood flow might be available, with a discrete, 1D flow and species transport model. This may be classified as a hybrid formulation according to the definition of Vilanova et al. 5 A similar approach for mass transport in tissue has been pursued in the so-called composite smeared finite element method developed by Kojic et al. [35][36][37] In the latter two contributions, the authors have coupled fluid flow and species transport between 1D capillaries with a homogenized or smeared representation of the capillary bed and tissue domain. By means of several examples, they have demonstrated that a homogenized and a resolved representation of the capillary bed yield comparable results for blood flow and passive scalar transport. However, one restriction of their approach is that the 1D and the 3D domain cannot be discretized independently. This is a major drawback of the method because it only allows to couple spatially conforming discretizations. If a more complex microvascular geometry is considered, the meshing procedure may become quite intricate. This effort can be considerably reduced by allowing arbitrary 1D and 3D configurations. Such embedded multiscale approaches in the context of diffusion-reaction scenarios have been theoretically studied by D'Angelo and Quarteroni 38,39 and further applied to drug delivery to tumors 40,41 and nanoparticle-mediated hyperthermia cancer treatment. 42,43 To enforce the coupling between the embedded and the homogenized vasculature with independent 1D and 3D discretizations, we have developed two line-based penalty methods, namely a Gauss-point-to-segment (GPTS) and a mortar penalty (MP) scheme. In addition, a large deformation approach is adopted throughout allowing also deformations of the embedded fluid network. The main purpose for which we have developed our nonconforming coupling method is vascular tumor growth. However, it is not restricted to this specific problem but may be applied to any mass transport phenomenon involving a coupling between 1D and 3D domains, particularly if a part of the domain can be considered in a homogenized sense because the full morphology of smaller scales is not of interest. Such problems might be other biological mass transport phenomena involving the microcirculation but also geomechanical applications. Moreover, a coupling between 1D slender objects embedded into an enclosing 3D structure as appearing in many engineering materials could be realized.
The remainder of this paper is structured as follows: in Section 2, we demonstrate how the embedded 1D fluid network can be incorporated into our vascular multiphase tumor growth model. Especially, the coupling between the preexisting vasculature and the NV formed through angiogenesis will be addressed. Two different strategies to enforce this coupling are introduced in Section 3. After space discretization with finite elements and time discretization with the one-step-theta scheme, the entire framework is solved with a monolithic coupling scheme. We show the principal applicability of the novel formulation to vascular tumor growth by means of several illustrative examples in Section 4. Section 5 contains a short summary of our findings and possible future extensions.

| INCORPORATION OF THE EMBEDDED 1D FLUID NETWORK INTO THE VASCULAR TUMOR GROWTH MODEL
Our mathematical model is based on distinguishing two domains: a 1D inclusion representing the preexisting capillary network described in Section 2.1 and a 3D continuum corresponding to the surrounding tissue including the NV, see Section 2.2. These two domains are fully coupled via transcapillary exchange terms but also through tissue deformation, as sketched in Figure 1.
Before going into detail about the governing equations of the model and their derivation, we want to present its general layout. For that, we have collected all entities (phases and species) and their interaction through mass exchange terms in Figure 2. The model comprises five phases: 10 The extracellular matrix (ECM) is a deformable, porous solid phase; three liquid phases, namely tumor cells (TCs), host cells (HCs), and the interstitial fluid (IF) flow in the pores of the ECM; and the vasculature, which we either resolve as an embedded 1D fluid network (preexisting vasculature

FIGURE 2
Schematic overview of the model and mass transfer relations (PRE)) or as a homogenized additional pore space inside the ECM (neovasculature (NV)). Blood, which is also modeled as a liquid, flows in the vasculature. The lymphatic system is not present as an additional phase, but the drainage of excessive fluid from the IF is taken into account through a suitable mass transfer relation. Furthermore, each phase can consist of multiple species. In general, we can specify an arbitrary number of species, but here, only the ones depicted in Figure 2B are considered. Oxygen is transported in the vasculature and may pass the blood vessel wall into the IF. It is then consumed by TCs and HCs. TCs can be either necrotic or viable. Therefore, the TC phase is composed of two species: necrotic tumor cells (NTCs) and living tumor cells (LTCs). The last species are socalled tumor angiogenic factors (TAF), which are produced during hypoxia by LTCs and diffuse in the IF to trigger angiogenesis. 13

| Governing equations for blood flow and species transport in 1D vasculature domain
We begin with blood flow and species transport in the 1D blood vessel domain Λ t corresponding to the preexisting vasculature. The computational effort for solving these equations can be tremendously reduced under the simplifying assumptions of Poiseuille flow. However, this 1D model is often also applied even when the simplifications are strictly speaking not valid, eg, when transcapillary flow is present but small compared to axial flow, see Nabil and Zunino 43 for a short discussion. Recent experimental data 44 suggests that the deformation of blood vessel networks in the brain subject to mechanical stress exerted by growing tumors may be significant. Therefore, we adopt a large deformation approach of the tissue, ie, the ECM as well as the embedded network. In the following, we will discuss our assumptions and their impact on the model formulation in an Arbitrary Lagrangian-Eulerian (ALE) setting. The main assumptions are as follows: 1. The deformation of the blood vessel network is completely described through the deformation of the underlying porous medium domain. 2. Blood vessel segments have negligible stiffness compared to the ECM. 3. The area of the blood vessel segments remains constant under deformation.
Our approach could also be extended to remove one or all of the above assumptions, but as they appear reasonable for our problems of interest, they are introduced to not artificially complicate the method. The first two assumptions are based on the fact that we only consider small capillaries. Then, the surrounding ECM acts as a scaffold which provides mechanical integrity to the enclosed capillary. 45 If the scaffolding material deforms, the blood vessel deforms equally. Furthermore, we neglect any contribution of these small capillaries to the stiffness of the surrounding ECM since the capillary wall does not consist of smooth muscle cells but only a thin basement membrane which provides structural integrity. 46 The third assumption is a major simplification of the formulation because lateral deformation of blood vessels induced through the deformation of the ECM is neglected. For example, axial elongation of a blood vessel will not evoke any radial constriction. However, we believe that this effect is insignificant compared with the much more substantial influence of blood vessel constriction and collapse during tumor growth. It is well-known that the growing tumor mass and high interstitial pressure lead to blood vessel compression or even collapse. This could easily be included through a suitable constitutive law 20,47 and will be a topic of future research. However, the main focus of this contribution lies on the development of a hybrid, alternative treatment of angiogenesis, which could serve as a less costly yet accurate alternative to classical fully resolving models.
The notation of the problem is depicted in Figure 1. We denote the 1D vasculature and the 3D tissue domain in reference configuration as Λ 0 and Ω 0 , respectively. Their counterparts in spatial configuration are termed Λ t and Ω t . Assumption 1 from above holds if a point X Λ ∈Λ 0 on the 1D vasculature reference configuration and a reference point X Ω ∈Ω 0 of the 3D domain, which are at the same location in reference configuration, share a common point x Ω ∈Ω t throughout the deformation. Hence, the superscripts for the distinct domains can be dropped. Furthermore, we refer to the arc-length coordinate of each vessel segment in reference configuration as S and in spatial configuration as s. The two are related through the deformation of the domain via with unit tangent vector in reference configuration t 0 and deformation gradient The corresponding unit tangent vector in spatial configuration follows as In the following, quantities such as pressures or velocities defined in the vasculature domain will be identified by

| Blood flow in the 1D vasculature domain
The mass conservation equation for an incompressible fluid in the vasculature domain may be written as 48 ∂A ∂t with blood vessel area A, area-averaged fluid velocity v b v , and blood density ρ b v . The considered time interval is denoted as t 0 ; t E ½ . Furthermore, a generic right-hand side term has been defined. In our notation, it expresses a possible mass transfer term between the 1D vasculature and the IF. For example, this term can be employed to model leakage of fluid over the blood vessel wall into the interstitium. All mass transfer terms are written in generic form in the model formulation. Their exact form will be given in Section 2.4, see also Remark 1. Assuming cylindrical blood vessels with constant radius R (Assumption 1) and the Poiseuille relation the balance of mass can be rewritten as In the Poiseuille Equation (5), we have included the projection of the skeleton velocity v s in current tangential segment direction t t in order to account for the underlying solid phase movement. Hence, this term is similar to the Darcy equation in a porous medium. The pressure gradient in the capillaries induces a flow relative to the solid phase movement in axial vessel direction. Blood viscosity μ b v is assumed constant in the following. In reality, however, blood behaves as a non-Newtonian fluid due to the presence of red blood cells. Therefore, the viscosity depends on blood vessel radius R and hematocrit, ie, the volume fraction of red blood cells within blood. Constitutive relationships for the apparent viscosity of blood based on experimental data 49,50 could easily be integrated.

| Species transport in the 1D vasculature domain
Blood flow in the vascular network can transport several species such as oxygen, anticancer drugs, or nanoparticles. These are advected by the flow and may cross the vessel walls into the IF. The mass balance of a species i dispersed in the blood vessel network with mass fraction ω ib v may be written as Such 1D diffusion-advection-reaction equations with diffusivity D ib v have been used previously in related work. [41][42][43] However, here the formulation is in terms of a species mass fraction ω ib v rather than a concentration which makes it easier to couple it to our multiphase model where species are generally identified by mass fractions. Again, the righthand side represents a generic species mass transfer term from species i in the IF to the considered species i in the preexisting vasculature b v. Note that we have omitted any intra-phase reaction terms as those are not present in the current model. We invoke the relation to transform the temporal derivative on the 1D domain to a derivative w.r.t. the reference configuration. As in (5), the solid velocity v s projected in current tangential segment direction t t is employed. Inserting into (7) and applying the product rule for the convective term yields an ALE equation Herein, the balance of mass of blood (4) under the assumption of constant area has been invoked. The convective term in ALE form can then be replaced by the Poiseuille law (5), which yields the final mass balance equation for species i Equations (6) and (10) are the governing equations for fluid flow and species transport in the preexisting vascular network. They are coupled in three ways to our multiphase tumor growth model for the surrounding 3D tissue. First, in a one-way coupling, the deformation of the ECM induces also deformation of the embedded blood vessels. Second, fluid flow and species transport in the preexisting vasculature and in the surrounding IF are two-way coupled via transvascular exchange over the blood vessel wall. Third, fluid flow and species transport in the preexisting vasculature and the NV formed during angiogenesis are coupled through suitable constraints, see Section 2.3.

| Governing equations of the vascular multiphase tumor growth model
The coupling with the embedded 1D fluid network requires some minor additions to our vascular multiphase tumor growth model, which we have introduced in Kremheller et al. 10 It is based on the original tumor growth model of Sciumé et al, 51 which has been derived on the basis of thermodynamically constrained averaging theory (TCAT). This intricate mathematical procedure guarantees physical, thermodynamical, and mathematical consistency between multiple scales. 52,53 Hence, at continuum scale, where the model is formulated, the macroscale variables are explicitly and precisely defined as averages of microscale variables. Over recent years, this TCAT tumor growth model has been enhanced to include three-phase flow, 54 a deformable ECM 55 and cell lysis and matrix deposition. 56 Microstructural and biomechanical features of the ECM have been studied in Santagiuliana et al. 57 As stated above, the vascular multiphase model comprises five phases: 10 The ECM, TCs, HCs, the IF, and the NV. The NV is modeled in a homogenized way as an additional pore space in the ECM, which develops during angiogenesis and in which blood flow takes place. The composition of a typical representative elementary volume at the microscale is sketched in Figure 3. Here, single phases and their interfaces are clearly discernible. However, at the macroscopic scale where our equations are valid, several phases α with corresponding volume fractions ε α are juxtaposed. For instance, single blood vessels of the NV are not modeled explicitly but in an averaged sense with volume fraction or capillary density ε v . Effectively, the formulation results in a double-porosity model with two separate porous networks, which are the pores of the ECM and the NV. In the following, the solid phase will be denoted with superscript s, the IF with l, HC with h, TC with t, and the NV with v. The sum of all volume fractions has to satisfy In addition, multiple species are present. In the TCAT framework, a generic species i which is dissolved in a specific phase α is described by its mass fraction ω i α in this phase. For example, oxygen (or nutrient) may be present in the NV with a mass fraction ω n v and in the IF with a mass fraction ω n l . As in previous versions of the model, the TCs partition into LTCs and NTCs. Therefore, NTCs constitute the second species with mass fraction ω N t . Then, the mass fraction of LTCs follows as ω L t ¼ 1 − ω N t , so that we do not explicitly have to solve the balance of mass of LTCs. The third species are the aforementioned TAF, which diffuse in the IF with mass fraction ω TAF l . The model equations for all phases and species are introduced in the following. In principle, additional types of cells such as cancer stem cells or cancerassociated fibroblasts could be introduced as additional phases or species into the model.

| Angiogenesis and blood flow in the neovasculature
Due to the vast consumption of oxygen by proliferating TCs, the tumor microenvironment can easily become hypoxic. Then, some tumors can form their own vasculature through angiogenesis, that is, the formation of blood vessels from a preexisting vasculature. This is achieved through the secretion of TAF, which triggers endothelial cell migration from the preexisting vasculature toward the tumor. Subsequently, these endothelial cells proliferate and sprout until a new vascular network, the NV, emerges, which provides oxygen and other nutrients to the tumor. As mentioned before, the NV exhibits abnormal and heterogeneous structure. Hence, we do not want to account for every single blood vessel but treat it in a homogenized or averaged sense. We assume that the volume fraction of the NV ε v follows an evolution equation which we have introduced in Kremheller et al. 10 This equation is an adaption of a very prominent formulation of angiogenesis originally proposed by Anderson and Chaplain. 13 It establishes a mathematical description of angiogenesis in a continuum sense. The main trigger of angiogenesis is the chemotactic response on a TAF gradient, which is modeled by the fourth term on the left-hand side of the equation with chemotactic coefficient χ ω TAF l , see also Equation (A7).
Random motility of endothelial cells is included as a diffusive term. We want to emphasize that the previous equation is not a mass balance equation but rather a constitutive equation which describes the evolution of the NV volume fraction due to endothelial cell migration. Previously, blood flow in the NV has been neglected by assuming constant blood pressure. 10 Now, also blood flow in the NV is modeled. This enables us to couple the model with the preexisting vasculature and study blood flow and species transport in the NV. It is crucial to include mass transport in the NV because of its importance during tumor progression: nutrients can be transported from preexisting into NV and to the tumor more efficiently than through diffusion in the IF such that the tumor has improved access to nutrients, which is why vascular tumor growth progresses much more rapidly. To model these phenomena, we introduce the balance of mass of blood in the NV in ALE form with blood velocity v v . Assuming constant blood density ρ v and a Darcy equation with (isotropic) NV permeability tensor k v and blood viscosity μ v , the equation may be rewritten as in terms of blood pressure in the NV p v . This equation constitutes a homogenized model for blood flow in the NV. In conclusion, our porous medium model for the NV comprises two equations: the first (12) governs angiogenesis, and the second one (14), blood flow. Therefore, the volume fraction of the NV ε v and its blood pressure p v are two primary variables of the vascular tumor growth model. However, in the current model, only the emergence of new vessels is considered in (12) but no vessel remodeling, eg, through blood vessel regression, collapse, or dilation. 16 In fully resolving vasculature models, these phenomena are usually taken into account quite naturally through suitable constitutive equations for the blood vessel radius. 20,24,26,27 In our homogenized treatment of the NV, vascular remodeling could also be incorporated as an additional effect in Equation (12), for instance, via a right-hand side term, which induces a decrease in vessel volume fraction based on the pressure exerted by the tumor on the NV. Similarly, remodeling effects could also be considered in the equation for blood flow (14) through a constitutive law for the permeability including effects of anisotropy.

| Tumor cells, host cells and the interstitial fluid
In the vascular model, we have defined two separate porous networks, see also Figure 3. The first one is the NV with blood flow. Cells and the IF occupy a second porous network, namely the space between the ECM fibers, that is, the pores of the ECM. Their volume fraction ε is given by From that, saturations of TCs, HCs, and the IF are readily obtained as which have to fulfill The mass balance equations for TCs, HCs, and the IF read as for α=t,h: for α=l: These equations are equivalent to the ones derived in Kremheller et al. 10 Herein, ψ α denotes the generic primary variable of phase α, which can be either its saturation S α , its pressure p α , or a pressure difference p αβ to another phase β. Furthermore, incompressible fluid and solid phases have been assumed. The equations above represent the balance of mass of TCs, HCs, and the IF. The first two (18) are simply the balance of mass of TCs and HCs, and the last one (19) can be obtained by summing up the mass balance equation of the three phases such that several terms cancel out due to the relation for saturations (17). The mass balance equations (18) and (19) have to be closed by an appropriate pressuresaturation relationship and the balance of mass of the solid phase in reference configuration, see Appendix A. In addition, a generalized Darcy law with isotropic permeability k and viscosity μ has been employed as in previous versions of the model. The right-hand side represents mass transfer between different phases.
The major modification due to the coupling with the 1D blood flow model is the addition of the mass transfer term scaled with the Dirac measure δ Λ t located on Λ t in (19). This term is the counterpart of the right hand side mass transfer in the balance of mass of blood (4) and expresses the coupling between flow in the embedded 1D domain and the surrounding tissue, typically through an outflow of fluid from the vascular network into the IF. In other words, the 1D problem is represented as a distributed Dirac source term embedded into the 3D problem. This singularity may lead to suboptimal convergence rates. 39,58 Recently, a different approach has been proposed 58 : Instead of prescribing the source term as a Dirac measure on the middle line of the inclusion, the coupling can be evaluated on the actual physical boundary, that is, the outer surface of the 1D vessels. Thus, the regularity of the solution is increased by lifting the dimension by one. An alternative would be to approximate the Dirac measure by smearing the term over a finite length. 26 Remark 1. The numerical discretization of the formulation proposed by D'Angelo and Quarteroni 38,39 is quite intricate since it involves evaluating the average value of quantities from the 3D domain w.r.t. a circle with the actual radius of the embedded domain. This local average of a generic quantity u for every point s of the embedded domain at a specific segment with radius R is defined via the integral Herein, n s; θ ð Þ is employed to describe the unit circle around Λ t (s) perpendicular to the segment direction t t . For practical problems with h>R, the integral can be approximated by as proposed by D'Angelo, 59 which is equivalent to our method.

| Species transport
For every species i dissolved in a phase α with mass fraction ω i α , the balance of mass reads as with and effective diffusivity D iα eff . Compared with the previous formulation, 10 we have now also included species transport in the NV, specifically of oxygen, see also Figure 3. In the aforementioned paper, the oxygen mass fraction in the NV has simply been assumed as constant. The first three mass transfer terms on the right-hand side are a generic mass transfer term from all species i in all phases κ to the considered species i in phase α, an intra-phase reaction term, and a term which arises when applying the product rule. The last two terms are due to the exchange of species from the embedded blood vessels to species in the IF, ie, these terms arise only for α=l in (22). Once more, they are represented as Dirac source terms δ Λ t on the 1D manifold Λ t . Species such as oxygen can pass the blood vessel wall into the IF. This exchange is included by mass transfer terms in the balance of species mass in the 1D vascular network (10) which is coupled to species mass transport in the 3D domain through the embedded multiscale method and the corresponding distribution of Dirac source terms in (22).

| Solid phase
We model the solid phase, ie, the ECM as the skeleton of the porous medium, where Terzaghi's principle of effective stress can be applied. Herein, the solid pressure is defined as a weighted average of the involved phases. 10 This formulation recovers the definition of the solid pressure as a weighted sum of TC, HC, and IF pressures with their saturations 51,60 if no NV is present. If angiogenesis occurs, that is, a second porous network is present, it is similar to the definition in multiple poroelastic network theory. 61 Pulled back into the material configuration Ω 0 , the balance of momentum of the solid phase can then be written as with the material divergence operator ∇ 0 . Herein, body and dynamic forces have been neglected.

| Coupling between 1D (resolved) and 3D (homogenized) representation of vasculature
The complexity of blood vessel networks in humans can be enormous. 36 Especially, the NV formed by tumors during angiogenesis is characterized by its chaotic structure and morphology accompanied with irregular flow patterns. [28][29][30] Hence, Kojic et al 36 have established an approach termed the composite smeared finite element method. Here, larger vessels are represented as 1D inclusions, while the capillary bed is smeared or homogenized. The parallels to our homogenized representation of the NV lie at hand. However, the two domains are coupled via so-called connectivity elements at spatially matching nodes of the two grids. 36 Therefore, the main benefit of the embedded multiscale method, which is the complete independence of 1D and 3D grids, is lost. In analogy to mesh tying in solid mechanics, we propose an alternative approach based on a nonconforming coupling between the embedded vasculature and the 3D domain, more precisely the homogenized vasculature. We define the constraint between the embedded and the surrounding domain as for coupling of blood flow or for coupling of species transport. Herein, a "gap" g between pressures or species mass fractions in the 1D domain w.r. t. the 3D domain has been introduced. It has to vanish in order to couple blood flow or species transport of the embedded blood vessel with the homogenized NV. Blood flow in the preexisting vasculature on Λ t is governed by Equation (6), while the corresponding equation for the NV in Ω t is Equation (14). The two domains are coupled via the constraint (27). The same holds for species transport in the preexisting vasculature (10) and the NV (22) with ε α =ε v . In Section 3.2, we will elaborate how to enforce this coupling between the two distinct domains with two penalty-based methods, a GPTS, and a MP approach.
Remark 2. The nonconforming coupling between preexisting and NV (27) is enforced along the axis of the preexisting vessels. This is motivated by the origination of angiogenesis laterally from a preexisting blood vessel through angiogenic sprouting followed by growth radially away from the pre-existing vessel and subsequent network formation. 62 In Section 4.2, we demonstrate how angiogenesis can be triggered from the preexisting vessel with our method. For alternative applications, where a vessel network originates from the tip of a major vessel, a formulation with a constraint on the tip of the resolved vessel and a homogenized continuation of a tree-like vessel network is also feasible.

| Constitutive mass transfer relations and oxygen transport model
In this section, we want to summarize the mass transfer relations sketched in Figure 2, which we employ to model vascular tumor growth. Several of the terms from the previous version of the multiphase model 10,51,54,55,63 have been reused and are listed in Appendix A. A summary of all relations has been collected in Table 1. Note that some terms which are of minor importance have been neglected, mostly the ones which would appear for species transport due to the application of the product rule in (10) and (22). For the five phases of the model, the terms for tumor growth M l→t growth , leakage of fluid from the NV into the IF M v→l leak , and drainage of excessive fluid by the lymph system M l→ly drainage have not been modified, see Appendix A. One addition for the multiphase model is that blood flow in the NV is now explicitly modeled through Equation (14). Thus, the leakage of fluid (A8) also appears as a mass transfer term in this equation. It is wellknown that the NV inside tumors is leaky, which might contribute to increased IF pressure inside tumors together with inefficient lymphatic drainage. 64,65 Furthermore, there is also leakage of fluid from the preexisting vasculature into the IF, which can be considered by the Starling equation with hydraulic conductivity L p;b v . Blood plasma leaks from the preexisting vasculature into the IF if the effective pressure is larger than the IF pressure p l . In the previous equation, the average osmotic reflection coefficient is denoted by ω osm and the osmotic pressures of the plasma and the IF by π blood and π l , respectively. Macaulay brackets · h i þ have been employed to allow only outflow of fluid. Please note that here and in the following mass transfer terms involving the preexisting vasculature, we have made use of the approximation (21).
For the three considered species, we have also adopted the mass transfer terms for consumption of oxygen by TCs M nl→t cons , the intra-phase reaction term ε t r Nt for necrosis, and the production of TAF by hypoxic TCs M TAFt→TAFl prod . The arithmetic expressions for these terms are again given in Appendix A. Consumption of oxygen by HCs has so far been neglected. Here, we include it with a similar term as the consumption of TCs (A12). The oxygen consumption by HCs is defined as Herein, ω nl env is the normal mass fraction of oxygen dissolved in the human plasma and γ nh 0 is a model constant related to the oxygen demand of HCs.
In Kremheller et al, 10 we have not taken oxygen transport in the NV into account. In this contribution, oxygen transport both in the preexisting and in the NV is included through the mass balance equations (10) and (22), respectively. However, in our TCAT framework, the primary variable is the mass fraction of oxygen instead of oxygen concentrations per unit volume of plasma or tissue C nα ½mlO 2 =ml or oxygen partial pressures P α oxy ½mmHg as commonly used for oxygen transport models. 27,[66][67][68][69] Hence, we need to convert mass fractions to partial pressures in order to reuse the mass transfer relations which are normally applied. The mass fraction of oxygen in the IF in terms of oxygen partial pressure in the IF P l oxy is given by Henry's law as where α l is the solubility of oxygen in the IF and ρ n and ρ l are the respective densities of oxygen and the IF. Oxygen transport in blood follows a more complex mechanism 27 since it can be either dissolved in the blood plasma or bound to hemoglobin. Hence, the mass fraction of oxygen in blood (either ω n v or ω nb v ) may be written as Tumor angiogenic factors Abbreviations: IF, interstitial fluid; NV, neovasculature; rhs, right-hand side.
As in (33), the first term represents the dissolved oxygen with effective solubility α v,eff , while the second one stands for the oxygen bound to hemoglobin (whose contribution to the total amount of oxygen is actually much larger). So, ω n v stands for the total mass fraction of oxygen which is dissolved in plasma and bound to hemoglobin. The same applies for ω nb v . Furthermore, H D is the discharge hematocrit, that is, the volume flux of red blood cells divided by total blood volume flux 50 and C nv 0 the concentration of oxygen at maximum saturation. We do not consider a heterogeneous distribution of hematocrit due to diverging bifurcations 50 as in other numerical models 27,68,70 but constant discharge hematocrit, see Table 4. For the binding of oxygen to hemoglobin, the Hill equation is typically applied. Here, the Hill exponent n and the partial pressure P v oxy; 50 at 50% oxygen saturation have been introduced. As stated above, mass fractions of oxygen are the primary variable of our oxygen transport model, so we actually need the inverses of (33) and (34) to get the oxygen partial pressure at a specific mass fraction. For oxygen in the IF (33), this is trivial; however, the relation for oxygen in blood (34) has to be inverted numerically by a small local Newton algorithm. For transvascular oxygen exchange from the preexisting vasculature into the IF, we then employ the relation in terms of partial pressures of oxygen which has been proposed by Welter et al 27 based on earlier works of Hellums et al. 71 Basically, the parameter b γ tv R ð Þ models the radial transport resistance of oxygen. Such a radial resistance 71 has also been employed elsewhere. [66][67][68] Here, a phenomenological fit for varying radii is reused, see appendix S1 of Welter et al. 27 For oxygen exchange from the NV into the IF, we use an equivalent relation We assume that oxygen mass exchange is proportional to the volume fraction of the NV ε v and its surface-to-volume ratio S=V ð Þ v . Macaulay brackets are employed to allow only oxygen transfer from the NV into the IF and not vice versa.

| NUMERICAL DISCRETIZATION AND COMPUTATIONAL FRAMEWORK
The framework developed for the vascular multiphase model 10 within our high-performance computing platform 72 has been enhanced with the embedded multiscale method. We employ the finite element method to discretize the system of equations in space and the one-step-theta scheme for time discretization. The space discretization of the multiphase system including species transport and structure deformation has been described in detail in Sciumé et al, 55 while the discretization of the 1D embedded domain is trivial. Therefore, we will focus on the interaction terms between the domains, namely the transcapillary mass transfer terms and the two proposed constraint enforcement strategies.

| Treatment of transcapillary exchange terms
We will only outline the discretization of the mass transfer relations between embedded vasculature and IF for the sake of brevity. A generic mass transfer term in our formulation may be written as for fluid mass transfer or for species mass transfer with function f · ð Þ. When brought to the left side, the contribution of the mass transfer term to the weak form of fluid flow resp. species transport in the embedded domain Λ t is given by while the contribution to the weak form of the equation for the surrounding tissue reads as Herein, test functions δφ b v ∈ H 1 0 Λ t ð Þand test functions δφ l ∈ H 1 0 Ω t ð Þhave been employed. 58 Only piecewise linear elements on Λ t and bi-, respectively, trilinear elements in Ω t are used in the following. Furthermore, ·; · ð Þ denotes the standard inner product on either Λ t or in Ω t . Spatially discretizing this term requires Gauss integration of products of shape functions along the length of the inclusion as described in Section 3.4.

| Constraint enforcement strategies
In the following, we show how the constraint (27) for coupling blood flow and species transport between preexisting and NV can be fulfilled through two different constraint enforcement strategies, a GPTS, and a penalized mortar-type method, which have been adapted from mesh tying and contact formulations in solid mechanics. These two methods perfectly fit into the finite element discretization of the transcapillary mass transfer terms as described above.

| Gauss-point-to-segment approach
Our first choice to fulfill the constraint (27) of equal pressures or species mass fractions is through a GPTS approach. For that, we define a penalty potential as The weak formulation of (43) follows as with corresponding test functions δφ b v and δφ v as above. The major benefit of this GPTS approach is that it can easily be combined with the embedded multiscale method by adding (44) to the weak form of the 1D-3D problem which is already two-way coupled through the mass transfer terms with similar form to (44), that is, Equations (41) and (42). Indeed, the coupling term may be interpreted as a mass exchange term between preexisting and NV with very large permeability ϵ GPTS such that pressures and species in the two domains immediately equalize. We have termed it GPTS approach here since evaluating (44) in its discretized form will ultimately require Gauss integration of the term along the vessel centerline, see also Section 3.4. The drawback of such GPTS methods is the choice of the penalty parameter ϵ GPTS . A too low value will not fulfill the constraint with sufficient accuracy, while a too high value will lead to an overconstrained problem, see also Section 4.1.

| Mortar approach with penalty regularization
An alternative constraint enforcement scheme, which does not suffer from the aforementioned shortcomings, is the Lagrange multiplier (LM) method whose contribution to the weak form can be written as Again, this allows interpreting the (scalar-valued) Lagrange multiplier λ defined on the embedded domain Λ t as a mass exchange term between embedded vasculature and NV, ie, Typically, a LM formulation transforms the resulting system of equations into a saddle point problem meaning that its solution will be a maximum with respect to the LMs and a minimum with respect to the primary variables. In the context of solid mechanics mesh tying and contact algorithms, such a formulation with a discretization of the LM field is usually termed mortar method. 73 The nodal Lagrange multipliers enter the system of equations as additional unknowns or are condensed out via a dual approach. 74,75 However, here we pursue a different strategy with a penalty regularization of the LM method similar to Yang et al, 76 which we term mortar penalty method in the following. The spatial discretization of (46) may be written as with the well-known mortar matrices and In the spatially discretized form of (47), we have introduced a nodal LM λ j as well as nodal variations of the primary variables δφ b v k on the discretized 1D domain Λ t,h and nodal variations of the primary variables δφ v l in the discretized 3D domain Ω t,h . The indices j and k identify nodes on the discretized 1D domain, while l is employed for nodes in the 3D domain. S and M denote the subsets of nodes on the 1D domain, respectively, in the 3D domain, which actually form the discrete interface Λ t,h . Shape functions in those domains are denoted with b N k and N l , respectively. Furthermore, the LM shape function at a node belonging to the 1D discretization is b Φ j . We assume that every node of the 1D domain carries a nodal LM with corresponding shape function, which has to be chosen from a suitable function space to ensure inf-sup stability. 73 In the following, we exclusively apply linear shape functions for both primary variables and Lagrange multiplier interpolation on the 1D domain, that is, b For the penalty regularization of the constraint, we define the discretized weighted pressure or species gap at node j of the discretized 1D domain as with As in Yang et al, 76 the scaling with the inverse of (51) guarantees consistent units of the weighted gap g j . It can be employed for a penalty regularization of the nodal Lagrange multiplier λ j as λ j ¼ ϵ MP ·g j (52) with (mortar) penalty parameter ϵ MP . This definition is then inserted into (47) to eliminate the LMs from the global problem. Again, the penalty parameter determines how accurately the constraint will be fulfilled. However, as we will show in Section 4.1, the MP method is not prone to overconstrainment as opposed to the GPTS scheme.
Finally, we want to emphasize that the dimensions of the mass transfer equations defined on the embedded domain and the surrounding tissue are different. The balance of mass of blood (6) and species (10) are written in terms of [length] 2 /[time], while the ones of the 3D tissue domains (14), (18), (19), and (22) are formulated in terms of 1/[time]. The coupling between the equations is achieved via the Dirac term δ Λ t , whose dimensions are actually 1/[length] 2 . 41 Hence, the units of the penalty parameters ϵ GPTS and ϵ MP are [length] 2 /[time·pressure] for the coupling of pressures and [length] 2 /[time] for the coupling of species, which makes their interpretation as permeabilities straightforward. From a modeling point of view, such finite permeabilities might even be advantageous. For instance, they could be employed to account for the partition of red blood cell flux between preexisting and NV. Red blood cells might preferentially follow the higher blood flow in the preexisting, developed vasculature such that hematocrit and correspondingly the oxygen mass fraction in the NV is lower. 50 Remark 3. Alternatively, the formulation of the mesh tying constraint (27) could also be written w.r.t. the averaged quantity from the 3D domain φ v according to (20). Then, the weak form for both schemes would have been equivalent to the formulation of Köppl et al 58 with a mass transfer term between preexisting and NV with large permeability living on the boundary of the inclusion.

| Discretized system of equations and computational algorithm
After space and time discretization, the nodal primary variables defined on the discretized embedded domain Λ t,h at that is, nodal pressures and nodal species mass fractions. Here, a total of b n spec species is transported in the preexisting vasculature, which is discretized with b n nodes nodes. For the surrounding tissue domain, the respective nodal primary variables are defined by d s nþ1 ∈ R n nodes ·n dim ; ψ t;h;l nþ1 ∈ R n nodes ·3 ; ϵ v nþ1 ∈ R n nodes ; p v nþ1 ∈ R n nodes ; ω nþ1 ∈ R n nodes ·n spec (54) as nodal displacements, nodal generic primary variables of TCs, HCs, and the IF, the nodal NV volume fraction, the nodal NV pressure, and nodal species mass fractions with n spec species in the multiphase model. Furthermore, n dim denotes the number of spatial dimensions and n nodes the number of nodes in the tissue domain. We can then define a coupled system of discrete nonlinear residuals R at time step n+1 as where the primary variable of each equation has been underlined. Equation (55) corresponds to the Hagen-Poiseuille equation for blood flow in the preexisting vasculature (6), while (56) constitutes the discretized form of species transport in the preexisting vasculature (10). Furthermore, the balance of linear momentum (26) is entering the system through the discrete residual (57). The multiphase system with the evolution equation for the NV volume fraction (12), the balance of mass of blood in the NV (14), and the balance of mass of TCs, HCs, and the IF (18) and (19) is represented by (58). Finally, species transport (22) in the multiphase system is described through (59). In our previous publication, 10 we have shown that a fully monolithic algorithm is the most efficient choice to solve the strongly coupled problem of the vascular multiphase model. Hence, we will exclusively apply a monolithic solution algorithm with a single Newton-Raphson loop per time step to the coupled problem (55) to (59). For that, the nonlinear system is fully linearized with the exception of the first two residuals w.r.t. to the solid phase displacement d s , ie, the deformation of the 1D network following the underlying porous medium. Neglecting this term did not show any significant influence on the convergence of the Newton scheme for typical deformations occuring during tumor growth in the model. A linear system of equations with 5×5 block structure corresponding to the five discrete residual blocks has to be solved for each Newton step. To apply a standard GMRES iterative solver, sophisticated preconditioners are required. For the TCAT tumor model, we reuse the ones developed by Verdugo and Wall, 77 more specifically the AMG-BGS variant, which proved more efficient when coupling the model with the embedded blood vessel network than the previously employed BGS-AMG scheme. 10

| Remarks on the implementation
In general, we study the embedding of 1D inclusions into 3D domains. However, for certain examples, we can also consider an interaction of 1D inclusions with 2D domains. After spatial discretization, our nonconforming 1D-2D/3D coupling requires the numerical integration of products of shape functions defined on the embedded discretized domain Λ t,h . These can be either shape functions on the 1D domain or the 2D/3D domain. Also, the MP approach requires a similar integration of LM shape functions and shape functions of the surrounding domain, comp. (48) and (49). We have chosen to apply a so-called segment-based integration to minimize the error of the numerical integration. 78 The general concept is sketched in Figure 4. The considered 1D element interacts with three elements of the 2D domain. Hence, the integration is split into three segments, where Gauss points are defined on the 1D domain. This avoids integration over weak discontinuities of shape functions across the (2D/3D) element boundaries. The proposed approach involves projecting a Gauss point of the 1D element into the 2D/3D domain to find its counterpart in the parameter space of the interacting 2D/3D element. If the surrounding tissue domain is discretized with a regular grid and birespectively trilinear shape functions, the mapping of Gauss points is trivial. However, if distorted elements are present, it might become nonlinear. This can especially occur during deformation where initially straight blood vessel segments become distorted. In that case, a local Newton algorithm is applied to map Gauss points from the 1D discretization into the 3D mesh.
Another important aspect is the treatment of the equations for blood flow (14) and species transport in the NV ((22) with ε α =ε v ). Since the NV only develops during the simulation, there are certain areas in the computational domain without any NV, ie, ε v =0. Naturally, the two aforementioned equations are not valid here but only in the part of the domain with ε v >0. In our implementation, we decide element-wise if the corresponding equations can be solved if the NV volume fraction inside the element from the last time step is bigger than a threshold value of ε v =0.01. Only in these elements the aforementioned equations are solved while in the rest of the domain they are not evaluated. This is equivalent with a no-flux boundary condition at the interface between the areas in which NV is already present and where angiogenesis has not yet occurred. Considering the underlying physical problem, this assumption states that there is no flux of species or fluid across the sprout tips during angiogenesis, which seems reasonable. The element-wise treatment is only a very crude approximation of the actual interface but showed to be sufficient here.
We want to conclude this section with some qualitative considerations about the computational cost of a fully resolved vasculature model as compared with a hybrid or homogenized approach. First of all, it is important to emphasize that our proposed coupling does not introduce any significant computational overhead into the formulation since the domains are already coupled through transcapillary exchange terms along the blood vessel network. In terms of computational cost of a homogenized formulation compared with a fully resolved blood vessel network, several factors play a role such as the relative cost of evaluating a 1D element as compared with a 3D homogenized element, the cost of the coupling method between the two distinct domains, and how many 1D elements can potentially be replaced by a homogenized description with a 3D element. It is hard to argue which effects dominate without performing a comparison which was not the main interest of this contribution. Nevertheless, one potential major computational advantage of a homogenized or hybrid approach persists, which is the fact that it could possibly eliminate the need of several simulation runs on the same setup, eg, Welter and Rieger 25,26 performed 15 simulation runs to obtain information about microvascular densities and other averaged quantities.

| NUMERICAL EXAMPLES
The primary aim of this paper is the formulation of the coupling between the TCAT tumor growth model and the embedded multiscale method, especially our proposed hybrid treatment of angiogenesis. A detailed validation with experimental data is beyond its scope. However, the following numerical examples have been designed to demonstrate the principal applicability of the novel model to vascular tumor growth.

| Comparison between GPTS and MP constraint enforcement
A validation and comparison of the GPTS and MP constraint enforcement strategies is performed in this section. For that, we have simplified the model by neglecting any deformability and solving only the equations which is a simplified version of (6), and which is a simplified version of (14), together with the constraint ie, the constraint equation (27) for pressures. Figure 5 depicts the computational domain which is used to solve these equations. A blood vessel is embedded into a porous block with dimensions 1×1×1. The straight blood vessel with radius  Figure 5, and on the upper right end, the pressure is set to zero. The 3D block only carries no-flux boundary conditions. It is discretized with a regular grid of 10×10×10 trilinear elements while the 1D domain is discretized with 23 equally spaced linear elements such that the two discretizations are nonconforming. Figure 5 illustrates the pressure distribution in both domains. For a more detailed investigation, we have plotted the pressures on the 1D domain p b v and in the 3D domain p v along the inclusion Λ 0,h in Figure 6 for the two different methods. Please note that the chosen penalty parameters are at least two orders larger than the characteristic order of magnitude defined by the permeability k v , which leads to good fulfilment of the constraint for all depicted values.
The advantage of the MP formulation compared with the GPTS approach becomes evident from these plots. For very large values of the penalty parameter, the latter formulation becomes overconstrained due to a lack of inf-sup stability. When increasing the penalty parameter, it converges toward a collocation method where the discrete constraint has to be fulfilled at every Gauss point. In this case, the number of discrete constraints is too high for the number of discrete degrees of freedom. 79 No convergence toward a physically reasonable result can be expected for ϵ GPTS →∞ due to the overconstraint. Indeed, if the penalty parameter is chosen even higher than the values shown in Figure 6A, a linear pressure drop from inflow to outflow emerges since this is the only solution which satisfies the constraint (62) at every Gauss point. By contrast, the MP method converges to the solution of the LM method. No visible difference in the result is present for penalty parameters larger than ϵ MP =0.1 (not shown here). Still, in the range of moderate penalty parameters, both methods perform with similar accuracy such that the GPTS approach seems applicable especially when a very high accuracy for fulfilling the constraints (27) for species and pressures is not needed. Indeed, from an implementational point of view, the GPTS approach can also be integrated more easily into the embedded multiscale method since the additional terms can be evaluated on element level together with the mass exchange terms. By contrast, the penalized LM method requires a global assembly of the mortar matrices and only then the LMs can be eliminated. Nevertheless, we will apply the MP method in the following due to its theoretical advantages outlined in this section.

| Two-dimensional growth of a tumor close to a preexisting blood vessel
We study the growth of an initially avascular tumor which is located close to a preexisting blood vessel from which angiogenesis can occur. By means of this simple case, we want to exemplify the coupling between the discrete preexisting vasculature and the NV and how angiogenesis can be triggered from the preexisting network in our model. The 2D domain with the embedded simple blood vessel network is sketched in Figure 7 where also initial and boundary conditions are listed. This setup corresponds to a vessel-ingrowth scenario where the whole tumor vasculature grows from outside into the tumor. 16 While this situation with a distance of the initial tumor spheroid of 0.4 mm to the nearest blood vessel is nonphysiological, the numerical example which has been designed similar to corresponding discrete models 13,18,22 enables us to methodologically investigate our proposed approach. For this case, we consider the full five-phase model as described in Section 2 including the proposed hybrid embedded/homogenized treatment of the vasculature. The primary variables for TCs, HCs, and IF are chosen as p th , p hl , and p l , that is, the pressure difference between TCs and HCs, the pressure difference between HCs and the IF, and the IF pressure. In addition, four species are present, namely oxygen in the IF, oxygen in the vasculature (preexisting and NV), TAF in the IF, and NTCs as a part of TCs. At the beginning of the simulation, the tumor covers a circular area with radius 0.05 mm. The outline of the 2D domain is fixed but the tissue inside is assumed to be deformable. At the inflow on the left end of the 1D inclusion, we set a Dirichlet boundary condition for the pressure in the preexisting vasculature p b v and the oxygen mass fraction in the preexisting vasculature ω nb v . At the two outflows of the vessel domain, the pressure is also fixed. Angiogenesis occurs from the preexisting blood vessel network which has a constant radius of R=0.015 mm. For that, we assume a boundary condition of ε v =0.1 along the 1D network. To trigger angiogenesis from the location of the preexisting blood vessel, this value is applied on all 2D elements which are "cut" by the blood vessels. The pressure in the NV p v and the oxygen mass fraction in the NV ω n v are coupled with the values in the 1D domain with the mortar penalty approach developed in Section 3.2.2. We employ a penalty parameter ϵ MP ¼ 1·10 −10 m 2 = Pas ð Þ for the coupling of pressures and ϵ MP ¼ 1·10 −5 m 2 =s for the coupling of species. Both domains are discretized in space completely independent from each other, which is the major advantage of our embedded coupling. The 2D domain is meshed with a regular grid of 210×140 bilinear elements. The 1D vasculature is represented by 805 linear elements. We consider the growth of the tumor over 24 days and employ a one-step-theta scheme with θ=0.5 and a time step size of Δt=900 s.
All parameters of the model are listed in Tables 2 to 6. They have been taken either from previous contributions on the TCAT model or from values reported in the literature. Some parameters have been estimated, which is highlighted through the footnotes in the corresponding tables. The ECM is modeled with a Neo-Hookean material law, whose parameters are summarized in Table 2.  The example of this section has been designed to imitate the classical model for the angiogenic switch. 62 If we did not allow angiogenesis to occur from the preexisting blood vessel, the initial tumor of this example would not grow any further since a steady state of proliferation and TC death due to the hypoxic conditions is reached. By contrast, in the case studied here, the onset of angiogenesis enables rapid tumor growth. The evolution of the LTC volume fraction is shown in Figure 8 as well as the NV volume fraction ε v which is visualized with white contour lines. The preexisting blood vessels are plotted in red. After 6 to 12 days, a slightly unsymmetrical growth toward the preexisting blood vessels can be observed. Oxygen is provided by the preexisting network and diffuses in the IF, which is why the tumor grows toward the blood vessel with higher nutrient availability. A necrotic region develops in the part of the tumor further away from the vessels. Its outline is shown in Figure 8 by the black contour line. Concurrently, the Effective solubility of oxygen in blood 27 Discharge hematocrit 27 Concentration of oxygen at maximum saturation 27 C nv Hill exponent 69 n 2.7 - Value for a radius of 0.015 mm according to the fit employed by Welter et al, 27 Appendix S1.
b Since we regard the neovasculature in a homogenized way, we assume that the average radius of the blood vessels in the neovasculature is 0.01 mm. The value for γ tv is then calculated with the fit employed by Welter et al, 27 Appendix S1. c This value is a simple approximation. It has been obtained from the value of oxygen diffusion in blood plasma of D p ¼ 2:75·10 −9 m 2 =s. 27 In our model, the mass fractions ω n v resp. ω nb v represent both the oxygen dissolved in plasma and bound to hemoglobin. At an oxygen partial pressure of P v oxy ¼ P v oxy; 50 ¼ 37 mmHg, the ratio of the mass fraction of dissolved oxygen to the mass fraction of the total oxygen present is approximately 0.01. Hence, we have scaled the value for D p 27 by this factor to only include diffusion of the oxygen dissolved in the plasma and not the one bound to hemoglobin. Therefore, our oxygen transport model is similar to the one of Beard and Bassingthwaighte 81 and Fang et al. 82 hypoxic TCs constantly produce TAF which triggers endothelial cell migration toward the tumor from the preexisting vasculature where the boundary condition on the NV volume fraction has been set. Please note that a small NV volume fraction is also present below the preexisting blood vessel due to the diffusive term in the formulation for Osmotic pressure difference 20 ω osm π blood − π l À Á 1333 Pa (31), A9) Hydraulic conductivity for transcapillary flow 83 Surface-to-volume ratio for transcapillary flow 83   In earlier publications, 10,63 a value of D TAFl =3.5·10 −4 was used, which is the nondimensional value employed by Anderson and Chaplain. 13 Here, we employ the correct dimensional value. Therefore, the production rate of TAF under hypoxia γ TAF production has been reduced accordingly compared to Kremheller et al. 10  angiogenesis (12). After 18 days, both angiogenesis and tumor growth have continued such that the NV has reached the region with a high proliferation of TCs. A characteristic bulge of the NV volume fraction toward the tumor starts to develop since the TAF concentration is highest there. The higher availability of oxygen enables rapid tumor progression in the interval between 18 and 24 days. After 24 days, the tumor has grown in a half-moon shape toward the preexisting blood vessels which is also visualized by the motion of its center of mass. In the period from 18 to 24 days, it moves with an average velocity of 0.013 mm/d in direction of the vessel. At the same time, the increased TC and IF pressure in the LTC region leads to a deformation of the ECM. Under our assumption that the blood vessel network completely follows the movement of the underlying ECM, this induces also a slight deformation of the initially straight blood vessels which are, in this case, slightly pushed away from the tumor.
To elucidate the coupling between preexisting vasculature and NV, pressures in the preexisting vasculature p b v and in the NV p v after 24 days are depicted in Figure 9. Clearly, the constraint of equal pressures on the embedded domain Λ t,h is fulfilled such that flow in preexisting and NV are coupled. As stated in Section 3.4, the balance of mass of blood in the NV (14) is only valid in the area the NV has already reached, which explains the shape of the pressure distribution p v in Figure 9 including the area below the blood vessel. Only those elements where we can actually solve the aforementioned equation are depicted. In the rest of the domain, the equation is not evaluated, which is equivalent to a no-flux boundary condition at the edge of the NV. A considerable pressure drop from the preexisting vasculature toward the edge of the NV can be observed which is due to the large leakage of fluid from the NV into the IF.
Also, species transport of oxygen in the preexisting and NV are coupled via the MP approach. The corresponding distributions are shown in Figure 10. Again, species transport of oxygen in the NV can only be solved in the portion of the domain where NV is present. As described in Section 3.4, the respective equation is only evaluated in elements with a threshold value of ε v >0.01 resulting in a no-flux boundary condition across the edge of the NV. The proposed approach worked well in capturing the shape of the NV domain in this case, but for different scenarios, the threshold value could be defined as a certain percentage of the degree of vascularization over the whole domain. Transcapillary exchange of oxygen from the preexisting and the NV into the IF cause the depicted oxygen distributions. In the tumor area, a lot of oxygen is required due to the vast oxygen consumption by proliferating TCs. At the beginning of tumor progression, oxygen is provided from the preexisting vasculature into the IF through transcapillary exchange. Diffusion in the IF is necessary to reach the tumor. During angiogenesis, the developing NV enables a more efficient transport of oxygen. Oxygen can now also be transported from the preexisting into the NV due to the coupling of species transport and then from the NV into the IF much closer to the site of the tumor. This behavior is augmented by the pressure drop in the NV, which develops due to the large leakage of fluid, such that oxygen is advected from the inflow through the vasculature. Still, the high oxygen demand leads to very low oxygen mass fractions in both the NV and the IF, see also Figure 11 (left). Also in the IF, most of the oxygen is present close to the preexisting vasculature while a hypoxic region emerges inside the tumor. In the right part of Figure 11, also the TAF mass fraction in the IF which is produced during hypoxia by LTCs and triggers angiogenesis is shown.
With this numerical example for 2D vascular tumor growth, we have demonstrated that our treatment of angiogenesis with a homogenized representation of the NV and a discrete representation of the preexisting vasculature is able to produce biologically and physically reasonable results. The information about the structure of the preexisting vasculature can be preserved, while not every single capillary segment of the NV has to be resolved as in most other hybrid vascular tumor growth models. At least qualitatively, the results are in good agreement with ingrowth simulations using

FIGURE 10
Mass fraction of oxygen in preexisting (left) and neovasculature (right) after 24 days discrete vasculature models. 18,22 It becomes evident how angiogenesis completely shifts the supply of the tumor with nutrients from diffusion in the IF to the more efficient transport in the NV. One drawback of the model is the difficulty to initiate sprouts from the preexisting vasculature. For now, we have set a Dirichlet boundary condition on the preexisting vasculature to trigger angiogenesis.
For this example, preexisting and NV are of comparable scale. Hence, our motivation for representing one as a 1D inclusion and the other in a homogenized way is not scale separation but rather structure/morphology and function of tumor NV, as described in the introduction. Due to these complexities and the associated heterogeneity, we believe that our homogenized approach for modeling the NV with a blood vessel density rather than single blood vessel segments is more suitable than a resolved approach. Naturally, quantities of interest such as microvascular density or average blood flow and species transport in the NV are available. Another possible application of the formulation apart from angiogenesis might be the modeling of larger tissue domains at the scale of a whole organ. Here, it might also be possible to model only the larger vessels patient-specifically as 1D inclusions and treat the smaller blood vessels, ie, the capillary bed in a homogenized sense as in the composite smeared finite element method. 36 In contrast to their approach, the major advantage of our formulation is that the 1D and the 3D grid can be completely independent.
Compared with the experimental results shown for brain tissue by Seano et al, 44 the deformation of the preexisting vasculature during tumor growth is actually quite small due to the chosen material parameters and boundary conditions. Hence, it might be justified to neglect the deformability of blood vessels in this case and simply evaluate the respective terms in reference configuration. Especially, the solid velocity term in (6) is insignificant since time scales during tumor growth are quite large.

| Three-dimensional growth of a tumor along a preexisting blood vessel network
The final example for the coupling of the TCAT tumor growth model with the embedded 1D fluid network is the growth of a tumor inside a 3D preexisting blood vessel network. In the previous section, we have considered a case inspired by the classical model of angiogenesis where the tumor first grows avascularly and then angiogenesis occurs from nearby vessels. However, a different mechanism of growth is possible for certain tumor types such as astrocytomas. 62,85,86 They can first acquire access to blood circulation by co-opting preexisting blood vessels and growing along them, which makes them a non-angiogenic but nevertheless well-vascularized tumor. Subsequently, the preexisting vasculature regresses such that a necrotic core inside the tumor evolves due to lack of nutrients. Only then angiogenesis at the tumor boundary is initiated to enable further tumor growth. 86 We will show how our model can capture this behavior.
For that, we consider the computational domain of Figure 12. The depicted blood vessel network has been obtained from R3230AC mammary carcinoma in rat dorsal skin flap preparation by Secomb et al 87 and has been employed in several other publications to study oxygen transport, 67 drug delivery, 40,41 and hyperthermia treatment. 42,43 The network topology is publicly available. 88 We have enlarged the enclosing 3D tissue domain in z-direction by 0.056 mm on the top and bottom of the domain as compared with the dimensions of the network geometry. The 3D domain is discretized uniformly in space with 55×52×30 trilinear elements and the 1D network with 8298 linear elements. Again, the meshes of both domains are completely independent. Note that the 1D network is discretized rather finely compared with the 3D domain to avoid instabilities due to convection-dominated oxygen transport in this example. Following Nabil et al, 42,43 we assume a uniform radius of R=0.00764 mm in the entire network. The growth process is simulated for a time period of 360 hours. For time discretization, the one-step-theta scheme with θ=0.5 and a time step size of Δt=1800 s is applied.
For this example, angiogenesis is not present because we only want to simulate the first growth stage along the blood vessel network as described above. For the TCAT multiphase system, this reduces the model to a three-phase model of TCs, HCs, and the IF along with two species, namely NTCs and oxygen in the IF. In the 1D embedded domain, we solve for blood pressure and oxygen mass fraction. Furthermore, deformations of the surrounding ECM and, hence, also the blood vessels are neglected. Consequently, only a simplified version of the model without angiogenesis and structural deformation, that is, without Equations (12), (14), and (26) is considered here. Nevertheless, the model is still fully coupled with the embedded vasculature via the appropriate exchange terms of Table 1. Only the terms for the NV, oxygen in the NV, and TAF are not considered because these phases and species are not taken into account. Moreover, no coupling between NV and preexisting blood vessels is performed since angiogenesis is not present. Boundary and initial conditions for the example are given in Figure 12. Initially, a spherical tumor with radius r 0 =0.03 mm around the point [0.25,0.25,0.11] is present inside the domain. Primary variables for TCs, HCs, the IF, and for necrotic TCs are fixed on the entire boundary of the cuboid Γ Ω . Nine open ends of the network are identified as inflows Γ Λ,i (denoted by red arrows) where pressure and mass fraction of oxygen in the 1D network are fixed. The remaining eight open ends are outflows Γ Λ,o (denoted by blue arrows) with fixed pressure. This approach has been proposed to prescribe the pressure drop along the network such that a physiologically reasonable blood velocity is obtained. [40][41][42] The parameters employed in this example are the same as given in Tables 2 to 6 apart from the coefficient for transvascular oxygen exchange and the coefficients for oxygen consumption. The factor for transvascular oxygen exchange 27 is simply calculated for the given radius as b γ tv ¼ 1:61·10 −8 m=ðmmHg sÞ. In addition, we have elevated the oxygen consumption rate w.r.t. the example from the previous section under the assumption that the tissue considered here is well-vascularized, which implies a heavy oxygen or nutrient demand. Therefore, the consumption by TCs has been increased to γ nt growth ¼ 9:6·10 −4 kg=ðm 3 sÞ and γ nt 0 ¼ 2:4·10 −3 kg=ðm 3 sÞ and for HCs to γ nh 0 ¼ 1·10 −3 kg=ðm 3 sÞ compared with the values of Table 5. Note that all the values applied for oxygen consumption lie in the physiological range of oxygen consumption rates of 0.01−0.3 min −1 (in terms of oxygen concentration). 89 The results for tumor growth over the considered time span are visualized in Figure 13. As stated above, we aim to demonstrate that our model can capture the initial blood vessel co-option by TCs during growth. However, we have not included the compression by excessive TC pressure or blood vessel regression. Nevertheless, the influence of these phenomena on the growth of the tumor can be imitated if the exchange terms between preexisting vasculature and IF, that is, the mass transfer terms for leakage and oxygen exchange are switched off once the tumor has reached a specific embedded 1D element. This approach has been applied to obtain the results of Figure 13. Initially, the tumor grows radially. Then, a necrotic core starts to appear due to lack of oxygen in its interior since nutrient diffusion from the outside to the inside of the tumor is a limiting factor. Nevertheless, tumor growth can continue along the vasculature FIGURE 12 Geometry of the example for 3D vascular tumor growth by co-option of the blood vessels because this is the region with the highest oxygen concentration in the IF which is provided by the vasculature there and consumed by TCs during growth.
The goal of this example was to highlight the capabilities of the nonconforming coupling to represent the complexity of in vivo vessel networks and their interaction with the TCAT multiphase tumor growth model. Here, the major advantage of nonmatching meshes for the 1D and the 3D domain becomes evident. Therefore, the proposed methodology for embedding arbitrary blood vessel networks into surrounding tissue could in principle also be applied for fully resolving tumor vasculature. However, we are aware that this is only the first step toward a more realistic description of in vivo tumor growth. For instance, we have excluded the second stage of angiogenesis on the tumor boundary occurring after blood vessel regression. In addition, blood vessel regression, compression by TCs, and adaption is not modeled satisfactorily yet.

| CONCLUSION
Capturing the complexity of cancer through mathematical models requires cutting-edge numerical algorithms. In this contribution, we have introduced a holistic framework which allows for studying vascular tumor growth. The preexisting blood vessels are represented as 1D inclusions in the surrounding 3D tissue domain. In the 1D vessels, blood flow and species transport as well as adequate mass exchange with the enclosing tissue are considered. A major novelty is the coupling between blood flow and species transport in a smeared or homogenized representation of the vasculature and the resolved (1D) part of the domain. This has been realized by two different penalty-based approaches, a GPTS, and a mortar scheme with penalty regularization, which allow for dissimilar meshes on both domains, the main advantage of the embedded multiscale method. In principle, both methods work but if a higher accuracy is necessary, the LM scheme should be the method of choice especially also because of its sound mathematical properties. Our nonconforming approach can not only be applied to vascular tumor growth but might have numerous applications in other mass transfer problems of biological interest as well, where a full resolution of the capillary bed is not necessary. The main problem for which we have applied the developed methods is vascular tumor growth. The complexity of the tumor NV formed during angiogenesis is enormous. Together with its structural and functional abnormalities, a homogenized treatment of blood flow and species transport lies at hand. Discrete or fully resolving angiogenesis models claim to provide more insight into the specific network structure. We question that this is necessary because quantitities of interest are rather the degree of vascularization as expressed by a (neo-)vasculature volume fraction, information about regions which are especially highly or badly vascularized or (time-and space-)averaged information about blood flow and species transport in the tumor micro-environment. Our model inherently offers these quantities while discrete models require a similar kind of homogenization via an average over an ensemble of simulations. Qualitatively, our results are comparable with earlier fully resolved models. Nevertheless, it might be interesting to compare results between our continuous representation of the NV with corresponding discrete models. An advantage of the latter ones is that vessel remodeling and regression can be treated quite naturally, while those significant phenomena are still missing in homogenized formulations. Our hybrid approach is a sensible compromise because the geometry of the preexisting vasculature is still contained in the model while we do not aim to fully resolve every capillary of the NV. While information about larger blood vessels and transport therein might be available through suitable imaging techniques which would enable a patient-specific model, the smallest vessels can usually only be characterized through their vascular density and average transport therein. To make the model applicable to larger scales with a dense capillary network, it might also be possible to model the smaller scales of the capillary bed and the NV with the homogenized approach and couple it to the larger, resolved vessels with the methods developed here.
We have shown the capabilities of the TCAT tumor growth model including the proposed hybrid embedded/homogenized approach for the study of vascular tumor growth in silico with three illustrative examples. The first one has been employed to study our coupling schemes. The second one is based on the classical model of an angiogenic switch where at first avascular growth is succeeded by angiogenesis from preexisting blood vessels. We have demonstrated how angiogenesis deregulates mass transport in the tumor micro-environment in favor of rapid tumor growth. The third example has been tailored to simulate co-option of preexisting vasculature by invasive tumor growth in well-vascularized tissue. However, these examples have not been validated with experimental data but prove only that we can phenomenologically model tumor growth under these circumstances. A detailed verification and validation with experimental data will be the main focus of further research. A logical next step is also the inclusion of drugs or nano-particles as additional species transported by the (neo-)vasculature and the IF to study drug delivery with the model.

ADDITIONAL CONSTITUTIVE EQUATIONS AND MASS TRANSFER TERMS
For the sake of completeness, we summarize all constitutive relations and mass transfer terms in this appendix. The volume fraction of the pores of the ECM (15) can be derived from the balance of mass of the solid phase in reference configuration as 10