# Modeling and simulation of vascular tumors
embedded in evolving capillary networks

Marvin Fritz1, Prashant K. Jha2, Tobias Köppl1,∗, J. Tinsley Oden2,
Andreas Wagner1, and Barbara Wohlmuth1,3

###### Abstract.

In this work, we present a coupled 3D-1D model of solid tumor growth within a
dynamically changing vascular network to facilitate realistic simulations of
angiogenesis. Additionally, the model includes erosion of the extracellular
matrix, interstitial flow, and coupled flow in blood vessels and tissue. We
employ continuum mixture theory with stochastic Cahn–Hilliard type phase-field
models of tumor growth. The interstitial flow is governed by a mesoscale
version of Darcy’s law. The flow in the blood vessels is controlled by
Poiseuille flow, and Starling’s law is applied to model the mass transfer in
and out of blood vessels. The evolution of the network of blood vessels is
orchestrated by the concentration of the tumor angiogenesis factors (TAFs);
blood vessels grow towards the increasing TAFs concentrations. This process is
not deterministic, allowing random growth of blood vessels and, therefore, due
to the coupling of nutrients in tissue and vessels, makes the growth of tumors
stochastic. We demonstrate the performance of the model by applying it to a
variety of scenarios. Numerical experiments illustrate the flexibility of the
model and its ability to generate satellite tumors. Simulations of the effects
of angiogenesis on tumor growth are presented as well as sample-independent
features of cancer.

###### Key words and phrases:

tumor growth, 3D-1D coupled blood flow models, angiogenesis, finite elements,
finite volume

###### 2020 Mathematics Subject Classification:

65M08, 65M60, 76S05, 76Z05, 92C17, 92C42.

∗Corresponding author

1Department of Mathematics, Technical University of Munich, Germany
2Oden Institute for Computational Engineering and Sciences, The University of
Texas at Austin, USA
3Department of Mathematics, University of Bergen, Allegaten 41, 5020 Bergen,
Norway

Keywords:  tumor growth, 3D-1D coupled blood flow models, angiogenesis, finite
elements, finite volume

##  1\. Introduction

In this work, we present new computational models and algorithms for
simulating and predicting a broad range of biological and physical phenomena
related to cancer at the tissue scale. We consider the growth of solid
vascular tumors inside living tissue containing a dynamically evolving
vasculature. One of the main goals of this work is to provide realistic
simulations of the vascular growth characterizing angiogenesis, whereby blood
vessels sprout and invade the domain of the solid tumor when prompted by
concentrations of proteins collectively 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 tumor cell concentrations are categorized as proliferative, hypoxic, or
necrotic. To capture the complex interaction of cell species and the evolving
interfaces between species, continuum mixture theory is used as a framework
for constructing mesoscale stochastic phase-field models of the Cahn–Hilliard
type. Other critical phenomena are also addressed by this class of models,
including the erosion of the extracellular matrix (ECM) due to concentrations
of matrix-degrading enzymes (MDEs, such as matrix metalloproteeinase and
urokinase plasminogen activators) that erode the ECM and permit the invasion
of tumor cells as a prelude to metastasis [25, 41].

The volume of an isolated colony of tumor cells will not generally grow beyond
approximately $1.0$ $\mathrm{m}\mathrm{m}$3 [31, 46, 42] unless sufficient
nutrients and oxygen are supplied for proliferation. To acquire such
nutrients, cancerous cells promote angiogenesis [8, 47]. Low levels of oxygen
and nutrient result in tumor cells entering the hypoxia phase during which
they remain dormant and release various proteins such as TAFs that promote the
proliferation of endothelial cells and new vessel formation. Similarly, low
oxygen levels can generate irregular invasive tumors governed by haptotaxis
[25, 41]. Because angiogenesis is one of the major processes through which
tumors grow, anti-angiogenic drugs that inhibit the formation of the new
vascular 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-angiogenic drugs.

Typically, the vasculature near the tumor core in the early stages of tumor
growth may not effectively supply nutrients to the tumor. The vasculature
evolves rapidly and, therefore, the vessel walls are not fully developed and
may be destroyed due to pressure (proliferation of tumor cells result in
higher pressure nearby), the pruning of vessels due to insufficient flow for a
sustained period and, vasculature adaptation and remodeling [45, 52, 50, 51,
59]. Highly interconnected and irregular vasculature with inefficient blood
vessels causes low blood flow rates to the tumor, making it possible that
therapeutic drugs miss the tumor mass altogether [59]. Shear and
circumferential stresses due to blood flow result in vascular adaptation
effects such as vessel radii adaptation, see [60, 38, 51]. All of these
phenomena are represented by the models described herein.

Earlier models taking into account angiogenesis include lattice-probabilistic
network models, see [3, 69, 37, 59, 60, 38, 45]. An overview of this class of
models is given in [15]. Another class of models referred to as agent-based
models has been proposed and extensively studied. There, the idea is to
introduce a phase-field for the tip endothelial cells that takes a value $1$
inside the vessel and $0$ outside and through the agents, which can move
anywhere in the simulation domain following certain rules, which can be
designed to trigger the sprouting of new vessels; see [36, 62, 64, 48]. These
models do not capture blood circulation in the vessel and, therefore, are
unable to be truly coupled to the tumor growth. In [66], a dimensionally
coupled model for drug delivery based on MRI data and a study of dosing
protocols is considered with drug flow in the vessels governed by algebraic
rules instead of PDEs. More recently, vasculature models involving a network
of straight cylindrical vessels supporting the 1D flow of nutrient, oxygen,
and therapeutic drugs and coupled to the 3D tissue domain by versions of the
Starling or Kedem–Katchalsky law have been presented; see [67, 68, 33].

We consider a class of 3D-1D vascular tumor models [21] that approximates the
flow within the blood vessels 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 in
the vessels and tissue, the blood vessels’ three-dimensional nature is
retained by approximating the vessels as 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-dimensional characterizations of vessel segments in the vascular
network embedded in three-dimensional domains of the tissue and the tumor
within the tissue and the assignment of mechanical models to this 3D-1D system
to depict interstitial flow and pressure fluctuations. Mathematical analysis
showing well-posedness and existence of weak solutions for the class of 3D-1D
model considered in this work is performed in a recent paper [21].

In our model, flow in vessels is governed by 1D Poiseuille law, whereas the
flow in tissue is derived by treating the tissue domain as a porous medium and
applying a version of Darcy’s law. The model consists of nutrients in the
tissue and vessels; nutrients in the vessels are governed by the 1D advection-
diffusion equation and advection-diffusion-reaction equation in the tissue.
Flow and nutrients in the tissue and vessel are coupled; we assume that vessel
walls are porous, resulting in the advection and diffusion-driven exchange of
nutrients and coupling of the extravascular and intravascular pressures. Some
aspects of the 1D model architecture and coupling of 3D and 1D models are
based on previous works, see [30, 33, 32]. The 3D tissue domain 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 theory as in [7, 12, 43, 36, 24, 13], and representations of
the principal mechanisms governing the development and evolution of cancer,
see, e.g., [36, 27]. Especially, we note the comprehensive developments of
diffuse-interface multispecies models presented in [65, 19], the ten species
models derived in [36], and the multispecies nonlocal models of adhesion and
promotes a tumor invasion due to ECM degradation described in [22].
Angiogenesis models embedded in models of hypoxic and cell growth are
presented in [36, 67, 68, 21]. Related models of extracellular matrix (ECM)
degradation due to matrix-degrading enzymes (MDEs) and subsequent tumor
invasion and metastasis are discussed in [22, 11, 10, 17]. Several of the
earlier mechanistic models of tumor growth focused on modeling the effects of
mechanical deformation and blood flow, and fluid pressure on tumor growth,
e.g., [1, 2, 49, 5, 6, 34].

A key new feature of the models proposed here is the dynamic growth/deletion
of the vascular network and full coupling between the dynamic network and
tumor system in the tissue microenvironment. In response to TAF generated by
nutrient-starved hypoxic cells, new vessels are formed. Due to the formation
of new vessels, the local conditions such as nutrient concentration changes
near the tumor, affecting TAF production and promoting a higher proliferation
of tumor cells. The rules by which the network grows, or existing vessels are
deleted due to insufficient flow and dormancy, are based on the experimentally
known causes of angiogenesis and are parameterized so that various aspects of
the network growth algorithm can be adjusted based on available experimental
data. By including the time-evolution of the larger vascular tissue domain and
the sprouting, growth, bifurcation, and pruning of the vascular network
orchestrated by a combination of blood supply and tumor-generated growth
factors, a more realistic depiction of tumor growth than 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 as the tissue and 1D network domain and the
equations governing various fields. The details associated with the vessel
network growth are presented in Section 3. Spatial and temporal discretization
and solver schemes for the 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 step are solved using a semi-implicit fixed point iteration scheme.
In Section 5, the model is applied to various situations, and several
simulation experiments are presented. For further details on our
implementation of the solver, we refer to the open-source code at
https://github.com/CancerModeling/Angiogenesis3D1D. Concluding comments are
given in Section 6.

##  2\. Mathematical Modeling

In this work, a colony of tumor cells in an open bounded domain
$\Omega\subset\mathbb{R}^{3}$, e.g., representing an organ, is considered. It
is supported by a system of macromolecules consisting of collagen, enzymes,
and various proteins, that constitute the extracellular matrix. We focus on
phenomenological characterizations to capture mesoscale and macroscale events.
Additionally, we consider a one-dimensional graph-like structure $\Lambda$
inside of $\Omega$ forming a microvascular network, see Figure 1.

Figure 1. Setup of the domain $\Omega$ with the 1D microvascular network
$\Lambda$ and a tumor mass, which is composed of its proliferative
($\phi_{P}$), hypoxic ($\phi_{H}$) and necrotic ($\phi_{N}$) phases.

The single edges of $\Lambda$ of vessel components are denoted by
$\Lambda_{i}$ such that $\Lambda$ is given by
$\Lambda=\bigcup_{i=1}^{N}\Lambda_{i}$ and each edge $\Lambda_{i}$,
$i\in\\{1,\dots,N\\}$, is parameterized by a corresponding curve parameter
$s_{i}$ such that

$\Lambda_{i}=\mathopen{}\mathclose{{}\left\\{\bm{x}\in\Omega\mathopen{}\mathclose{{}\left|\;\bm{x}=\Lambda_{i}(s_{i})=\bm{x}_{i,1}+s_{i}\cdot(\bm{x}_{i,2}-\bm{x}_{i,1}),\;s_{i}\in(0,1)}\right.}\right\\},$

where $\bm{x}_{i,1}\in\Omega$ and $\bm{x}_{i,2}\in\Omega$ mark the boundary
nodes of $\Lambda_{i}$, see Figure 2. For the total 1D network $\Lambda$, we
introduce a global curve parameter $s$, which is interpreted in the following
way: $s=s_{i}$, if $\bm{x}=\Lambda(s)=\Lambda_{i}(s_{i})$. At each value of
the curve parameter $s$, various 1D constituents exist, which interact with
their respective 3D counterpart in $\Omega$.

We introduce the surface $\Gamma$ of the microvascular network $\Lambda$ to
formulate the coupling between the 3D and 1D constituents in Subsection 2.2
and Subsection 2.4. For simplicity, it is assumed that the surface for a
single vessel is approximated by a cylinder with a constant radius, see Figure
2. The radius of a vessel that is associated with $\Lambda_{i}$ is given by
$R_{i}$ and the corresponding surface is denoted by $\Gamma_{i}$; i.e., we
have as the total surface $\Gamma=\bigcup_{i=1}^{N}\Gamma_{i}.$

Figure 2. Modeling of a single edge $\Lambda_{i}$ contained in the 1D graph-
like structure with boundary nodes $\bm{x}_{i,1}$ and $\bm{x}_{i,2}$. The
cylinder $\Gamma_{i}$ has a constant radius $R_{i}$.

###  2.1. Governing constituents

The principal dependent variables characterizing the growth and decline of the
tumor mass are taken to be a set of scalar-valued fields $\phi_{\alpha}$ with
values $\phi_{\alpha}(\bm{x},t)$ at a time $t\in[0,T]$ and point
$\bm{x}\in\Omega\subset\mathbb{R}^{3}$, representing the volume fractions of
constituents in the space-time domain $\Omega\times[0,T]$. The primary feature
of our model of tumor growth is the application of the framework of continuum
mixture theory in which multiple mechanical and chemical species can exist at
a point $\bm{x}\in\Omega$ at time $t>0$. Therefore, for a medium with
${N_{\alpha}}\in\mathbb{N}$ interacting constituents, the volume fraction of
each species $\phi_{\alpha}$, $\alpha\in\\{1,\dots,{N_{\alpha}}\\}$, is
represented by a field $\phi_{\alpha}$ with the value
$\phi_{\alpha}(\bm{x},t)$ at $(\bm{x},t)$ and the property
$\sum_{\alpha}\phi_{\alpha}(\bm{x},t)=1$.

We separate the tumor volume fraction $\phi_{T}=\phi_{T}(\bm{x},t)$ into the
sum of three phases $\phi_{T}=\phi_{P}+\phi_{H}+\phi_{N}$, where
$\phi_{P}=\phi_{P}(\bm{x},t)$ is the volume fraction of proliferative cells,
$\phi_{H}=\phi_{H}(\bm{x},t)$ that of hypoxic cells, and
$\phi_{N}=\phi_{N}(\bm{x},t)$ is the volume fraction of necrotic cells, see
Figure 1. Proliferative cells have a high probability of mitosis, i.e.,
division into twin cells, and to produce growth of the tumor mass. Hypoxic
cells are those tumor cells which are deprived of sufficient nutrient to
become or remain proliferative. Lastly, necrotic cells have died due to the
lack of nutrients.

The nutrient concentration and the tumor angiogenesis factor (TAF) over
$\Omega\times[0,T]$ are represented by scalar fields
$\phi_{\sigma}=\phi_{\sigma}(\bm{x},t)$ and $\phi_{TAF}=\phi_{TAF}(\bm{x},t)$,
respectively. The tumor cells response to hypoxia, i.e., when $\phi_{\sigma}$
is below a certain threshold, is the production of an enzyme that increases
cell mobility and activates the secretion of angiogenesis promoting factors
characterized by $\phi_{{TAF}}$. As a particular case of TAFs, we consider the
vascular endothelial growth factor (VEGF), which promotes sprouting of
endothelial cells forming the tubular structure of blood vessels, which grow
into new vessels and supply nutrients to the hypoxic volume fraction
$\phi_{H}$.

Another consequence of hypoxia is the release of matrix-degrading enzymes
(MDEs), e.g., urokinase plasminogen and matrix metalloproteinases, by the
hypoxic cells. We denote the volume fraction of the MDEs by
$\phi_{{MDE}}=\phi_{{MDE}}(\bm{x},t)$. The primary feature of the MDEs is the
erosion of the extracellular matrix, whose volume fraction is denoted by
$\phi_{{ECM}}=\phi_{{ECM}}(\bm{x},t)$. Consequently, the erosion of the ECM
produces room for the invasion of tumor cells, which increases $\phi_{T}$ in
the ECM domain and therefore, raises the likelihood of metastasis. Below a
certain level necrosis occurs and cells die, entering the necrotic phase
$\phi_{N}$. Tumor cells may also die naturally, in a process which is called
apoptosis.

Within the one-dimensional network $\Lambda$, we introduce the constituents
$\phi_{v}=\phi_{v}(s,t)$ and $v_{v}=v_{v}(s,t)$, which represent the one-
dimensional counterparts of the local nutrient concentration $\phi_{\sigma}$
and the volume-averaged velocity $v$. Additionally, we consider the pressures
$p_{v}=p_{v}(s,t)$ and $p=p(\bm{x},t)$ in the network and tissue domain,
respectively. In summary, we refer to the table below for the primary
variables and constituents of the model.

Constituents
---
${\bm{\phi}}$ |  Vector of all 3D species volume fractions
$\phi_{\alpha}$ |  Volume fraction of 3D species $\alpha\in\mathcal{A}=\\{P,H,N,\sigma,MDE,TAF,ECM\\}$
$\mu_{\beta}$ |  Chemical potential, $\beta\in\\{P,H\\}$
$\phi_{v}$ |  Volume fraction of nutrients in 1D network $\Lambda$
Flow model
$v$ |  Convective velocity in tissue domain $\Omega$
$p$ |  Pressure in tissue domain $\Omega$
$p_{v}$ |  Pressure in 1D network domain $\Lambda$
$v_{v}$ |  Velocity of interstitial flow in tissue domain $\Lambda$
Functions
$\Psi$ |  Double-well potential, see Eq. 2.2
$m_{\alpha}$ |  Mobility function, see Eq. 2.1
$S_{\alpha}$ |  Mass source, see Eq. 2.5
$W_{\alpha}$ |  Wiener process, see Eq. 2.3
$J_{\sigma v}$ |  Mass source density of nutrient due to 1D network Eq. 2.6

###  2.2. Three-dimensional model governing the tumor constituents

The evolution of the constituents $\phi_{\alpha}$ must obey the balance laws
of continuum mixture theory (e.g., see [36, 23]). Assuming constant and equal
mass densities of the constituents, the mass balance equations for the mixture
read as follows:

$\partial_{t}\phi_{\alpha}+\text{div}(\phi_{\alpha}v_{\alpha})=-\text{div}J_{\alpha}(\bm{\phi})+S_{\alpha}(\bm{\phi}),$

where $v_{\alpha}$ is the cell velocity of the $\alpha$-th constituent, and
$S_{\alpha}$ describes a mass source term that may depend on all species
$\bm{\phi}=(\phi_{P},\phi_{H},\phi_{N},\phi_{\sigma},\phi_{{MDE}},\phi_{{TAF}},\phi_{{ECM}})$.
Moreover, $J_{\alpha}$ represents the mass flux of the $\alpha$-th constituent
and is given by:

$J_{\alpha}(\bm{\phi})=-m_{\alpha}(\bm{\phi})\nabla\mu_{\alpha},$

where $\mu_{\alpha}$ denotes the chemical potential of the $\alpha$-th
species, and $m_{\alpha}$ is its corresponding mobility function. Generally,
the mobilities may depend on many species, but in this work we consider the
following cases,

(2.1) $\displaystyle m_{\alpha}(\bm{\phi})$
$\displaystyle=M_{\alpha}\phi_{\alpha}^{2}(1-\phi_{T})^{2}I_{d},$
$\displaystyle\alpha\in\\{P,H\\},$ $\displaystyle m_{\beta}(\bm{\phi})$
$\displaystyle=M_{\beta}I_{d},$
$\displaystyle\beta\in\\{\sigma,{MDE},{TAF}\\},$

