Angiogenic Factors produced by Hypoxic Cells are a leading driver of Anastomoses in Sprouting Angiogenesis–a computational study

Angiogenesis - the growth of new blood vessels from a pre-existing vasculature - is key in both physiological processes and on several pathological scenarios such as cancer progression or diabetic retinopathy. For the new vascular networks to be functional, it is required that the growing sprouts merge either with an existing functional mature vessel or with another growing sprout. This process is called anastomosis. We present a systematic 2D and 3D computational study of vessel growth in a tissue to address the capability of angiogenic factor gradients to drive anastomosis formation. We consider that these growth factors are produced only by tissue cells in hypoxia, i.e. until nearby vessels merge and become capable of carrying blood and irrigating their vicinity. We demonstrate that this increased production of angiogenic factors by hypoxic cells is able to promote vessel anastomoses events in both 2D and 3D. The simulations also verify that the morphology of these networks has an increased resilience toward variations in the endothelial cell’s proliferation and chemotactic response. The distribution of tissue cells and the concentration of the growth factors they produce are the major factors in determining the final morphology of the network.


S1. Angiogenesis Dynamics
In the present model, the vessels are described by an order parameter φ(r, t) and the cells in hypoxia by a second order parameter ψ(r, t), following [1,2]. The three phases of the system (the extra cellular matrix (ECM), the vessel cells and the tissue cells) have specific values for the order parameters φ and ψ (see Fig. 1 of the main document). The tissue cells will not alter their shape and the modelling is focused on the evolution of the capillary.
According to the phase-field formalism, we start from a phenomenological free-energy, F [φ, ψ], given by [2][3][4] In this free energy we include the surface tensions of the different interfaces in the system (proportional to ε and to ε h ), and a phenomenological potential which guarantees that the minimum energy is obtained in the three phases of the system: the vessel (with φ = 1 and ψ = −1), the ECM (with φ = −1 and ψ = −1), and the tissue cells (with φ = −1 and ψ = 1). The parameters ε and ε h in (1) are the widths of the interfaces and are proportional to their surface tensions. The last term of the potential V (φ, ψ) guarantees that there is an increase in energy when the vessel overlaps the tissue cells, and γ represents the energy cost for this overlap; we set this parameter high enough in order to avoid appreciable overlap. Due to the large value of γ in our simulations, the system only has interfaces between ECM and the vessels and between the tissue cells and the ECM.
The evolution of the vessel is then obtained by calculating the system's free energy functional derivative with respect to φ, according to [1]: This equation becomes Equation (1) of the main text once the proliferation term is replaced. The angiogenic factor in the model is modelled as a diffusive field T , that is consumed by the endothelial cells (ECs) with the rate α T : where Θ is the Heaviside step function.
A new tip cell agent is introduced at a site where there is a minimum angiogenic factor concentration T c , a minimum angiogenic factor gradient G m , if the site for the candidate is located at a distance larger than 4R c from all other tip cell agents (due to Delta-Notch signalling) and if it is at a minimum distance of R c from the ECM. All the activated tip cells migrate by chemotaxis, and their velocity is proportional to the angiogenic factor gradient: where v is the agent velocity, χ is the chemotactic response, G M is the gradient modulus for which the maximum velocity is attained, and G ≡ |∇T |. An tip cell is not able to enter inside a tissue cell (independently if it is in hypoxia or not): when the distance between the tip cell and a tissue cell is lower than 2R c , then the radial component of the velocity is set to zero and the tip cell moves in the axial direction around the tissue cell.
To link the phase-field approach described in Equation (3) to the agent-based model we set the order parameter value inside the tip cell agent to a constant value φ c that depends on the tip cell velocity |v| and on the local ECs proliferation rate α p (T ) [1]: To improve the computational efficiency, we do not include angiogenic factor consumption inside the tip cells. As the tip cell is very small, the VEGF gradient is not significantly altered in such a small distance, and we opted to simplify the computation by measuring the gradient at a single point: the cell center. Preforming a single measurement of the VEGF gradient saves time (specially in 3D) that would be used in calculating several interpolations to obtain the VEGF concentrations at the cell boundary and averaging over all of the differences of concentrations. Since the amount of VEGF consumed by the tip cells would be very small when compared to the VEGF consumed by all the cells in the network, this choice does not change in a significant way the total amount of VEGF in the system, it is able to obtain the correct gradient, and it is faster. The blood is able to flow through the neo-vessels, once there are anastomoses in the vascular network. In order to calculate the regions irrigated in the tissue, we introduce a simple system to identify the vessels where blood is flowing. We do not model explicitly the oxygen diffusion in the simulation because the oxygen's consumption rate is sufficiently high that it reaches the stationary state much faster than any other relevant process in the simulation [5]). We start by extracting the medial lines of the vascular structure, thus finding the bifurcations sites and the length of every vessel. The extraction method used for this purpose is detailed in [6] and consists in the application of a set of topological masks that determines if a point is a simple point (marked to be removed) or an essential point (part of the final skeleton). After some iterations the skeleton of the blood vessel network is obtained (as seen in Fig. 1

).
A B Figure 1. The column A shows the vessel network before the thinning algorithm from two different viewing angles and the column B shows the skeleton after the application of the algorithm.
Afterwards we calculate the blood flow rate in every vessel using the Hagen-Poiseuille law: where ∆p is the difference in pressure between the beginning and end of every vessel, L the vessel length, R the vessel radius, Q is the flow rate, and K is a constant proportional to the inverse of the blood viscosity. The value of K was set to K = 1, since we are only interested in identifying the vessels with Q = 0, and not the precise value of the flow rate. In the simulation there is a single input node (at the top of the simulation box) and a single exit node (at the bottom of the simulation box).

