Modeling and simulation of vascular tumors embedded in evolving capillary networks
Marvin Fritz, Prashant K. Jha, Tobias Köppl, J. Tinsley Oden, Andreas Wagner, Barbara Wohlmuth
MMODELING AND SIMULATION OF VASCULAR TUMORSEMBEDDED IN EVOLVING CAPILLARY NETWORKS
MARVIN FRITZ , PRASHANT K. JHA , TOBIAS K ¨OPPL , ∗ , J. TINSLEY ODEN ,ANDREAS WAGNER , AND BARBARA WOHLMUTH , Department of Mathematics, Technical University of Munich, Germany Oden Institute for Computational Engineering and Sciences, The University of Texas at Austin, USA Department of Mathematics, University of Bergen, Allegaten 41, 5020 Bergen, Norway
Abstract.
In this work, we present a coupled 3D-1D model of tumor growth within a dynamically changingvascular network to facilitate realistic simulations of angiogenesis. Additionally, the model includes ECMerosion, interstitial flow, and coupled flow in vessels and tissue. We employ continuum mixture theorywith stochastic Cahn–Hilliard type phase-field models of tumor growth. The interstitial flow is governedby a mesoscale version of Darcy’s law. The flow in the blood vessels is controlled by Poiseuille flow, andStarling’s law is applied to model the mass transfer in and out of blood vessels. The evolution of thenetwork of blood vessels is orchestrated by the concentration of the tumor angiogenesis factor (TAF) bygrowing towards increasing TAF concentration. The process is not deterministic, allowing random growthof blood vessels and, therefore, due to the coupling of nutrients in tissue and vessels, stochastic tumorgrowth. We demonstrate the performance of the model by applying it to a variety of scenarios. Numericalexperiments illustrate the flexibility of the model and its ability to generate satellite tumors. Simulations ofthe effects of angiogenesis on tumor growth are presented as well as sample-independent features of cancer.
Keywords: tumor growth, 3D-1D coupled blood flow models, angiogenesis, finite elements, finite volume1.
Introduction
In this work, we present new computational models and algorithms for simulating and predicting a broadrange of biological and physical phenomena related to cancer at the tissue scale. We consider the growth ofvascular tumors inside living tissue containing a dynamically evolving vasculature. One of the main goalsof this work is to provide realistic simulations of the vascular growth characterizing angiogenesis, wherebyblood vessels sprout and invade the domain of the tumor mass when prompted by concentrations of proteinscollectively referred to as tumor angiogenesis factors (TAFs); these proteins are produced by nutrients-starved cancerous cells. The tumor growth is necessarily depicted by a multispecies model in which tumorcell concentrations are categorized as proliferative, hypoxic, or necrotic. To capture the complex interactionof cell species and the evolving interfaces between species, continuum mixture theory is used as a frameworkfor constructing mesoscale stochastic phase-field models of the Cahn–Hilliard type. Other critical phenomenaare also addressed by this class of models, including the erosion of the extracellular matrix (ECM) due toconcentrations of matrix-degenerative enzymes (MDEs) and the invasion of tumor cells through the ECMas a prelude to metastasis.The volume of an isolated colony of tumor cells will not generally grow beyond approximately the 1mm cube [28,39,43] unless sufficient nutrients and oxygen are supplied for proliferation. To acquire such nutrients,cancerous cells promote angiogenesis [8, 44]. Low levels of oxygen and nutrient result in tumor cells enteringthe hypoxia phase during which they remain dormant and release various proteins such as TAFs that promotethe proliferation of tip endothelial cells and new vessel formation, and MDEs that erode the ECM allowingthe invasion of tumor cells into the surrounding tissue and blood vessels. Because angiogenesis is one of the E-mail address : { fritzm, koepplto, wagneran, wohlmuth } @ma.tum.de, [email protected], [email protected] .2020 Mathematics Subject Classification.
Key words and phrases. tumor growth, 3D-1D coupled blood flow models, angiogenesis, finite elements, finite volume. ∗ Corresponding author. a r X i v : . [ q - b i o . T O ] J a n M. FRITZ, P. K. JHA, T. K ¨OPPL, J. T. ODEN, A. WAGNER, AND B. WOHLMUTH major processes through which tumors grow, anti-angiogenic drugs that inhibit the formation of the newvascular structure are often identified as one of the approaches to delay or arrest the growth of cancer. Thus,a realistic model of angiogenesis is of critical importance for studying the effectiveness of anti-angiogenicdrugs.Typically, the vasculature near the tumor core in the early stages of tumor growth may not effectivelysupply nutrients to the tumor. The vasculature evolves rapidly and, therefore, the vessel walls are notfully developed and may be destroyed due to pressure (proliferation of tumor cells result in higher pressurenearby), the pruning of vessels due to insufficient flow for a sustained period and, vasculature adaptationand remodeling [42, 47–49, 56]. Highly interconnected and irregular vasculature with inefficient blood vesselscauses low blood flow rates to the tumor, making it possible that therapeutic drugs miss the tumor massaltogether [56]. Shear and circumferential stresses due to blood flow result in vascular adaptation effects suchas vessel radii adaptation, see [35, 48, 57]. All of these phenomena are represented by the models describedherein.Earlier models taking into account angiogenesis include lattice-probabilistic network models, see [3, 34,35, 42, 56, 57, 66]. An overview of this class of models is given in [14]. Another class of models referredto as agent-based models has been proposed and extensively studied. There, the idea is to introduce aphase-field for the tip endothelial cells that takes a value 1 inside the vessel and 0 outside and through theagents, which can move anywhere in the simulation domain following certain rules, which can be designed totrigger the sprouting of new vessels; see [33, 45, 59, 61]. These models do not capture blood circulation in thevessel and, therefore, are unable to be truly coupled to the tumor growth. In [63], a dimensionally coupledmodel for drug delivery based on MRI data and a study of dosing protocols is considered with drug flowin the vessels governed by algebraic rules instead of PDEs. More recently, vasculature models involving anetwork of straight cylindrical vessels supporting the 1D flow of nutrient, oxygen, and therapeutic drugs andcoupled to the 3D tissue domain by versions of the Starling or Kedem–Katchalsky law have been presented;see [30, 64, 65].We consider a class of 3D-1D vascular tumor models [19] that approximates the flow within the bloodvessels by one-dimensional flow based on the Poiseuille law; effectively reducing the flow in the three-dimensional vessels to the flow in a network of one-dimensional vessel segments. While coupling the flow inthe vessels and tissue, the blood vessels’ three-dimensional nature is retained by approximating the vesselsas a network of straight-cylinders and applying the fluid exchange at the walls of cylindrical segments.From a mathematical and computational point of view, a complicating factor is the use of one-dimensionalcharacterizations of vessel segments in the vascular network embedded in three-dimensional domains of thetissue and the tumor within the tissue and the assignment of mechanical models to this 3D-1D system todepict interstitial flow and pressure fluctuations. Mathematical analysis showing well-posedness and existenceof weak solutions for the class of 3D-1D model considered in this work is performed in a recent paper [19].In our model, flow in vessels is governed by 1D Poiseuille law, whereas the flow in tissue is derived bytreating the tissue domain as a porous medium and applying a version of Darcy’s law. The model consistsof nutrients in the tissue and vessels; nutrients in the vessels are governed by the 1D advection-diffusionequation and advection-diffusion-reaction equation in the tissue. Flow and nutrients in the tissue and vesselare coupled; we assume that vessel walls are porous, resulting in the advection and diffusion-driven exchangeof nutrients and coupling of the extravascular and intravascular pressures. Some aspects of the 1D modelarchitecture and coupling of 3D and 1D models are based on previous works, see [27, 29, 30]. The 3D tissuedomain includes, in addition to the nutrients, ECM, tumor species such as proliferative, hypoxic, necrotic,and diffusive molecules such as TAF and MDE.As noted earlier, the 3D tumor model is derived from the balance laws of continuum mixture theoryas in [7, 12, 13, 22, 33, 40], and representations of the principal mechanisms governing the development andevolution of cancer, see, e.g., [24,33]. Especially, we note the comprehensive developments of diffuse-interfacemultispecies models presented in [18,62], the ten species models derived in [33], and the multispecies nonlocalmodels of adhesion and promotes a tumor invasion due to ECM degradation described in [20]. Angiogenesismodels embedded in models of hypoxic and cell growth are presented in [19, 33, 64, 65]. Related models ofextracellular matrix (ECM) degradation due to matrix-degenerative enzymes (MDEs) released by hypoxiccell concentrations and subsequent tumor invasion and metastasis are discussed in [10, 11, 16, 20]. Several ofthe earlier mechanistic models of tumor growth focused on modeling the effects of mechanical deformationand blood flow, and fluid pressure on tumor growth, e.g., [1, 2, 5, 6, 31, 46].
ODELING AND SIMULATION OF VASCULAR TUMORS 3
A key new feature of the models proposed here is the dynamic growth/deletion of the vascular networkand full coupling between the dynamic network and tumor system in the tissue microenvironment. Inresponse to TAF generated by nutrient-starved hypoxic cells, new vessels are formed. Due to the formationof new vessels, the local conditions such as nutrient concentration changes near the tumor, affecting TAFproduction and promoting a higher proliferation of tumor cells. The rules by which the network grows, orexisting vessels are deleted due to insufficient flow and dormancy, are based on the experimentally knowncauses of angiogenesis and are parameterized so that various aspects of the network growth algorithm canbe adjusted based on available experimental data. By including the time-evolution of the larger vasculartissue domain and the sprouting, growth, bifurcation, and pruning of the vascular network orchestrated by acombination of blood supply and tumor-generated growth factors, a more realistic depiction of tumor growththan the more common isolated-tumor (avascular) models is obtained.This article is organized as follows: In Section 2, we introduce various components of the model, such asthe tissue and 1D network domain and the equations governing various fields. The details associated with thevessel network growth are presented in Section 3. Spatial and temporal discretization and solver schemes forthe highly nonlinear coupled system of equations are discussed in Section 4. We apply mixed finite volumes,finite element approximations to the model equations. The systems of equations arising in each time stepare solved using a semi-implicit fixed point iteration scheme. In Section 5, the model is applied to varioussituations, and several simulation experiments are presented. Concluding comments are given in Section 6.2.
Mathematical Modeling
In this work, a colony of tumor cells in an open bounded domain Ω ⊂ R , e.g., representing an organ,is considered. It is supported by a system of macromolecules consisting of collagen, enzymes, and variousproteins, that constitute the extracellular matrix. We focus on phenomenological characterizations to capturemesoscale and macroscale events. Additionally, we consider a one-dimensional graph-like structure Λ insideof Ω forming a microvascular network, see Figure 1. Ω φ H φ P φ N Λ EpitheliumBasement membraneStroma
Figure 1.
Setup of the domain Ω with the 1D microvascular network Λ and a tumor mass,which is composed of its proliferative ( φ P ), hypoxic ( φ H ) and necrotic ( φ N ) phases.The single edges of Λ of vessel components are denoted by Λ i such that Λ is given by Λ = (cid:83) Ni =1 Λ i andeach edge Λ i , i ∈ { , . . . , N } , is parameterized by a corresponding curve parameter s i such thatΛ i = { x ∈ Ω | x = Λ i ( s i ) = x i, + s i · ( x i, − x i, ) , s i ∈ (0 , } , where x i, ∈ Ω and x i, ∈ Ω mark the boundary nodes of Λ i , see Figure 2. For the total 1D network Λ, weintroduce a global curve parameter s , which is interpreted in the following way: s = s i , if x = Λ( s ) = Λ i ( s i ).At each value of the curve parameter s , various 1D constituents exist, which interact with their respective3D counterpart in Ω. M. FRITZ, P. K. JHA, T. K ¨OPPL, J. T. ODEN, A. WAGNER, AND B. WOHLMUTH
We introduce the surface Γ of the microvascular network Λ to formulate the coupling between the 3D and1D constituents in Subsection 2.2 and Subsection 2.4. For simplicity, it is assumed that the surface for asingle vessel is approximated by a cylinder with a constant radius, see Figure 2. The radius of a vessel thatis associated with Λ i is given by R i and the corresponding surface is denoted by Γ i ; i.e., we have as the totalsurface Γ = (cid:83) Ni =1 Γ i . Λ i Γ i x i, x i, s i ∂ B R i ( s i ) R i Figure 2.
Modeling of a single edge Λ i of the 1D graph-like structure with boundary nodes x i, and x i, by approximating it by a cylinder Γ i with constant radius R i .2.1. Governing constituents.
The principal dependent variables characterizing the growth and decline ofthe tumor mass are taken to be a set of scalar-valued fields φ α with values φ α ( x , t ) at a time t ∈ [0 , T ] andpoint x ∈ Ω ⊂ R , representing the volume fractions of constituents in the space-time cylinder Ω × [0 , T ].The primary feature of our model of tumor growth is the application of the framework of continuum mixturetheory in which multiple mechanical and chemical species can exist at a point x ∈ Ω at time t >
0. Therefore,for a media with N ∈ N interacting constituents, the volume fraction of each species φ α , α ∈ { , . . . , N } , isrepresented by a field φ α with the value φ α ( x , t ) at ( x , t ) with the property (cid:80) α φ α ( x , t ) = 1.We separate the tumor volume fraction φ T = φ T ( x , t ) into the sum of three phases φ T = φ P + φ H + φ N ,where φ P = φ P ( x , t ) is the volume fraction of proliferative cells, φ H = φ H ( x , t ) that of hypoxic cells, and φ N = φ N ( x , t ) is the volume fraction of necrotic cells, see Figure 1. Proliferative cells have a high probabilityof mitosis, i.e., division into twin cells, and to produce growth of the tumor mass. Hypoxic cells are thosetumor cells which are deprived of sufficient nutrient to become or remain proliferative. Lastly, necrotic cellshave died due to the lack of nutrients.The nutrient concentration and the tumor angiogenesis factor (TAF) over Ω × [0 , T ] are represented byscalar fields φ σ = φ σ ( x , t ) and φ TAF = φ TAF ( x , t ), respectively. The tumor cells response to hypoxia,i.e., when φ σ is below a certain threshold, is the production of an enzyme that increases cell mobility andactivates the secretion of angiogenesis promoting factors characterized by φ TAF . As a particular case ofTAFs, we consider the vascular endothelial growth factor (VEGF), which promotes sprouting of endothelialcells forming the tubular structure of blood vessels, which grow into new vessels and supply nutrients to thehypoxic volume fraction φ H .Another consequence of hypoxia is the release of matrix-degenerative enzymes (MDEs), e.g., urokinaseplasminogen and matrix metalloproteinases, by the hypoxic cells. We denote the volume fraction of theMDEs by φ MDE = φ MDE ( x , t ). The primary feature of the MDEs is the erosion of the extracellular matrix,whose density is denoted by φ ECM = φ ECM ( x , t ). Consequently, the erosion of the ECM produces room forthe invasion of tumor cells, which increases φ T in the ECM domain and therefore, raises the likelihood ofmetastasis. Below a certain level necrosis occurs and cells die, entering the necrotic phase φ N . Tumor cellsmay also die naturally, in a process which is called apoptosis.Within the one-dimensional network Λ, we introduce the constituents φ v = φ v ( s, t ) and v v = v v ( s, t ),which represent the one-dimensional counterparts of the local nutrient concentration φ σ and the volume-averaged velocity v . Additionally, we consider the pressures p v = p v ( s, t ) and p = p ( x , t ) in the network andtissue domain, respectively.2.2. Three-dimensional model governing the tumor constituents.
The evolution of the constituents φ α must obey the balance laws of continuum mixture theory (e.g., see [21,33]). Assuming constant and equalmass densities of the constituents, the mass balance equations for the mixture read as follows: ∂ t φ α + div( φ α v α ) = − div J α ( φ ) + S α ( φ ) , where v α is the cell velocity of the α -th constituent, and S α describes a mass source term that may dependon all species φ = ( φ P , φ H , φ N , φ σ , φ MDE , φ
TAF , φ
ECM ). Moreover, J α represents the mass flux of the α -thconstituent and is given by: J α ( φ ) = − m α ( φ ) ∇ µ α , ODELING AND SIMULATION OF VASCULAR TUMORS 5 where µ α denotes the chemical potential of the α -th species, and m α is its corresponding mobility function.Generally, the mobilities may depend on many species, but in this work we consider the following cases, m α ( φ ) = M α φ α (1 − φ T ) I d , α ∈ { P, H } ,m β ( φ ) = M β I d , β ∈ { σ, MDE, TAF } , where M α are mobility constants, and I d is the ( d × d )-dimensional identity matrix. For the remaining species φ N and φ ECM , we choose m N = m ECM = 0 in accordance to the non-diffusivity of the necrotic cells andthe ECM; see [38]. Following [25, 33, 62], we define the chemical potential µ α as the first variation (Gˆateauxderivative) with respect to φ α of the Ginzburg–Landau–Helmholtz free energy functional, E ( φ ) = (cid:90) Ω (cid:110) Ψ( φ P , φ H , φ N ) + (cid:88) α ∈{ P,H } ε α |∇ φ α | + (cid:88) β ∈RD D β φ β − ( χ c φ σ + χ h φ ECM ) (cid:88) α ∈{ P,H } φ α (cid:111) d x , where CH = { P, H, N } and RD = { σ, MDE, TAF , ECM } . Here, Ψ is a double-well potential, such asΨ( φ P , φ H , φ N ) = (cid:80) α ∈{ T,H,P } C Ψ α φ α (1 − φ α ) with appropriate prefactors C Ψ α , χ c and χ h are the chemo-taxis and haptotaxis parameters, respectively, see [26, 58]. The quantity, ε α is a parameter associated withthe interfacial thickness separating the different cell species. We assume a volume-averaged velocity v for theproliferative cells, hypoxic cells, and the nutrients concentration. This assumption is regarded as reasonablewhenever cells are tightly packed. Following these assumptions and conventions, we arrive at the followingsystem of equations governing the model:(2.1) ∂ t φ P + div( φ P v ) = div( m P ( φ ) ∇ µ P ) + S P ( φ ) + G P ( φ P , φ H , φ N ) ˙ W P ,µ P = ∂ φ P Ψ( φ P , φ H , φ N ) − ε P ∆ φ P − χ c φ σ − χ h φ ECM ,∂ t φ H + div( φ H v ) = div( m H ( φ ) ∇ µ H ) + S H ( φ ) + G H ( φ P , φ H , φ N ) ˙ W H ,µ H = ∂ φ H Ψ( φ P , φ H , φ N ) − ε H ∆ φ H − χ c φ σ − χ h φ ECM ,∂ t φ N = S N ( φ ) ,∂ t φ σ + div( φ σ v ) = div( m σ ( φ )( D σ ∇ φ σ − χ c ∇ ( φ P + φ H )) + S σ ( φ ) + J σv ( φ σ , p, Π Γ φ v , Π Γ p v ) δ Γ ,∂ t φ MDE + div( φ MDE v ) = div( m MDE ( φ ) D MDE ∇ φ MDE ) + S MDE ( φ ) ,∂ t φ TAF + div( φ TAF v ) = div( m TAF ( φ ) D TAF ∇ φ TAF ) + S TAF ( φ ) ,∂ t φ ECM = S ECM ( φ ) , − div( K ∇ p ) = J pv ( p, Π Γ p v ) δ Γ − div( KS p ( φ , µ P , µ H )) ,v = − K ( ∇ p − S p ( φ , µ P , µ H )) , in the space-time domain Ω × (0 , T ) and we supplement the system with homogeneous Neumann boundaryconditions. Here J pv and J σv are the fluid flux and nutrient flux as described in Subsection 2.3. We considerthe following choices of the coupling source functions; see [19],(2.2) S P ( φ ) = λ pro P φ σ φ P (1 − φ T ) − λ deg P φ P − λ PH H ( σ PH − φ σ ) φ P + λ HP H ( φ σ − σ HP ) φ H ,S H ( φ ) = λ pro H φ σ φ H (1 − φ T ) − λ deg H φ H + λ PH H ( σ PH − φ σ ) φ P − λ HP H ( φ σ − σ HP ) φ H − λ HN H ( σ HN − φ σ ) φ H ,S N ( φ ) = λ HN H ( σ HN − φ σ ) φ H ,S ECM ( φ ) = − λ deg ECM φ ECM φ MDE + λ pro ECM φ σ (1 − φ ECM ) H ( φ ECM − φ pro ECM ) ,S σ ( φ ) = − λ pro P φ σ φ P − λ pro H φ σ φ H + λ deg P φ P + λ deg H φ H − λ pro ECM φ σ (1 − φ ECM ) H ( φ ECM − φ pro ECM )+ λ deg ECM φ ECM φ MDE ,S MDE ( φ ) = − λ deg MDE φ MDE + λ pro MDE ( φ P + φ H ) φ ECM σ HP σ HP + φ σ (1 − φ MDE ) − λ deg ECM φ ECM φ MDE ,S TAF ( φ ) = λ pro TAF (1 − φ TAF ) φ H H ( φ H − φ H P ) − λ deg TAF φ TAF ,S p ( φ , µ P , µ H ) = ( µ P + χ c φ σ + χ h φ ECM ) ∇ φ P + ( µ H + χ c φ σ + χ h φ ECM ) ∇ φ H . Here, λ pro α and λ deg α denote the proliferation and degradation rate of the α -th species, respectively, λ αβ thetransition rate from the α -th to the β -th volume fraction, σ αβ the corresponding nutrient threshold for the M. FRITZ, P. K. JHA, T. K ¨OPPL, J. T. ODEN, A. WAGNER, AND B. WOHLMUTH transition, and H is the Heaviside step function. Further, φ pro ECM denotes the threshold level for the ECMdensity in order to activate the production of ECM fibers. The symbol Π Γ denotes the projection of the 1Dvascular functions onto Γ.Additionally, we add randomness in the equations governing φ P and φ H in the form of G P ˙ W P and G H ˙ W H , respectively. Here, W P and W H are independent cylindrical Wiener processes on L (Ω). Moreover,the functions G α , α ∈ { P, H } , are given by(2.3) G α ( φ P , φ H , φ N ) = ω α H (( φ α − φ ωα )(1 − φ α − φ ωα )) H (( φ T − φ ωT )(1 − φ T − φ ωT )) . The parameters ω P and ω H denote the intensity of the noise. Moreover, the Heaviside functions with theparameters φ ωP , φ ωH and φ ωT control the regions of randomness, e.g., here restricted to subdomains at theinterfaces of the phase-fields. We refer to reference [41] for an application of the stochastic Cahn–Hilliardequation to tumor modeling where the noise can represent environmental disruptions due to therapeutictreatment and to [4] for discussion of the numerical approximation of stochastic Cahn–Hilliard equation.2.3. Interaction between the 3D and 1D model.
We apply the Kedem–Katchalsky law [23] to quantifythe flux of nutrients across the vessel surface; i.e., J σv in (2.1) is given by(2.4) J σv ( φ σ , p, φ v , p v ) = (1 − r σ ) J pv ( p, p v ) φ vσ + L σ ( φ v − φ σ ) , where J pv denotes the nutrient flux, which is caused by the flux of blood plasma from the vessels into thetissue or vice versa. Further, J pv is governed by the Starling law [53], i.e., J pv ( p, p v ) = L p ( p v − p ) where p denotes an averaged pressure over the circumference of cylinder cross-sections and is computed in thefollowing way: For each point s i on the curve Λ i , we consider the circle ∂B R i ( s i ) of radius R i , which isperpendicular to Λ i ; see Figure 2. Thus, the tissue pressure p is averaged with respect to ∂B R i ( s i ), p ( s i ) = 12 πR i (cid:90) ∂B Ri ( s i ) p | Γ d S. The part J pv φ vσ in the Kedem–Katchalsky law (2.4) is weighted by a factor 1 − r σ , r σ being a reflectionparameter, introduced to account for the permeability of the vessel wall with respect to the nutrients. Thevalue of φ vσ is set to φ v for p v ≥ p and to φ σ otherwise. The second term on the right of the law (2.4) is aFickian type law, accounting for the tendency of the nutrients to balance their concentration levels, wherethe permeability of the vessel wall is represented by the parameter L σ .The interaction between the vascular network and the tissue occur at the vessel surface Γ, and thus, weconcentrate the flux J σv by means of the Dirac measure δ Γ ; i.e., we define δ Γ ( ϕ ) = (cid:90) Γ ϕ | Γ d S, for a sufficiently smooth test function ϕ with compact support. Moreover, we introduce the projection Π Γ ofthe 1D quantities onto the cylinder Γ via extending its function values Π Γ φ v ( s ) = φ v ( s i ) for all s ∈ ∂B R i ( s i ).2.4. One-dimensional model for transport in the vascular network.
The one-dimensional vesselvariables φ v and p v represent averages across cross-section of the blood vessels. Thus, the one-dimensionalvariables φ v and p v on a 1D vessel Λ i , i ∈ { , . . . , N } , depend only on s i . See also [30] for more detailsrelated to the derivation of the 1D pipe flow and transport models. With these conventions, the 1D modelequations for flow and transport on Λ i are given by(2.5) ∂ t φ v + ∂ s i ( v v φ v ) = ∂ s i ( m v ( φ v ) D v ∂ s i φ v ) − πR i J σv ( φ σ , p, φ v , p v ) ,R i π ∂ s i ( K v,i ∂ s i p v ) = 2 πR i J pv ( p, p v ) . Here, we have introduced the permeability K v,i = R i /µ bl of the i -th vessel with µ bl being the viscosityof blood. We assign µ bl a constant value, i.e., non-Newtonian behavior of blood is not considered. Thediffusivity parameter D v is set to the same value as D σ . The blood velocity v v is given by the Darcyequation v v = − K v,i ∂ s i p v . In order to interconnect p v and φ v on Λ i at the inner networks nodes on the intersections x ∈ ∂ Λ i \ ∂ Λ , we require continuity conditions on p v and φ v . Moreover, we enforce conservation of mass to obtain aphysically relevant solution. To formulate these coupling conditions in a mathematical way, we define foreach bifurcation point x an index set N ( x ) = { i ∈ { , . . . , N } | x ∈ ∂ Λ i } . ODELING AND SIMULATION OF VASCULAR TUMORS 7
We state the following continuity and mass conservation conditions at an inner node x ∈ ∂ Λ i : p v | Λ i ( x ) − p v | Λ j ( x ) = 0 , for all j ∈ N ( x ) \{ i } ,φ v | Λ i ( x ) − φ v | Λ j ( x ) = 0 , for all j ∈ N ( x ) \{ i } , (cid:88) j ∈ N ( x ) − R j π µ bl ∂p v ∂s j (cid:12)(cid:12)(cid:12) Λ j ( x ) = 0 , (cid:88) j ∈ N ( x ) (cid:16) v v φ v − m v ( φ v ) D v ∂φ v ∂s j (cid:17)(cid:12)(cid:12)(cid:12) Λ j ( x ) = 0 . Angiogenesis: Network Growth Algorithm
As noted earlier, angiogenesis is triggered by an increased TAF concentration φ TAF around the pre-existingblood vessels. After the TAF molecules are emitted by the dying hypoxic tumor cells, they move throughthe tissue matrix and may encounter sensor ligands on the vessel surfaces. If the TAF concentration islarge enough at the vessel surfaces, an increased number of sensors in the vessel wall are activated anda reproduction of endothelial cells forming the vessel walls is initiated. As a result, the affected vesselscan elongate, resulting in two different kinds of vessel elongations or growth. In medical literature, see,e.g., [15, 51], this process is referred to as apical growth and sprouting of vessels. The term apical growthis derived from the term apex denoting the tip of a blood vessel, i.e., apical growth is the type of growthoccurring at the tip of a vessel. On the other hand, the sprouting of new vessels results in the formationof new vessels at other places on the vessel surface. In order to increase or decrease the flow of blood andnutrients through the vessels, it is observed that the newly formed blood vessels can adapt their vessel radiiwhich is caused, e.g., by an increased wall shear stress at the inner side of the vessel walls. Combining thesemechanisms, an increased supply of nutrients for both the healthy and cancerous tissue can be achieved suchthat the tumor can continue to grow.In the following, we describe how an angiogenesis step can be realized within our mathematical model.It is assumed that in such a step the apical growth is considered first and then the sprouting of new vesselsis simulated. At the end of an angiogenesis step, the radii of the vessels are adapted to regulate the bloodflow. The 1D network that is updated during an angiogenesis step is denoted by Λ old .3.1.
Apical growth.
Since the apical growth occurs only at the tips of the blood vessels, we consider allthe boundary nodes x of the network Λ old contained in the inner part of Ω, i.e., x ∈ ∂ Λ old and x / ∈ ∂ Ω.Moreover, we assume that x is contained in the segment Λ i ⊂ Λ old . At x , the value of the TAF concentrationis denoted by φ TAF ( x ). If this value exceeds a certain threshold T h
TAF : φ TAF ( x ) ≥ T h
TAF , the tip of thecorresponding vessel is considered as a candidate for growth.There are two types of growth that are allowed to occur at the apex of a vessel: either the vessel canfurther elongate or it can bifurcate. In order to decide which event occurs, a probabilistic method is used.According to [54] and the references therein, the ratio r i = l i /R i of the vessel Λ i follows a log-normaldistribution:(3.1) p b ( r ) ∼ LN ( r, µ r , σ r ) = 1 r (cid:112) πσ r exp (cid:18) − (log r − µ r ) σ r (cid:19) . The parameters µ r and σ r represent the mean value and standard deviation of the probability distribution p b ,respectively. Using the cumulative distribution function of p b , we decide whether a bifurcation is consideredor not. This means that a bifurcation at x ∈ ∂ Λ i ∪ ∂ Λ is formed with a probability of:(3.2) P b ( r ) = Φ (cid:18) ln r − µ r σ r (cid:19) = 12 + 12 erf (cid:18) ln r − µ r (cid:112) σ r (cid:19) , where Φ denotes the standard normal cumulative distribution function and x (cid:55)→ erf( x ) the Gaussian errorfunction. We refer to Figure 3 for the illustration of an exemplary vessel, which bifurcates. Moreover, wedepict the distribution of the ratio l i /R i according to (3.1), the radii of its bifurcations, see (3.4) below, andthe probability of the occurrence of a bifurcation, see (3.2). M. FRITZ, P. K. JHA, T. K ¨OPPL, J. T. ODEN, A. WAGNER, AND B. WOHLMUTH . . . . l/R p b ( l / R ) distribution of l i /R i .
51 ratio l/R P b ( l / R ) probability of bifucation R min R c R max ,
000 radius R N t ( R , R c , R c / ) distribution of R i,b j l i R i R i,b R i,b (a) (b) (c) Figure 3.
Given a vessel with length l i and radius R i , we plot the probability of theoccurrence of a bifurcation (red curve in figure (b)), the ratio of its length over the radius(blue curve in figure (a)), and the distribution of the radii of the sproutings (green curvein figure (c)); we choose R i = 1 . · − , R c = 2 − R i according to (3.4), µ r = 1, σ r = 0 . R min = 9 · − , R max = 3 . · − according to Table 2, R max = max { R max , R i } = R i .If a single vessel is formed at x , the direction of growth d g is based on the TAF concentration: d g ( x ) = ∇ φ TAF ( x ) (cid:107)∇ φ TAF ( x ) (cid:107) + λ g d i (cid:107) d i (cid:107) , where (cid:107) · (cid:107) denotes the Euclidean norm. The vector d i = x i, − x i, is the orientation of the vessel Λ i , andthe value λ g ∈ (0 ,
1] represents a regularization parameter that can be used to circumvent the formation ofsharp bendings and corners. This is necessary if the TAF gradient at x encloses an acute angle with d i . Theradius R i (cid:48) of the new vessel Λ i (cid:48) is taken over from Λ i i.e. R i (cid:48) = R i . Having the radius R i (cid:48) at hand, we use(3.2) to determine the length l i (cid:48) of Λ i (cid:48) . Before Λ i (cid:48) is incorporated into the network Λ old , we check, whetherit intersects another vessel in the network. If this is the case, Λ i (cid:48) is not added to Λ old .In the case of bifurcations, we have to choose the radii, orientations and lengths of the new branches b and b . The radii of the new branches are computed based on a Murray-type law. It relates the radius R i of the father vessel to the radius R i,b of branch b and the radius R i,b of branch b as follows [37]:(3.3) R γi = R γi,b + R γi,b , where γ denotes the bifurcation exponent. According to [54], γ can vary between 2 . .
5. In additionto (3.3), we require an additional equation to determine the radii of the branches. Towards this end, it isassumed that R b follows a truncated Gaussian normal distribution:(3.4) R c = 2 − γ R i , R b k ∼ N t ( R, µ = R c , σ = R c / , k ∈ { , } , ODELING AND SIMULATION OF VASCULAR TUMORS 9 which is set to zero outside of the interval [ R min , R max ] with R max = max { R max , R i } ; we refer to Table 2for a choice of parameters for R min and R max . Additionally, the radius of the parent vessel acts as a naturalbound for the radius of its bifurcations.The selection of R b k is motivated as follows: Using the radius R i of Λ i , we compute the expected radius R c resulting from Murray’s law for a symmetric bifurcation ( R b = R b ). Here, R c is used as a mean value fora Gaussian normal distribution, with a small standard derivation. This yields bifurcations that are slightlydeviating from a symmetric bifurcation which is in accordance with Murray’s law. Having R b and R b athand, we compute the corresponding lengths l b and l b as in the case of a single vessel.We refer to Figure 4 for the distribution of the radii of the bifurcating vessels. We note that the ideal caseis a symmetric bifurcation, that means both radii which correspond to the mean. Further, we also depicttwo asymmetric cases where the radii deviate from the mean. R min R c R max , ,
200 radius R N t ( R , R c , R c / ) distribution of radii(a) (b) (c) Figure 4.
Distribution of the radii of the bifurcating vessels, choosing R i = 0 . R c =2 − R i . Examples of bifurcations with different radii are given, R = 1 . · − (case (a)), R = R c (case (b)), R = 1 . · − (case (c)).The creation of a bifurcation is accomplished by specifying the orientations of the two branches. At first,we define the plane in which the bifurcation is contained. The normal vector n p of this plane is given by thecross product of the vessel orientation d i and the growth direction d g from the non-bifurcating case: n p ( x ) = d i × d g (cid:107) d i × d g (cid:107) . The exact location of the plane is determined such that the vessel Λ i is contained in this plane. Furtherconstraints for the bifurcation configuration are related to the bifurcation angles. In [36,37], it is shown howoptimality principles like minimum work and minimum energy dissipation can be utilized to derive formulasrelating the radii of the branches to the branching angles α (1) i and α (2) i :(3.5) cos (cid:0) α (1) i (cid:1) = R i + R b − R b · R i R b and cos (cid:0) α (2) i (cid:1) = R i + R b − R b · R i R b . The value α ( k ) i denotes the bifurcation angle between branch k ∈ { , } and the father vessel. Rotatingthe vector d g at x around the axis defined by n p ( x ) counterclockwise by α (1) i + α (2) i , we obtain two newgrowth directions d b = d g and d b . These vectors are used to define the main axes of the two cylindersrepresenting the two branches. This choice of the growth directions can be considered as a compromisebetween the optimality principles provided by [36, 37] and the tendency of the network to adapt its growthdirection to the nutrient demand of the surrounding tissue. At the end of the apical growth phase, we obtaina 1D network denoted by Λ ap . Sprouting of new vessels.
In the second phase of the angiogenesis process, we examine each vesselor segment Λ i ⊂ Λ ap . As ligands has been already mentioned, the sprouting of inner vessels is triggeredby TAF molecules touching some sensor ligands in the vessel walls. Therefore, we determine for the middleregion of each segment, i.e., Λ i ( s i ) ⊂ Λ i , s i ∈ (0 . , .
75) at which place an averaged TAF concentration φ TAF attains its maximum φ (max) TAF . As in the previous section φ TAF is determined by means of an integralexpression: φ TAF ( s i ) = 12 πR i (cid:90) ∂B Ri ( s i ) φ TAF ( x ) d S, s i ∈ (0 . , . . We consider only the parameters s i ∈ (0 . , . i . Furthermore, boundary edges are not considered, and we demand that the edges shouldhave a minimal length l min to avoid the formation of tiny vessels. If φ (max) TAF is larger than
T h
TAF , we attacha new vessel Λ i (cid:48) at x . As in the case of apical growth, the local TAF gradient is considered as the preferredgrowth direction of the new vessel: d g ( x ) = ∇ φ TAF ( x ) (cid:107)∇ φ TAF ( x ) (cid:107) . In order to prevent Λ i (cid:48) from growing in the direction of Λ i , we demand that d g draws an angle of at least π . The new radius R i (cid:48) is computed as follows:˜ R i = Ψ R i , ˜ R i (cid:48) = ( ˜ R i − R i ) γ = (Ψ − γ R i , R i (cid:48) = (cid:40) R i (cid:48) ∼ U (1 . · R min , ˜ R i (cid:48) ) if 1 . · R min < ˜ R i (cid:48) R min otherwise.Ψ > R min denotes the minimal radius of a blood vessel, and U stands for theuniform distribution, i.e., new segment radius R i (cid:48) is chosen from the interval [ R min , ˜ R i (cid:48) ] based on a uniformdistribution. For a given radius ˜ R i (cid:48) , the new length l i (cid:48) of Λ i (cid:48) is determined by means of (3.1).Finally, three new vessels Λ i , Λ i and Λ i (cid:48) are added to the network Λ ap . As in the case of apical growth,we test whether a new vessel intersects an existing vessel, before we incorporate Λ i (cid:48) into Λ ap . In addition,we check whether a terminal vessel, i.e., a vessel that is part of ∂ Λ ap can be linked to another vessel. For thispurpose, the distance of the point x b ∈ ∂ Λ ap ∪ ∂ Λ i to its neighboring network nodes that are not directlylinked to x b is computed. If the distance is below a certain threshold dist link , the corresponding networknode is considered as a candidate to be linked with x b . If x b is part of an artery or the high pressure regionof Λ ap , we link it preferably with a candidate at minimal distance and whose pressure is in the low pressureregion (venous part). If x b is part of a vein, the roles are switched.3.3. Adaption of the vessel radii.
In the final phase of the angiogenesis step, we iterate over all thevessels Λ i ⊂ Λ sp and compute for each vessel the wall shear stress τ w by: τ w,i = 4 . µ bl πR i | Q i | , Q i = − K v,i R i π ∆ p v,i l i , where ∆ p v,i is the pressure drop along Λ i . By means of the 1D wall shear stress, the wall shear stressstimulus for the vessel adaption is given by [57]: S WSS ,i = log( τ w,i + τ ref ) .τ ref is a constant that is included to avoid a singular behavior at lower wall shear stresses [47]. Followingthe theoretical model presented in [55], the change with respect to the radius ∆ R i over a time step size ∆ t ,is proportional to this stimulus and R i :∆ R i = k WSS · ∆ t · R i · S WSS ,i and the new vessel R new ,i is given by: R new ,i = R i + ∆ R i . k W SS is a proportionality constant. If R new ,i ∈ [ R min , . · R i ] , we update the vessel radius of Λ i . Otherwise, if R i < R min and ∂ Λ i ∪ ∂ Λ sp (cid:54) = ∅ , the vessel is removed fromthe network. Finally, we obtain a new network Λ new that can be updated within the following angiogenesisstep. ODELING AND SIMULATION OF VASCULAR TUMORS 11 Numerical Discretization
With our mathematical models for tumor growth, blood flow and nutrient transport as well as angiogenesisprocesses presented in previous sections, we now turn our attention to numerical solution strategies. Towardthis end, let us consider a time step n given by the interval [ t n , t n +1 ], with ∆ t = t n +1 − t n . At the beginningof a time step n , we decide, whether an angiogenesis process has to be simulated or not. As examples ofrelevant simulations, we consider an angiogenesis process after each third time step. If angiogenesis has tobe taken into account, we follow the steps described in Section 3. Given the 1D network Λ at the timepoint t n , we first apply the algorithm for the apical growth. Afterwards, the sprouting of new vessels andthe adaption of the vessel radii is simulated. Finally, we obtain a new network Λ new for the new time point t n +1 . If the simulation of angiogenesis is omitted in the respective time step n , Λ is directly used for thesimulation of the tumor growth as well as blood flow and nutrient transport, see Figure 5. A n g i o g e n e s i s Simulate angiogenesis?Input: 1D network Apical growthSprouting growthAdapt radiiNew networkAdapt boundary conditions T u m o r g r o w t h a n d b l oo d fl o w Y e s N o Solve the 3D-1D blood fl ow model Solve the 3D-1D nutrient transport equationsSolve the tumor growth equations Figure 5.
Simulation steps within a single time step.For the time discretization of the 3D model equations in Section 2, the semi-implicit Euler method is usedi.e. we keep the linear terms implicit and the nonlinear terms explicit with respect to time. Discretizingthe model equations in space, standard conforming trilinear Q p ) and nutrienttransport ( φ σ ) are solved by means of cell centered finite volume methods. The computational mesh is givenby a union of cubes having an edge length of h D .We use finite elements to approximate the higher-order Cahn–Hilliard type equations as well advection-reaction-diffusion equations corresponding to the species φ T AF , φ
ECM , φ
MDE in (2.1). In order to ensuremass conservation for both flow and nutrient transport in the interstitial domain, finite volume schemes aretaken into account, since they are locally mass conservative.In order to solve the 3D-1D coupled system, such as pressure ( p v , p ), the iterative Block-Gauss-Seidelmethod is used, i.e., in each iteration, we first solve the equation system for the 1D system. Then theupdated 1D solution is used to solve the equation system derived from the 3D problem. We stop theiteration when the change in the current and previous iteration solution is within small tolerance. At eachtime step, we first solve the ( p v , p ) coupled system. Afterwards the ( φ v , φ σ ) coupled system is solved. Next,we solve the remaining equations in the 3D system. This is summarized in Algorithm 1. In the remainderof this section, the discretizations of the 1D and 3D systems are outlined. VGM discretization of the 1D PDEs.
It remains to specify the numerical solution techniques forthe 1D network equations. The time integration is based on the implicit Euler method. For the spatial dis-cretization of the 1D equations, the Vascular Graph Method (VGM) is considered. This method correspondsto a node centered finite volume method [50, 60]. We then briefly describe this numerical method as wellas the discretization of the terms arising in the context of the 3D-1D coupling. We restrict ourselves to thepressure equations.As mentioned in Section 2, the 1D network is given by a graph-like structure, consisting of edges Λ i ⊂ Λand network nodes x i ∈ Λ. In a first step, we assign to each network node x i an unknown for the pressurethat is denoted by p v,i . Let us assume that the edges containing x i are given by Λ i , . . . , Λ i N and itsmidpoints by m i , . . . , m i N , see Figure 6. Exchange with tissue
Figure 6.
Notation for the Vascular Graph Method.On each edge, Λ k ∈ { Λ i , . . . , Λ i N } , we consider the following PDE for the pressure; see also (2.5). Forconvenience, the curve parameter is simply denoted by s . − R k π ∂ s ( K v,k ∂ s p v ) = − πR k L p ( p v − p ) . Next, we establish for the node x i a mass balance equation taking the fluxes across the cylinders Z i l intoaccount. Z i l is a cylinder having the edge Λ i l as a rotation axis and the radius R i l . Furthermore its top andbottom facets are located at m i l and x i , respectively (see Figure 6). The corresponding curve parametersare denoted by s ( x i ) and s ( x i l ) , l ∈ { , . . . , N } . Accordingly, the mass balance equation reads as follows: − N (cid:88) l =1 (cid:90) s ( m il ) s ( x i ) R i l π ∂ s ( K v,i l ∂ s p v ) d s = − πL p N (cid:88) l =1 (cid:90) s ( m il ) s ( x i ) R i l ( p v − p ) d s. Integration by parts yields: − N (cid:88) l =1 R i l πK v,i l ∂ s p v | s ( m il ) + N (cid:88) l =1 R i l πK v,i l ∂ s p v | s ( x i ) = − πL p n (cid:88) l =1 R i l (cid:90) s ( m il ) s ( x i ) p v − p d s. Approximating the derivatives by central finite differences and using the mass conservation equation (seeSubsection 2.4): N (cid:88) l =1 R i l πK v,i l ∂ s p v | s ( x i ) = 0 , at an inner node x i , it follows that N (cid:88) l =1 R i l πK v,i l p v,i − p v,i l l i l = − πL p N (cid:88) l =1 R i l (cid:90) s ( m il ) s ( x i ) p v − p d s, ODELING AND SIMULATION OF VASCULAR TUMORS 13 where l i l denotes the length of the edge Λ i l . Denoting the mantle surface of Z i l by S i l , we have: N (cid:88) l =1 R i l πK v,i l l i l ( p v,i − p v,i l ) = − L p N (cid:88) l =1 | S i l | p v,i + L p N (cid:88) l =1 (cid:90) S il p d S. Computing the integrals (cid:82) S il p d S , we introduce the decomposition of Ω into M finite volume cells CV k :Ω = (cid:83) Mk =1 CV k . The pressure unknown assigned to CV k is given by p k . Using this notation, one obtains: (cid:90) S il p d S = (cid:88) CV k ∩ S il (cid:54) = ∅ (cid:90) CV k ∩ S il p d S ≈ | S i l | (cid:88) CV k ∩ S il (cid:54) = ∅ | CV k ∩ S i l || S i l | (cid:124) (cid:123)(cid:122) (cid:125) =: w kil p k = | S i l | (cid:88) CV k ∩ S il (cid:54) = ∅ w ki l p k . We note that one has to guarantee that the relation (cid:80) CV k ∩ S il (cid:54) = ∅ w ki l = 1 holds. Otherwise, a consistentmass exchange between the vascular system and the tissue can not be enforced. All in all, we obtain a linearsystem of equations for computing the pressure values. On closer examination, it can be noted that thismatrix is composed of four blocks, i.e., two coupling blocks as wells as a block for the 1D diffusion term andthe 3D diffusion term. As said earlier, at each time step, we decouple the 1D and 3D pressure equationsand use a Block-Gauss-Seidel iteration to solve the two systems until the 3D pressure is converged. Thediscretization of the nutrient equation is exerted in a similar manner, where the main difference consists inadding an upwinding procedure for the convective term. At each time step, the nutrient equations are alsosolved using a Block-Gauss-Seidel iteration.4.1.1. Initial and boundary conditions for the 1D PDEs.
Since we use a transient PDE to simulate thetransport of nutrients, we require an initial condition for the variable φ v . In doing so, a threshold R T for theradii is introduced in order to distinguish between arteries and veins. If the radius of a certain vessel is below R T , the vessel is considered as an artery and otherwise as a vein. In case of an artery, we set φ v ( t = 0) = 1and in case of a vein φ v ( t = 0) = 0 is used. When the network starts to grow, initial values for the newlycreated vessels have to be provided. If a new vessel is created due to sprouting growth, we consider thevessel or edge to which the new vessel is attached. At the point of the given vessel, where the new vessel oredge is added, φ v is interpolated linearly. For this purpose the two values of φ v located at the nodes of theexisting vessel are used. The interpolated value is assigned to both nodes of the newly created vessel.When apical growth takes place and a new vessel is added to x ∈ ∂ Λ. In this case, we consider φ v ( x , t )for a time point t and assign it to the newly created node, since we assume that no flow boundary conditionsare enforced. With respect to the boundary conditions for the 1D pressure PDE, the following distinctionof cases is made: • Dirichlet boundary for p v if x ∈ ∂ Λ ∩ ∂ Ω . In this case, we set: p v ( x ) = p v,D ( x ), where p v,D is a givenDirichlet value at x . Numerically, we enforce this boundary condition by setting in the correspondingline of the matrix all entries to zero except for the entry on the diagonal which is fixed to the value 1.Additionally, the corresponding component of the right hand side vector contains the Dirichlet value p D .Let us assume that y is the neighbor of x on the edge Λ . The other edges adjacent to y are denoted byΛ , . . . , Λ N . Then the balance equation for y has to be adapted as follows to account for the Dirichletboundary condition p v,D : R πK v, l ( p v ( y ) − p v,D ) + N (cid:88) j =2 R j πK v,j l j ( p v ( y ) − p v,j )= − L p (cid:12)(cid:12)(cid:12) ˜ S (cid:12)(cid:12)(cid:12) p v ( y ) + L p (cid:88) CV k ∩ ˜ S (cid:54) = ∅ w k p k − L p N (cid:88) j =2 | S j | p v ( y ) + L p | S j | (cid:88) CV k ∩ S j (cid:54) = ∅ ,j> w kj p k , where ˜ S is the mantle surface of the cylinder covering the whole edge Λ . • Homogeneous Neumann boundary for p v if x ∈ ∂ Λ ∩ Ω . Let x ∈ ∂ Λ i ∩ ∂ Λ ∩ Ω, then we set − R i π K v,i ∂ s p v (cid:12)(cid:12) x = 0, resulting in the following discretization: R i π K v,i p v ( x ) − p v ( y ) l i = − L p | S i | p v,i + L p | S i | (cid:88) CV k ∩ S i (cid:54) = ∅ w ki l p k , where y ∈ ∂ Λ i ∩ Λ and l i is the length of the edge Λ i .Summarizing, we consider for the pressure in the network Dirichlet boundaries at the boundary of the 3Ddomain Ω and homogeneous Neumann boundary conditions in the inner part of Ω. For the nutrients, theimplementation of boundary conditions is more challenging, since an upwinding procedure has to be takeninto account. • Dirichlet boundary for φ v if x ∈ ∂ Λ i ∩ ∂ Λ ∩ ∂ Ω and v v ( x ) ≈ − K v,i p v ( y ) − p v,D ( x ) l i > . In this case, we set: φ v ( x , t ) = φ v,D ( x ), where φ v,D is a given Dirichlet value at x . The numericalimplementation can be exerted analogously to the case of the pressure p v . • Homogeneous Neumann boundary for φ v if x ∈ ∂ Λ ∩ Ω . Let x ∈ ∂ Λ i ∩ ∂ Λ ∩ Ω and we set( v v φ v − D v ∂ s i φ v ) | x = 0 , resulting in the following discretization: l i φ v ( x , t + ∆ t ) − φ v ( x , t )∆ t + v v φ v | m i − D v φ v ( y , t + ∆ t ) − φ v ( x , t + ∆ t ) l i = − πR i (cid:34) (1 − r σ ) (cid:90) s ( m i ) s ( x ) J pv ( p, p v ) · φ vσ ( s, t + ∆ t ) d s + L σ (cid:90) s ( m i ) s ( x ) φ v ( s, t + ∆ t ) − φ σ ( s, t + ∆ t ) d s (cid:35) , where y ∈ ∂ Λ i ∩ Λ. The integrals modeling the exchange terms are discretized as in the case of thepressure equations. • Upwinding boundary for φ v if x ∈ ∂ Λ ∩ ∂ Λ i ∩ ∂ Ω and v v ( m i ) ≈ − K v,i p v ( y ) − p v ( x ) l i ≤ . Here, we obtain the following semi-discrete equation: l i φ v ( x , t + ∆ t ) − φ v ( x , t )∆ t + v v | m i φ v ( y , t + ∆ t ) − v v | m i φ v ( x , t + ∆ t )= − πR i (cid:34) (1 − r σ ) (cid:90) s ( m i ) s ( x ) J pv ( p, p v ) · φ vσ ( s, t + ∆ t ) d s + L σ (cid:90) s ( m i ) s ( x ) φ v ( s, t + ∆ t ) − φ σ ( s, t + ∆ t ) d s (cid:35) , where y ∈ ∂ Λ i ∩ Λ.4.2.
Discretization of the 3D PDEs.
Suppose φ α n , µ α n , p n , p v n , v n denote the various fields at time t n .Let V h be a subspace of H (Ω , R ) consisting of continuous piecewise trilinear functions on a uniform meshΩ h . We consider φ α n ∈ V h for α ∈ { P, H, N, T AF, ECM, M DE } , µ P n , µ H n ∈ V h , and v h ∈ [ V h ] . The testfunctions are denoted by ˜ φ ∈ V h for species in V h , ˜ µ ∈ V h for chemical potentials µ P , µ H , and ˜ v ∈ [ V h ] forvelocity.Given a time step n and solutions φ α n , µ α n , p n , p v n , v n , we are interested in the solution at the next timestep. For the 3D-1D coupled pressure ( p v n +1 , p n +1 ), as mentioned earlier, we utilize a block Gauss-Seideliteration, where the discretization of the 1D equation is discussed in Subsection 4.1 and discretization of the3D equation using finite-volume scheme is provided in (4.3). Similarly, ( φ v n +1 , φ σ n +1 ) is solved using a blockGauss-Siedel iteration with the discretization of the 1D equation along the lines of the discretization of the1D pressure equation and discretization of the 3D equation provided in (4.4). We then solve the proliferative, ODELING AND SIMULATION OF VASCULAR TUMORS 15 hypoxic, necrotic, MDE, ECM, and TAF systems sequentially. Once we have pressure p n +1 , we compute thevelocity v n +1 ∈ [ V h ] using the weak form:(4.1) ( v n +1 , ˜ v ) = ( − K ( ∇ p n +1 − S p n ) , ˜ v ) , ∀ ˜ v ∈ [ V h ] , where S p n = S p ( φ n , µ P n , µ H n ), see (2.2). For the advection terms, using the fact that ∇ p · n = 0 on ∂ Ω andso v n +1 · n = 0 on ∂ Ω, we can write (cid:16) ∇ · ( φ α v n +1 ) , ˜ φ (cid:17) = − (cid:16) φ α v n +1 , ∇ ˜ φ (cid:17) , ∀ ˜ φ ∈ V h , ∀ α ∈ { P, H, T AF, M DE } . (4.2)In what follows, we consider the expression on the right hand side in above equation for the advection terms.For fields φ a , a ∈ { P, H, N, σ, T AF, ECM, M DE } , and chemical potentials µ P , µ H , we assume homoge-neous Neumann boundary condition on ∂ Ω. Next, we provide the discretization of the 3D fields. • Pressure.
Let CV ∈ Ω h denotes the typical finite volume cell and σ ∈ ∂CV face of a cell CV . Let( p kv n +1 , p kn +1 ) denote the pressures at k th iteration and time t n +1 . Suppose we have solved for p k +1 v n +1 following Subsection 4.1. To solve p k +1 n +1 , we consider, for all CV ∈ Ω h , − (cid:88) σ ∈ ∂CV (cid:90) σ K ∇ p k +1 n +1 · n d S = − (cid:88) σ ∈ ∂CV (cid:90) σ KS p n · n d S + (cid:90) Γ ∩ CV J pv ( p k +1 n +1 , Π Γ p k +1 v n +1 ) d S = − (cid:88) σ ∈ ∂CV (cid:90) σ KS p n · n d S + N (cid:88) i =1 (cid:90) Γ i ∩ CV J pv ( p k +1 n +1 , Π Γ i p k +1 v n +1 ) d S. (4.3)Above follows by the integration of the pressure equation in (2.1) over CV and using the divergence theo-rem. Here J pv ( p, p v ) = L p ( p v − p ), Γ = ∪ Ni =1 Γ i is the total vascular surface, and Π Γ i ( p v ) is the projectionof the 1D pressure defined on the centerline Λ i onto the surface of the cylinder, Γ i . • Nutrients.
Suppose we have solved for φ k +1 v n +1 . To solve φ k +1 σ n +1 , we consider, for all CV , (cid:90) CV φ k +1 σ n +1 − φ σ n ∆ t d V + (cid:88) σ ∈ ∂CV (cid:90) σ φ k +1 σ n +1 ˆ v n +1 · n d S − (cid:88) σ ∈ ∂CV (cid:90) σ m σ ( φ n ) (cid:16) D σ ∇ φ k +1 σ n +1 (cid:17) · n d S + (cid:90) CV λ pro P φ P n φ k +1 σ n +1 d V + (cid:90) CV λ pro H φ H n φ k +1 σ n +1 d V + (cid:90) CV λ pro ECM (1 − φ ECM n ) H ( φ ECM n − φ pro ECM ) φ k +1 σ n +1 d V = (cid:90) CV λ deg P φ P n + λ deg H φ H n d V + (cid:90) CV λ deg ECM φ ECM n φ MDEn d V − (cid:88) σ ∈ ∂CV (cid:90) σ m σ ( φ n ) χ c ∇ ( φ P n + φ H n ) · n d S + N (cid:88) i =1 (cid:90) Γ i ∩ CV J σv ( φ k +1 σ n +1 , p k +1 n +1 , Π Γ i φ k +1 v n +1 , Π Γ i p k +1 v n +1 ) d S, (4.4)where J σv is given by (2.4). Noting that the velocity isˆ v n +1 = − K ∇ p n +1 + KS p n , we divide the advection term, for σ ∈ ∂CV , into two parts: (cid:90) σ φ k +1 σ n +1 ˆ v n +1 · n d S = − K (cid:90) σ φ k +1 σ n +1 ∇ p n +1 · n d S + K (cid:90) σ φ k +1 σ n +1 S p n · n d S. (4.5)For the first term we apply the upwinding scheme. For the second term, we perform quadrature approx-imation to compute the integral over the face σ . • Proliferative.
For a general double-well potential Ψ( φ P , φ H , φ N ) = (cid:80) a ∈{ T,P,H } C Ψ a φ a (1 − φ a ) , weconsider the convex-concave splitting, see [17], as followsΨ( φ P , φ H , φ N ) = (cid:88) a ∈{ T,P,H } C Ψ a φ a + (cid:88) a ∈{ T,P,H } C Ψ a ( φ a − φ a − φ a ) . (4.6) This results in ∂ φ P Ψ( φ P , φ H , φ N ) = (cid:88) a ∈{ T,P } C Ψ a φ a + (cid:88) a ∈{ T,P } C Ψ a φ a (4 φ a − φ a − . (4.7)The expression for ∂ φ H Ψ can be derived analogously. In our implementation, φ P , φ H , φ N are the mainstate variables and φ T is computed using φ T = φ P + φ H + φ N . Let the mobility ¯ m P n at the current stepbe given by ¯ m P n = M P (cid:2) ( φ P n ) + (1 − φ T n ) + (cid:3) , (4.8)where for a field f , ( f ) + is the projection onto [0 ,
1] given by( f ) + = f if f ∈ [0 , , f ≤ , f ≥ . (4.9) We solve for φ P n +1 , µ P n +1 using the weak forms below (cid:18) φ P n +1 − φ P n ∆ t , ˜ φ (cid:19) − (cid:16) φ P n +1 v n +1 , ∇ ˜ φ (cid:17) + (cid:16) ¯ m P n ∇ µ P n +1 , ∇ ˜ φ (cid:17) − (cid:16) λ pro P φ σ n +1 (1 − φ T n ) + φ P n +1 , ˜ φ (cid:17) + (cid:16) λ deg P φ P n +1 , ˜ φ (cid:17) = (cid:16) λ HP H ( φ σ n +1 − σ HP )( φ H n ) + , ˜ φ (cid:17) − (cid:16) λ PH H ( σ PH − φ σ n +1 )( φ P n ) + , ˜ φ (cid:17) + 1∆ t (cid:18) G P n (cid:90) t n +1 t n d W P , ˜ φ (cid:19) (4.10)and (cid:0) µ P n +1 , ˜ µ (cid:1) − (cid:0) C Ψ T + C Ψ P ) φ P n +1 , ˜ µ (cid:1) − (cid:0) (cid:15) P ∇ φ P n +1 , ∇ ˜ µ (cid:1) = (cid:0) C Ψ T φ T n (4 φ T n − φ T n − , ˜ µ (cid:1) + (cid:0) C Ψ P φ P n (4 φ P n − φ P n − , ˜ µ (cid:1) + (3 C Ψ T ( φ H n + φ N n ) , ˜ µ ) − (cid:0) χ c φ σ n +1 + χ h φ ECM n , ˜ µ (cid:1) , (4.11)where ( · ) + is the projection to [0 ,
1] defined in (4.9), G P n = G P ( φ P n , φ H n , φ N n ) is given by (2.3), and W P is the cylindrical Wiener process. We discuss the computation of stochastic term in more detail inSubsection 4.2.1. • Hypoxic.
Let the mobility ¯ m H n be given by¯ m H n = M H (cid:2) ( φ H n ) + (1 − φ T n ) + (cid:3) . (4.12)To solve for φ H n +1 , µ H n +1 , we consider (cid:18) φ H n +1 − φ H n ∆ t , ˜ φ (cid:19) − (cid:16) φ H n +1 v n +1 , ∇ ˜ φ (cid:17) + (cid:16) ¯ m H n ∇ µ H n +1 , ∇ ˜ φ (cid:17) − (cid:16) λ pro H φ σ n +1 (1 − φ T n ) + φ H n +1 , ˜ φ (cid:17) + (cid:16) λ deg H φ H n +1 , ˜ φ (cid:17) = (cid:16) λ P H H ( σ P H − φ σ n +1 )( φ P n ) + , ˜ φ (cid:17) − (cid:16) λ HP H ( φ σ n +1 − σ HP )( φ H n ) + , ˜ φ (cid:17) − (cid:16) λ HN H ( σ HN − φ σ n +1 )( φ H n ) + , ˜ φ (cid:17) + 1∆ t (cid:18) G H n (cid:90) t n +1 t n d W H , ˜ φ (cid:19) (4.13)and (cid:0) µ H n +1 , ˜ µ (cid:1) − (cid:0) C Ψ T + C Ψ H ) φ H n +1 , ˜ µ (cid:1) − (cid:0) (cid:15) H ∇ φ H n +1 , ∇ ˜ µ (cid:1) = (cid:0) C Ψ T φ T n (4 φ T n − φ T n − , ˜ µ (cid:1) + (cid:0) C Ψ H φ H n (4 φ H n − φ H n − , ˜ µ (cid:1) + (3 C Ψ T ( φ P n + φ N n ) , ˜ µ ) − (cid:0) χ c φ σ n +1 + χ h φ ECM n , ˜ µ (cid:1) , (4.14)where G H n = G H ( φ P n , φ H n , φ N n ) is given by (2.3) and W H is the cylindrical Wiener process. ODELING AND SIMULATION OF VASCULAR TUMORS 17 • Necrotic. (cid:18) φ N n +1 − φ N n ∆ t , ˜ φ (cid:19) = (cid:16) λ HN H ( σ HN − φ σ n +1 )( φ H n ) + , ˜ φ (cid:17) . (4.15) • MDE. (cid:18) φ MDEn +1 − φ MDEn ∆ t , ˜ φ (cid:19) − (cid:16) φ MDEn +1 v n +1 , ∇ ˜ φ (cid:17) + (cid:16) m MDE ( φ n ) D MDE ∇ φ MDEn +1 , ∇ ˜ φ (cid:17) + (cid:16) λ deg MDE φ MDEn +1 , ˜ φ (cid:17) + (cid:18) λ pro MDE ( φ P n + φ H n ) φ ECM n σ HP σ HP + φ σ n +1 φ MDEn +1 , ˜ φ (cid:19) = (cid:18) λ pro MDE ( φ P n + φ H n ) φ ECM n σ HP σ HP + φ σ n +1 , ˜ φ (cid:19) − (cid:16) λ deg ECM φ ECM n φ MDEn , ˜ φ (cid:17) . (4.16) • ECM. (cid:18) φ ECM n +1 − φ ECM n ∆ t , ˜ φ (cid:19) = (cid:16) λ pro ECM φ σ n +1 H ( φ ECM n − φ pro ECM )(1 − φ ECM n ) , ˜ φ (cid:17) − (cid:16) λ deg ECM φ MDEn .φ ECM n , ˜ φ (cid:17) (4.17) • TAF. (cid:18) φ TAF n +1 − φ TAF n ∆ t , ˜ φ (cid:19) − (cid:16) φ TAF n +1 v n +1 , ∇ ˜ φ (cid:17) + (cid:16) D T AF ∇ φ TAF n +1 , ∇ ˜ φ (cid:17) + (cid:16) λ pro TAF φ H n +1 φ TAF n +1 H ( φ H n +1 − φ H P ) , ˜ φ (cid:17) = (cid:16) λ pro TAF φ H n +1 H ( φ H n +1 − φ H P ) , ˜ φ (cid:17) − (cid:16) λ deg TAF φ TAF n , ˜ φ (cid:17) . (4.18)The steps followed in solving the coupled system of equations including the angiogenesis step are sum-marized in Algorithm 1. Algorithm 1:
The 3D-1D tumor growth model solver with angiogenesis step Input :Model parameters, φ α , v , ∆ t, T, TOL for α ∈ A := { P, H, N, σ, TAF , MDE, ECM } Output : φ α n , µ P n , µ H n , p n , v n , p v n , φ σv b for all n n = 0, t = 0 while t ≤ T do φ α n = φ α n +1 ∀ α ∈ A , µ P n = µ P n +1 , µ H n = µ H n +1 , p v n = p v n +1 , φ v n = φ v n +1 if apply angiogenesis(n) == True then apply angiogensis model described in Section 3 update
1D systems if the network is changed end solve coupled linear system ( p v n +1 , p n +1 ) using block Gauss-Seidel iteration and Subsection 4.1 solve coupled linear system ( φ v n +1 , φ σ n +1 ) using block Gauss-Seidel iteration and Subsection 4.1 solve velocity v n +1 using (4.1) solve (cid:0) φ P n +1 , µ P n +1 (cid:1) using the semi-implicit scheme in (4.10) and (4.11) solve (cid:0) φ H n +1 , µ H n +1 (cid:1) using the semi-implicit scheme in (4.13) and (4.14) solve φ N n +1 using the semi-implicit scheme in (4.15) solve φ MDEn +1 using the semi-implicit scheme in (4.16) solve φ ECM n +1 using the semi-implicit scheme in (4.17) solve φ TAF n +1 using the semi-implicit scheme in (4.18) n (cid:55)→ n + 1, t (cid:55)→ t + ∆ t end Stochastic component of the system.
Generally, the cylindrical Wiener processes W α , α ∈ { P, H } , on L (Ω) with Ω = (0 , can be written as W α ( t )( x ) = ∞ (cid:88) i,j,k =1 η αijk ( t ) cos( iπx /L ) cos( jπx /L ) cos( kπx /L ) (cid:124) (cid:123)(cid:122) (cid:125) =: e ijk , where x = ( x , x , x ), L = diam(Ω), { e ijk } is an orthonormal basis in L (Ω), and { η αijk } i,j,k ∈ N is a familyof real-valued, independent, and identically distributed Brownian motions. Following [4, 9], we approximatethe term involving the Wiener process in the fully discretized system as follows(4.19) 1∆ t (cid:18)(cid:90) t n +1 t n d W α ( t ) , ξ (cid:19) L ≈ t (cid:88) i,j,k,i + j + k In this section, we apply the models described in Sections 2 and 3 and use the numerical discretizationsteps discussed in Section 4. We consider examples that showcase the effects of angiogenesis on the tumorgrowth. For this purpose, the model parameters and the basic setting for our simulations are introduced inSubsection 5.1. In the base setting, we consider two vessels, one representing an artery and the other a vein,and introduce an initially spherical tumor core. Based on this setting, tumor growth is simulated first withoutconsidering angiogenesis, i.e., the growth algorithm from Section 3 is not applied. Afterwards, we repeatthe same simulation including the angiogenesis effects and study the differences between the correspondingsimulations results. Initial tumor ArteryVein Figure 7. Initial setting with boundary conditions for a first numerical experiment. Pres-sure at top plane ( z = 0) and bottom ( z = 2) ends of the artery are 3000 and 2000respectively. Similarly, the pressure at top and bottom ends of the vein are fixed at 1100and 1600, respectively. The high pressure end of the artery is the inlet end, and there weassign the nutrients value to 1. The high pressure end of the vein is also the inlet, and herewe assign the nutrients value to 0. At the remaining ends, we apply the upwinding schemefor solving the nutrient equation.We then consider a scenario consisting of a tumor core surrounded by a small capillary network. Weobtain the network from source . First, we rescale the network so that it fits into the domain Ω = (0 , . https://physiology.arizona.edu/sites/default/files/brain99.txt ODELING AND SIMULATION OF VASCULAR TUMORS 19 The vessel radii are rescaled such that the maximal and minimal vessel radius is given non-dimensionally by0 . . C Ψ T φ T (1 − φ T ) , where C Ψ T is a constant. Since the model involves stochastic PDEs as well as stochastic angiogenesis,we employ Monte-Carlo approximation based on samples of the probability distributions characterizing thewhite noise terms, using 10 samples for the case without angiogenesis and 50 samples for the case withangiogenesis. We use the samples to report statistics of quantity of interests such as total tumor volume,vessel volume density, etc. Table 1. List of parameters and their values for the numerical simulations. Unlisted pa-rameters are set to zero. φ ωα , ω α , and J α are parameters related to Wiener process, see (2.3)and (4.19). Parameter Value Parameter Value Parameter Value λ pro P λ pro H λ deg P , λ deg H λ pro ECM λ pro MDE , λ deg MDE λ deg ECM λ PH λ HP λ HN σ PH σ HP σ HN M P M H C Ψ T ε P ε H λ pro TAF D TAF M TAF L p − D σ M σ K − D MDE M MDE L σ . D v . µ bl r σ . J α , α ∈ { P, H } φ ωα , α ∈ { P, H, T } . ω α , α ∈ { P, H } . Table 2. List of parameters and their values for the growth algorithm and numericaldiscretization Parameter Value Function T h TAF . · − Threshold for the TAF concentration (sprouting) µ r . σ r . λ g . γ . R min . · − Minimal vessel radius l min . 13 Minimal vessel length for which sprouting is activated R max . 035 Maximal vessel radius R T . 05 Threshold for the radius to distinguish between arteries and veins for t = 0 Ψ . 05 Sprouting parameterdist link . 08 Maximal distance at which a terminal vessel is linked to the network τ ref . k WSS . t . h D . h D . 25 Mesh size of the initial 1D grid∆ t net t Angiogenesis (network update) time interval5.1. Setup and model parameters for the two vessel problem. As a computational domain Ω, wechoose a cube, Ω = (0 , . Within Ω two different vessels are introduced: an artery and a vein; see Figure 7. The radius of the vein R v is given by R v = 0 . R a is set to R a = 0 . . , . , 0) and ends at (0 . , . , . , . , 0) and ends at(1 . , . , φ v , we choose φ v = 1 in the artery and φ v = 0 in the vein. The initial value forthe nutrient variable φ σ in the tissue matrix is given by φ σ = 0 . 5. In order to define the initial conditionsfor the tumor, we consider a ball B T of radius r c = 0 . x c = (1 . , . , . B T , thetotal tumor volume fraction φ T smoothly decays from 1 in the center to 0 on the boundary of the ball: φ T ( x , t = 0) = exp (cid:18) − − ( | x − x c | /r c ) (cid:19) , if | x − x c | < r c , , otherwise . (5.1)Thereby, the necrotic and hypoxic volume fractions, φ N and φ H , are set initially to zero. In the other partsof the domain, all the volume fractions for the tumor species are set to 0 at t = 0. In Table 1, the parametersfor the model equations in Section 2 are listed and Table 2 contains the parameters for the growth algorithmdescribed in Section 3 as well as the discretization parameters. In particular, the parameters for the stochasticdistributions are listed, which determine the radii and vessel lengths, the probability of bifurcations, and thesprouting probability of new vessels.5.2. Tumor growth without angiogenesis. The simulation results for tumor growth without angiogenesisare summarized in Figure 8. For t = 8, the tumor cell distribution within the plane perpendicular to the z -axis at z = 1 . φ T = φ P + φ H + φ N , φ N , aswell as the nutrients φ σ are presented. It can be observed that the primary tumor is enlarged and smallsatellites are formed in the vicinity of the main tumor. The distribution of the necrotic cells indicates thatthe main tumor consists mostly of necrotic cells, while the hypoxic and proliferative cells gather around thenutrient-rich blood vessels. This means that the tumor cells can migrate against the flow from the arterytowards the vein. Apparently, the chemical potential caused by the nutrient source dominates the interstitialflow.These observations are consistent with simulation results and measurements discussed, e.g., in [18,32,52].In [32, 52] a tumor is introduced into a mouse. At the same time, anti-VEGF agents were injected into themouse, such that the sprouting of new vessels growing towards the tumor is prevented. This process leadsto the formation of satellites located in the vicinity of the primary tumor as well as the accumulation oftumor cells at nutrient-rich vessels and cells. Furthermore, the primary tumor stops growing and forms alarge necrotic core.In Figure 9, the L -norms of the tumor species over time are presented for the case when angiogenesiswas inactive and when it was active. While the profiles for different species in two cases look similar, thetotal tumor is about 70% times higher when angiogenesis is active. Diffusivity D σ = 3 is large enough,and therefore, the nutrients originating from the nutrient-rich vessels diffuse quickly throughout the domain.The result is that in the case of angiogenesis, the overall nutrient concentration is higher compared to thecase without angiogenesis, while the spatial variation of the nutrient is the same in the two; and hence thegrowth of the tumor in the two cases are similar except that the tumor grows more rapidly in the caseof angiogenesis. We will see in our second example, where D σ = 0 . 05 is much smaller, that the nutrientconcentration is higher near the nutrient rich vessels and tumor growth is more concentrated near theseregions, see Figure 13.In summary, one can conclude that without angiogenesis a tumor can grow to a certain extent, before theprimary tumor starts to collapse, i.e., a large necrotic core is formed. However, this does not mean that thetumor cells are entirely removed from the healthy tissue. If there is a source of nutrients such as an arteryclose by, transporting nutrient rich blood, a portion of tumor cells can survive by migrating towards theneighboring nutrient source. ODELING AND SIMULATION OF VASCULAR TUMORS 21 Main tumor Satellites Main tumorSatellitesNecrotic kernel v Figure 8. Top left: Tumor growing in a mouse that is treated with anti-VEGF agents.As a consequence tumor satellites in the vicinity of the main tumor can be detected. Imagetaken from [52], with permission from Elsevier. Top right: φ T presented in a plane at z = 1 perpendicular to the z -axis. As seen in the medical experiments the formation ofsatellites and the accumulation of tumor cells at the nutrient-rich artery are reproduced inthe simulations. Bottom left: Distribution of necrotic cells ( φ N ). It can be seen that themain tumor consists of a large necrotic kernel. Bottom right: Distribution of nutrients( φ σ ).5.3. Tumor growth with angiogenesis. As in the previous subsection, we compute the L -norms of thetumor species φ T at time t = 8. However, since several stochastic processes are involved in the networkgrowth and also Wiener processes appear in the proliferative and hypoxic cell mass balances, several datasets have to be created in order to rule out statistical variations. In this context, the issue arises as to howmany data sets are needed to obtain a representative value. In order to investigate this, we compute forevery sample i , the L -norm of the tumor species, denoted by φ α L ,i , α ∈ { T, P, H, N } . Additionally, thevolume of the blood vessel network V i is computed. For each data set with i samples, we compute the meanvalues: mean i ( V ) = 1 i i (cid:88) j =1 V j , mean (cid:0) φ α L ,i (cid:1) = 1 i i (cid:88) j =1 φ α L ,j , α ∈ { T, P, H, N } . In Figure 10, the mean values mean i ( V ) and mean (cid:0) φ α L ,i (cid:1) , α ∈ { T, P, H, N } , are shown. From the plots wesee that the mean of || φ T || L stabilizes after about 25 samples. For the vessel volume, fluctuations in thesample means reduce with increasing sample and get small after 30 samples. While the results in Figure 10show that mean of the QoIs stabilizes with increasing sample and converge to some number, the trajectory in Time L - no r m Without angiogenesis THNP Time L - no r m With angiogenesis THNP Figure 9. Quantities of interests (QoIs) related to tumor species over a time interval [0 , µ α ( t ) − σ α ( t ) , µ α ( t ) + σ α ( t )) ∩ [0 , ∞ ), for t ∈ [0 , T ], where µ α ( t ) , σ α ( t ) are the mean and standard deviations of QoI α ∈ {(cid:107) φ T (cid:107) L , (cid:107) φ P (cid:107) L , (cid:107) φ H (cid:107) L , (cid:107) φ N (cid:107) L } at time t . The variations in the QoIs for the non-angiogenesis case are very small. The mean of the total tumor volume fraction (cid:107) φ T (cid:107) L atthe final time for the angiogenesis case is about 1.7 times that of the non-angiogenesis case.the figure could change with change in sample values. For example, if we shuffle the samples and recomputethe quantities in Figure 10, various curves may look different. Number of datasets M ean L - no r m THNP Number of datasets M ean v e ss e l v o l u m e Vessel volume Figure 10. The mean values for the L -norm of the tumor cell volume fractions φ T , φ P , φ H , φ N and the volume of the blood vessel network at time t = 8 from increas-ing number of samples. Results correspond to the two-vessel setting. The mean of the totaltumor volume fraction appears to be stable after about 28 samples. The mean of the vesselvolume shows smaller fluctuations as the number of samples in the data set grows. ODELING AND SIMULATION OF VASCULAR TUMORS 23 As mentioned earlier, Figure 9 presents the L -norms of tumor species. For the angiogenesis simulations,we compute the mean and standard deviation using 50 samples. We see that the total tumor volume fractionvaries from sample to sample, as expected. Both the hypoxic and the tumor cells show an exponential growthafter t ≈ 4; see Figure 9. After decreasing until t ≈ 3, the proliferative cell mass grows from t ≈ t = 0.24 t = 0.64t = 3.20 t = 5.60 Figure 11. Growth of the tumor and the network at four different time points t ∈ { . , . , . , . } and one specific sample. We show contour φ T = 0 . φ T = 0 . t = 0 . 24 (top-left figure), we observe that new vessels originate from the artery and movetowards the hypoxic core; the directions of these vessels being based on the gradient of TAF with randomperturbations. At t = 0 . 64 (top-right), we see a large number of new vessels formed as predicted by theangiogenesis algorithm. However, at time t = 3 . t = 3 . t = 5 . z = 1 plane along with the nutrient distribution in the vessels inFigure 12. The plot corresponds to time t = 8. The plots corresponding to the necrotic species show thatthe necrotic concentration is typically higher away from the nutrient-rich vessels. From the hypoxic plot,we see that it is higher near the vessels, and this is explained by the fact that as soon as the proliferationof new tumor cells takes place, due to nutrient concentration below the proliferative-to-hypoxic transitionthreshold, these newly added proliferative tumor cells convert to hypoxic cells. Further transition to necrotic Tumor NecroticHypoxic Proliferative Figure 12. Tumor cell distributions in the plane z = 1 with angiogenesis for one sample. Top left: φ T . Top right: φ N . Bottom left: φ H . Bottom right: φ P .cells would take place if the nutrients are even below the hypoxic-to-necrotic transition threshold. This isalso consistent with the increase concentration of the proliferative cells near the outer tumor-healthy cellphase interface.5.4. Angiogenesis and tumor growth for the capillary network. Returning to (5.1), let us considera smooth spherical tumor core with center at x c = (1 . , . , . 7) and radius r c = 0 . , . The initial blood vessel network and boundary conditions for pressure and nutrient on the networkis described in Figure 15.In the simulation, the vessels are uniformly refined two times. We fix φ a = 0 for a ∈ { H, N, TAF } and φ σ = 0 . t = 0. The domain is discretized uniformly with a mesh size h D = 0 . t = 0 . φ v = φ v in = 1 andpressure p v = p in = 8000 . p v = p out = 1000 . t = 0, p v and φ v at internal network nodes are set to zero. Furthermore, we set L σ = 0 . D σ = 0 . D T AF = 0 . T h TAF = 0 . µ r = 1 . k WSS = 0 . 45, and ∆ t net = 10∆ t . All the other parameter valuesremain unchanged, see Table 1 and Table 2.We first compare the tumor volume with and without angiogenesis; see Figure 13. The results are similarto the two-vessel setting. They show that the overall tumor growth is higher with angiogenesis as expected.We also observe that proliferative cells start to grow rapidly at t ≈ . t ≈ . 75 without angiogenesis. Production of necrotic cells is higher in the non-angiogenesis case until time t ≈ 5. Compared to Figure 9 for the two-vessel setting, the variations in the tumor species related QoIs are ODELING AND SIMULATION OF VASCULAR TUMORS 25 Time L - no r m Without angiogenesis THNP Time L - no r m With angiogenesis THNP Figure 13. Quantities of interests (QoIs) related to tumor species over a time interval[0 , 8] for the capillary network setting. Similar to the two-vessels simulation, we computethe mean and standard deviation using 50 and 10 samples for the case with and withoutangiogenesis, respectively. We refer to Figure 9 for more details on the plots. As in the caseof two-vessels setting, the variation in the QoIs are much smaller for the non-angiogenesiscase. The mean of the total tumor volume fraction (cid:107) φ T (cid:107) L for the angiogenesis case is about1.62 times that of the non-angiogenesis case. Number of datasets M ean L - no r m THNP Number of datasets M ean v e ss e l v o l u m e Vessel volume Figure 14. The mean values for the L -norm of the tumor cell volume fractions φ T , φ P , φ H , φ N and the volume of the blood vessel network at time t = 8. Results cor-respond to the capillary network setting. The mean of the total tumor volume fractionstabilizes with small samples. This agrees with Figure 13 that shows that the variations in L -norm QoIs are overall smaller. The mean of the vessel volume shows some change whensamples are small and stabilizes as the size of data set grows.much smaller in Figure 13. This may be due to the fact that diffusivity of the nutrients in the latter case ismuch smaller, and that L σ is also smaller in the later case resulting in a smaller exchange of nutrients. Next,we plot the mean of QoIs as we increase the size of data in data sets in Figure 14; results show that mean of tumor species related QoIs is stable and can be computed accurately using fewer samples. This relatesto the fact that we see smaller variation in the L -norm of the tumor species φ T in Figure 13. The mean ofvessel volume shows some variations for smaller data sets and the variations get smaller later on; still thevariations are very small and contained in range [0 . , . t = 8 for three samples. In Figure 17, we compare the the tumor species at time t = 5 . 12 for the angiogenesis and non-angiogenesis case. As in the two-vessel case, the tumor starts to growfaster after it is vascularized. Apparently, the tumor cells tend to migrate towards the nutrient rich partof the computational domain despite the fact that they have to move against the direction of flow which isinduced by the pressure gradient. Not surprisingly, the volume fraction of the necrotic cells is larger in thepart that is facing away from the nutrient rich part and related to the whole tumor it remains relativelysmall.It is interesting to observe that, as in the two-vessel case, the contour plot of φ T for φ T = 0 . Initial Setting t = 1.20t = 5.60t = 3.04 vv vv Figure 15. Top left: Spherical tumor core (Contour plot for φ T = 0 . 9) at x c =(1 . , . , . 7) with radius r c = 0 . φ v = φ v in = 1 . p v = p in = 8000 is prescribed as a Dirichlet boundary condition. Top right: Formation offirst sprouts at t = 1 . 20 growing towards the tumor core. Bottom left: Around t = 3 . 04a complex vascular network is formed and the tumor starts to grow towards the nutrients. Bottom right: At t = 5 . 60 the tumor is significantly enlarged and creates satellites nearthe nutrient vessels. ODELING AND SIMULATION OF VASCULAR TUMORS 27 Figure 16. Plot of vessel network at t = 8 from three samples for the capillary networksetting. Tumor (with angiogenesis) Tumor (without angiogenesis)Necrotic (with angiogenesis) Necrotic (without angiogenesis)Hypoxic (with angiogenesis) Hypoxic (without angiogenesis) Figure 17. Distribution of tumor cells for t = 5 . 12. The simulation results for tumorgrowth supported by angiogenesis are shown on the left hand side, while the results fortumor growth without angiogenesis are presented on the right. The necrotic ( φ N ) andhypoxic ( φ H ) volume fractions are visualized in the z -plane at z = 0 . 7. For φ T , in bothcases a contour plot for φ T = 0 . Summary and outlook In this work, we presented a stochastic model for tumor growth characterized by a coupled system ofnonlinear PDEs of Cahn-Hilliard types coupled with a model of an evolving vascular network. To accountfor the effects of angiogenesis, a 3D-1D coupled PDE system is developed to simulate flow and nutrienttransport within both the network and porous tissue. In this model, the blood vessel network is given by a1D graph-like structure coupling the flow and transport processes in the network and tissue. Furthermore,the model facilitates the handling of a growing network structure which is crucial for the simulation ofangiogenesis.The angiogenesis process in our model is simulated by an iterative algorithm starting from a single arteryand a vein or a given network. The blood vessel network supplying the tumor employs Murray’s law todetermine the radii at a bifurcation of network capillaries. The choice of the other radii, lengths as well asthe probability for bifurcations is governed by stochastic algorithms. The direction of growth of the vesselsis given by the gradient of the local TAF concentration.The model can reproduce qualitative behaviors of tumor growth. As in experiments with mice, oursimulations indicate that is suffering a lack of nutrient supply the tumor forms a large necrotic kernel andsatellites located preferably at nutrient rich vessels. If tumor growth is combined with angiogenesis, thenan excessive exponential growth of the tumor mass can be observed in agreement with existing literature onangiogenesis.We believe that our model can serve as a starting point for important predictive simulations of cancertherapy. Clearly, much experimental data is needed to calibrate model parameters of both the tumor modeland the vascular model. In our work, the adaption of the vessel radii is related to the wall shear stress.However, other stimuli can influence the vessel radii that could be included, such as metabolic hematocrit-related stimulus, which may also lead to a significant pruning and restructuring of the network.Another important effect that occurs for certain types of tumor is that when the tumor starts to growis the deformation of healthy tissue as well as the blood vessel networks. To capture such phenomena,mechanical deformations may have to be added and a coupling concept accounting for the deformation ofthe networks has to be designed. Having such a model in hand, the effect of therapies by drugs preventingthe sprouting of new vessels and the uptake of metastases could be studied. There are, of course, aspectsof the complex numerical algorithms employed in the 1D-3D models that could be modified and expanded.We leave these issues and other extensions to be explored in future work. Acknowledgements The authors gratefully acknowledge the support of the Deutsche Forschungsgemeinschaft (DFG) throughTUM International Graduate School of Science and Engineering (IGSSE), GSC 81. The work of PKJ andJTO was supported by the U.S. Department of Energy, Office of Science, Office of Advanced ScientificComputing Research, Mathematical Multifaceted Integrated Capability Centers (MMICCS), under AwardNumber DE-SC0019303. References [1] D. Ambrosi, F. Bussolino, and L. Preziosi , A review of vasculogenesis models , Journal of Theoretical Medicine, 6(2005), pp. 1–19.[2] D. Ambrosi and L. Preziosi , Cell adhesion mechanisms and stress relaxation in the mechanics of tumours , Biomechanicsand Modeling in Mechanobiology, 8 (2009), pp. 397–413.[3] A. R. Anderson and M. A. J. Chaplain , Continuous and discrete mathematical models of tumor-induced angiogenesis ,Bulletin of Mathematical Biology, 60 (1998), pp. 857–899.[4] D. Antonopoulou, L. Banas, R. N¨urnberg, and A. Prohl , Numerical approximation of the stochastic cahn-hilliardequation near the sharp interface limit , arXiv preprint arXiv:1905.11050, (2019).[5] N. Bellomo, N. Li, and P. K. Maini , On the foundations of cancer modelling: Selected topics, speculations, and per-spectives , Mathematical Models and Methods in Applied Sciences, 18 (2008), pp. 593–646.[6] N. Bellomo and L. Preziosi , Modelling and mathematical problems related to tumor evolution and its interaction withthe immune system , Mathematical and Computer Modelling, 32 (2000), pp. 413–452.[7] H. Byrne and L. Preziosi , Modelling solid tumour growth using the theory of mixtures , Mathematical Medicine andBiology: A Journal of the IMA, 20 (2003), pp. 341–366.[8] P. Carmeliet and R. K. Jain , Molecular mechanisms and clinical applications of angiogenesis , Nature, 473 (2011),pp. 298–307. ODELING AND SIMULATION OF VASCULAR TUMORS 29 [9] S. Chai, Y. Cao, Y. Zou, and W. Zhao , Conforming finite element methods for the stochastic cahn–hilliard–cook equation ,Applied Numerical Mathematics, 124 (2018), pp. 44–56.[10] M. A. Chaplain, M. Lachowicz, Z. Szyma´nska, and D. Wrzosek , Mathematical modelling of cancer invasion: Theimportance of cell-cell adhesion and cell-matrix adhesion , Mathematical Models and Methods in Applied Sciences, 21(2011), pp. 719–743.[11] M. A. Chaplain and G. Lolas , Mathematical modelling of cancer cell invasion of tissue: The role of the urokinaseplasminogen activation system , Mathematical Models and Methods in Applied Sciences, 15 (2005), pp. 1685–1734.[12] V. Cristini, X. Li, J. S. Lowengrub, and S. M. Wise , Nonlinear simulations of solid tumor growth using a mixturemodel: invasion and branching , Journal of Mathematical Biology, 58:723 (2009).[13] V. Cristini and J. Lowengrub , Multiscale Modeling of Cancer: An Integrated Experimental and Mathematical ModelingApproach , Cambridge University Press, 2010.[14] M. Dorraki, A. Fouladzadeh, A. Allison, C. S. Bonder, and D. Abbott , Angiogenic networks in tumors—insightsvia mathematical modeling , IEEE Access, 8 (2020), pp. 43215–43228.[15] H. M. Eilken and R. H. Adams , Dynamics of endothelial cell behavior in sprouting angiogenesis , Current Opinion in CellBiology, 22 (2010), pp. 617–625.[16] C. Engwer, C. Stinner, and C. Surulescu , On a structured multiscale model for acid-mediated tumor invasion: Theeffects of adhesion and proliferation , Mathematical Models and Methods in Applied Sciences, 27 (2017), pp. 1355–1390.[17] D. J. Eyre , Unconditionally gradient stable time marching the cahn-hilliard equation , in Materials Research SocietySymposium Proceedings, vol. 529, Materials Research Society, 1998, pp. 39–46.[18] H. B. Frieboes, F. Jin, Y.-L. Chuang, S. M. Wise, J. S. Lowengrub, and V. Cristini , Three-dimensional multispeciesnonlinear tumor growth – II: Tumor invasion and angiogenesis , Journal of Theoretical Biology, 264 (2010), pp. 1254–1278.[19] M. Fritz, P. K. Jha, T. K¨oppl, J. T. Oden, and B. Wohlmuth , Analysis of a new multispecies tumor growth modelcoupling 3D phase-fields with a 1D vascular network , arXiv preprint arXiv:2006.10477 [math.AP], (2020).[20] M. Fritz, E. Lima, V. Nikolic, J. T. Oden, and B. Wohlmuth , Local and nonlocal phase-field models of tumor growthand invasion due to ECM degradation , Mathematical Models and Methods in Applied Sciences, 29 (2019), pp. 2433–2468.[21] M. Fritz, E. Lima, J. T. Oden, and B. Wohlmuth , On the unsteady Darcy–Forchheimer–Brinkman equation in localand nonlocal tumor growth models , Mathematical Models and Methods in Applied Sciences, 29 (2019), pp. 1691–1731.[22] H. Garcke, K. F. Lam, R. N¨urnberg, and E. Sitka , A multiphase Cahn–Hilliard–Darcy model for tumour growth withnecrosis , Mathematical Models and Methods in Applied Sciences, 28 (2018), pp. 525–577.[23] B. Ginzburg and A. Katchalsky , The frictional coefficients of the flows of non-electrolytes through artificial membranes ,The Journal of General Physiology, 47 (1963), pp. 403–418.[24] D. Hanahan and R. A. Weinberg , Hallmarks of cancer: The next generation , Cell, 144 (2011), pp. 646–674.[25] A. Hawkins-Daarud, K. G. van der Zee, and J. T. Oden , Numerical simulation of a thermodynamically consistentfour-species tumor growth model , International Journal for Numerical Methods in Biomedical Engineering, 28 (2012),pp. 3–24.[26] T. Hillen, K. J. Painter, and M. Winkler , Convergence of a cancer invasion model to a logistic chemotaxis model ,Mathematical Models and Methods in Applied Sciences, 23 (2013), pp. 165–198.[27] E. Hodneland, X. Hu, and J. Nordbotten , Well-posedness, discretization and preconditioners for a class of models formixed-dimensional problems with high dimensional gap , arXiv preprint arXiv:2006.12273, (2020).[28] L. Holmgren, M. S. O’Reilly, and J. Folkman , Dormancy of micrometastases: balanced proliferation and apoptosis inthe presence of angiogenesis suppression , Nature medicine, 1 (1995), pp. 149–153.[29] T. Koch, M. Schneider, R. Helmig, and P. Jenny , Modeling tissue perfusion in terms of 1d-3d embedded mixed-dimension coupled problems with distributed sources , Journal of Computational Physics, 410:100050 (2020).[30] T. K¨oppl, E. Vidotto, and B. Wohlmuth , A 3D-1D coupled blood flow and oxygen transport model to generate mi-crovascular networks , International Journal for Numerical Methods in Biomedical Engineering, e3386 (2020).[31] P. Koumoutsakos, I. Pivkin, and F. Milde , The fluid mechanics of cancer and its therapy , Annual review of fluidmechanics, 45 (2013), pp. 325–355.[32] P. Kunkel, U. Ulbricht, P. Bohlen, M. Brockmann, R. Fillbrandt, D. Stavrou, M. Westphal, and K. Lamszus , Inhibition of glioma angiogenesis and growth in vivo by systemic treatment with a monoclonal antibody against vascularendothelial growth factor receptor-2 , Cancer research, 61 (2001), pp. 6624–6628.[33] E. Lima, J. T. Oden, and R. Almeida , A hybrid ten-species phase-field model of tumor growth , Mathematical Modelsand Methods in Applied Sciences, 24 (2014), pp. 2569–2599.[34] S. R. McDougall, A. Anderson, M. Chaplain, and J. Sherratt , Mathematical modelling of flow through vascularnetworks: implications for tumour-induced angiogenesis and chemotherapy strategies , Bulletin of mathematical biology,64 (2002), pp. 673–702.[35] S. R. McDougall, A. R. Anderson, and M. A. Chaplain , Mathematical modelling of dynamic adaptive tumour-inducedangiogenesis: clinical implications and therapeutic targeting strategies , Journal of theoretical biology, 241 (2006), pp. 564–589.[36] C. Murray , The physiological principle of minimum work applied to the angle of branching of arteries , The Journal ofgeneral physiology, 9 (1926), pp. 835–841.[37] , The physiological principle of minimum work: I. the vascular system and the cost of blood volume , Proceedings ofthe National Academy of Sciences of the United States of America, 12 (1926), pp. 207–214.[38] N. Nargis and R. Aldredge , Effects of matrix metalloproteinase on tumour growth and morphology via haptotaxis , J.Bioengineer. & Biomedical Sci., 6:1000207 (2016). [39] N. Nishida, H. Yano, T. Nishida, T. Kamura, and M. Kojiro , Angiogenesis in cancer , Vascular health and riskmanagement, 2 (2006), pp. 213–219.[40] J. T. Oden, E. Lima, R. C. Almeida, Y. Feng, M. N. Rylander, D. Fuentes, D. Faghihi, M. M. Rahman, M. DeWitt,M. Gadde, et al. , Toward predictive multiscale modeling of vascular tumor growth , Archives of Computational Methodsin Engineering, 23 (2016), pp. 735–779.[41] C. Orrieri, E. Rocca, and L. Scarpa , Optimal control of stochastic phase-field models related to tumor growth , ESAIM:Control, Optimisation and Calculus of Variations, to appear (2020).[42] M. R. Owen, T. Alarc´on, P. K. Maini, and H. M. Byrne , Angiogenesis and vascular remodelling in normal andcancerous tissues , Journal of mathematical biology, 58:689 (2009).[43] S. Parangi, M. O’Reilly, G. Christofori, L. Holmgren, J. Grosfeld, J. Folkman, and D. Hanahan , Antiangiogenictherapy of transgenic mice impairs de novo tumor growth , Proceedings of the National Academy of Sciences, 93 (1996),pp. 2002–2007.[44] C. Patsch, L. Challet-Meylan, E. C. Thoma, E. Urich, T. Heckel, J. F. O’Sullivan, S. J. Grainger, F. G. Kapp,L. Sun, K. Christensen, et al. , Generation of vascular endothelial and smooth muscle cells from human pluripotentstem cells , Nature cell biology, 17 (2015), pp. 994–1003.[45] C. M. Phillips, E. A. Lima, R. T. Woodall, A. Brock, and T. E. Yankeelov , A hybrid model of tumor growth andangiogenesis: In silico experiments , Plos one, 15 (2020), p. e0231137.[46] L. Preziosi , Cancer modelling and simulation , CRC Press, 2003.[47] A. Pries, B. Reglin, and T. Secomb , Structural adaptation of vascular networks: role of the pressure response , Hyper-tension, 38 (2001), pp. 1476–1479.[48] A. Pries, B. Reglin, and T. W. Secomb , Structural adaptation of microvascular networks: functional roles of adaptiveresponses , American Journal of Physiology-Heart and Circulatory Physiology, 281 (2001), pp. H1015–H1025.[49] A. Pries, T. Secomb, and P. Gaehtgens , Structural adaptation and stability of microvascular networks: theory andsimulations , American Journal of Physiology-Heart and Circulatory Physiology, 275 (1998), pp. H349–H360.[50] J. Reichold, M. Stampanoni, A. L. Keller, A. Buck, P. Jenny, and B. Weber , Vascular graph model to simulate thecerebral blood flow in realistic vascular networks , Journal of Cerebral Blood Flow & Metabolism, 29 (2009), pp. 1429–1443.[51] D. Ribatti and E. Crivellato , “sprouting angiogenesis”, a reappraisal , Developmental biology, 372 (2012), pp. 157–165.[52] J. Rubenstein, J. Kim, T. Ozawa, M. Zhang, M. Westphal, D. Deen, and M. Shuman , Anti-VEGF antibody treatmentof glioblastoma prolongs survival but results in increased vascular cooption , Neoplasia, 2 (2000), pp. 306–314.[53] E. P. Salathe and K.-N. An , A mathematical analysis of fluid movement across capillary walls , Microvascular research,11 (1976), pp. 1–23.[54] M. Schneider, J. Reichold, B. Weber, G. Sz´ekely, and S. Hirsch , Tissue metabolism driven arterial tree generation ,Medical image analysis, 16 (2012), pp. 1397–1414.[55] T. Secomb, J. Alberding, R. Hsu, M. Dewhirst, and A. Pries , Angiogenesis: an adaptive dynamic biological patterningproblem , PLoS computational biology, 9:e1002983 (2013).[56] A. Stephanou, S. R. McDougall, A. R. Anderson, and M. A. Chaplain , Mathematical modelling of flow in 2d and3d vascular networks: applications to anti-angiogenic and chemotherapeutic drug strategies , Mathematical and ComputerModelling, 41 (2005), pp. 1137–1156.[57] A. St´ephanou, S. R. McDougall, A. R. Anderson, and M. A. Chaplain , Mathematical modelling of the influenceof blood rheological properties upon adaptative tumour-induced angiogenesis , Mathematical and Computer Modelling, 44(2006), pp. 96–123.[58] Y. Tao and M. Winkler , A chemotaxis-haptotaxis model: the roles of nonlinear diffusion and logistic source , SIAMJournal on Mathematical Analysis, 43 (2011), pp. 685–704.[59] R. D. Travasso, E. C. Poir´e, M. Castro, J. C. Rodrguez-Manzaneque, and A. Hern´andez-Machado , Tumorangiogenesis and vascular patterning: a mathematical model , PloS one, 6:e19989 (2011).[60] E. Vidotto, T. Koch, T. K¨oppl, R. Helmig, and B. Wohlmuth , Hybrid models for simulating blood flow in microvas-cular networks , Multiscale Modeling & Simulation, 17 (2019), pp. 1076–1102.[61] G. Vilanova, M. Bur´es, I. Colominas, and H. Gomez , Computational modelling suggests complex interactions betweeninterstitial flow and tumour angiogenesis , Journal of The Royal Society Interface, 15:20180415 (2018).[62] S. M. Wise, J. S. Lowengrub, H. B. Frieboes, and V. Cristini , Three-dimensional multispecies nonlinear tumor growth– I: Model and numerical method , Journal of Theoretical Biology, 253 (2008), pp. 524–543.[63] C. Wu, D. A. Hormuth, T. A. Oliver, F. Pineda, G. Lorenzo, G. S. Karczmar, R. D. Moser, and T. E. Yankeelov , Patient-specific characterization of breast cancer hemodynamics using image-guided computational fluid dynamics , IEEETransactions on Medical Imaging, (2020).[64] J. Xu, G. Vilanova, and H. Gomez , A mathematical model coupling tumor growth and angiogenesis , PloS one, 11 (2016).[65] , Full-scale, three-dimensional simulation of early-stage tumor growth: The onset of malignancy , Computer Methodsin Applied Mechanics and Engineering, 314 (2017), pp. 126–146.[66]