where $M_{\alpha}$ are mobility constants, and $I_{d}$ is the $(d\times
d)$-dimensional identity matrix. For the remaining species $\phi_{N}$ and
$\phi_{{ECM}}$, we choose $m_{N}=m_{{ECM}}=0$ in accordance to the non-
diffusivity of the necrotic cells and the ECM; see [41]. Following [28, 36,
65], we define the chemical potential $\mu_{\alpha}$ as the first variation
(Gâteaux derivative) with respect to $\phi_{\alpha}$ of the
Ginzburg–Landau–Helmholtz free energy functional $\mathcal{E}(\bm{\phi})$. The
free energy in this work is designed to capture the following key effects:

  * •

Phase change in tumor species $\phi_{T},\phi_{P},\phi_{H}$. For example,
$\phi_{T}$ can change (conditions permitting) from a healthy phase
$\phi_{T}=0$ to a tumor phase $\phi_{T}=1$. This is typically achieved by
introducing a double-well potential

(2.2) $
\Psi=\Psi(\phi_{T},\phi_{P},\phi_{H})=\sum_{\alpha\in\\{T,P,H\\}}C_{\Psi_{\alpha}}\phi_{\alpha}^{2}(1-\phi_{\alpha})^{2}$

to the free energy, where $C_{\Psi_{\alpha}}$, $\alpha\in\\{T,P,H\\}$, are
constants. In addition to phase separation between healthy and cancer phases
(using the energy term $C_{\Psi_{T}}\phi_{T}^{2}(1-\phi_{T})^{2}$), we have
also introduced energy terms that promote phase separation between
proliferative and non-proliferative and hypoxic and non-hypoxic phases. It is
possible to consider different forms of the double-well potential [21],
however, in this work we will consider $\Psi$ in Eq. 2.2 with
$C_{\Psi_{P}}=C_{\Psi_{H}}=0$, see Section 5.

  * •

Promote phase separation between two phases of species
$\phi_{T},\phi_{P},\phi_{H}$. For example, a model could exhibit phase values
at $\bm{x}$ between, say, $\phi_{\alpha}=0$ and $\phi_{\alpha}=1$, with a
change in gradient, $\nabla\phi_{\alpha}$, at the interface of these phases.
Such changes are manifested as surface energy terms in the form of penalties
on the magnitude of $\nabla\phi_{\alpha}$ of the form

$\frac{\varepsilon_{\alpha}^{2}}{2}|\nabla\phi_{\alpha}|^{2},$

where $\varepsilon_{\alpha}$ controls the thickness of the phase interface.

  * •

Diffusion driven mobilities of species
$\phi_{\sigma},\phi_{{TAF}},\phi_{{MDE}}$. These effects are captured by
adding the diffusive energies

$\frac{D_{\beta}}{2}\phi_{\beta}^{2},$

where $D_{\beta}$, $\beta\in\\{\sigma,{TAF},{MDE}\\}$, are diffusion
coefficients.

  * •

Chemotaxis and haptotaxis effects. Chemotaxis represents a movement of cells
towards a gradient of nutrients (i.e., along the direction of increasing
nutrients). Similar to chemotaxis, the tumor cells show a tendency to move
along the ECM gradient, and this phenomenon is referred to as haptotaxis.
These effects are incorporated via the terms [29, 61]

$-(\chi_{c}\phi_{\sigma}+\chi_{h}\phi_{{ECM}})\sum_{\alpha\in\\{P,H\\}}\phi_{\alpha},$

where $\chi_{c},\chi_{h}$ are chemotaxis and haptotaxis coefficients,
respectively. In the above energy terms, we exclude necrotic cells to be
consistent with our assumption that necrotic cells are immobile.

Combining these effects, the free energy takes the form

$\mathcal{E}(\bm{\phi})=\int_{\Omega}\Big{\\{}\Psi(\phi_{P},\phi_{H},\phi_{N})+\sum_{\alpha\in\\{P,H\\}}\frac{\varepsilon_{\alpha}^{2}}{2}|\nabla\phi_{\alpha}|^{2}+\sum_{\beta\in{\mathcal{R}\mathcal{D}}}\frac{D_{\beta}}{2}\phi_{\beta}^{2}-(\chi_{c}\phi_{\sigma}+\chi_{h}\phi_{{ECM}})\sum_{\alpha\in\\{P,H\\}}\phi_{\alpha}\Big{\\}}\text{
d}\bm{x},$

where ${\mathcal{R}\mathcal{D}}=\\{\sigma,{MDE},{TAF},{ECM}\\}$ is the set of
species driven by reaction-diffusion type equations. We assume a volume-
averaged velocity $v$ for the proliferative cells, hypoxic cells, and the
nutrients concentration. This assumption is regarded as reasonable whenever
cells are tightly packed.

In thin subdomains at the interfaces of the phase fields, stochastic
variations of the phase concentrations are possible. The variations in these
regions of random behavior are bounded by noise parameters
$\phi_{\alpha}^{\omega}$ and noise intensity $\omega_{\alpha}$; the variations
(along with the noise intensity) in $\phi_{\alpha}$, $\alpha\in\\{P,H\\}$, are
restricted to interface regions using function $G_{\alpha}$ given by

(2.3)
$G_{\alpha}(\phi_{P},\phi_{H},\phi_{N})=\omega_{\alpha}\mathcal{H}((\phi_{\alpha}-\phi_{\alpha}^{\omega})(1-\phi_{\alpha}-\phi_{\alpha}^{\omega}))\mathcal{H}((\phi_{T}-\phi_{T}^{\omega})(1-\phi_{T}-\phi_{T}^{\omega})).$

Here, $\mathcal{H}$ denotes the Heaviside step function. Typically, the
randomness in the evolution of species near the interface is incorporated in
the model in the form of cylindrical Wiener process on $L^{2}(\Omega)$, see
[14, 44, 4]; we add $G_{P}\dot{W}_{P}$ and $G_{H}\dot{W}_{H}$ to the mass
balance equation for $\phi_{P}$ and $\phi_{H}$. To keep the mass balance
equations in standard form, we slightly abuse the standard notation and use
notation $\dot{W}_{\alpha}$ such that
$\dot{W}_{\alpha}\textup{d}t=\textup{d}W_{\alpha}$. Further details on Wiener
processes $W_{\alpha}$ and numerical discretization are provided in Section 4.

Following these assumptions and conventions, we arrive at the following system
of equations governing the model:

(2.4) $\displaystyle\partial_{t}\phi_{P}+\textup{div}(\phi_{P}v)$
$\displaystyle=\textup{div}(m_{P}(\bm{\phi})\nabla\mu_{P})+S_{P}(\bm{\phi})+G_{P}(\phi_{P},\phi_{H},\phi_{N})\dot{W}_{P},$
$\displaystyle\mu_{P}$
$\displaystyle=\partial_{\phi_{P}}\Psi(\phi_{P},\phi_{H},\phi_{N})-\varepsilon^{2}_{P}\Delta\phi_{P}-\chi_{c}\phi_{\sigma}-\chi_{h}\phi_{{ECM}},$
$\displaystyle\partial_{t}\phi_{H}+\textup{div}(\phi_{H}v)$
$\displaystyle=\textup{div}(m_{H}(\bm{\phi})\nabla\mu_{H})+S_{H}(\bm{\phi})+G_{H}(\phi_{P},\phi_{H},\phi_{N})\dot{W}_{H},$
$\displaystyle\mu_{H}$
$\displaystyle=\partial_{\phi_{H}}\Psi(\phi_{P},\phi_{H},\phi_{N})-\varepsilon^{2}_{H}\Delta\phi_{H}-\chi_{c}\phi_{\sigma}-\chi_{h}\phi_{{ECM}},$
$\displaystyle\partial_{t}\phi_{N}$ $\displaystyle=S_{N}(\bm{\phi}),$
$\displaystyle\partial_{t}\phi_{\sigma}+\textup{div}(\phi_{\sigma}v)$
$\displaystyle=\textup{div}(m_{\sigma}(\bm{\phi})(D_{\sigma}\nabla\phi_{\sigma}\\!-\\!\chi_{c}\nabla(\phi_{P}+\phi_{H}))+S_{\sigma}(\bm{\phi})+J_{\sigma
v}(\phi_{\sigma},p,\Pi_{\Gamma}\phi_{v},\Pi_{\Gamma}p_{v})\delta_{\Gamma},$
$\displaystyle\partial_{t}\phi_{{MDE}}+\textup{div}(\phi_{{MDE}}v)$
$\displaystyle=\textup{div}(m_{{MDE}}(\bm{\phi})D_{{MDE}}\nabla\phi_{{MDE}})+S_{{MDE}}(\bm{\phi}),$
$\displaystyle\partial_{t}\phi_{{TAF}}+\textup{div}(\phi_{{TAF}}v)$
$\displaystyle=\textup{div}(m_{{TAF}}(\bm{\phi})D_{{TAF}}\nabla\phi_{{TAF}})+S_{{TAF}}(\bm{\phi}),$
$\displaystyle\partial_{t}\phi_{{ECM}}$ $\displaystyle=S_{{ECM}}(\bm{\phi}),$
$\displaystyle-\textup{div}(K\nabla p)$
$\displaystyle=J_{pv}(p,\Pi_{\Gamma}p_{v})\delta_{\Gamma}-\textup{div}(KS_{p}(\bm{\phi},\mu_{P},\mu_{H})),$
$\displaystyle v$ $\displaystyle=-K(\nabla
p-S_{p}(\bm{\phi},\mu_{P},\mu_{H})),$

in the space-time domain $\Omega\times(0,T)$ and we supplement the system with
homogeneous Neumann boundary conditions. In the above set of governing
equations, the velocity $v$ is given by modified Darcy’s law, where $K$
denotes the hydraulic conductivity. The source term $S_{p}$ (defined below)
represents a form of the elastic Korteweg force, e.g., see [20], and includes
a correction of the chemical potential by the haptotaxis and chemotaxis
adhesion terms following [24]. Here $J_{pv}$ and $J_{\sigma v}$ are the fluid
flux and nutrient flux as described in Subsection 2.3. We consider the
following choices of the coupling source functions; see [21],

(2.5) $\displaystyle S_{P}(\bm{\phi})$
$\displaystyle=\lambda^{\\!\textup{pro}}_{P}\phi_{\sigma}\phi_{P}(1-\phi_{T})-\lambda^{\\!\textup{deg}}_{P}\phi_{P}-\lambda_{P\\!H}\mathcal{H}(\sigma_{P\\!H}-\phi_{\sigma})\phi_{P}+\lambda_{H\\!P}\mathcal{H}(\phi_{\sigma}-\sigma_{H\\!P})\phi_{H},$
$\displaystyle S_{H}(\bm{\phi})$
$\displaystyle=\lambda^{\\!\textup{pro}}_{H}\phi_{\sigma}\phi_{H}(1-\phi_{T})-\lambda^{\\!\textup{deg}}_{H}\phi_{H}+\lambda_{P\\!H}\mathcal{H}(\sigma_{P\\!H}-\phi_{\sigma})\phi_{P}-\lambda_{H\\!P}\mathcal{H}(\phi_{\sigma}-\sigma_{H\\!P})\phi_{H}$
$\displaystyle\qquad-\lambda_{H\\!N}\mathcal{H}(\sigma_{H\\!N}-\phi_{\sigma})\phi_{H},$
$\displaystyle S_{N}(\bm{\phi})$
$\displaystyle=\lambda_{H\\!N}\mathcal{H}(\sigma_{H\\!N}-\phi_{\sigma})\phi_{H},$
$\displaystyle S_{{ECM}}(\bm{\phi})$
$\displaystyle=-\lambda^{\\!\textup{deg}}_{{ECM}}\phi_{{ECM}}\phi_{{MDE}}+\lambda^{\\!\textup{pro}}_{{ECM}}\phi_{\sigma}(1-\phi_{{ECM}})\mathcal{H}(\phi_{{ECM}}-\phi^{\text{pro}}_{{ECM}}),$
$\displaystyle S_{\sigma}(\bm{\phi})$
$\displaystyle=-\lambda^{\\!\textup{pro}}_{P}\phi_{\sigma}\phi_{P}-\lambda^{\\!\textup{pro}}_{H}\phi_{\sigma}\phi_{H}+\lambda^{\\!\textup{deg}}_{P}\phi_{P}+\lambda^{\\!\textup{deg}}_{H}\phi_{H}-\lambda^{\\!\textup{pro}}_{{ECM}}\phi_{\sigma}(1-\phi_{{ECM}})\mathcal{H}(\phi_{{ECM}}-\phi^{\text{pro}}_{{ECM}})$
$\displaystyle\qquad+\lambda^{\\!\textup{deg}}_{{ECM}}\phi_{{ECM}}\phi_{{MDE}},$
$\displaystyle S_{{MDE}}(\bm{\phi})$
$\displaystyle=-\lambda^{\\!\textup{deg}}_{{MDE}}\phi_{{MDE}}+\lambda^{\\!\textup{pro}}_{{MDE}}(\phi_{P}+\phi_{H})\phi_{{ECM}}\frac{\sigma_{H\\!P}}{\sigma_{H\\!P}+\phi_{\sigma}}(1-\phi_{{MDE}})-\lambda^{\\!\textup{deg}}_{{ECM}}\phi_{{ECM}}\phi_{{MDE}},$
$\displaystyle S_{{TAF}}(\bm{\phi})$
$\displaystyle=\lambda^{\\!\textup{pro}}_{{TAF}}(1-\phi_{{TAF}})\phi_{H}\mathcal{H}(\phi_{H}-\phi_{H_{P}})-\lambda_{TAF}^{\deg}\phi_{TAF},$
$\displaystyle S_{p}(\bm{\phi},\mu_{P},\mu_{H})$
$\displaystyle=(\mu_{P}+\chi_{c}\phi_{\sigma}+\chi_{h}\phi_{{ECM}})\nabla\phi_{P}+(\mu_{H}+\chi_{c}\phi_{\sigma}+\chi_{h}\phi_{{ECM}})\nabla\phi_{H}.$

Here, $\lambda^{\\!\textup{pro}}_{\alpha}$ and
$\lambda^{\\!\textup{deg}}_{\alpha}$ denote the proliferation and degradation
rate of the $\alpha$-th species, respectively, $\lambda_{\alpha\beta}$ the
transition rate from the $\alpha$-th to the $\beta$-th volume fraction,
$\sigma_{\alpha\beta}$ the corresponding nutrient threshold for the
transition, and $\mathcal{H}$ is the Heaviside step function. Further,
$\phi^{\text{pro}}_{{ECM}}$ denotes the threshold level for the ECM density in
order to activate the production of ECM fibers. Moreover, we introduce the
projection $\Pi_{\Gamma}$ of the 1D quantities onto the cylinder $\Gamma$ via
extending its function values $\Pi_{\Gamma}\phi_{v}(s)=\phi_{v}(s_{i})$ for
all $s\in\partial B_{R_{i}}(s_{i})$.

###  2.3. Interaction between the 3D and 1D model

We apply the Kedem–Katchalsky law [26] to quantify the flux of nutrients
across the vessel surface; i.e., $J_{\sigma v}$ in Eq. 2.4 is given by

(2.6) $J_{\sigma
v}(\overline{\phi}_{\sigma},\overline{p},\phi_{v},p_{v})=(1-r_{\sigma})J_{pv}(\overline{p},p_{v})\phi_{\sigma}^{v}+L_{\sigma}(\phi_{v}-\overline{\phi}_{\sigma}),$

where $J_{pv}$ denotes the flux, which is caused by the flux of blood plasma
from the vessels into the tissue or vice versa. Further, $J_{pv}$ is governed
by Starling’s law [56], i.e.,
$J_{pv}(\overline{p},p_{v})=L_{p}(p_{v}-\overline{p})$ where $\overline{p}$
denotes an averaged pressure over the circumference of cylinder cross-sections
and is computed in the following way: For each point $s_{i}$ on the curve
$\Lambda_{i}$, we consider the circle $\partial B_{R_{i}}(s_{i})$ of radius
$R_{i}$, which is perpendicular to $\Lambda_{i}$; see Figure 2. Thus, the
tissue pressure $p$ is averaged with respect to $\partial B_{R_{i}}(s_{i})$,

$\overline{p}(s_{i})=\frac{1}{2\pi R_{i}}\int_{\partial
B_{R_{i}}(s_{i})}p|_{\Gamma}\,\textup{d}S.$

The part $J_{pv}\phi_{\sigma}^{v}$ in the Kedem–Katchalsky law Eq. 2.6 is
weighted by a factor $1-r_{\sigma}$, $r_{\sigma}$ being a reflection
parameter, introduced to account for the permeability of the vessel wall with
respect to the nutrients. The value of $\phi_{\sigma}^{v}$ is set to
$\phi_{v}$ for $p_{v}\geq\overline{p}$ and to $\overline{\phi}_{\sigma}$
otherwise. The second term on the right hand side of Eq. 2.6 is a Fickian type
law, accounting for the tendency of the nutrients to balance their
concentration levels, where the permeability of the vessel wall is represented
by the parameter $L_{\sigma}$.

The interaction between the vascular network and the tissue occur at the
vessel surface $\Gamma$, and thus, we concentrate the flux $J_{\sigma v}$ by
means of the Dirac measure $\delta_{\Gamma}$; i.e., we define

$\delta_{\Gamma}(\varphi)=\int_{\Gamma}\varphi|_{\Gamma}\,\textup{d}S,$

for a sufficiently smooth test function $\varphi$ with compact support.

###  2.4. One-dimensional model for transport in the vascular network

The one-dimensional vessel variables $\phi_{v}$ and $p_{v}$ represent averages
across cross-section of the blood vessels. Thus, the one-dimensional variables
$\phi_{v}$ and $p_{v}$ on a 1D vessel $\Lambda_{i}$, $i\in\\{1,\dots,N\\}$,
depend only on $s_{i}$. See also [33] for more details related to the
derivation of the 1D pipe flow and transport models. With these conventions,
the 1D model equations for flow and transport on $\Lambda_{i}$ are given by

(2.7) $\displaystyle\partial_{t}\phi_{v}+\partial_{s_{i}}(v_{v}\phi_{v})$
$\displaystyle=\partial_{s_{i}}(m_{v}(\phi_{v})D_{v}\partial_{s_{i}}\phi_{v})-2\pi
R_{i}J_{\sigma v}(\overline{\phi}_{\sigma},\overline{p},\phi_{v},p_{v}),$
$\displaystyle R_{i}^{2}\pi\;\partial_{s_{i}}(K_{v,i}\;\partial_{s_{i}}p_{v})$
$\displaystyle=2\pi R_{i}J_{pv}(\overline{p},p_{v}).$

Here, we have introduced the permeability
$K_{v,i}=\tfrac{1}{8}R_{i}^{2}/\mu_{bl}$ of the $i$-th vessel with $\mu_{bl}$
being the viscosity of blood. We assign $\mu_{bl}$ a constant value, i.e.,
non-Newtonian behavior of blood is not considered. The diffusivity parameter
$D_{v}$ is set to the same value as $D_{\sigma}$. The blood velocity $v_{v}$
is given by the Darcy equation $v_{v}=-K_{v,i}\partial_{s_{i}}p_{v}.$

In order to interconnect $p_{v}$ and $\phi_{v}$ on $\Lambda_{i}$ at the inner
networks nodes on the intersections
$\bm{x}\in\partial\Lambda_{i}\setminus\partial\Lambda,$ we require continuity
conditions on $p_{v}$ and $\phi_{v}$. Moreover, we enforce conservation of
mass to obtain a physically relevant solution. To formulate these coupling
conditions in a mathematical way, we define for each bifurcation point
$\bm{x}$ an index set

$N(\bm{x})=\mathopen{}\mathclose{{}\left\\{\mathopen{}\mathclose{{}\left.i\in\mathopen{}\mathclose{{}\left\\{1,\ldots,N}\right\\}\;}\right|\;\bm{x}\in\partial\Lambda_{i}}\right\\}.$

We state the following continuity and mass conservation conditions at an inner
node $\bm{x}\in\partial\Lambda_{i}$:

$\displaystyle p_{v}|_{\Lambda_{i}}(\bm{x})-p_{v}|_{\Lambda_{j}}(\bm{x})$
$\displaystyle=0,\quad\text{ for all }j\in N(\bm{x})\backslash\\{i\\},$
$\displaystyle\phi_{v}|_{\Lambda_{i}}(\bm{x})-\phi_{v}|_{\Lambda_{j}}(\bm{x})$
$\displaystyle=0,\quad\text{ for all }j\in N(\bm{x})\backslash\\{i\\},$
$\displaystyle\sum_{j\in
N(\bm{x})}-\frac{R_{j}^{4}\pi}{8\mu_{\mathrm{bl}}}\frac{\partial
p_{v}}{\partial s_{j}}\Big{|}_{\Lambda_{j}}(\bm{x})$ $\displaystyle=0,$
$\displaystyle\sum_{j\in
N(\bm{x})}\Big{(}v_{v}\phi_{v}-m_{v}(\phi_{v})D_{v}\frac{\partial\phi_{v}}{\partial
s_{j}}\Big{)}\Big{|}_{\Lambda_{j}}(\bm{x})$ $\displaystyle=0.$

##  3\. Angiogenesis: Network Growth Algorithm

As noted earlier, angiogenesis is triggered by an increased TAF concentration
$\phi_{{TAF}}$ around the pre-existing blood vessels. After the TAF molecules
are emitted by the dying hypoxic tumor cells, they move through the tissue
matrix and may encounter sensor ligands on the vessel surfaces. If the TAF
concentration is large enough at the vessel surfaces, an increased number of
sensors in the vessel wall are activated and a reproduction of endothelial
cells forming the vessel walls is initiated. As a result, the affected vessels
can elongate, resulting in two different kinds of vessel elongations or
growth. In medical literature, see, e.g., [54, 16], this process is referred
to as apical growth and sprouting of vessels. The term apical growth is
derived from the term apex denoting the tip of a blood vessel, i.e., apical
growth is the type of growth occurring at the tip of a vessel. On the other
hand, the sprouting of new vessels results in the formation of new vessels at
other places on the vessel surface. In order to increase or decrease the flow
of blood and nutrients through the vessels, it is observed that the newly
formed blood vessels can adapt their vessel radii which is caused, e.g., by an
increased wall shear stress at the inner side of the vessel walls. Combining
these mechanisms, an increased supply of nutrients for both the healthy and
cancerous tissue can be achieved such that 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 vessels is simulated. At the
end of an angiogenesis step, the radii of the vessels are adapted to regulate
the blood flow. The 1D network that is updated during an angiogenesis step is
denoted by $\Lambda_{{\text{old}}}$.

###  3.1. Apical growth

Since the apical growth occurs only at the tips of the blood vessels, we
consider all the boundary nodes $\bm{x}$ of the network
$\Lambda_{{\text{old}}}$ contained in the inner part of $\Omega$, i.e.,
$\bm{x}\in\partial\Lambda_{{\text{old}}}$ and $\bm{x}\notin\partial\Omega$.
Moreover, we assume that $\bm{x}$ is contained in the segment
$\Lambda_{i}\subset\Lambda_{{\text{old}}}$. At $\bm{x}$, the value of the TAF
concentration is denoted by $\phi_{{TAF}}(\bm{x})$. If this value exceeds a
certain threshold $Th_{{TAF}}$: $\phi_{{TAF}}(\bm{x})\geq Th_{{TAF}}$, the tip
of the corresponding 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 can further elongate or it can bifurcate. In order
to decide which event occurs, a probabilistic method is used. According to
[57] and the references therein, the ratio $r_{i}=l_{i}/R_{i}$ of the vessel
$\Lambda_{i}$ follows a log-normal distribution:

(3.1)
$p_{b}(r)\sim\mathcal{L}\mathcal{N}(r,\mu_{r},\sigma_{r})=\frac{1}{r\sqrt{2\pi\sigma_{r}^{2}}}\exp\bigg{(}-\frac{({\ln
r}-\mu_{r})^{2}}{2\sigma_{r}^{2}}\bigg{)}.$

The parameters $\mu_{r}$ and $\sigma_{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 considered or not. This means that a bifurcation at
$\bm{x}\in\partial\Lambda_{i}\cup\partial\Lambda$ is formed with a probability
of:

(3.2) $P_{b}(r)=\Phi\bigg{(}\frac{\ln
r-\mu_{r}}{\sigma_{r}}\bigg{)}=\frac{1}{2}+\frac{1}{2}\text{erf}\bigg{(}\frac{\ln
r-\mu_{r}}{\sqrt{2\sigma_{r}^{2}}}\bigg{)},$

where $\Phi$ denotes the standard normal cumulative distribution function and
$x\mapsto\text{erf}(x)$ the Gaussian error function. We refer to Figure 3 for
the illustration of an exemplary vessel, which bifurcates. Moreover, we depict
the distribution of the ratio $l_{i}/R_{i}$ according to Eq. 3.1, the radii of
its bifurcations, see Eq. 3.5 below, and the probability of the occurrence of
a bifurcation, see Eq. 3.2.

Figure 3. Given a vessel with length $l_{i}$ and radius $R_{i}$, we plot the
probability of the occurrence 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 curve in figure (c)); we
choose $R_{i}=1.5\cdot 10^{-2}$, $R_{c}=2^{-\frac{1}{3}}R_{i}$ according to
Eq. 3.5, $\mu_{r}=1$, $\sigma_{r}=0.2$, $R_{\min}=9\cdot 10^{-3}$,
$R_{\max}=3.5\cdot 10^{-2}$ according to Table 2,
$\overline{R}_{\max}=\max\\{R_{\max},R_{i}\\}=R_{i}$.

If a single vessel is formed at $\bm{x}$, the direction of growth
$\mathbf{d}_{g}$ is based on the TAF concentration:

(3.3)
$\mathbf{d}_{g}(\bm{x})=\frac{\nabla\phi_{{TAF}}(\bm{x})}{\mathopen{}\mathclose{{}\left\|\nabla\phi_{{TAF}}(\bm{x})}\right\|}+\lambda_{g}\frac{\mathbf{d}_{i}}{\mathopen{}\mathclose{{}\left\|\mathbf{d}_{i}}\right\|},$

where $\|\cdot\|$ denotes the Euclidean norm. The vector
$\mathbf{d}_{i}=\bm{x}_{i,2}-\bm{x}_{i,1}$ is the orientation of the vessel
$\Lambda_{i}$, and the value
$\lambda_{g}\in\mathopen{}\mathclose{{}\left(0,1}\right]$ represents a
regularization parameter that can be used to circumvent the formation of sharp
bendings and corners. This is necessary if the TAF gradient at $\bm{x}$
encloses an acute angle with $\mathbf{d}_{i}$. The radius $R_{i^{\prime}}$ of
the new vessel $\Lambda_{i^{\prime}}$ is taken over from $\Lambda_{i}$ i.e.
$R_{i^{\prime}}=R_{i}$. Having the radius $R_{i^{\prime}}$ at hand, we use
(3.2) to determine the length $l_{i^{\prime}}$ of $\Lambda_{i^{\prime}}$.
Before $\Lambda_{i^{\prime}}$ is incorporated into the network
$\Lambda_{{\text{old}}}$, we check, whether it intersects another vessel in
the network. If this is the case, $\Lambda_{i^{\prime}}$ is not added to
$\Lambda_{{\text{old}}}$. In order to test whether a new vessel intersects an
existing vessel that is not directly connected, we compute the distance
between the centerlines of the new vessel and the existing vessel. If this
distance is smaller than the sum of the radii for any of the existing vessels,
the new vessel is considered too close to existing vessels, and, therefore,
the new vessel is not inserted into the network.

In the case of bifurcations, we have to choose the radii, orientations and
lengths of the new branches $b_{1}$ and $b_{2}$. 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_{1}}$ of branch $b_{1}$ and the radius
$R_{i,b_{2}}$ of branch $b_{2}$ as follows [40]:

(3.4) $R_{i}^{\gamma}=R_{i,b_{1}}^{\gamma}+R_{i,b_{2}}^{\gamma},$

where $\gamma$ denotes the bifurcation exponent. According to [57], $\gamma$
can vary between $2.0$ and $3.5$. In addition to (3.4), we require an
additional equation to determine the radii of the branches. Towards this end,
it is assumed that $R_{b_{1}}$ follows a truncated Gaussian normal
distribution:

(3.5) $R_{c}=2^{-\frac{1}{\gamma}}R_{i},\;\qquad
R_{b_{k}}\sim\mathcal{N}^{t}(R,\mu=R_{c},\sigma=R_{c}/35),\;\qquad
k\in\mathopen{}\mathclose{{}\left\\{1,2}\right\\},$

which is set to zero outside of the interval $[R_{\min},\overline{R}_{\max}]$
with $\overline{R}_{\max}=\max\\{R_{\max},R_{i}\\}$; we refer to Table 2 for a
choice of parameters for $R_{\min}$ and $R_{\max}$. Additionally, the radius
of the parent vessel acts as a natural bound for the radius of its
bifurcations.

The selection of $R_{b_{k}}$ is motivated as follows: Using the radius $R_{i}$
of $\Lambda_{i}$, we compute the expected radius $R_{c}$ resulting from
Murray’s law for a symmetric bifurcation $(R_{b_{1}}=R_{b_{2}})$. Here,
$R_{c}$ is used as a mean value for a Gaussian normal distribution, with a
small standard deviation. This yields bifurcations that are slightly deviating
from a symmetric bifurcation which is in accordance with Murray’s law. Having
$R_{b_{1}}$ and $R_{b_{2}}$ at hand, we compute the corresponding lengths
$l_{b_{1}}$ and $l_{b_{2}}$ 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 case is a symmetric bifurcation, that means
both radii which correspond to the mean. Further, we also depict two
asymmetric cases where the radii deviate from the mean.