S2. Model Parameters
The parameters for the model were obtained, when possible, from experimental data available on the literature and are presented in Table 1.

S3. Computational Methods
The simulations of the angiogenesis model were performed on parallel task farm on the computer cluster of the Centre of Physics of the University of Coimbra (see example of simulation in Figure 2). The gradients and laplacians were calculated using the finite differences over a discrete mesh with dimensions: 200 × 200 × 100 for the 3D simulation and 300 × 300 for the 2D, with the lattice parameter a = 1. In the 3D simulations were employed periodic boundary conditions on the radial direction and Neumann boundary conditions on the axial direction (no flux condition). The 2D simulations were performed with no flux conditions along both x and y directions. The initial conditions contain in one main vessel placed at the center (border) of the 3D (2D) simulation box and randomly distributed hypoxic cells at a distance larger than 25 µm of the main vessel (about 80 cells in the 2D simulations and 240 cells in the 3D simulations). The order parameter ψ(r, t) was integrated in time before the start of the simulation in order to determine the diffuse interface between the hypoxic cells and the ECM. The process used to obtain the pressure at every bifurcation (needed to calculate the blood flow rate from equation (7)), leads to a large sparse matrix which requires the application of an efficient method to invert matrices. The method of choice was the PLU decomposition with partial pivoting, implemented in the ScaLAPACK [12].

S4. Variable chemotaxis response
The chemotactic response of the tip cells depends on the amount of VEGF receptor 2 (VEGFR2) molecules at their cytoplasm membrane. However, it is known that the expression of VEGFR2 is lower in hypoxic environments [13,14], leading to a lower chemotactic response. Here we address if this fact alters the conclusions drawn regarding how the numbers of branches and anastomoses depend on the maximum tip cell velocity. With this aim, we run the 2D system for both rules, but with the following altered tip cell velocity: where In our model, when the tissue cells are not irrigated they will produce VEGF, and therefore we will use the concentration of VEGF as a proxy to quantify hypoxia. In this way, the new term χ T (T ) that is included in this equation for the tip cell velocity implies that the chemotactic response decreases with hypoxia (See Figure 3). For simulations with Rules 1 and 2, we have varied the value of χ and quantified the density of anastomoses and branches as well as the branches' diameter (see Figure 6). We observed even more striking results than in Figure 4 of the main text. Namely, Rule 1 leads to a clear increase in the number of branches and anastomoses, but at the same time, to much more stable vascular networks (regarding branches and anastomoses densities as well as vessel diameters). For Rule 2, as the tip cell velocity increases, the number of branches increases until a maximum is reached, and the density of anastomoses follows, qualitatively, the same non-monotonic behaviour. For Rule 1 the number of branches and anastomoses are significantly higher and their relative variations with tip cell velocity are very small. As in Figure 4 of the main text, the vessel diameter decreases with the tip cell velocity, but now it is clear that Rule 1 leads to a smaller variation of the vessel diameter. We conclude that the results described in the main text are robust to the inclusion of variable chemotactic response.

S5. Growth factor degradation
In the model in the main text only the VEGF uptake by the vessel cells was considered and not the VEGF degradation. However, the conclusions drawn are not altered when we add VEGF degradation. The degradation rate that is observed in vivo is of the order of λ deg = 1.0 × 10 −4 s −1 [15]. This degradation rate, together with the VEGF diffusion constant, leads to a typical VEGF diffusion length of 30 µm, slightly larger than the oxygen diffusion length.
With VEGF degradation, equation (2) of the main document is replaced by We simulated the vascular growth as a function of the endothelial cell proliferation rate for Rule 1 and for Rule 2. The corresponding branches and anastomoses densities are plotted in Figure 5.
These results were obtained with the same parameters as in the main text, except for the value of T s which was doubled (when introducing degradation in the model the total supply of VEGF decreased, and we increased T s so that the VEGF concentration at the vessels was still large enough to induce sprouting).
As for the case without VEGF degradation, we observe that Rule 1 leads to vascular morphologies much more independent of EC proliferation and with higher number of anastomoses.

S6. Variable hypoxic cell density
In order to address the question of how the number of vessels depends on the number of cells in hypoxia, we have performed simulations where the chemotactic velocity and the proliferation rate were set as constants, and the density of hypoxic cells was varied. Below we show, for 2D, that the number of branches and the number of anastomoses increase linearly with the number of cells in hypoxia for Rule 1. Again, we observe a stronger trend with Rule 1 for the number of vessels to follow the tissue demand: more cells in hypoxia lead to a higher vessel density. S7. Varying the consumption rate α T We run 2D simulations to predict how the morphologies of 2D vascular networks depend on EC proliferation rate for different values of VEGF consumption rate α T . In particular, we measured the anastomoses density for α T values one third lower and 50% higher than the value used in the main simulations in the manuscript. The results are presented in Figure 7.
A B C Figure 7. Anastomoses density as a function of the proliferation rate α p for three values of α T : The values of α T in panels A, B and C are respectively α T = 2α T,0 /3 = 2.64 × 10 −3 s −1 , α T = α T,0 = 4.0 × 10 −3 s −1 and α T = 3α T,0 /2 = 6.0 × 10 −3 s −1 , where α T,0 is the consumption rate used in the main document. We observe that for the three consumption rates, Rule 1 always leads to more anastomoses and to more resilient networks than with Rule 2.
For the three VEGF consumption rates, we observe that Rule 1 leads to vascular morphologies with a higher number of anastomoses than Rule 2. Moreover, similarly to what is observed in the main document, the anastomoses density in the networks that grow under Rule 1 has a much smaller relative variation with EC proliferation rate than the anastomoses density in the networks that grow under Rule 2. Therefore, the results observed and the conclusions drawn in the main document are robust to variations in the VEGF consumption rate.