Figure 4. Distribution of the radii of the bifurcating vessels, choosing
$R_{i}=0.015$, $R_{c}=2^{-\frac{1}{3}}R_{i}$. Examples of bifurcations with
different radii are given, $R=1.08\cdot 10^{-2}$ (case (a)), $R=R_{c}$ (case
(b)), $R=1.25\cdot 10^{-2}$ (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 $\mathbf{n}_{p}$ of this plane is given by the
cross product of the vessel orientation $\mathbf{d}_{i}$ and the growth
direction $\mathbf{d}_{g}$ from the non-bifurcating case:

(3.6)
$\mathbf{n}_{p}(\bm{x})=\frac{\mathbf{d}_{i}\times\mathbf{d}_{g}}{\mathopen{}\mathclose{{}\left\|\mathbf{d}_{i}\times\mathbf{d}_{g}}\right\|}.$

The exact location of the plane is determined such that the vessel
$\Lambda_{i}$ is contained in this plane. Further constraints for the
bifurcation configuration are related to the bifurcation angles. In [40, 39],
it is shown how optimality principles like minimum work and minimum energy
dissipation can be utilized to derive formulas relating the radii of the
branches to the branching angles $\alpha_{i}^{(1)}$ and $\alpha_{i}^{(2)}$:

(3.7)
$\cos\big{(}\alpha_{i}^{(1)}\big{)}=\frac{R_{i}^{4}+R_{b_{1}}^{4}-R_{b_{2}}^{4}}{2\cdot
R_{i}^{2}R_{b_{1}}^{2}}\;\text{ and
}\;\cos\big{(}\alpha_{i}^{(2)}\big{)}=\frac{R_{i}^{4}+R_{b_{2}}^{4}-R_{b_{1}}^{4}}{2\cdot
R_{i}^{2}R_{b_{2}}^{2}}.$

The value $\alpha_{i}^{(k)}$ denotes the bifurcation angle between branch
$k\in\mathopen{}\mathclose{{}\left\\{1,2}\right\\}$ and the father vessel.
Rotating the vector $\mathbf{d}_{g}$ at $\mathbf{x}$ around the axis defined
by $\mathbf{n}_{p}(\bm{x})$ counterclockwise by
$\alpha_{i}^{(1)}+\alpha_{i}^{(2)}$, we obtain two new growth directions
$\mathbf{d}_{b_{1}}=\mathbf{d}_{g}$ and $\mathbf{d}_{b_{2}}$. These vectors
are used to define the main axes of the two cylinders representing the two
branches. This choice of the growth directions can be considered as a
compromise between the optimality principles provided by [40, 39] and the
tendency of the network to adapt its growth direction to the nutrient demand
of the surrounding tissue. At the end of the apical growth phase, we obtain a
1D network denoted by $\Lambda_{\text{ap}}$.

1 Input: Network $\Lambda_{old}$, Output: New network $\Lambda_{ap}$

2 for  _each $\mathbf{x}\in\partial\Lambda\cap\Omega$_ do

3 Compute the TAF concentration at $\mathbf{x}$: $\phi_{{TAF}}(\bm{x})$;

4 Consider the TAF threshold $Th_{{TAF}}$;

5 if  _$\phi_{{TAF}}(\bm{x})\geq Th_{{TAF}}$_ then

6 Consider the edge $\Lambda_{i}$ containing $\mathbf{x}$ i.e.
$\mathbf{x}\in\partial\Lambda_{i}\cap\partial\Lambda$;

7 $\Lambda_{i}$ has the orientation $\mathbf{d}_{i}$, the radius $R_{i}$

8 Compute the gradient $\nabla\phi_{{TAF}}(\bm{x})$;

9 Compute the new growth direction $\mathbf{d}_{g}$ using (3.3);

10 Compute the probability
$P_{b}\mathopen{}\mathclose{{}\left(\mathbf{x}}\right)$ given by (3.2);

11 Form a bifurcation with probability
$P_{b}\mathopen{}\mathclose{{}\left(\mathbf{x}}\right)$;

12 if  _a bifurcation is formed_ then

13 Determine the radii of the new branches $R_{b_{1}}$ and $R_{b_{2}}$
according to (3.4) and (3.5);

14 Compute the bifurcation angels $\alpha_{i}^{(1)}$ and $\alpha_{i}^{(2)}$
according to (3.7);

15 Rotate $\mathbf{d}_{g}(\bm{x})$ by the angle
$\alpha_{i}^{(1)}+\alpha_{i}^{(2)}$ counterclockwise around the rotation axis
defined by the vector $\mathbf{n}_{p}(\bm{x})$ (computed using (3.6)) to
obtain a second growth direction $\mathbf{d}_{b_{2}}(\bm{x})$;

16 Determine the ratios $r_{b_{1}}$ and $r_{b_{2}}$ according to the
probability distribution (3.1);

17 Construct new edges $\Lambda_{b_{1}}$ and $\Lambda_{b_{2}}$ having the
radii $R_{b_{1}}$ and $R_{b_{2}}$,

18 the lengths $l_{b_{1}}=r_{b_{1}}\cdot R_{b_{1}}$ and
$l_{b_{2}}=r_{b_{2}}\cdot R_{b_{2}}$ as well as the orientations

19
$\mathbf{d}_{b_{1}}=\mathbf{d}_{g}\mathopen{}\mathclose{{}\left(\mathbf{x}}\right)$
and $\mathbf{d}_{b_{2}}$;

20 if  _$\Lambda_{b_{1}}$ and $\Lambda_{b_{2}}$ are not intersecting and
$R_{b_{1}},R_{b_{2}}\in\mathopen{}\mathclose{{}\left[R_{\min},\overline{R}_{\max}}\right]$_
then

21 Add $\Lambda_{b_{1}}$ and $\Lambda_{b_{2}}$ to $\Lambda_{i}$ at the node
$\mathbf{x}$;

22

23 end if

24 else

25 Continue;

26

27 end if

28

29 end if

30 else

31 The radius for the new edge is set to $R_{i}$;

32 Determine the ratio $r_{i}$ according to the probability distribution
(3.1);

33 Construct a new edge $\Lambda_{i^{\prime}}$ having the radius $R_{i}$,

34 the length $l_{i^{\prime}}=r_{i}\cdot R_{i}$ and the orientation
$\mathbf{d}_{g}\mathopen{}\mathclose{{}\left(\mathbf{x}}\right)$;

35 Check whether $\Lambda_{i^{\prime}}$ intersects;

36 if  _$\Lambda_{i^{\prime}}$ is not intersecting and
$R_{i}\in\mathopen{}\mathclose{{}\left[R_{\min},\overline{R}_{\max}}\right]$_
then

37 Add $\Lambda_{i^{\prime}}$ to $\Lambda_{i}$ at the node $\mathbf{x}$;

38

39 end if

40 else

41 Continue;

42

43 end if

44

45 end if

46

47 end if

48 else

49 Continue;

50 end if

51

52 end for

Algorithm 1 Apical growth algorithm

###  3.2. Sprouting of new vessels

In the second phase of the angiogenesis process, we examine each vessel or
segment $\Lambda_{i}\subset\Lambda_{\text{ap}}$. As ligands has been already
mentioned, the sprouting of inner vessels is triggered by TAF molecules
touching some sensor ligands in the vessel walls. Therefore, we determine for
the middle region of each segment, i.e.,
$\Lambda_{i}(s_{i})\subset\Lambda_{i},\;s_{i}\in(0.25,0.75)$ at which place an
averaged TAF concentration $\overline{\phi}_{{TAF}}$ attains its maximum
$\overline{\phi}_{{TAF}}^{(\max)}$. As in the previous section
$\overline{\phi}_{{TAF}}$ is determined by means of an integral expression:

$\overline{\phi}_{{TAF}}(s_{i})=\frac{1}{2\pi R_{i}}\int_{\partial
B_{R_{i}}(s_{i})}\phi_{{TAF}}(\bm{x})\,\textup{d}S,\;s_{i}\in(0.25,0.75).$

We consider only the parameters $s_{i}\in(0.25,0.75)$, since we want to avoid
a sprouting of new vessels at the boundaries of $\Lambda_{i}$. Furthermore,
boundary edges are not considered, and we demand that the edges should have a
minimal length $l_{\min}$ to avoid the formation of tiny vessels. If
$\overline{\phi}_{{TAF}}^{(\max)}$ is larger than $Th_{{TAF}}$, we attach a
new vessel $\Lambda_{i^{\prime}}$ at $\bm{x}$. As in the case of apical
growth, the local TAF gradient is considered as the preferred growth direction
of the new vessel:

$\mathbf{d}_{g}(\bm{x})=\frac{\nabla\phi_{{TAF}}(\bm{x})}{\mathopen{}\mathclose{{}\left\|\nabla\phi_{{TAF}}(\bm{x})}\right\|}.$

In order to prevent $\Lambda_{i^{\prime}}$ from growing in the direction of
$\Lambda_{i}$, we demand that $\mathbf{d}_{g}$ draws an angle of at least
$\frac{10}{180}\pi$. The new radius $R_{i^{\prime}}$ is computed as follows:

$\tilde{R}_{i}={\zeta}R_{i},\;\tilde{R}_{i^{\prime}}=(\tilde{R}_{i}-R_{i})^{\frac{1}{\gamma}}=({\zeta}-1)^{\frac{1}{\gamma}}R_{i},\;R_{i^{\prime}}=\begin{cases}R_{i^{\prime}}\sim\mathcal{U}(1.25\cdot
R_{\min},\tilde{R}_{i^{\prime}})\text{ if }1.25\cdot
R_{\min}<\tilde{R}_{i^{\prime}}\\\ R_{\min}\text{ otherwise.}\end{cases}$

${\zeta}>1$ is a fixed parameter, $R_{\min}$ denotes the minimal radius of a
blood vessel, and $\mathcal{U}$ stands for the uniform distribution, i.e., new
segment radius $R_{i^{\prime}}$ is chosen from the interval
$[R_{\min},\tilde{R}_{i^{\prime}}]$ based on a uniform distribution. For a
given radius $\tilde{R}_{i^{\prime}}$, the new length $l_{i^{\prime}}$ of
$\Lambda_{i^{\prime}}$ is determined by means of (3.1).

Finally, three new vessels $\Lambda_{i_{1}}$, $\Lambda_{i_{2}}$ and
$\Lambda_{i^{\prime}}$ are added to the network $\Lambda_{\text{ap}}$. As in
the case of apical growth, we test whether a new vessel intersects an existing
vessel, before we incorporate $\Lambda_{i^{\prime}}$ into
$\Lambda_{\text{ap}}$. In addition, we check whether a terminal vessel, i.e.,
a vessel that is part of $\partial\Lambda_{\text{ap}}$ can be linked to
another vessel. For this purpose, the distance of the point
$\bm{x}_{b}\in\partial\Lambda_{\text{ap}}\cup\partial\Lambda_{i}$ to its
neighboring network nodes that are not directly linked to $\bm{x}_{b}$ is
computed. If the distance is below a certain threshold
$\text{dist}_{\text{link}}$, the corresponding network node is considered as a
candidate to be linked with $\bm{x}_{b}$. If $\bm{x}_{b}$ is part of an artery
or the high pressure region of $\Lambda_{\text{ap}}$, we link it preferably
with a candidate at minimal distance and whose pressure is in the low pressure
region (venous part). If $\bm{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 the vessels
$\Lambda_{i}\subset\Lambda_{\text{sp}}$ and compute for each vessel the wall
shear stress $\mathbf{\tau}_{w}$ by:

$\mathbf{\tau}_{w,i}=\frac{4.0\;\mu_{\text{bl}}}{\pi
R_{i}^{3}}\mathopen{}\mathclose{{}\left|Q_{i}}\right|,\;Q_{i}=-K_{v,i}\frac{R_{i}^{2}\pi\Delta
p_{v,i}}{l_{i}},$

where $\Delta p_{v,i}$ is the pressure drop along $\Lambda_{i}$. By means of
the 1D wall shear stress, the wall shear stress stimulus for the vessel
adaption is given by [60]:

(3.8) $S_{{\text{WSS}},i}={\ln}(\mathbf{\tau}_{w,i}+\tau_{\text{ref}}).$

Here, $\tau_{\text{ref}}$ is a constant that is included to avoid a singular
behavior at lower wall shear stresses [50]. Following the model for radius
adaptation in [58], the change in radius $\Delta R_{i}$ over a time step
$\Delta t$ is assumed to be proportional to the stimulus $S_{{\text{WSS}},i}$
and current radius $R_{i}$:

(3.9) $\Delta R_{i}=\mathopen{}\mathclose{{}\left(k_{{\text{WSS}}}\cdot
S_{{\text{WSS}},i}-k_{s}}\right)\cdot\Delta t\cdot R_{i},$

where $k_{s}$ is a constant that controls the natural shrinking tendency of
the blood vessel and $k_{{\text{WSS}}}$ a proportionality constant that
controls the effect of stimulus $S_{{\text{WSS}},i}$. Once we have $\Delta
R_{i}$, we can compute the updated radius of vessels using
$R_{\textup{new},i}=R_{i}+\Delta R_{i}$. If

$R_{\textup{new},i}\in\mathopen{}\mathclose{{}\left[R_{\min},1.25\cdot
R_{i}}\right],$

where $R_{\min}$ is some fixed constant, we update the vessel radius of
$\Lambda_{i}$. Otherwise, if $R_{i}<R_{\min}$ and
$\partial\Lambda_{i}\cup\partial\Lambda_{\text{sp}}\neq\emptyset$, the vessel
is removed from the network. Finally, after following the procedure discussed
in this section, we obtain a new network $\Lambda_{\textup{new}}$.

##  4\. Numerical Discretization

With our mathematical models for tumor growth, blood flow and nutrient
transport as well as angiogenesis processes presented in previous sections, we
now turn our attention to numerical solution strategies. Toward this end, let
us consider a time step $n$ given by the interval
$\mathopen{}\mathclose{{}\left[t_{n},t_{n+1}}\right]$, with $\Delta
t=t_{n+1}-t_{n}$. At the beginning of a time step $n$, we decide whether an
angiogenesis process has to be simulated or not. As examples of relevant
simulations, we consider an angiogenesis process after each third time step.
If angiogenesis has to be taken into account, we follow the steps described in
Section 3. Given the 1D network $\Lambda$ at the time point $t_{n}$, we first
apply the algorithm for the apical growth. Afterwards, the sprouting of new
vessels and the adaption of the vessel radii is simulated. Finally, we obtain
a new network $\Lambda_{\textup{new}}$ for the new time point $t_{n+1}$. If
the simulation of angiogenesis is omitted in the respective time step $n$,
$\Lambda$ is directly used for the simulation of the tumor growth as well as
blood flow and nutrient transport, see Figure 5.

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 used i.e. we keep the linear terms implicit and the
nonlinear terms explicit with respect to time. Discretizing the model
equations in space, standard conforming trilinear $Q1$ finite elements are
employed for the partial differential equations governing the tumor growth
(2.4), whereas the PDEs for pressure $(p)$ and nutrient transport
$(\phi_{\sigma})$ are solved by means of cell centered finite volume methods.
The computational mesh is given by a union of cubes having an edge length of
$h_{3D}$.

We use finite elements to approximate the higher-order Cahn–Hilliard type
equations as well advection-reaction-diffusion equations corresponding to the
species $\phi_{TAF},\phi_{ECM},\phi_{MDE}$ in Eq. 2.4. In order to ensure mass
conservation for both flow and nutrient transport in the interstitial domain,
finite volume schemes are taken 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-Seidel method is used, i.e., in each iteration, we first
solve the equation system for the 1D system. Then the updated 1D solution is
used to solve the equation system derived from the 3D problem. We stop the
iteration when the change in the current and previous iteration solution is
within a small tolerance. At each time step, we first solve the $(p_{v},p)$
coupled system. Afterwards the $(\phi_{v},\phi_{\sigma})$ coupled system is
solved. Next, we solve the remaining equations in the 3D system. This is
summarized in Algorithm 2. In the remainder of this section, the
discretizations of the 1D and 3D systems are outlined.

###  4.1. VGM discretization of the 1D PDEs

It remains to specify the numerical solution techniques for the 1D network
equations. The time integration is based on the implicit Euler method. For the
spatial discretization of the 1D equations, the Vascular Graph Method (VGM) is
considered. This method corresponds to a node centered finite volume method
[53, 63]. We then briefly describe this numerical method as well as the
discretization of the terms arising in the context of the 3D-1D coupling. We
restrict ourselves to the pressure equations.

As mentioned in Section 2, the 1D network is given by a graph-like structure,
consisting of edges $\Lambda_{i}\subset\Lambda$ and network nodes
$\bm{x}_{i}\in\Lambda$. In a first step, we assign to each network node
$\bm{x}_{i}$ an unknown for the pressure that is denoted by $p_{v,i}$. Let us
assume that the edges containing $\bm{x}_{i}$ are given by
$\Lambda_{i_{1}},\ldots,\Lambda_{i_{N}}$ and its midpoints by
$\mathbf{m}_{i_{1}},\ldots,\mathbf{m}_{i_{N}}$, see Figure 6.

Figure 6.  Notation for the Vascular Graph Method.

On each edge,
$\Lambda_{k}\in\mathopen{}\mathclose{{}\left\\{\Lambda_{i_{1}},\ldots,\Lambda_{i_{N}}}\right\\}$,
we consider the following PDE for the pressure; see also (2.7). For
convenience, the curve parameter is simply denoted by $s$.

$-R_{k}^{2}\pi\;\partial_{s}(K_{v,k}\;\partial_{s}p_{v})=-2\pi
R_{k}L_{p}(p_{v}-\overline{p}).$

Next, we establish for the node $\bm{x}_{i}$ a mass balance equation taking
the fluxes across the cylinders $Z_{i_{l}}$ into account. $Z_{i_{l}}$ is a
cylinder having the edge $\Lambda_{i_{l}}$ as a rotation axis and the radius
$R_{i_{l}}$. Furthermore its top and bottom facets are located at
$\mathbf{m}_{i_{l}}$ and $\bm{x}_{i}$, respectively (see Figure 6). The
corresponding curve parameters are denoted by $s(\bm{x}_{i})$ and
$s(\bm{x}_{i_{l}}),\;l\in\mathopen{}\mathclose{{}\left\\{1,\ldots,N}\right\\}$.
Accordingly, the mass balance equation reads as follows:

$-\sum_{l=1}^{N}\int_{s(\bm{x}_{i})}^{s(\mathbf{m}_{i_{l}})}R_{i_{l}}^{2}\pi\;\partial_{s}(K_{v,{i_{l}}}\;\partial_{s}p_{v})\,\textup{d}s=-2\pi
L_{p}\sum_{l=1}^{N}\int_{s(\bm{x}_{i})}^{s(\mathbf{m}_{i_{l}})}R_{i_{l}}(p_{v}-\overline{p})\,\textup{d}s.$

Integration yields:

$-\sum_{l=1}^{N}R_{i_{l}}^{2}\pi
K_{v,{i_{l}}}\mathopen{}\mathclose{{}\left.\partial_{s}p_{v}}\right|_{s(\mathbf{m}_{i_{l}})}+\sum_{l=1}^{N}R_{i_{l}}^{2}\pi
K_{v,{i_{l}}}\mathopen{}\mathclose{{}\left.\partial_{s}p_{v}}\right|_{s(\bm{x}_{i})}=-2\pi
L_{p}\sum_{l=1}^{n}\int_{s(\bm{x}_{i})}^{s(\mathbf{m}_{i_{l}})}{R_{i_{l}}(p_{v}-\overline{p})}\,\textup{d}s.$

Approximating the derivatives by central finite differences and using the mass
conservation equation (see Subsection 2.4):

$\sum_{l=1}^{N}R_{i_{l}}^{2}\pi
K_{v,{i_{l}}}\mathopen{}\mathclose{{}\left.\partial_{s}p_{v}}\right|_{s(\bm{x}_{i})}=0,$

at an inner node $\bm{x}_{i}$, it follows that

$\sum_{l=1}^{N}R_{i_{l}}^{2}\pi
K_{v,{i_{l}}}\;\frac{p_{v,i}-p_{v,i_{l}}}{l_{i_{l}}}=-2\pi
L_{p}\sum_{l=1}^{N}\int_{s(\bm{x}_{i})}^{s(\mathbf{m}_{i_{l}})}{R_{i_{l}}(p_{v}-\overline{p})}\,\textup{d}s,$

where $l_{i_{l}}$ denotes the length of the edge $\Lambda_{i_{l}}$. Denoting
the mantle surface of $Z_{i_{l}}$ by $S_{i_{l}}$, we have:

$\sum_{l=1}^{N}\frac{R_{i_{l}}^{2}\pi
K_{v,{i_{l}}}}{l_{i_{l}}}\;(p_{v,i}-p_{v,i_{l}})=-L_{p}\sum_{l=1}^{N}\mathopen{}\mathclose{{}\left|S_{i_{l}}}\right|p_{v,i}+L_{p}\sum_{l=1}^{N}\int_{S_{i_{l}}}p\,\textup{d}S.$

Computing the integrals $\int_{S_{i_{l}}}p\,\textup{d}S$, we introduce the
decomposition of $\Omega$ into $M$ finite volume cells ${CV}_{k}$:
$\Omega=\bigcup_{k=1}^{M}{CV}_{k}$. The pressure unknown assigned to
${CV}_{k}$ is given by $p_{k}$. Using this notation, one obtains:

$\int_{S_{i_{l}}}p\,\textup{d}S=\sum_{{CV}_{k}\cap
S_{i_{l}}\neq\emptyset}\int_{{CV}_{k}\cap
S_{i_{l}}}p\,\textup{d}S\approx\mathopen{}\mathclose{{}\left|S_{i_{l}}}\right|\sum_{{CV}_{k}\cap
S_{i_{l}}\neq\emptyset}\underbrace{\frac{\mathopen{}\mathclose{{}\left|{CV}_{k}\cap
S_{i_{l}}}\right|}{\mathopen{}\mathclose{{}\left|S_{i_{l}}}\right|}}_{=:w_{ki_{l}}}p_{k}=\mathopen{}\mathclose{{}\left|S_{i_{l}}}\right|\sum_{{CV}_{k}\cap
S_{i_{l}}\neq\emptyset}w_{ki_{l}}\;p_{k}.$

In order to estimate the weights $w_{ki_{l}}$ we discretize the mantle surface
$S_{i_{l}}$ by $N_{s}$ nodes. For our simulations, we used $N_{s}=400$ nodes.
$S_{i_{l}}$ intersects some finite volume cells $CV_{k}$. The number of nodes
contained in $CV_{k}$ is denoted by $N_{ki_{l}}$. Using these definitions, the
weights $w_{ki_{l}}$ are computed as follows: $w_{ki_{l}}=N_{ki_{l}}/N_{s}$.
As an example, consider Figure 7 below, where we show the discretization of
the surface of a cylinder. We note that one has to guarantee that the relation
$\sum_{{CV}_{k}\cap S_{i_{l}}\neq\emptyset}w_{ki_{l}}=1$ holds. Otherwise, a
consistent mass exchange between the vascular system and the tissue could not
be enforced. All in all, we obtain a linear system of equations for computing
the pressure values.

Figure 7.  Typical discretization of the cylinder surface (left). Cross
section through a mesh composed of finite volume cells $CV_{k}$ and a cylinder
with the mantle surface $S_{i_{l}}$. $S_{i_{l}}$ is discretized by $N_{s}$
nodes, which are contained in different finite volume cells. Nodes belonging
to different cells are colored differently. The number of nodes contained in
$CV_{k}$ is denoted by $N_{ki_{l}}$.

On closer examination, it can be noted that the corresponding matrix is
composed of four blocks, i.e., two coupling blocks as wells as a block for the
1D diffusion term and the 3D diffusion term. As said earlier, at each time
step, we decouple the 1D and 3D pressure equations and use a Block-Gauss-
Seidel iteration to solve the two systems until the 3D pressure is converged.
The discretization of the nutrient equation is exerted in a similar manner,
where the main difference consists in adding an upwinding procedure for the
convective term. At each time step, the nutrient equations are also solved
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 the transport of nutrients, we
require an initial condition for the variable $\phi_{v}$. In doing so, a
threshold $R_{T}$ for the radii 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 $\phi_{v}(t=0)=1$ and in case of a vein $\phi_{v}(t=0)=0$ is
used. When the network starts to grow, initial values for the newly created
vessels have to be provided. If a new vessel is created due to sprouting
growth, we consider the vessel or edge to which the new vessel is attached. At
the point of the given vessel, where the new vessel or edge is added,
$\phi_{v}$ is interpolated linearly. For this purpose the two values of
$\phi_{v}$ located at the nodes of the existing vessel are used. The
interpolated value is assigned to both nodes of the newly created vessel.

When apical growth takes place, a new vessel is added to
$\bm{x}\in\partial\Lambda$. In this case, we consider $\phi_{v}(\bm{x},t)$ for
a time point $t$ and assign it to the newly created node, since we assume that
no flow boundary conditions are enforced. With respect to the boundary
conditions for the 1D pressure PDE, the following distinction of cases is
made:

  * •

Dirichlet boundary for $p_{v}$ if
$\bm{x}\in\partial\Lambda\cap\partial\Omega$. In this case, we set:
$p_{v}(\bm{x})=p_{v,D}(\bm{x})$, where $p_{v,D}$ is a given Dirichlet value at
$\bm{x}$. Numerically, we enforce this boundary condition by setting in the
corresponding line 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 $\bm{y}$ is the neighbor of $\bm{x}$ on the edge
$\Lambda_{1}$. The other edges adjacent to $\bm{y}$ are denoted by
$\Lambda_{2},\ldots,\Lambda_{N}$. Then the balance equation for $\bm{y}$ has
to be adapted as follows to account for the Dirichlet boundary condition
$p_{v,D}$:

$\displaystyle\frac{R_{1}^{2}\pi
K_{v,{1}}}{l_{1}}\;(p_{v}\mathopen{}\mathclose{{}\left(\bm{y}}\right)-p_{v,D})+\sum_{j=2}^{N}\frac{R_{j}^{2}\pi
K_{v,{j}}}{l_{j}}\;(p_{v}\mathopen{}\mathclose{{}\left(\bm{y}}\right)-p_{v,j})$
$\displaystyle\quad=-L_{p}\mathopen{}\mathclose{{}\left|\tilde{S}_{1}}\right|p_{v}\mathopen{}\mathclose{{}\left(\bm{y}}\right)+L_{p}\sum_{{CV}_{k}\cap\tilde{S}_{1}\neq\emptyset}w_{k1}\;p_{k}-L_{p}\sum_{j=2}^{N}\mathopen{}\mathclose{{}\left|S_{j}}\right|p_{v}\mathopen{}\mathclose{{}\left(\bm{y}}\right)+L_{p}\mathopen{}\mathclose{{}\left|S_{j}}\right|\sum_{{CV}_{k}\cap
S_{j}\neq\emptyset,j>1}w_{kj}\;p_{k},$

where $\tilde{S}_{1}$ is the mantle surface of the cylinder covering the whole
edge $\Lambda_{1}$.

  * •

Homogeneous Neumann boundary for $p_{v}$ if
$\bm{x}\in\partial\Lambda\cap\Omega$. Let
$\bm{x}\in\partial\Lambda_{i}\cap\partial\Lambda\cap\Omega$, then we set
$\mathopen{}\mathclose{{}\left.-R_{i}^{2}\pi\;K_{v,i}\partial_{s}p_{v}}\right|_{\bm{x}}=0$,
resulting in the following discretization:

$R_{i}^{2}\pi\;K_{v,i}\frac{p_{v}(\bm{x})-p_{v}(\bm{y})}{l_{i}}=-L_{p}\mathopen{}\mathclose{{}\left|S_{i_{1}}}\right|p_{v,i}+L_{p}\mathopen{}\mathclose{{}\left|S_{i_{1}}}\right|\sum_{{CV}_{k}\cap
S_{i_{1}}\neq\emptyset}w_{ki_{l}}\;p_{k},$

where $\bm{y}\in\partial\Lambda_{i}\cap\Lambda$ and $l_{i}$ is the length of
the edge $\Lambda_{i}$.

Summarizing, we consider for the pressure in the network Dirichlet boundaries
at the boundary of the 3D domain $\Omega$ and homogeneous Neumann boundary
conditions in the inner part of $\Omega$. For the nutrients, the
implementation of boundary conditions is more challenging, since an upwinding
procedure has to be taken into account.

  * •

Dirichlet boundary for $\phi_{v}$ if
$\bm{x}\in\partial\Lambda_{i}\cap\partial\Lambda\cap\partial\Omega$ and

$v_{v}(\bm{x})\approx-
K_{v,i}\frac{p_{v}\mathopen{}\mathclose{{}\left(\bm{y}}\right)-p_{v,D}(\bm{x})}{l_{i}}>0.$

In this case, we set: $\phi_{v}(\bm{x},t)=\phi_{v,D}(\bm{x})$, where
$\phi_{v,D}$ is a given Dirichlet value at $\bm{x}$. The numerical
implementation can be exerted analogously to the case of the pressure $p_{v}$.

  * •

Homogeneous Neumann boundary for $\phi_{v}$ if
$\bm{x}\in\partial\Lambda\cap\Omega$. Let
$\bm{x}\in\partial\Lambda_{i}\cap\partial\Lambda\cap\Omega$ and we set

$(\mathopen{}\mathclose{{}\left.v_{v}\phi_{v}-D_{v}\partial_{s_{i}}\phi_{v})}\right|_{\bm{x}}=0,$

resulting in the following discretization:

$\displaystyle\frac{l_{i}}{2}\frac{\phi_{v}(\bm{x},t+\Delta
t)-\phi_{v}(\bm{x},t)}{\Delta
t}+\mathopen{}\mathclose{{}\left.v_{v}\phi_{v}}\right|_{\mathbf{m}_{i}}-D_{v}\frac{\phi_{v}(\bm{y},t+\Delta
t)-\phi_{v}(\bm{x},t+\Delta t)}{l_{i}}$ $\displaystyle\quad=-2\pi
R_{i}\mathopen{}\mathclose{{}\left[(1-r_{\sigma})\int_{s(\bm{x})}^{s(\mathbf{m}_{i})}J_{pv}(\overline{p},p_{v})\cdot\phi_{\sigma}^{v}(s,t+\Delta
t)\;\textup{d}s+L_{\sigma}\int_{s(\bm{x})}^{s(\mathbf{m}_{i})}\phi_{v}(s,t+\Delta
t)-\overline{\phi}_{\sigma}(s,t+\Delta t)\,\textup{d}s}\right],$

where $\bm{y}\in\partial\Lambda_{i}\cap\Lambda$. The integrals modeling the
exchange terms are discretized as in the case of the pressure equations.

  * •

Upwinding boundary for $\phi_{v}$ if
$\bm{x}\in\partial\Lambda\cap\partial\Lambda_{i}\cap\partial\Omega$ and

$v_{v}(\mathbf{m}_{i})\approx-
K_{v,i}\frac{p_{v}\mathopen{}\mathclose{{}\left(\bm{y}}\right)-p_{v}\mathopen{}\mathclose{{}\left(\mathbf{x}}\right)}{l_{i}}\leq
0.$

Here, we obtain the following semi-discrete equation:

$\displaystyle\frac{l_{i}}{2}\frac{\phi_{v}(\bm{x},t+\Delta
t)-\phi_{v}(\bm{x},t)}{\Delta
t}+\mathopen{}\mathclose{{}\left.v_{v}}\right|_{\mathbf{m}_{i}}\phi_{v}(\bm{y},t+\Delta
t)-\mathopen{}\mathclose{{}\left.v_{v}}\right|_{\mathbf{m}_{i}}\phi_{v}(\bm{x},t+\Delta
t)$ $\displaystyle\quad=-2\pi
R_{i}\mathopen{}\mathclose{{}\left[(1-r_{\sigma})\int_{s(\bm{x})}^{s(\mathbf{m}_{i})}J_{pv}(\overline{p},p_{v})\cdot\phi_{\sigma}^{v}(s,t+\Delta
t)\;\textup{d}s+L_{\sigma}\int_{s(\bm{x})}^{s(\mathbf{m}_{i})}\phi_{v}(s,t+\Delta
t)-\overline{\phi}_{\sigma}(s,t+\Delta t)\,\textup{d}s}\right],$

where $\bm{y}\in\partial\Lambda_{i}\cap\Lambda$.

###  4.2. Discretization of the 3D PDEs

Suppose that $\phi_{\alpha_{n}},\mu_{\alpha_{n}},p_{n},p_{v_{n}},\bm{v}_{n}$
denote the various fields at time $t_{n}$. Let $V_{h}$ be a subspace of
$H^{1}(\Omega,\mathbb{R})$ consisting of continuous piecewise trilinear
functions on a uniform mesh $\Omega_{h}$. We consider $\phi_{\alpha_{n}}\in
V_{h}$ for $\alpha\in\\{P,H,N,TAF,ECM,MDE\\}$, $\mu_{P_{n}},\mu_{H_{n}}\in
V_{h}$, and $\bm{v}_{h}\in[V_{h}]^{3}$. The test functions are denoted by
$\tilde{\phi}\in V_{h}$ for species in $V_{h}$, $\tilde{\mu}\in V_{h}$ for
chemical potentials $\mu_{P},\mu_{H}$, and $\tilde{\bm{v}}\in[V_{h}]^{3}$ for
velocity.

Given a time step $n$ and solutions
$\phi_{\alpha_{n}},\mu_{\alpha_{n}},p_{n},p_{v_{n}},\bm{v}_{n}$, we are
interested in the solution at the next time step. For the 3D-1D coupled
pressure $(p_{v_{n+1}},p_{n+1})$, as mentioned earlier, we utilize a block
Gauss-Seidel iteration, where the discretization of the 1D equation is
discussed in Subsection 4.1 and discretization of the 3D equation using
finite-volume scheme is provided in • ‣ Subsection 4.2. Similarly,
$(\phi_{v_{n+1}},\phi_{\sigma_{n+1}})$ is solved using a block Gauss-Seidel
iteration with the discretization of the 1D equation along the lines of the
discretization of the 1D pressure equation and discretization of the 3D
equation provided in • ‣ Subsection 4.2. We then solve the proliferative,
hypoxic, necrotic, MDE, ECM, and TAF systems sequentially. Once we have
pressure $p_{n+1}$, we compute the velocity $\bm{v}_{n+1}\in[V_{h}]^{3}$ using
the weak form:

(4.1) $(\bm{v}_{n+1},\tilde{\bm{v}})=(-K(\nabla
p_{n+1}-S_{p_{n}}),\tilde{\bm{v}}),\qquad\forall\tilde{\bm{v}}\in[V_{h}]^{3},$

where $S_{p_{n}}=S_{p}({\bm{\phi}}_{n},\mu_{P_{n}},\mu_{H_{n}})$, see Eq. 2.5.
For the advection terms, using the fact that $\nabla p\cdot\bm{n}=0$ on
$\partial\Omega$ and so $\bm{v}_{n+1}\cdot\bm{n}=0$ on $\partial\Omega$, we
can write

(4.2)
$\displaystyle\mathopen{}\mathclose{{}\left({\nabla\cdot(\phi_{\alpha}\bm{v}_{n+1})},{\tilde{\phi}}}\right)=-\mathopen{}\mathclose{{}\left({\phi_{\alpha}\bm{v}_{n+1}},{\nabla\tilde{\phi}}}\right),\qquad\forall\tilde{\phi}\in
V_{h},\;\forall\alpha\in\\{P,H,TAF,MDE\\}.$

In what follows, we consider the expression on the right hand side in above
equation for the advection terms.

For fields $\phi_{a}$, $a\in\\{P,H,N,\sigma,TAF,ECM,MDE\\}$, and chemical
potentials $\mu_{P},\mu_{H}$, we assume homogeneous Neumann boundary condition
on $\partial\Omega$. Next, we describe the discretization of the scalar fields
in the 3D model.

  * •

Pressure. Let ${CV}\in\Omega_{h}$ denote the typical finite volume cell and
$\sigma\in\partial{CV}$ face of a cell ${CV}$. Let
$(p_{v_{n+1}}^{k},p_{{n+1}}^{k})$ denote the pressures at $k^{\text{th}}$
iteration and time $t_{n+1}$. Suppose we have solved for $p_{v_{n+1}}^{{k+1}}$
following Subsection 4.1. To solve $p_{n+1}^{{k+1}}$, we consider, for all
${CV}\in\Omega_{h}$,

$\displaystyle-\sum_{\sigma\in\partial{CV}}\int_{\sigma}{K\nabla
p_{n+1}^{{k+1}}}\cdot\bm{n}\,\textup{d}S$
$\displaystyle=-\sum_{\sigma\in\partial{CV}}\int_{\sigma}{KS_{p_{n}}}\cdot\bm{n}\,\textup{d}S+\int_{\Gamma\cap{CV}}J_{pv}(p_{n+1}^{{k+1}},\Pi_{\Gamma}p_{v_{n+1}}^{{k+1}})\,\textup{d}S$
(4.3)
$\displaystyle=-\sum_{\sigma\in\partial{CV}}\int_{\sigma}{KS_{p_{n}}}\cdot\bm{n}\,\textup{d}S+\sum_{i=1}^{N}\int_{\Gamma_{i}\cap{CV}}J_{pv}(p_{n+1}^{{k+1}},\Pi_{\Gamma_{i}}p_{v_{n+1}}^{{k+1}})\,\textup{d}S.$

Above follows by the integration of the pressure equation in Eq. 2.4 over
${CV}$ and using the divergence theorem. Here
$J_{pv}(p,p_{v})=L_{p}(p_{v}-p)$, $\Gamma=\cup_{i=1}^{N}\Gamma_{i}$ is the
total vascular surface, and $\Pi_{\Gamma_{i}}(p_{v})$ is the projection of the
1D pressure defined on the centerline $\Lambda_{i}$ onto the surface of the
cylinder, $\Gamma_{i}$.

  * •

Nutrients. Suppose we have solved for $\phi_{v_{n+1}}^{{k+1}}$. To solve
$\phi_{{\sigma}_{n+1}}^{{k+1}}$, we consider, for all ${CV}$,

$\displaystyle\int_{{CV}}{\frac{\phi_{\sigma_{n+1}}^{{k+1}}-\phi_{\sigma_{n}}}{\Delta
t}}\,\textup{d}V+\sum_{\sigma\in\partial{CV}}\int_{\sigma}{\phi_{\sigma_{n+1}}^{{k+1}}\hat{\bm{v}}_{n+1}}\cdot\bm{n}\,\textup{d}S-\sum_{\sigma\in\partial{CV}}\int_{\sigma}{m_{\sigma}(\bm{\phi}_{n})\mathopen{}\mathclose{{}\left(D_{\sigma}\nabla\phi_{\sigma_{n+1}}^{{k+1}}}\right)}\cdot\bm{n}\,\textup{d}S$
$\displaystyle\quad+\int_{{CV}}{\lambda^{\\!\textup{pro}}_{P}\phi_{P_{n}}\phi_{\sigma_{n+1}}^{{k+1}}}\,\textup{d}V+\int_{{CV}}{\lambda^{\\!\textup{pro}}_{H}\phi_{H_{n}}\phi_{\sigma_{n+1}}^{{k+1}}}\,\textup{d}V$
$\displaystyle\quad+\int_{{CV}}{\lambda^{\\!\textup{pro}}_{ECM}(1-{\phi_{{ECM}}}_{n})\mathcal{H}({\phi_{{ECM}}}_{n}-\phi^{\textmd{pro}}_{{ECM}})\phi_{\sigma_{n+1}}^{{k+1}}}\,\textup{d}V$
$\displaystyle=\int_{{CV}}{\lambda^{\\!\textup{deg}}_{P}\phi_{P_{n}}+\lambda^{\\!\textup{deg}}_{H}\phi_{H_{n}}}\,\textup{d}V+\int_{{CV}}{\lambda^{\\!\textup{deg}}_{{ECM}}{\phi_{{ECM}}}_{n}{\phi_{{MDE}}}_{n}}\,\textup{d}V$
(4.4)
$\displaystyle-\sum_{\sigma\in\partial{CV}}\int_{\sigma}{m_{\sigma}(\bm{\phi}_{n})\chi_{c}\nabla\mathopen{}\mathclose{{}\left(\phi_{P_{n}}+\phi_{H_{n}}}\right)}\cdot\bm{n}\,\textup{d}S+\sum_{i=1}^{N}\int_{\Gamma_{i}\cap{CV}}J_{\sigma
v}(\phi_{\sigma_{n+1}}^{{k+1}},p_{n+1}^{{k+1}},\Pi_{\Gamma_{i}}\phi_{{v}_{n+1}}^{{k+1}},\Pi_{\Gamma_{i}}p_{v_{n+1}}^{{k+1}})\,\textup{d}S,$

where $J_{\sigma v}$ is given by Eq. 2.6. Noting that the velocity is

$\hat{\bm{v}}_{n+1}=-K\nabla p_{n+1}+KS_{p_{n}},$

we divide the advection term, for $\sigma\in\partial{CV}$, into two parts:

(4.5)
$\displaystyle\int_{\sigma}{\phi_{\sigma_{n+1}}^{{k+1}}\hat{\bm{v}}_{n+1}}\cdot\bm{n}\,\textup{d}S$
$\displaystyle=-K\int_{\sigma}{\phi_{\sigma_{n+1}}^{{k+1}}\nabla
p_{n+1}}\cdot\bm{n}\,\textup{d}S+K\int_{\sigma}{\phi_{\sigma_{n+1}}^{{k+1}}S_{p_{n}}}\cdot\bm{n}\,\textup{d}S.$

For the first term we apply the upwinding scheme. For the second term, we
perform quadrature approximation to compute the integral over the face
$\sigma$.

  * •

Proliferative. For a general double-well potential
$\Psi(\phi_{P},\phi_{H},\phi_{N})=\sum_{a\in\\{T,P,H\\}}C_{\Psi_{a}}\phi_{a}^{2}(1-\phi_{a})^{2}$,
we consider the convex-concave splitting, see [18], as follows

(4.6)
$\displaystyle\Psi(\phi_{P},\phi_{H},\phi_{N})=\sum_{a\in\\{T,P,H\\}}\frac{3}{2}C_{\Psi_{a}}\phi_{a}^{2}+\sum_{a\in\\{T,P,H\\}}C_{\Psi_{a}}(\phi_{a}^{4}-2\phi_{a}^{3}-\frac{1}{2}\phi_{a}^{2}).$

This results in

(4.7)
$\displaystyle\partial_{\phi_{P}}\Psi(\phi_{P},\phi_{H},\phi_{N})=\sum_{a\in\\{T,P\\}}3C_{\Psi_{a}}\phi_{a}+\sum_{a\in\\{T,P\\}}C_{\Psi_{a}}\phi_{a}(4\phi_{a}^{2}-6\phi_{a}-1).$

The expression for $\partial_{\phi_{H}}\Psi$ can be derived analogously. In
our implementation, $\phi_{P},\phi_{H},\phi_{N}$ are the main state variables
and $\phi_{T}$ is computed using $\phi_{T}=\phi_{P}+\phi_{H}+\phi_{N}$. Let
the mobility $\bar{m}_{P_{n}}$ at the current step be given by

(4.8)
$\displaystyle\bar{m}_{P_{n}}=M_{P}\mathopen{}\mathclose{{}\left[(\phi_{P_{n}})^{+}(1-\phi_{T_{n}})^{+}}\right]^{2},$

where for a field $f$, $\mathopen{}\mathclose{{}\left(f}\right)^{+}$ is the
projection onto $[0,1]$ given by

(4.9) $\displaystyle\mathopen{}\mathclose{{}\left(f}\right)^{+}$
$\displaystyle=\begin{cases}f\qquad\text{if }f\in[0,1],\\\ 0\qquad\text{if
}f\leq 0,\\\ 1\qquad\text{if }f\geq 1.\end{cases}$

We solve for $\phi_{P_{n+1}},\mu_{P_{n+1}}$ using the weak forms below

$\displaystyle\mathopen{}\mathclose{{}\left({\frac{\phi_{P_{n+1}}-\phi_{P_{n}}}{\Delta
t}},{\tilde{\phi}}}\right)-\mathopen{}\mathclose{{}\left({\phi_{P_{n+1}}\bm{v}_{n+1}},{\nabla\tilde{\phi}}}\right)+\mathopen{}\mathclose{{}\left({\bar{m}_{P_{n}}\nabla\mu_{P_{n+1}}},{\nabla\tilde{\phi}}}\right)$
$\displaystyle\quad-\mathopen{}\mathclose{{}\left({\lambda^{\\!\textup{pro}}_{P}\phi_{\sigma_{n+1}}(1-\phi_{T_{n}})^{+}\phi_{P_{n+1}}},{\tilde{\phi}}}\right)+\mathopen{}\mathclose{{}\left({\lambda^{\\!\textup{deg}}_{P}\phi_{P_{n+1}}},{\tilde{\phi}}}\right)$
$\displaystyle=\mathopen{}\mathclose{{}\left({\lambda_{HP}\mathcal{H}(\phi_{\sigma_{n+1}}-\sigma_{HP})\mathopen{}\mathclose{{}\left(\phi_{H_{n}}}\right)^{+}},{\tilde{\phi}}}\right)-\mathopen{}\mathclose{{}\left({\lambda_{P\\!H}\mathcal{H}(\sigma_{P\\!H}-\phi_{\sigma_{n+1}})\mathopen{}\mathclose{{}\left(\phi_{P_{n}}}\right)^{+}},{\tilde{\phi}}}\right)$
(4.10) $\displaystyle\quad+\frac{1}{\Delta
t}\mathopen{}\mathclose{{}\left({G_{P_{n}}\int_{t_{n}}^{t_{n+1}}\textup{d}{W}_{P}},{\tilde{\phi}}}\right)$

and

$\displaystyle\mathopen{}\mathclose{{}\left({\mu_{P_{n+1}}},{\tilde{\mu}}}\right)-\mathopen{}\mathclose{{}\left({3(C_{\Psi_{T}}+C_{\Psi_{P}})\phi_{P_{n+1}}},{\tilde{\mu}}}\right)-\mathopen{}\mathclose{{}\left({\epsilon_{P}^{2}\nabla\phi_{P_{n+1}}},{\nabla\tilde{\mu}}}\right)$
$\displaystyle=\mathopen{}\mathclose{{}\left({C_{\Psi_{T}}\phi_{T_{n}}(4\phi_{T_{n}}^{2}-6\phi_{T_{n}}-1)},{\tilde{\mu}}}\right)+\mathopen{}\mathclose{{}\left({C_{\Psi_{P}}\phi_{P_{n}}(4\phi_{P_{n}}^{2}-6\phi_{P_{n}}-1)},{\tilde{\mu}}}\right)$
(4.11)
$\displaystyle\quad+\mathopen{}\mathclose{{}\left({3C_{\Psi_{T}}(\phi_{H_{n}}+\phi_{N_{n}})},{\tilde{\mu}}}\right)-\mathopen{}\mathclose{{}\left({\chi_{c}\phi_{\sigma_{n+1}}+\chi_{h}{\phi_{{ECM}}}_{n}},{\tilde{\mu}}}\right),$

where $(\cdot)^{+}$ is the projection to $[0,1]$ defined in Eq. 4.9,
$G_{P_{n}}=G_{P}(\phi_{P_{n}},\phi_{H_{n}},\phi_{N_{n}})$ is given by Eq. 2.3,
and $W_{P}$ is the cylindrical Wiener process. We discuss the computation of
stochastic term in more detail in Subsection 4.2.1.

  * •

Hypoxic. Let the mobility $\bar{m}_{H_{n}}$ be given by

(4.12)
$\displaystyle\bar{m}_{H_{n}}=M_{H}\mathopen{}\mathclose{{}\left[(\phi_{H_{n}})^{+}(1-\phi_{T_{n}})^{+}}\right]^{2}.$

To solve for $\phi_{H_{n+1}},\mu_{H_{n+1}}$, we consider

$\displaystyle\mathopen{}\mathclose{{}\left({\frac{\phi_{H_{n+1}}-\phi_{H_{n}}}{\Delta
t}},{\tilde{\phi}}}\right)-\mathopen{}\mathclose{{}\left({\phi_{H_{n+1}}\bm{v}_{n+1}},{\nabla\tilde{\phi}}}\right)+\mathopen{}\mathclose{{}\left({\bar{m}_{H_{n}}\nabla\mu_{H_{n+1}}},{\nabla\tilde{\phi}}}\right)$
$\displaystyle\quad-\mathopen{}\mathclose{{}\left({\lambda^{\\!\textup{pro}}_{H}\phi_{\sigma_{n+1}}(1-\phi_{T_{n}})^{+}\phi_{H_{n+1}}},{\tilde{\phi}}}\right)+\mathopen{}\mathclose{{}\left({\lambda^{\\!\textup{deg}}_{H}\phi_{H_{n+1}}},{\tilde{\phi}}}\right)$
$\displaystyle=\mathopen{}\mathclose{{}\left({\lambda_{PH}\mathcal{H}(\sigma_{PH}-\phi_{\sigma_{n+1}})\mathopen{}\mathclose{{}\left(\phi_{P_{n}}}\right)^{+}},{\tilde{\phi}}}\right)-\mathopen{}\mathclose{{}\left({\lambda_{H\\!P}\mathcal{H}(\phi_{\sigma_{n+1}}-\sigma_{H\\!P})\mathopen{}\mathclose{{}\left(\phi_{H_{n}}}\right)^{+}},{\tilde{\phi}}}\right)$
(4.13)
$\displaystyle\quad-\mathopen{}\mathclose{{}\left({\lambda_{H\\!N}\mathcal{H}(\sigma_{H\\!N}-\phi_{\sigma_{n+1}})\mathopen{}\mathclose{{}\left(\phi_{H_{n}}}\right)^{+}},{\tilde{\phi}}}\right)+\frac{1}{\Delta
t}\mathopen{}\mathclose{{}\left({G_{H_{n}}\int_{t_{n}}^{t_{n+1}}\textup{d}{W}_{H}},{\tilde{\phi}}}\right)$

and

$\displaystyle\mathopen{}\mathclose{{}\left({\mu_{H_{n+1}}},{\tilde{\mu}}}\right)-\mathopen{}\mathclose{{}\left({3(C_{\Psi_{T}}+C_{\Psi_{H}})\phi_{H_{n+1}}},{\tilde{\mu}}}\right)-\mathopen{}\mathclose{{}\left({\epsilon_{H}^{2}\nabla\phi_{H_{n+1}}},{\nabla\tilde{\mu}}}\right)$
$\displaystyle=\mathopen{}\mathclose{{}\left({C_{\Psi_{T}}\phi_{T_{n}}(4\phi_{T_{n}}^{2}-6\phi_{T_{n}}-1)},{\tilde{\mu}}}\right)+\mathopen{}\mathclose{{}\left({C_{\Psi_{H}}\phi_{H_{n}}(4\phi_{H_{n}}^{2}-6\phi_{H_{n}}-1)},{\tilde{\mu}}}\right)$
(4.14)
$\displaystyle\quad+\mathopen{}\mathclose{{}\left({3C_{\Psi_{T}}(\phi_{P_{n}}+\phi_{N_{n}})},{\tilde{\mu}}}\right)-\mathopen{}\mathclose{{}\left({\chi_{c}\phi_{\sigma_{n+1}}+\chi_{h}{\phi_{{ECM}}}_{n}},{\tilde{\mu}}}\right),$

where $G_{H_{n}}=G_{H}(\phi_{P_{n}},\phi_{H_{n}},\phi_{N_{n}})$ is given by
Eq. 2.3 and $W_{H}$ is the cylindrical Wiener process.

  * •

Necrotic.

(4.15)
$\displaystyle\mathopen{}\mathclose{{}\left({\frac{\phi_{N_{n+1}}-\phi_{N_{n}}}{\Delta
t}},{\tilde{\phi}}}\right)=\mathopen{}\mathclose{{}\left({\lambda_{HN}\mathcal{H}(\sigma_{HN}-\phi_{\sigma_{n+1}})\mathopen{}\mathclose{{}\left(\phi_{H_{n}}}\right)^{+}},{\tilde{\phi}}}\right).$

  * •

MDE.

$\displaystyle\mathopen{}\mathclose{{}\left({\frac{{\phi_{{MDE}}}_{n+1}-{\phi_{{MDE}}}_{n}}{\Delta
t}},{\tilde{\phi}}}\right)-\mathopen{}\mathclose{{}\left({{\phi_{{MDE}}}_{n+1}\bm{v}_{n+1}},{\nabla\tilde{\phi}}}\right)+\mathopen{}\mathclose{{}\left({m_{{MDE}}(\bm{\phi}_{n})D_{{MDE}}\nabla{\phi_{{MDE}}}_{n+1}},{\nabla\tilde{\phi}}}\right)$
$\displaystyle\quad+\mathopen{}\mathclose{{}\left({\lambda^{\\!\textup{deg}}_{{MDE}}{\phi_{{MDE}}}_{n+1}},{\tilde{\phi}}}\right)+\mathopen{}\mathclose{{}\left({\lambda^{\\!\textup{pro}}_{{MDE}}(\phi_{P_{n}}+\phi_{H_{n}}){\phi_{{ECM}}}_{n}\frac{\sigma_{H\\!P}}{\sigma_{H\\!P}+\phi_{\sigma_{n+1}}}{\phi_{{MDE}}}_{n+1}},{\tilde{\phi}}}\right)$
$\displaystyle=\mathopen{}\mathclose{{}\left({\lambda^{\\!\textup{pro}}_{{MDE}}(\phi_{P_{n}}+\phi_{H_{n}}){\phi_{{ECM}}}_{n}\frac{\sigma_{H\\!P}}{\sigma_{H\\!P}+\phi_{\sigma_{n+1}}}},{\tilde{\phi}}}\right)$
(4.16)
$\displaystyle\quad-\mathopen{}\mathclose{{}\left({\lambda^{\\!\textup{deg}}_{ECM}{\phi_{{ECM}}}_{n}{\phi_{{MDE}}}_{n}},{\tilde{\phi}}}\right).$

  * •

ECM.

(4.17)
$\displaystyle\mathopen{}\mathclose{{}\left({\frac{{\phi_{{ECM}}}_{n+1}-{\phi_{{ECM}}}_{n}}{\Delta
t}},{\tilde{\phi}}}\right)$
$\displaystyle=\mathopen{}\mathclose{{}\left({\lambda^{\\!\textup{pro}}_{ECM}\phi_{\sigma_{n+1}}\mathcal{H}({\phi_{{ECM}}}_{n}-\phi^{\text{pro}}_{{ECM}})(1-{\phi_{{ECM}}}_{n})},{\tilde{\phi}}}\right)-\mathopen{}\mathclose{{}\left({\lambda^{\\!\textup{deg}}_{ECM}{\phi_{{MDE}}}_{n}.{\phi_{{ECM}}}_{n}},{\tilde{\phi}}}\right)$

  * •

TAF.

$\displaystyle\mathopen{}\mathclose{{}\left({\frac{{\phi_{{TAF}}}_{n+1}-{\phi_{{TAF}}}_{n}}{\Delta
t}},{\tilde{\phi}}}\right)-\mathopen{}\mathclose{{}\left({{\phi_{{TAF}}}_{n+1}\bm{v}_{n+1}},{\nabla\tilde{\phi}}}\right)+\mathopen{}\mathclose{{}\left({D_{TAF}\nabla{\phi_{{TAF}}}_{n+1}},{\nabla\tilde{\phi}}}\right)$
$\displaystyle\quad+\mathopen{}\mathclose{{}\left({\lambda^{\\!\textup{pro}}_{{TAF}}\phi_{H_{n+1}}{\phi_{{TAF}}}_{n+1}\mathcal{H}(\phi_{H_{n+1}}-\phi_{H_{P}})},{\tilde{\phi}}}\right)$
(4.18)
$\displaystyle=\mathopen{}\mathclose{{}\left({\lambda^{\\!\textup{pro}}_{{TAF}}\phi_{H_{n+1}}\mathcal{H}(\phi_{H_{n+1}}-\phi_{H_{P}})},{\tilde{\phi}}}\right)-\mathopen{}\mathclose{{}\left({\lambda^{\\!\textup{deg}}_{{TAF}}\phi_{{TAF}_{n}}},{\tilde{\phi}}}\right).$

The steps followed in solving the coupled system of equations including the
angiogenesis step are summarized in Algorithm 2.

1 Input:Model parameters, $\phi_{\alpha_{0}},v_{0},\Delta t,T,\text{TOL}$ for
$\alpha\in\mathcal{A}:=\\{P,H,N,\sigma,{TAF},{MDE},{ECM}\\}$

2 Output:
$\phi_{\alpha_{n}},\mu_{P_{n}},\mu_{H_{n}},p_{n},\bm{v}_{n},p_{v_{n}},\phi_{{\sigma
v}_{b}}$ for all $n$

3 $n=0$, $t=0$

4 while  _$t\leq T$ _ do

5 $\phi_{\alpha_{n}}=\phi_{\alpha_{n+1}}$ $\forall\alpha\in\mathcal{A}$,
$\mu_{P_{n}}=\mu_{P_{n+1}}$, $\mu_{H_{n}}=\mu_{H_{n+1}}$,
$p_{v_{n}}=p_{v_{n+1}}$, $\phi_{{v}_{n}}=\phi_{{v}_{n+1}}$

6 if  _$\text{apply\\_angiogenesis(n)}==\text{True}$_ then

7 apply angiogensis model described in Section 3

8 update 1D systems if the network is changed

9 end if

10

11 solve coupled linear system $(p_{v_{n+1}},p_{n+1})$ using block Gauss-
Seidel iteration and Subsection 4.1

12 solve coupled linear system $(\phi_{v_{n+1}},\phi_{\sigma_{n+1}})$ using
block Gauss-Seidel iteration and Subsection 4.1

13 solve velocity $\bm{v}_{n+1}$ using Eq. 4.1

14 solve $\mathopen{}\mathclose{{}\left(\phi_{P_{n+1}},\mu_{P_{n+1}}}\right)$
using the semi-implicit scheme in • ‣ Subsection 4.2 and • ‣ Subsection 4.2

15 solve $\mathopen{}\mathclose{{}\left(\phi_{H_{n+1}},\mu_{H_{n+1}}}\right)$
using the semi-implicit scheme in • ‣ Subsection 4.2 and • ‣ Subsection 4.2

16 solve $\phi_{N_{n+1}}$ using the semi-implicit scheme in Eq. 4.15

17 solve ${\phi_{{MDE}}}_{n+1}$ using the semi-implicit scheme in • ‣
Subsection 4.2

18 solve ${\phi_{{ECM}}}_{n+1}$ using the semi-implicit scheme in Eq. 4.17

19 solve ${\phi_{{TAF}}}_{n+1}$ using the semi-implicit scheme in • ‣
Subsection 4.2

20 $n\mapsto n+1$, $t\mapsto t+\Delta t$

21 end while

Algorithm 2 The 3D-1D tumor growth model solver with angiogenesis step

######  Remark 1.

If we ignore the advection and reaction terms of the given system and set
$\chi_{c}=0$ we can show that our algorithm is unconditionally gradient
stable.

This is due to the fact that if we freeze the field $\phi_{H_{n}}$, in both
the convex and concave part of our double-well potential and solve Eq. • ‣ 4.2
for $\phi_{P_{n+1}}$, then due to the given convex-concave splitting we get

$\mathcal{E}(\phi_{\sigma_{n+1}},\phi_{P_{n+1}},\phi_{H_{n}},\phi_{N_{n}},\phi_{MDE_{n}},\phi_{ECM_{n}},\phi_{TAF_{n}})-\mathcal{E}(\phi_{\sigma_{n+1}},\phi_{P_{n}},\phi_{H_{n}},\phi_{N_{n}},\phi_{MDE_{n}},\phi_{ECM_{n}},\phi_{TAF_{n}})\leq
0.$

Similarly, if we now freeze $\phi_{P_{n+1}}$ in both parts of the potential
and solve • ‣ 4.2 for $\phi_{H_{n+1}}$ we get from the unconditional gradient
stability of the sub scheme that

$\mathcal{E}(\phi_{\sigma_{n+1}},\phi_{P_{n+1}},\phi_{H_{n+1}},\phi_{N_{n}},\phi_{MDE_{n}},\phi_{ECM_{n}},\phi_{TAF_{n}})-\mathcal{E}(\phi_{\sigma_{n+1}},\phi_{P_{n+1}},\phi_{H_{n}},\phi_{N_{n}},\phi_{MDE_{n}},\phi_{ECM_{n}},\phi_{TAF_{n}})\leq
0.$

This observation extends to arbitrarily large systems of Cahn-Hilliard
equations and since Eqs. (• ‣ 4.2), (• ‣ 4.2), (• ‣ 4.2) can be considered as
very simple Cahn-Hilliard equations it also extends to them. Finally we note
that $\phi_{N_{n+1}}=\phi_{N_{n}}$ and $\phi_{ECM_{n+1}}=\phi_{ECM_{n}}$ holds
trivially without source terms. Hence, using a telescope sum over all the
energy decrements due to solving Eqs. (• ‣ 4.2), (• ‣ 4.2), (• ‣ 4.2), (4.15),
(• ‣ 4.2), (4.17) and (• ‣ 4.2) we get

$\mathcal{E}(\mathbf{\bm{\phi}}_{n+1})-\mathcal{E}(\mathbf{\bm{\phi}}_{n})\leq
0$

independent of our time-step size, which provides a strong motivation for the
stability of our algorithm.

With the stochastic terms we can generate tumor mass and hence $\mathcal{E}$
does not have to decrease. And even though the reaction terms all add up to
zero, this does not necessarily mean that $\mathcal{E}$ has to decrease, since
they are not part of our gradient-flow. For arbitrary initial states we
therefore cannot expect that $\frac{d}{dt}\mathcal{E}\leq 0$ holds even for
our continuous system.

####  4.2.1. Stochastic component of the system

Generally, the cylindrical Wiener processes $W_{\alpha}$,
$\alpha\in\\{P,H\\}$, on $L^{2}(\Omega)$ with $\Omega=(0,2)^{3}$ can be
written as

$W_{\alpha}(t)(\bm{x})=\sum_{i,j,k=1}^{\infty}\eta^{\alpha}_{ijk}(t)\underbrace{\cos(i\pi
x_{1}/L)\cos(j\pi x_{2}/L)\cos(k\pi x_{3}/L)}_{=:e_{ijk}},$

where $\bm{x}=(x_{1},x_{2},x_{3})$, $L$ is the edge length of the cubed domain
$\Omega$, $\\{e_{ijk}\\}$ form the orthonormal basis of $L^{2}(\Omega)$, and
$\\{\eta^{\alpha}_{ijk}\\}_{i,j,k\in\mathbb{N}}$ is a family of real-valued,
independent, and identically distributed Brownian motions. Following [9, 4],
we approximate the term involving the Wiener process in the fully discretized
system as follows

(4.19) $\frac{1}{\Delta
t}\mathopen{}\mathclose{{}\left(\int_{t_{n}}^{t_{n+1}}\textup{d}W_{\alpha}(t),\xi}\right)_{L^{2}}\approx\frac{1}{\Delta
t}\sum_{\begin{subarray}{c}i,j,k,\\\
i+j+k<I_{\alpha}\end{subarray}}\eta^{\alpha}_{ijk}(e_{ijk},\xi)_{L^{2}},$

where $\xi\in V_{h}$ is a test function,
$\eta^{\alpha}_{ijk}\sim\mathcal{N}(0,\Delta t)$ are independent Gaussians,
and $I_{\alpha}$ controls the number of basis functions.

##  5\. Numerical simulations

In this section, we apply the models described in Sections 2 and 3 and use the
numerical discretization steps discussed in Section 4. We consider examples
that showcase the effects of angiogenesis on the tumor growth. For this
purpose, the model parameters and the basic setting for our simulations are
introduced in Table 2. 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
without considering angiogenesis, i.e., the growth algorithm from Section 3 is
not applied. Afterwards, we repeat the same simulation including the
angiogenesis effects and study the differences between the corresponding
simulations results.

Figure 8.  Initial setting with boundary conditions for a first numerical
experiment. Pressure at top plane ($z=0$) and bottom ($z=2$) ends of the
artery are $3000$ and $2000$ respectively. Similarly, the pressure at top and
bottom ends of the vein are fixed at $1100$ and $1600$, respectively. The high
pressure end of the artery is the inlet end, and there we assign the nutrients
value to $1$. The high pressure end of the vein is also the inlet, and here we
assign the nutrients value to $0$. At the remaining ends, we apply the
upwinding scheme for solving the nutrient equation.

We then consider a scenario consisting of a tumor core surrounded by a small
capillary network. We obtain the network from
source111https://physiology.arizona.edu/sites/default/files/brain99.txt.
First, we rescale the network so that it fits into the domain
$\Omega=(0,2)^{3}$. The vessel radii are rescaled such that the maximal and
minimal vessel radius is given non-dimensionally by $0.05606$ and $0.025$,
respectively. In all of the simulations, we consider the double-well potential
of the form:

$\Psi=C_{\Psi_{T}}\phi_{T}^{2}(1-\phi_{T})^{2},$

where $C_{\Psi_{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 the white noise
terms, using 10 samples for the case without angiogenesis and 50 samples for
the case with angiogenesis. We use the samples to report statistics of
quantity of interests such as total tumor volume, vessel volume density, etc.

###  5.1. Setup and model parameters for the two vessel problem

Table 1. List of parameters and their values for the numerical simulations. Unlisted parameters are set to zero. $\phi_{\alpha}^{\omega}$, $\omega_{\alpha}$, and $I_{\alpha}$ are parameters related to Wiener process, see Eq. 2.3 and Eq. 4.19. Parameter | Value | Parameter | Value | Parameter | Value
---|---|---|---|---|---
$\lambda^{\\!\textup{pro}}_{P}$ | 5 | $\lambda^{\\!\textup{pro}}_{H}$ | 0.5 | $\lambda^{\\!\textup{deg}}_{P},\lambda^{\\!\textup{deg}}_{H}$ | 0.005
$\lambda^{\\!\textup{pro}}_{{ECM}}$ | 0.01 | $\lambda^{\\!\textup{pro}}_{{MDE}},\lambda^{\\!\textup{deg}}_{{MDE}}$ | 1 | $\lambda^{\\!\textup{deg}}_{{ECM}}$ | 5
$\lambda_{P\\!H}$ | 1 | $\lambda_{H\\!P}$ | 1 | $\lambda_{H\\!N}$ | 1
$\sigma_{P\\!H}$ | 0.55 | $\sigma_{H\\!P}$ | 0.65 | $\sigma_{H\\!N}$ | 0.44
$M_{P}$ | 50 | $M_{H}$ | 25 | $C_{\Psi_{T}}$ | 0.03
$\varepsilon_{P}$ | 0.005 | $\varepsilon_{H}$ | 0.005 | $\lambda^{\\!\textup{pro}}_{{TAF}}$ | 10
$D_{{TAF}}$ | 0.5 | $M_{{TAF}}$ | 1 | $L_{p}$ | $10^{-7}$
$D_{\sigma}$ | 3 | $M_{\sigma}$ | 1 | $K$ | $10^{-9}$
$D_{{MDE}}$ | 0.5 | $M_{{MDE}}$ | 1 | $L_{\sigma}$ | $4.5$
$D_{v}$ | $0.1$ | $\mu_{\text{bl}}$ | 1 | $r_{\sigma}$ | $0.95$
$I_{\alpha}$, $\alpha\in\\{P,H\\}$ | $17$ |  $\phi_{\alpha}^{\omega}$, $\alpha\in\\{P,H,T\\}$ | $0.1$ |  $\omega_{\alpha}$, $\alpha\in\\{P,H\\}$ | $0.0025$
Table 2. List of parameters and their values for the growth algorithm and numerical discretization Parameter | Value | Function
---|---|---
$Th_{{TAF}}$ | $7.5\cdot 10^{-3}$ | Threshold for the TAF concentration (sprouting)
$\mu_{r}$ | $1.0$ | Mean value for the log-normal distribution (ratio radius/vessel length)
$\sigma_{r}$ | $0.2$ | Standard dev. for the log-normal distribution (ratio radius/vessel length)
$\lambda_{g}$ | $1.0$ | Regularization parameter to avoid bendings and sharp corners
$\gamma$ | $3.0$ | Murray parameter determining the radii at a bifurcation
$R_{\min}$ | $9.0\cdot 10^{-3}$ | Minimal vessel radius
$l_{\min}$ | $0.13$ | Minimal vessel length for which sprouting is activated
$R_{\max}$ | $0.035$ | Maximal vessel radius
$R_{T}$ | $0.05$ | Threshold for the radius to distinguish between arteries and veins for $t=0$
${\zeta}$ | $1.05$ | Sprouting parameter
$\text{dist}_{\text{link}}$ | $0.08$ | Maximal distance at which a terminal vessel is linked to the network
$\tau_{\text{ref}}$ | $0.02$ | Lower bound for the wall shear stress
$k_{{\text{WSS}}}$ | $0.4$ | Proportional constant (wall shear stress)
$k_{s}$ | $0.14$ | Shrinking parameter
$\Delta t$ | $0.0025$ | Time step size
$h_{3D}$ | $0.0364$ | Mesh size of the 3D grid
$h_{1D}$ | $0.25$ | Mesh size of the initial 1D grid
$\Delta t_{net}$ | $2\Delta t$ | Angiogenesis (network update) time interval

As a computational domain $\Omega$, we choose a cube,
$\Omega=\mathopen{}\mathclose{{}\left(0,2}\right)^{3}$. Within $\Omega$ two
different vessels are introduced: an artery and a vein; see Figure 8. The
radius of the vein $R_{v}$ is given by $R_{v}=0.0625$, and the radius of the
artery $R_{a}$ is set to $R_{a}=0.047$. The centerlines of both vessels are
given by straight lines. In case of the artery, the centerline starts at
$\mathopen{}\mathclose{{}\left(0.1875,0.1875,0}\right)$ and ends at
$\mathopen{}\mathclose{{}\left(0.1875,0.1875,2}\right)$, whereas the vein
starts at $\mathopen{}\mathclose{{}\left(1.8125,1.8125,0}\right)$ and ends at
$\mathopen{}\mathclose{{}\left(1.8125,1.8125,2}\right)$.

At the boundaries of the vessels, we choose Dirichlet boundaries for the
pressure, see also Figure 8. These boundary conditions imply that the artery
provides nutrients for the tissue block $\Omega$, while the vein will take up
nutrients. For the nutrients in the blood vessels mixed boundary conditions
are considered, as depicted in Figure 8.

As initial conditions for $\phi_{v}$, we choose $\phi_{v}=1$ in the artery and
$\phi_{v}=0$ in the vein. The initial value for the nutrient variable
$\phi_{\sigma}$ in the tissue matrix is given by $\phi_{\sigma}=0.5$. In order
to define the initial conditions for the tumor, we consider a ball $B_{T}$ of
radius $r_{c}=0.3$ around the center
$\bm{x}_{c}=\mathopen{}\mathclose{{}\left(1.0,0.8,1.0}\right)$. Within
$B_{T}$, the total tumor volume fraction $\phi_{T}$ smoothly decays from $1$
in the center to $0$ on the boundary of the ball:

(5.1) $\displaystyle\phi_{T}(\bm{x},t=0)$
$\displaystyle=\begin{cases}\begin{aligned}
&\exp\mathopen{}\mathclose{{}\left(1-\frac{1}{1-(|\bm{x}-\bm{x}_{c}|/r_{c})^{4}}}\right),&&\text{if
}|\bm{x}-\bm{x}_{c}|<r_{c},\\\ &0,&&\text{otherwise}.\end{aligned}\end{cases}$

Thereby, the necrotic and hypoxic volume fractions, $\phi_{N}$ and $\phi_{H}$,
are set initially to zero. In the other parts of the domain, all the volume
fractions for the tumor species are set to $0$ at $t=0$. In Table 1, the
parameters for the model equations in Section 2 are listed and Table 2
contains the parameters for the growth algorithm described in Section 3 as
well as the discretization parameters. In particular, the parameters for the
stochastic distributions are listed, which determine the radii and vessel
lengths, the probability of bifurcations, and the sprouting probability of new
vessels.

###  5.2. Robustness of the 3D-1D solver

To ascertain the accuracy and robustness of the proposed solver, we performed
several studies where we changed mesh size and time steps and found that the
solver is robust and the size of time step and mesh size employed in the
studies in sections below balance the computational cost and numerical
accuracy pretty well. To strengthen the claims, we consider a two-vessel setup
described above without the stochasticity and network growth. We run the
simulations using four different time steps $\Delta t_{i}=0.01/2^{i-1}$,
$i=1,2,3,4$, and compute the rate of convergence of quantity of interests such
as $L^{2}$ or $L^{1}$ norms of tumor species and nutrients. In Figure 9, we
plot the $L^{2}$ norm of tumor species and nutrients for different time steps.
We see that the difference between the curves for different $\Delta t$ is very
small. Let $Q_{i}(t)$ denote the quantity of interest ($L^{2}$ norm of
species) at time $t$ for $\Delta t_{i}$. We can approximately compute the rate
of convergence of $Q$ using the formula:

$\displaystyle
r(t)=\frac{\log(|Q_{1}(t)-Q_{4}(t)|)-\log(|Q_{2}(t)-Q_{4}(t)|)}{\log(\Delta
t_{1})-\log(\Delta t_{2})}.$

For $Q(t)=\mathopen{}\mathclose{{}\left\|\phi_{T}(t)}\right\|_{L^{2}}$, we
found $r(1)=0.894,r(2)=1.03,r(3)=1.025,r(4)=0.531,r(5)=1.692$.

We also remark that the proposed solver, see Algorithm 2, does not involve
nonlinear iterations to compute the
$\phi_{P},\phi_{H},\phi_{N},\phi_{{MDE}},\phi_{{ECM}}$ solutions at current
time step. We compared the results of current solver and the solver involving
nonlinear iterations and observed that the solver with nonlinear iterations
still required us to consider small time steps. Also, the error in solution
from two solvers decreases with mesh refinement and smaller time steps. These
observations motivated us to use the proposed solver for all numerical tests
in the sections below.

Figure 9. Plot of the $L^{2}$ norm of various species using four different
time steps.

###  5.3. Tumor growth without angiogenesis

Figure 10.  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. Image taken from [55], with permission from Elsevier. Top
right: $\phi_{T}$ presented in a plane at $z=1$ perpendicular to the $z$-axis.
As seen in the medical experiments the formation of satellites and the
accumulation of tumor cells at the nutrient-rich artery are reproduced in the
simulations. Bottom left: Distribution of necrotic cells ($\phi_{N}$). It can
be seen that the main tumor consists of a large necrotic kernel. Bottom right:
Distribution of nutrients ($\phi_{\sigma}$).

The simulation results for tumor growth without angiogenesis are summarized in
Figure 10. For $t=8$, the tumor cell distribution within the plane
perpendicular to the $z$-axis at $z=1.0$ is shown. In three subfigures, the
volume fraction variables $\phi_{T}=\phi_{P}+\phi_{H}+\phi_{N}$, $\phi_{N}$,
as well as the nutrients $\phi_{\sigma}$ are presented. It can be observed
that the primary tumor is enlarged and small satellites are formed in the
vicinity of the main tumor. The distribution of the necrotic cells indicates
that the main tumor consists mostly of necrotic cells, while the hypoxic and
proliferative cells gather around the nutrient-rich blood vessels. This means
that the tumor cells can migrate against the flow from the artery towards the
vein. Apparently, the chemical potential caused by the nutrient source
dominates the interstitial flow.

These observations are consistent with simulation results and measurements
discussed, e.g., in [19, 35, 55]. In [35, 55] a tumor is introduced into a
mouse. At the same time, anti-VEGF agents were injected into the mouse, such
that the sprouting of new vessels growing towards the tumor is prevented. This
process leads to the formation of satellites located in the vicinity of the
primary tumor as well as the accumulation of tumor cells at nutrient-rich
vessels and cells. Furthermore, the primary tumor stops growing and forms a
large necrotic core.

Figure 11.  Quantities of interests (QoIs) related to tumor species over a
time interval $[0,8]$ for the two-vessels setting. For the case without
angiogenesis, the mean QoIs computed using 10 samples is shown. For the
angiogenesis case, we compute the mean and standard deviation from 50 samples.
The solid line shows the mean QoI as a function of time. The thick layer
around the solid line corresponds to interval
$(\mu_{\alpha}(t)-\sigma_{\alpha}(t),\mu_{\alpha}(t)+\sigma_{\alpha}(t))\cap[0,\infty)$,
for $t\in[0,T]$, where $\mu_{\alpha}(t),\sigma_{\alpha}(t)$ are the mean and
standard deviations of QoI
$\alpha\in\\{\mathopen{}\mathclose{{}\left\|\phi_{T}}\right\|_{L^{1}},\mathopen{}\mathclose{{}\left\|\phi_{P}}\right\|_{L^{1}},\mathopen{}\mathclose{{}\left\|\phi_{H}}\right\|_{L^{1}},\mathopen{}\mathclose{{}\left\|\phi_{N}}\right\|_{L^{1}}\\}$
at time $t$. The variations in the QoIs for the non-angiogenesis case are very
small. The mean of the total tumor volume fraction
$\mathopen{}\mathclose{{}\left\|\phi_{T}}\right\|_{L^{1}}$ at the final time
for the angiogenesis case is about 1.7 times that of the non-angiogenesis
case.

In Figure 11, the $L^{1}$-norms of the tumor species over time are presented
for the case when angiogenesis was inactive and when it was active. While the
profiles for different species in two cases look similar, the total tumor is
about 70$\%$ times higher when angiogenesis is active. Diffusivity
$D_{\sigma}=3$ is large enough, and therefore, the nutrients originating from
the nutrient-rich vessels diffuse quickly throughout the domain.

In summary, one can conclude that without angiogenesis a tumor can grow to a
certain extent, before the primary tumor starts to collapse, i.e., a large
necrotic core is formed. However, this does not mean that the tumor cells are
entirely removed from the healthy tissue. If there is a source of nutrients
such as an artery close by, transporting nutrient rich blood, a portion of
tumor cells can survive by migrating towards the neighboring nutrient source.

###  5.4. Tumor growth with angiogenesis

As in the previous subsection, we compute the $L^{1}$-norms of the tumor
species $\phi_{T}$ at time $t=8$. However, since several stochastic processes
are involved in the network growth and also Wiener processes appear in the
proliferative and hypoxic cell mass balances, several data sets have to be
created in order to rule out statistical variations. In this context, the
issue arises as to how many data sets are needed to obtain a representative
value. In order to investigate this, we compute for every sample $i$, the
$L^{1}$-norm of the tumor species, denoted by
$\phi_{\alpha_{L_{1},i}},\;\alpha\in\mathopen{}\mathclose{{}\left\\{T,P,H,N}\right\\}$.
Additionally, the volume of the blood vessel network $V_{i}$ is computed. For
each data set with $i$ samples, we compute the mean values:

$\textup{mean}_{i}\mathopen{}\mathclose{{}\left(V}\right)=\frac{1}{i}\sum_{j=1}^{i}V_{j},\;\qquad\textup{mean}\mathopen{}\mathclose{{}\left(\phi_{\alpha_{L_{1},i}}}\right)=\frac{1}{i}\sum_{j=1}^{i}\phi_{\alpha_{L_{1},j}},\;\qquad\alpha\in\mathopen{}\mathclose{{}\left\\{T,P,H,N}\right\\}.$

In Figure 12, the mean values
$\textup{mean}_{i}\mathopen{}\mathclose{{}\left(V}\right)$ and
$\textup{mean}\mathopen{}\mathclose{{}\left(\phi_{\alpha_{L_{1},i}}}\right)$,
$\alpha\in\\{T,P,H,N\\}$, are shown. From the plots we see that the mean of
$||\phi_{T}||_{L^{1}}$ stabilizes after about 25 samples. For the vessel
volume, fluctuations in the sample means reduce with increasing sample and get
small after 30 samples. While the results in Figure 12 show that mean of the
QoIs stabilizes with increasing sample and converge to some number, the
trajectory in the figure could change with change in sample values. For
example, if we shuffle the samples and recompute the quantities in Figure 12,
various curves may look different.

Figure 12.  The mean values for the $L^{1}$-norm of the tumor cell volume
fractions $\phi_{T},\phi_{P},\phi_{H},\phi_{N}$ and the volume of the blood
vessel network at time $t=8$ from increasing number of samples. Results
correspond to the two-vessel setting. The mean of the total tumor volume
fraction appears to be stable after about 28 samples. The mean of the vessel
volume shows smaller fluctuations as the number of samples in the data set
grows.

As mentioned earlier, Figure 11 presents the $L^{1}$-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 fraction varies from
sample to sample, as expected. Both the hypoxic and the tumor cells show an
exponential growth after $t\approx 4$; see Figure 11. After decreasing until
$t\approx 3$, the proliferative cell mass grows from $t\approx 4$ onward. The
result is that in the case of angiogenesis, the overall nutrient concentration
is higher compared to the case without angiogenesis, while the spatial
variation of the nutrient is the same in the two; and hence the growth of the
tumor in the two cases are similar except that the tumor grows more rapidly in
the case of angiogenesis. We will see in our second example, where
$D_{\sigma}=0.05$ is much smaller, that the nutrient concentration is higher
near the nutrient rich vessels and tumor growth is more concentrated near
these regions, see Figure 18.

Figure 13.  Growth of the tumor and the network at four different time points
$t\in\mathopen{}\mathclose{{}\left\\{0.24,0.64,3.20,5.60}\right\\}$ and one
specific sample. We show contour $\phi_{T}=0.9$. Figure 14.  Tumor cell
distributions in the plane $z=1$ with angiogenesis for one sample. Top left:
$\phi_{T}$. Top right: $\phi_{N}$. Bottom left: $\phi_{H}$. Bottom right:
$\phi_{P}$.

In Figure 13, we show the evolving network together with the contour plot
$\phi_{T}=0.9$ of the total tumor species. At time $t=0.24$ (top-left figure),
we observe that new vessels originate from the artery and move towards the
hypoxic core; the directions of these vessels being based on the gradient of
TAF with random perturbations. At $t=0.64$ (top-right), we see a large number
of new vessels formed as predicted by the angiogenesis algorithm. However, at
time $t=3.2$ (bottom-left), vessels adapt and due to lower flow rates in some
newly created vessels, some vessels are gradually removed, and thus the number
of vessels decreases. Comparing $t=3.2$ and $t=5.6$ (bottom-right), we see
that the network has stabilized and little has changed in this time window.
From Figure 13, we can also summarize that the tumor growth is directed
towards the nutrient-rich vessels.

Next, we plot the tumor species at the $z=1$ plane along with the nutrient
distribution in the vessels in Figure 14. The plot corresponds to time $t=8$.
The plots corresponding to the necrotic species show that the necrotic
concentration is typically higher away from the nutrient-rich vessels. From
the hypoxic plot, we see that it is higher near these vessels, and this is
explained by the fact that as soon as the proliferation of new tumor cells
takes place, due to nutrient concentration below the proliferative-to-hypoxic
transition threshold, these newly added proliferative tumor cells convert to
hypoxic cells. Further transition to necrotic cells would take place if the
nutrients are even below the hypoxic-to-necrotic transition threshold. This is
also consistent with the increase concentration of the proliferative cells
near the outer tumor-healthy cell phase interface.

###  5.5. Sensitivity of the growth parameters

Figure 15.  Results of a parametric study in which certain growth parameters
$\gamma$ and $k_{{\text{WSS}}}$ are varied to measure their impact on the
total network length and volume. For these studies, we considered a coarse
mesh for the 3D domain with $16^{3}$ uniform cells and time step $\Delta
t=0.01$. The network update time step was fixed to $\Delta t_{net}=2\Delta t$.
Figure 16.  Network structure for different values of $\gamma$ at time
$t=2.4$. Left: $\gamma=2$. Middle: $\gamma=3$. Right: $\gamma=4$.  Figure 17.
Network structure for different values of $k_{{\text{WSS}}}$ at time $t=6$.
Left: $k_{WSS}=4$. Middle: $k_{WSS}=0.4$. Right: $k_{WSS}=0.001$.

In Figure 15, we present the results of a parametric study designed to test
the robustness of the vascular network model to changes in the values of the
parameters $\gamma$ and $k_{WSS}$. It is observed that changes in these
parameters can produce significant changes in the network structure for given
values of the other model parameters.

The parameter $\gamma$, for example, appears in Murray’s law (3.4), and
controls the radii of network branches, with increase in $\gamma$ leading to
larger radii of bifurcating vessels. Such larger radii vessels have a higher
probability for connecting with neighboring vessels so as to increase the flow
and to continue to evolve; for high $\gamma$, the networks are more dense and
the total network length is higher (see Figure 16, right). Conversely, small
values of $\gamma$ promote thin network segments with lower probability for
connecting to neighboring vessels, see Figure 16, left.

The change in vessel radius due to the vessel wall shear stress is
proportional to the constant $k_{{\text{WSS}}}$ and stimulus
$S_{{\text{WSS}},i}$, see (3.8) and (3.9). Further, the vessels shrink
naturally and this effect is controlled by constant $k_{s}$ (higher $k_{s}$
means radius decay is higher). In our study, we varied the values of parameter
$k_{{\text{WSS}}}$ and found that, for a large $k_{{\text{WSS}}}$, radii of
sprouting vessels decrease, and new sprouts are removed in their early stage
of growth before they could join the nearby vessels, see Figure 17 left. As a
result of new sprouts getting removed in early stages, the total network
length and the vessel volume stay constant with respect to time with constant
values very close to the initial values. For a very small $k_{WSS}$ with
$k_{s}$ being fixed, the radii during the early phase of simulations do not
change much. But in the later phases of the simulation, the radii begin to
decay and even with large flow rate (which means large wall shear stress);
their decay is unavoidable as the term $k_{{\text{WSS}}}S_{{\text{WSS}},i}$ is
small as $k_{{\text{WSS}}}$ is small and can not counter the effects of
$k_{s}$. In summary, in the long run the radii of vessels decrease with time,
see Figure 17 right.We also observed that for values of $k_{{\text{WSS}}}$
within certain bounds, its impact on the network morphology is low. However,
when $k_{{\text{WSS}}}$ is outside the bound, some care is required so that
vessel radii do not tend to zero with time.

###  5.6. Angiogenesis and tumor growth for the capillary network

Returning to (5.1), let us consider a smooth spherical tumor core with center
at $\bm{x}_{c}=(1.3,0.9,0.7)$ and radius $r_{c}=0.3$ in the domain
$\Omega=(0,2)^{3}$. The initial blood vessel network and boundary conditions
for pressure and nutrient on the network is described in Figure 20.

In the simulation, the vessels are uniformly refined two times. We fix
$\phi_{a}=0$ for $a\in\\{H,N,{TAF}\\}$ and $\phi_{\sigma}=0.5$ at $t=0$. The
domain is discretized uniformly with a mesh size $h_{3D}=0.0364$ and the time
step size is $\Delta t=0.0025$. We identify four inlet ends (see Figure 20) at
which the unit nutrient $\phi_{v}=\phi_{v_{in}}=1$ and pressure
$p_{v}=p_{in}=8000.0$ is prescribed as the Dirichlet boundary condition. At
the remaining boundaries, we prescribe the pressure $p_{v}=p_{out}=1000.0$ and
apply an upwinding scheme for the nutrients.

At $t=0$, $p_{v}$ and $\phi_{v}$ at internal network nodes are set to zero.
Furthermore, we set $L_{\sigma}=0.5$, $D_{\sigma}=0.05$, $D_{TAF}=0.1$,
$Th_{{TAF}}=0.0075$, $\mu_{r}=1.5$, $k_{{\text{WSS}}}=0.45$, and $\Delta
t_{net}=10\Delta t$. All the other parameter values remain unchanged, see
Table 1 and Table 2.

Figure 18.  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 compute the mean and standard deviation using 50 and 10
samples for the case with and without angiogenesis, respectively. We refer to
Figure 11 for more details on the plots. As in the case of two-vessels
setting, the variation in the QoIs are much smaller for the non-angiogenesis
case. The mean of the total tumor volume fraction
$\mathopen{}\mathclose{{}\left\|\phi_{T}}\right\|_{L^{1}}$ for the
angiogenesis case is about 1.62 times that of the non-angiogenesis case.
Figure 19.  The mean values for the $L^{1}$-norm of the tumor cell volume
fractions $\phi_{T},\phi_{P},\phi_{H},\phi_{N}$ and the volume of the blood
vessel network at time $t=8$. Results correspond to the capillary network
setting. The mean of the total tumor volume fraction stabilizes with small
samples. This agrees with Figure 18 that shows that the variations in
$L^{1}$-norm QoIs are overall smaller. The mean of the vessel volume shows
some change when samples are small and stabilizes as the size of data set
grows.

We first compare the tumor volume with and without angiogenesis; see Figure
18. The results are similar to 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\approx 3.5$ with
angiogenesis as compared to $t\approx 4.75$ without angiogenesis. Production
of necrotic cells is higher in the non-angiogenesis case until time $t\approx
5$. Compared to Figure 11 for the two-vessel setting, the variations in the
tumor species related QoIs are much smaller in Figure 18. This may be due to
the fact that diffusivity of the nutrients in the latter case is much smaller,
and that $L_{\sigma}$ 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 19; results show that mean of tumor species
related QoIs is stable and can be computed accurately using fewer samples.
This relates to the fact that we see smaller variation in the $L^{1}$-norm of
the tumor species $\phi_{T}$ in Figure 18. The mean of vessel volume shows
some variations for smaller data sets and the variations get smaller later on;
still the variations are very small and contained in range $[0.117,0.121]$.

In Figures 20 and 22, some results for the capillary network are summarized,
where in Figure 20 the growing network is shown and in the Figure 21 the
vessel network as a result of angiogenesis is shown at the final simulation
time $t=8$ for three samples. In Figure 22, 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 grow faster after it is vascularized.
Apparently, the tumor cells tend to migrate towards the nutrient rich part of
the computational domain despite the fact that they have to move against the
direction of flow which is induced by the pressure gradient. Not surprisingly,
the volume fraction of the necrotic cells is larger in the part that is facing
away from the nutrient rich part and related to the whole tumor it remains
relatively small.

It is interesting to observe that, as in the two-vessel case, the contour plot
of $\phi_{T}$ for $\phi_{T}=0.9$ exhibits a secondary structure, while in the
simulation without angiogenesis, this effect cannot be seen. Moreover, as in
the two-vessel experiment, the tumor contains a large necrotic kernel if there
is no angiogenesis, indicating that the tumor has almost died. This simulation
portrays once again that angiogenesis can play a crucial role in the evolution
of tumor growth.

Figure 20.  Top left: Spherical tumor core
$\mathopen{}\mathclose{{}\left(\text{Contour plot for }\phi_{T}=0.9}\right)$
at $\bm{x}_{c}=(1.3,0.9,0.7)$ with radius $r_{c}=0.3$ surrounded by a network
of vessels. We identify four inlet ends (red cross-sections) at which the unit
nutrient $\phi_{v}=\phi_{v_{in}}=1.0$ and pressure $p_{v}=p_{in}=8000$ is
prescribed as a Dirichlet boundary condition. Top right: Formation of first
sprouts at $t=1.20$ growing towards the tumor core. Bottom left: Around
$t=3.04$ a 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 near the nutrient vessels. Figure 21.  Plot of
vessel network at $t=8$ from three samples for the capillary network setting.
Figure 22.  Distribution of tumor cells for $t=5.12$. The simulation results
for tumor growth supported by angiogenesis are shown on the left hand side,
while the results for tumor growth without angiogenesis are presented on the
right. The necrotic $\mathopen{}\mathclose{{}\left(\phi_{N}}\right)$ and
hypoxic $\mathopen{}\mathclose{{}\left(\phi_{H}}\right)$ volume fractions are
visualized in the $z$-plane at $z=0.7$. For $\phi_{T}$, in both cases a
contour plot for $\phi_{T}=0.9$ is shown.

##  6\. Summary and outlook

In this work, we presented a stochastic model for tumor growth characterized
by a coupled system of nonlinear PDEs of Cahn-Hilliard type coupled with a
model of an evolving vascular network. A 3D-1D model is developed to simulate
flow and nutrient transport within both the network and porous tissue so as to
depict the phenomena of angiogenesis. In this model, the blood vessel network
is given by a 1D 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 with bifurcation of growing vessels
which is crucial for the simulation of angiogenesis. The angiogenesis process
is simulated by an iterative algorithm starting from a single artery and a
vein or a given network. The blood vessel network supplying the tumor employs
Murray’s law to determine the radii at a bifurcation of network capillaries.
The choice of radii and lengths of new vessels as well as bifurcations are
governed by stochastic algorithms. The direction of growth of the vessels is
determined by the gradient of the local TAF concentration.

We demonstrate that the model is capable of simulating the development of
satellite of tumor concentrations in nutrient rich vessels near necrotic cores
in agreement with some experimental observations. Also, as expected, rapid
growth of solid tumor mass accompanies increased supply of nutrient through
angiogenesis.

We believe that our model can serve as a starting point for important
predictive simulations of cancer therapy; in particular the effect of anti-
angiogenic drugs could be studied using models of these types. However, models
of these kind require experimental data such as MR imaging data that inform
the vasculature in the tissue as well as the parameters in the tumor growth
model and vasculature flow model. In the present work, the adaption of the
vessel radii is related to the wall shear stress. However, other effects 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. Further work on refining and improving
computational algorithms is also needed, such as the development of efficient
distributed solvers for the 3D-1D systems. We hope to address these issues and
other extensions in future work.

## Acknowledgements

The authors gratefully acknowledge the support of the Deutsche
Forschungsgemeinschaft (DFG) through TUM International Graduate School of
Science and Engineering (IGSSE), GSC 81. The work of PKJ and JTO was supported
by the U.S. Department of Energy, Office of Science, Office of Advanced
Scientific Computing Research, Mathematical Multifaceted Integrated Capability
Centers (MMICCS), under Award Number 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, Biomechanics and 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, Ĺ. Baňas, R. Nürnberg, and A. Prohl, Numerical approximation of the stochastic Cahn–Hilliard equation near the sharp interface limit, Numerische Mathematik, 147 (2021), pp. 505–551.
  * [5] N. Bellomo, N. Li, and P. K. Maini, On the foundations of cancer modelling: Selected topics, speculations, and perspectives, 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 with the 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 and Biology: 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.
  * [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ńska, and D. Wrzosek, Mathematical modelling of cancer invasion: The importance 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 urokinase plasminogen 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 mixture model: 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 Modeling Approach, Cambridge University Press, 2010.
  * [14] G. Da Prato and A. Debussche, Stochastic Cahn–Hilliard equation, Nonlinear Analysis: Theory, Methods & Applications, 26 (1996), pp. 241–263.
  * [15] M. Dorraki, A. Fouladzadeh, A. Allison, C. S. Bonder, and D. Abbott, Angiogenic networks in tumors—insights via mathematical modeling, IEEE Access, 8 (2020), pp. 43215–43228.
  * [16] H. M. Eilken and R. H. Adams, Dynamics of endothelial cell behavior in sprouting angiogenesis, Current Opinion in Cell Biology, 22 (2010), pp. 617–625.
  * [17] C. Engwer, C. Stinner, and C. Surulescu, On a structured multiscale model for acid-mediated tumor invasion: The effects of adhesion and proliferation, Mathematical Models and Methods in Applied Sciences, 27 (2017), pp. 1355–1390.
  * [18] D. J. Eyre, Unconditionally gradient stable time marching the Cahn–Hilliard equation, in Materials Research Society Symposium Proceedings, vol. 529, Materials Research Society, 1998, pp. 39–46.
  * [19] H. B. Frieboes, F. Jin, Y.-L. Chuang, S. M. Wise, J. S. Lowengrub, and V. Cristini, Three-dimensional multispecies nonlinear tumor growth – II: Tumor invasion and angiogenesis, Journal of Theoretical Biology, 264 (2010), pp. 1254–1278.
  * [20] S. Frigeri, K. F. Lam, E. Rocca, and G. Schimperna, On a multi-species Cahn–Hilliard–Darcy tumor grwoth model with singular potentials, Communications in Mathematical Sciences, 16 (2018), pp. 821–856.
  * [21] M. Fritz, P. K. Jha, T. Köppl, J. T. Oden, and B. Wohlmuth, Analysis of a new multispecies tumor growth model coupling 3D phase-fields with a 1D vascular network, Nonlinear Analysis: Real World Applications, 61 (2021), p. 103331.
  * [22] M. Fritz, E. Lima, V. Nikolic, J. T. Oden, and B. Wohlmuth, Local and nonlocal phase-field models of tumor growth and invasion due to ECM degradation, Mathematical Models and Methods in Applied Sciences, 29 (2019), pp. 2433–2468.
  * [23] M. Fritz, E. Lima, J. T. Oden, and B. Wohlmuth, On the unsteady Darcy–Forchheimer–Brinkman equation in local and nonlocal tumor growth models, Mathematical Models and Methods in Applied Sciences, 29 (2019), pp. 1691–1731.
  * [24] H. Garcke, K. F. Lam, R. Nürnberg, and E. Sitka, A multiphase Cahn–Hilliard–Darcy model for tumour growth with necrosis, Mathematical Models and Methods in Applied Sciences, 28 (2018), pp. 525–577.
  * [25] A. Gerisch and M. Chaplain, Mathematical modelling of cancer cell invasion of tissue: Local and non-local models and the effect of adhesion, Journal of Theoretical Biology, 250 (2008), pp. 684–704.
  * [26] 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.
  * [27] D. Hanahan and R. A. Weinberg, Hallmarks of cancer: The next generation, Cell, 144 (2011), pp. 646–674.
  * [28] A. Hawkins-Daarud, K. G. van der Zee, and J. T. Oden, Numerical simulation of a thermodynamically consistent four-species tumor growth model, International Journal for Numerical Methods in Biomedical Engineering, 28 (2012), pp. 3–24.
  * [29] 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.
  * [30] E. Hodneland, X. Hu, and J. Nordbotten, Well-posedness, discretization and preconditioners for a class of models for mixed-dimensional problems with high dimensional gap, arXiv preprint arXiv:2006.12273, (2020).
  * [31] L. Holmgren, M. S. O’Reilly, and J. Folkman, Dormancy of micrometastases: balanced proliferation and apoptosis in the presence of angiogenesis suppression, Nature Medicine, 1 (1995), pp. 149–153.
  * [32] 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).
  * [33] T. Köppl, E. Vidotto, and B. Wohlmuth, A 3D-1D coupled blood flow and oxygen transport model to generate microvascular networks, International Journal for Numerical Methods in Biomedical Engineering, e3386 (2020).
  * [34] P. Koumoutsakos, I. Pivkin, and F. Milde, The fluid mechanics of cancer and its therapy, Annual Review of Fluid Mechanics, 45 (2013), pp. 325–355.