A simple agent-based model to simulate 3D tumor-induced angiogenesis considering the evolution of the hypoxic conditions of the cells

In tumors, angiogenesis (conformation of a new vasculature from another primal one) is produced with the releasing of tumor angiogenic factors from hypoxic cells. These angiogenic substances are distributed around the tumor micro-environment by diffusion. When they reach the primal blood vessel bed, the sprouting and branching of a new micro-vascular network is produced. These new capillaries will supply oxygen to cells so that their hypoxic state is overcome. In this work, a new and simple 3D agent-based model to simulate tumor-induced angiogenesis is presented. In this approach, the evolution of the hypoxic conditions in cells along the related conformation of the new micro-vessels is considered. The importance that the relative position of the primal vasculature and tumor structure takes in the final distribution of the new micro-vasculature has also been addressed. The diffusion of angiogenic factors and oxygen has been modelled at the targets by numerical convolution superposition of the analytical solution from the sources. Qualitative and quantitative results show the importance of tip endothelial cells in overcoming hypoxic conditions in cells at early stages of angiogenesis. At final stages, anastomosis plays an important role in the reduction of hypoxia in cells.


Introduction
Whereas the proposed approach may be employed in many types of tumors, the motivation of our work is to understand Glioblastoma Multiforme, where angiogenesis plays a special key role, through an efficient procedure which can progressively incorporate the variables and complex physiological interactions occurring in tumors. Therefore, the parameters in the simulations are, when possible, obtained for Glioblastoma.

Glioblastoma multiforme (GBM)
The nervous system is a complex part of the human body that regulates the efficient and proper functioning of the organs and cells in the body. It is composed of a peripherical and a central nervous system (CNS) consisting of brain, spinal cord and retina [33]. Brain diseases are considered one of most important problems in public health according to the World Health Organization (WHO) [16]. Malignant brain tumors are one of these concerns, not only because of their impact in the function of the nervous system, but also because they are one the most aggressive types of cancer [64]. Glioblastoma multiforme (GBM) is a grade IV neoplasia [55] of the glial cells and constitutes 82% of cases of malignant glioma (25% of all brain tumors) with 17, 000 new cases every year. The 5-year survival rate is only of 4.7% [27].
Nervous tissue consists in an extracellular matrix (ECM), blood vessels, a cerebrospinal fluid and neurons 1 3 and glial cells. There are three types of glial cells: astrocytes (feeding neurons from vessels through their cytoplasm), oligodendrocytes (producing myelin) and microgliocytes (phagocytes for waste and immune system) [33]. GBM usually develops de novo from stem cells (GS-cells) which can be found as subpopulations in GBM [36,53,75]. GS-cells are highly adaptative to acidic stress [46] and hypoxia [8,44] (which both promote conversion in a cancer stem cell phenotype [53]) and proved highly resistant, by activation of DNA damage response, to both radio and chemotherapy in vitro and in vivo [7,53]. Oxygen supply by vessels to GBM cells has a relevant role in its evolution. For example, hyperbaric oxygen promotes GBM proliferation and chemosensitization [90]. Very high vascular endothelial growth factor (VEGF) concentration values, ten times the values found in other brain tumors, and more than 100 times the values found in normal brain, are also found in GBM [78], promoting not only vascular growth but also microvascular permeability [67].

The role of angiogenic structures in GBM
During the last decades several attempts to seek treatments have been developed in order to fight against this kind of tumors. In line wih the mentioned observations, among tumors, gliomas have one of the highest degrees of neovascularization [62]. Then, some of these treatments have been focused in the angiogenesis process which is also typical of solid malignant tumors in general [30,32]; although it has been surprising that antiangiogenic therapies (like bevacuzimab) have to date limited success [39,88]. A possible cause is the specific structure of the neovascularization due to the improperly accelerated angiogenetic process, which results highly unstructured, glumeruloid and pseudosarcomatous, leaky, and with disruptive blood brain barrier (BBB) [62]. For each patient, the observed specific vascular structure from preoperative magnetic resonance images of GBM has been highly correlated with survival expectancies [48]. Furthermore, the role of the vascular structure in drug delivery and the development of promising GBM treatments has been emphasized [15], where disruption of the BBB and the high intersticial fluid pressure (IFP) are key players, which influence in the overall drug perfusion and the treatment success is known, but which accurate role determination needs research efforts. Therefore, it is important to understand the angiogenetic process and be capable of simulating it under different scenarios to provide much needed tools, for example, to take advantage of the BBB disruption window and predict drug perfusion [15], or investigate possible antiangiogenic/cytostatic combined therapies [62].

The process of angiogenesis
Angiogenesis is the process of the generation of new vessels from sprouting and branching the pre-existing vascular bed [20]. The basic scheme of sprouting angiogenesis is [56]: 1. Tip cell selection: an endothelial cell (EC) from the parent vessel changes its phenotype to become a migratory leading cell (tip endothelial cell, TEC) 2. Sprout extension: the tip cell migrates along the chemotactic path, thanks to the ECM degradation, so that a new capillary is created. 3. Lumen formation: connection of the luminal space of the sprout with the parent vessel.
Although human vasculature is normally quiescent in adults, with an equilibrium between pro-and anti-angiogenic agents, an appropriate stimulus ("angiogenic switch") can propitiate the growth of new capillaries [43] by favouring pro-angiogenic factors, which trigger the production of a new vasculature in few days [11,14]. As mentioned, this neovasculature has a specific structure clearly different from the primal one [62]. The most important group of pro-angiogenic factors are called VEGF (Vascular Endothelial Growth Factor) and the archetypic growth factor is VEGF-A [5]. The new vasculature induced by tumors entails new ways to release nutrients, drugs but also new ways to disseminate tumor cells in metastatic processes [63]. Solid Tumors cannot grow more than 1-2 mm without being provided by an additional amount of nutrients, such as oxygen, given by an appropriate vascular supply [49]. This is due to the fact that some tumoral cells are under hypoxia when tumors reach a size greater than the typical characteristic length of diffusion of oxygen ( ∼ 100 μ m) [72]. Hypoxia stimulates the migration of tumoral cells towards regions with higher concentrations of oxygen [25]. If hypoxia persists, tumoral cells release VEGF. When VEGF reaches some of the pre-existing blood vessels, ECs that line these vessels are activated so that one of them (the TEC) changes its phenotype from quiescent to migratory [45]. The adjacent endothelial cells are blocked from adopting a tip cell role in a process called lateral inhibition [51]. These ECs change their phenotype from quiescent to proliferative becoming stalk cells which are continuously being divided, due to the VEGF stimulation, allowing the conformation of the new capillary [38,86]. Its growth is conducted by some filaments deployed by TEC called filopodia that detect the levels of angiogenic factors and guide the formation of the new vessel towards the greatest VEGF gradient [76]. After the formation of the first new capillaries, a remodeling process can be produced when some vessels grow at the same time that the incipient vascular bed suffers some regression [9]. Regression is produced when a capillary is close enough to another one so that the oxygen locally released by the latter avoids VEGF concentration at TEC and therefore stops its growth [63]. Regrowth can also be promoted by the vascular scaffold of the basal membrane [57].
As well as filopodia guide TECs towards the gradient of VEGF, they are also able to detect closer filopodia from another TECs, or other vessels, promoting their merge and closing the blood flow circuit. This merging process is known as anastomosis [31,56,59]. The filopodia size is comprised between 3 and 7 μ m which means that when a TEC is near this distance, vessels anastomose [58]. The direction of growth has also an important role in anastomosis since filopodia deploy in circular sectors with an amplitude of around 60 • [83]. When contact has been produced, vascular endothelial-cadherin (VE-cadherin) definitely joins the two capillaries and stabilizes the joint by the formation of a luminal space between the anastomosed capillaries. Anastomosis stops by the loss of filopodial structures and when TEC recovers its quiescent phenotype [89]. In the case of tumoral vasculature, and specifically in high-grade astrocytic tumors (GBM), the new capillary bed presents a higher density of vessels irregularly shaped and structured, with possible dead-ends and poorer circulation than the regular capillary bed [73]. The extracellular matrix in tumors may obstruct these processes of local amplification, resulting in tumoral vasculatures that do not fit the typical pattern of compact and well-organized capillary networks [37]. Moreover, specifically, glioma cells adhere to cerebral blood vessels moving along them, invading the unaffected brain and lifting astrocyte endfeet [15,24,29,92], disrupting the BBB and possibly facilitating the infiltration of glioma cells [24]. Then, focus on "tumor remodelling" has been placed as a key strategy to improve chemoterapeutic drug delivery in GBM [15].
Hence, the experimental research of angiogenesis is crucial [47] since it would allow to design strategies to develop suitable combined treatments, to avoid metastatic cells dissemination, or to deplete oxygen supply to tumoral glial cells. Studying this complex multi-stage process in vivo is a difficult task [22], because just understanding isolated effects or predicting effect combinations in an equilibriumbased and structure-based system needs multiple parametric studies. At this point, the computational simulation of the angiogenesis process may help to understand this complex stage of tumors. However, the formulation of mathematical models capable of reproducing angiogenesis is a challenging task because the computational simulation of this complex process requires a high computational efficiency which can be achieved only from an adequate primal conceptual design.

Models for angiogenesis
Although the first approaches to model angiogenesis were based on very simplified equations [79], most of the mathematical models have been formulated through differential equations considering the tissue as a continuum medium. For example, Baish et al [6] developed an innovative approach based on percolation models that reproduced the tumor vasculature with a rule-based algorithm. Schugart et al. proposed a mathematical model in this framework for wound angiogenesis [71]. More recently, Cai et al formulated an extensive model that coupled glioblastoma growth, preexisting vessel co-option, angiogenesis and blood perfusion [18]. They considered successive stages that included a simulation of three-dimensional O 2 transportation [19]. Valero et al. [82] modeled wound angiogenesis in the skin coupling angiogenesis and wound contraction. Authors included biological, chemical and mechanical aspects in the modeling. Giverso and Ciarletta also formulated a chemo-mechanical model for tumoral angiogenesis considering sprouting as a consequence of a surface instability [40].
To complement the differential formulation and overcome some of the limitations, some hybrid models have also been developed coupling continuum and discrete formulations. Anderson and Chaplain [4] solved the differential equation that governed the endothelial cell density by finite differences and introduced a stochastic movement rule for the cells, extended later to 3D [21]. Recently, Suzuki et al. [77] complemented the Anderson and Chaplain model accounting for chemotaxis of TECs throughout the VEGF gradient produced by ECM degradation. Vilanova et al. [84,85] combined a phase-field method, to model the capillaries dynamics, with discrete agent-based models, to account for hypoxic cells and TECs. They also studied growth, regression and regrowth of tumor-induced angiogenesis using this approach [87]. Travasso et al. [81] also used the phase-field method to develop a model to study the chemotactic respose of endothelial cells and their proliferation in tumorinduced angiogenesis process. Some other multi-phase models that combine continuum and discrete approaches have been developed and applied to simulate angiogenesis in the growth of spheroids and melanoma [69], or in skeletal muscle [54]. Tumors and the angiogenesis process have also been modelled through purely agent-based approaches; see review [17]. In particular, Alarcón et al have extensively used cellular automata models [2] for tumor growth. This group also extended their approach to angiogenesis [65], also including different effects at different scales [66], and studying the influence on tumor growth of vessel normalization with respect to antiangiogenic and cytotoxic drugs [3].
There are relevant models available in the literature that graphically reproduce 3D angiogenesis processes, see for example [6,12,21,54,63,74,85]. However, they are based in mathematical and algorithmic formulations that involve large computational costs to obtain relevant simulations, an aspect which hinder their scalability to real problems. Furthermore, most of them are not concerned with tumor-induced angiogenesis or complement the graphic representation of angiogenesis process with suitable quantitative metrics in order to validate them [22]. As far as authors know, none of the available models have described how hypoxic conditions of cells change with the development of the new microvasculature (a very important aspect since it is the trigger of the process). We neither have found works addressing the structure (e.g. necrotic core) and relative position of the avascular tumor with respect the pre-existing vasculature (vital in possible drug perfusion studies). Thus, simple and efficient models capable of addressing new issues in angiogenesis with accuracy and with high computational efficiency (important to address real patient-specific analyses) are needed.
In this work, a simple, efficient and fully discrete 3D agent-based model is formulated to simulate the process of sprouting angiogenesis and anastomosis in GBM. It is based on simple algebraic operations grounded in the fundamental solution of the diffusion problem (Fick's law [52]) so that no differential equations or high-order tensor formulations are solved or approximated numerically. Both the total field of VEGF around the vascular network and the oxygen concentration in the cells are calculated by superposition of the contribution from the different sources for both species. Thus, the model formulation is highly efficient and only relies in a few parameters with clear physical meaning.
The model was implemented in Julia Language [13] to be integrated aposteriori in a global model of GBM being developed by the authors. Results are visualized in Paraview [1]. Additionally, a quantitative analysis has been performed, where different metrics have been measured throughout all the angiogenesis process, because several quantitative parameters are needed for a sufficient characterization and comparison of the angiogenesis process [50].
Results show that the reduction of the hypoxic conditions in cells is dependent on the TECs density at initial stages, whereas anastomosis plays an important role at final phases. The performed simulations show the importance of the structure and relative position of the initial vasculature with respect to the tumor in the final configuration in the development of the new micro-vascular network. As previously commented, the vascular structure of GBM is highly correlated to survival predictions, and its modelling will also be important in the development of therapies, both to guarantee proper perfusion of drugs and to design and understand the effects of combined therapies.

Theoretical framework
In this work, GBM cells and TECs are considered as individual agents. The former are static in the space (no migration process is considered in this work for GBM) whereas the latter migrate throughout the direction of the maximum VEGF gradient. GBM cells behave as sources of VEGF in hypoxic conditions and, along other cells, as sinks of the oxygen released by the vascular bed. Blood vessels are divided in small segments, segments hereinafter, so that they are simultaneously receptors of VEGF from cells and source of O 2 . For both O 2 and VEGF, Fick's law is a good approach to model their processes of diffusion [52]. Other models as the 3D O 2 transport model of Yan Cai et al [19] also adopted this assumption.

Diffusion equation solution
Fick's law represents the diffusion of matter in a porous media. Let J = −D∇C be the flux, and C(x, t) the concentration of some kind of matter, such as oxygen or VEGF. Fick's law is [34] where x = (x, y, z) is the spatial location, t is the time, D(x, t) is the diffusivity tensor, ∇C is the gradient of C and ∇ ⋅ (•) is the divergence of (•) . If we consider homogeneous isotropic media (which is an obvious and typical first order simplification), D(x, t) = DI , where D is the diffusivity constant and I is the identity tensor. Since this is a linear differential equation, we can apply the superposition principle. GBM tumors are much smaller than the brain, so we can also consider the solution in infinite media. Then spherical coordinates r = (r, , ) are handy in our problem, where a typical distance (e.g. to the source) r will be the relevant length magnitude. Assuming spherical symmetry and the aforementioned conditions, for 3D, the solution is [34] where √ 4Dt =∶ l is the characteristic length (or diffusion length) and m∕( 3 is the mean concentration in a characteristic l 3 volume. Ĉ (r, t) is the concentration change from the initial state due to an instantaneous point-source of strength m = 1.
If we consider now n s source functions m i (t) (probably constant within some time interval t 1 < t < t 2 ) at locations x i , the change of concentration at a given point x and time t is obtained by superposition as [42] (1) (2) is the convolution of the functions. For the case of a continuous source, when the substance is liberated under a continuous rate, say (t) as expected in the problem at hand, we can integrate Eq. (2), so taking reference as i = 0 (for notational simplicity, we omit now the convolution from different sources) which for a constant , defining the adimensionalized dis- we have that erf( ) → 0 and we recover the steady state solution, which is of the form C(r) = A∕r , where A∕r 0 is the steady-state concentration at r 0 ; i.e. A = ∕(4 D) and C 0 = ∕(4 Dr 0 ) , and with erfc(•) being the complementary error function and C ∞ (r) ∶= C 0 r 0 ∕r is the target steady state at location r, i.e. at t → ∞ . Note that C ∞ (r 0 ) = C 0 . Note also that the continuous rate mass constant is = C 0 r 0 (4 D) , i.e. it is proportional to the concentration.
Of particular interest below is the time derivative of the solution (6) and the gradient ∇C where r is the unit direction of the spherical variable r. Note that for t → ∞ , the gradient is constant, i.e.
where C ∞ (r)∕r is a constant at each r.

Diffusion and consumption of oxygen
Each unit segment length of vessel may be considered as a point source of oxygen. Applying the principle of superposition, see Eq. (3), the overall supply of oxygen is then computed by adding the supply from every unit length; Eq. (6) as where r i = |x − x i | is the distance between the point and the source, and i is the starting time for the source. The tissue surrounding the vessels, including cells of any type, is a sink of oxygen. We can assume that the consumption of oxygen is homogeneous, but depending on the concentration. The treatment of the distribution of oxygen in tissues is performed in terms of a characteristic length of diffusion, which is a fitting parameter, as l O 2 (∼ √ 4Dt) under an assumed steady (source-sink equilibrated) solution for the diffusion equation. The length l O 2 is the quantity measured and given in the literature. This assumes a steady-state balanced solution of the form given above, in which √ 4Dt is replaced by l O 2 , and the adimensionalized constant is now = r∕l O 2 . Then, if the solution for the oxygen concentration is of the form the flux is independent of time, namely J = −D∇C O 2 . Of course, in the steady solution, the concentration is also independent of time, so to comply with balance of mass, the divergence of the flux (net change in a control volume) must be the consumption by the cells in the tissue. This implies that-cf. Eq. (7) or equivalently, we can write the steady-state solution in rate form by leaving aside the only dependence on r in the exponential (the rest are constants) where Ċ O 2 (r) is the steady (rate) consumption of oxygen (units [(kg∕s)∕ m 3 ] ) in cells at a distance r of the source and ) is defined as the value at the vessel point r → 0 . Equation (13) just relates the steady-state source and consumption (availability) variations at a given distance in terms of the characteristic diffusion length; or alternatively at a relative distance = r∕l O 2 . We also assume that cancer cells occupy the location where previously normal cells were present, and that the consumption profile does not change; indeed the need for extra local consumption is what triggers the VEGF emission and generation of vessels towards such location to increase local oxygen availability and consumption.

Methodology
The proposed model can be implemented for any given initial arrangement of vasculature and hypoxic cells, with the same parameters and equations. Note, since we are not interested in transitory effects at times characteristic of diffusion, at any time interval we consider that the emission, flow and consumption of O 2 is stable, so the time increments may be relatively large, of the order of the events at cells. In essence, we consider valid a separation of time-scales, and to our problem only the larger one is of interest. Hence, the total process of angiogenesis will be divided in time steps so that at each time t the concentration of VEGF in the segments and the oxygen rate in the cells are calculated and updated through Eqs. (6) and (11)-or the rate form Eq. (13), respectively. From these data, modification of the VEGF emission at cells is assessed and the possibility of generating new vasculature and anastomosis is also determined. If the necessary requirements are met, new capillaries will grow and anastomosis will be performed. The vessels are all divided into segments. Each one of those is treated independently thanks to the decomposition capability of the formulation developed for the calculation of the VEGF and oxygen fields. It allows us to decouple the calculations in each term of the summations, leading to a calculation discretized by the cells and the vessel segments. We will name the set with all the vessel segments v i by , and the set with all the cells c i by ℂ . Within the set of the vessel segments, we will distinguish two types: the ones with a TEC, which are connected only in one endpoint to another vessel, and the ones without a TEC, which has both endpoints connected to other vessels.
The vessel segments with a TEC are the ones that will potentially grow by adding a new segment with a calculated orientation at the endpoint with the TEC. The only condition imposed in the growing is that the VEGF gradient needs to be positive up front, otherwise the vessel will not grow.
The vessel segments with a TEC can end as well in an anastomosis. Anastomosis happens when the TEC of the segment s i is close enough to another vessel segment s j . Then, the endpoint of s i with the TEC is merged with the endpoint of the vessel segment s j . The anastomosis threshold distance is determined by the parameter d a .
New branches can be created from any vessel segment without a TEC. To test this possibility, the VEGF concentration in the segments is evaluated, and if a certain segment has a local maximum of concentration compared with the neighbouring ones within a distance, then a new branch is created through the insertion of a new vessel with an endpoint in the vessel evaluated and a TEC in the other. The direction of this new segment will be calculated with the gradient of VEGF-see Eq. (16) below. In all the cases where a new segment is inserted, the length of this segment is the growing rate of the vessel g r multiplied by the time step of the calculation, g r Δt and weight constant (see below). Note that the actual value of VEGF concentration is not relevant, only relative ones are of interest in our model, because the growth rate is given independently by the parameter g r . Further model refinements could consider the influence of the concentration in g r . However, the amount from different cells, at different locations, may change the VEGF field and its singular points.
To calculate the VEGF concentration and gradient at every vessel segment s i , a loop is done over all the cells of the set ℂ with Equations (6) and (8). And to calculate the oxygen concentration at the position of every cell c i , the addition is carried out over all the cells of with Eq. (11). With computational efficiency in mind, both tasks may be performed simultaneously through a nested loop.
The model only has eight parameters, all with physical meaning and constant over the calculation. Those are given in Table 1.

Sprouting process: determination of the VEGF concentration and sprouting direction
For computing the function C V ∞ijs (r ij ) corresponding to the VEGF concentration from cell j at distance r ij from segment i, which is C ∞ (r) in Eq. (6) ). Then, we posit that, for every segment i, for cell j at step s where ⟨•⟩ is the Maccaulay bracket and C V 0 is the maximum VEGF concentration emission by a cell at r 0 (note that the actual relevant value is the strength C V 0 r 0 ). Then, when the target is met (so C Although there are other possibilities, for simplicity we considered C Table 1. Note that the VEGF concentration is only measured at the vessel segments to decide whether sprouting is initiated or not, whereas the oxygen rate is measured at the cells, to compute the scaling factor for the VEGF emission if the rate is insufficient. The treatment of O 2 is in rate (or steady) form of concentration because all the tissue is consuming oxygen, whereas the VEGF is in concentration because VEGF is only affecting the vessels, through the process mentioned in the Introduction. Note that the amount of VEGF arriving at a vessel segment depends on the VEGF emission by the different hypoxic cells j and the distance from the hypoxic cells to the segment i. For long times, the exponential (or erfc) contribution may be omitted (assumed ≈ 1 ), and the time stepping is just used to compute the effects of different events, so it is related to the cell-life cycle.
Consider now a vessel segment i. The gradient of VEGF at t → ∞ , if there are no further events, is-see Eqs. (9) and (14) which is constant, and renamed for short as g i . The unitary direction is denoted by ĝ i . Then, the VEGF acts over the wall of the vessel. The effective sprout direction is given by the projection of the gradient direction in the plane transverse to the vessel, i.e.
where v i is the segment direction and s i is the possible weighted sprout direction. Note that if the gradient is in the direction of the vessel, s i = 0 , so no sprouting takes place. Gradients perpendicular to the vessel have higher chance to promote sprouting. The TECs proliferation is restricted to one unit per segment. To simulate lateral inhibition, the concentration of VEGF in one segment is compared to the ones of its vicinity within a number of positions at its left and right within a characteristic length. If the level of concentration at the segment considered is higher than that at the neighbouring ones, a new vessel will be generated at this position. A number of segments within a characteristic length l inh are inhibited in order to avoid branching in subsequent steps (since there is not reliable data, we chose to inhibit the two neighbouring segments). The branch growth is produced following the gradient of VEGF, direction ĝ i . It is analytically calculated employing the time-dependent Eq. (8) using the approach explained for the equilibrated solution of Eq. (15).

Vessels growth
The gradient of VEGF is calculated also at the new vessel segments created in the previous steps, except in the previously anastomosed capillaries (see below). Once sprouted, for a nonanastomized vessel, the subsequent growth of the vessel is in the direction of the gradient, and the actual effective weighted rate is given by the filopodia at the front (projection in the direction of the vessel). Hence, in this case where the modulus acts as a scaling factor for g r .

Anastomosis assessement
To complete the ingredients of the simulation, a calculation is needed to assess if the tips are close enough to anastomose. The determination of which vessel is the closest one to a given tip may be assessed through the O 2 flux addends from other vessels, which can be computed just as if the TEC were a cell. If the addend is larger from a segment than that corresponding to the minimum distance of anastomosis, d a (an input parameter), they anastomose.

Parameters of the model and standard cases
The parameters that govern this model and the values adopted for the aforementioned cases are shown in Table 2.
c and D constitute VEGF diffusion parameters. On the other hand, C O 2 0 and l O 2 are the parameters corresponding to the oxygen diffusion. The process of anastomosis relies on d a that is the critical distance below which anastomosis is performed. Finally, the grow rate g r Δt is the length that blood vessels are enlarged after every time step.
The parameters have been estimated considering the literature. l O 2 has been estimated considering the porosity of the capillaries at different scales [70,72], and d a has been considered between 5 and 100 μ m, which are the range of typical lengths of filopodia [38,58]. Furthermore, the grow rates g r are consistent with the reference rates given by Brem and Folkman [14]. This value is also consistent with values adopted in other models like the one developed by Tong and Yuan for angiogenesis in the cornea [80].
The diffusivity D of VEGF in blood is given by Mac Gabhann et al [35] (note that the diffusivity of O 2 is immersed in l O 2 ). Pressure values of oxygen in GBM Tumors have been measured in the approximate range of 10-20 mmHg peritumoral and 2.5-15 mmHg intratumoral [10]. Brain tissue pressure pO 2 is given in values of about 12 mmHg for an arterial pressure of pO 2 = 100 mmHg [10]. Hypoxia levels are typically considered about 50% of those of normoxia levels. Cellular pO 2 levels below 1 mmHg are considered hypoxic levels [10].
The concentrations C V 0 and C O 2 0 and ranged according to the values of Salmon et al. [68], and the value C O 2 c is arbitrarily estimated within this range. On the other hand, VEGF concentrations used in studies are in the order of 6.6 ng/ml, e.g. [67], or [41] with a value of 10 ng/ml. Similar values have been found in GBM Tumors, which are one order of magnitude larger than values found in other brain Tumors [78]. Other works regarding to the VEGF influence show that there is a low variation of the global behaviour for VEGF concentrations between 10 and 1000 ng/ml [91]. Among this range, for this study the product C V 0 r 0 was chosen as 20 mm ng/ml, and can be trated as a single parameter of the model. It involves r 0 , which is the radial distance from the center of each cellular representative unit where C V 0 is achieved. In this case, due to the number of cells present in a glioblastoma, and an individual cellular size of 14 μ m [61], each point in the model represents a group of those cells, not necessarily with a radius r 0 . We have seen in the model that for Eqs. (15) and (6) the pair of the multiplication C V 0 r 0 can be considered constant for the different scales because the growing rate is fixed as g r for each case, and only the gradient direction, as an unitary vector, is required in the model.
The values of those parameters can be found within the ranges reported in the literature, but which is more significant in the model is the proportion among then, since the growing rate of the vessels is fixed and defined as g r . Only the growing direction and position of the local maximums are determined by the field of VEGF.
To test the scalability of the model proposed, there are validations at two scales, at mili and micro levels. The model is scalable with the same equations and algorithms by only changing the characteristic parameters of the equations. As shown in Table 2, in both models the concentration parameters are the same, since the key mechanisms and environment are equal and to the model, only relative quantities are relevant, but the parameters of the oxygen distribution and growing rate of the vessels and capillaries are different [14]. This multiscale versatility of the model will allow us to reproduce both the angiogenesis initiation in an incipient tumor, and the final vasculature in a big tumor. For this latter target in this paper the cells are static and the tumor is fully formed, but when mixed with an ABM model, the full process of growing and continuous dynamic remodelling will be modelled more realistically.
We present five cases with the same parameters but different distribution of the hypoxic cells, and a comparison with the continuous model of Vilanova et al. [87]. Those are summarized in Table 3, where the shape of the hypoxic cells in the tumor and the computational cost of the simulations (elapsed time) are also shown.
The computational cost is the elapsed time of calculation of the case studied in a HP ProBook 440 G7 Intel Core i5-10210U CPU 1.60 GHz-2.11 GHz and 8 GB of RAM. The code can be parallelized, but such parallelization is not implemented in the results shown.

Results
This section is focused in the validation of the model and the exploration of its capabilities. In the first example, we compare the results of our model with the results of Vilanova et al. [87] who developed a 2D continuous phase-field model. The subsequent cases are for different 3D tumoral configurations to show the behaviour of the model. Our purpose in those cases is to show the capabilities of the model at different scales, shapes and tumoral configurations, and assess whether the results are reasonable and according to typical qualitative observations.

Case 1: comparison with a continuous phase-field formulation
Vilanova et al. [87] proposed a 2D continuous model for angiogenesis. For continuous Finite Element formulations there is an important computational difference between 2D and 3D formulations. In this work, since our model is 3D, the tumor is modelled with a penny shape (i.e. like a coin, so the problem basically lies in a plane) aligned with the initial vessel. The geometry is the same as the one used by Vilanova et al. [87], with a circular primary vessel and a tumor with a radius of approximately 315 μ m. We have modelled the tumor with 4000 groups of hypoxic cells randomly distributed within the volume domain, each one represented by one point emitting VEGF and consuming oxygen. The number of selected points ensured statistically consistent quantitative measured results. The minimum distance between the border of the tumor and the primary vessel is of 0.06 μ m. The parameters of the model are the ones listed as small capillary in Table 2.
The graph of Fig. 1 shows a quantitative comparison between the topological results of both models given in Fig. 2. The final values, shape and the peak of TECs are very close, as well as the time when all the cells are oxygenated and the angiogenesis process ends. Figure 2 shows the temporal evolution of the proposed model, and the comparison of the vasculature developed at 2 days with the result of Vilanova et al. [87], which is the point when both models calculate the absence of hypoxic cells.
The snapshots of Fig. 2 are the initial penetration of the vessels in the tumor (0.34 days), the growing of the number of TECs (0.7 days), the growing of bifurcations (1.12 days), the peak of the number of TECs (1.4 days) and the final configuration (2 days). From both figures, it can be deduced that the model of Vilanova et al produces many more TECS initially which thereafter produce bifurcations. However, the final numbers are of the same order both in TECS and bifurcations. Nonetheless, the proposed model has more tertiary bifurcations and these are apparently also more anastomosed. Regarding the final configuration shown in Fig. 2, in Villanova et al the primary vessels cross the tumor and seem to continue the same path, whereas in our model the vessels stick to the contour of the tumor. Indeed, this is a tendency observed in the simulations; see that there are branches in the results from Vilanova et al that grow without ever touching the tumor, an observation which is not present in our simulation.  1 3 To further analyze the soundness of the predictions given by our model, we have performed for the simulation of Vilanova et al some parametric studies. Figure 3 shows the behaviour of the vasculature with different tumoral shapes, considering as a reference the model of Fig. 2 with the corresponding model parameters.
First, we analyze the possibility of randomness in the distribution or increasing tumor population by randomly activating the cells. The case 1 shown in Fig. 3 is a case where the start of the VEGF production of the cells is random from 0 to 1.4 days. The result at 2 days is to be compared with that given for the same time in Fig. 2. Whereas the morphology is different because of the stochastic nature and different activation times, they are very similar in qualitative distribution.
In cases 2, 3 and 4, all cells are activated at the same time, but we study different tumor geometries for the same model parameters as in Fig. 2. In case 2, we deactivate the cells in a zone in the proximity of the primary vessel. It is seen that since there are now several local maxima of VEGF concentration in the primary vessel, two families of vessels grow initially from different zones closest to the tumor, and soon another third family grows in the middle of the previous ones because of the trayectories of outward initial ones, which do not give sufficient suply to the central zone until later. In case 3, another zone in the center is deactivated, and the simulation shows that the central branch of case 2 does not appear until later, because there are fewer cells in hypoxia to be feed by that branch. Note also that the zones without hypoxic cells are not vascularized by the simulations. Finally, in case 4 we study the possibility of two hypoxic groups, one closer to the main vessel and one further away. It is seen in the simulations of Fig. 3 that, as expected, the simulation shows a primary vascularization of the closest tumor and then a new branch appears to oxygenate the furthest one. That is the only case when a vessel grows outside the tumour: when there is another VEGF production centre outside.
In Fig. 4, we study the influence of the model parameters. We have taken the parameters from estimations obtained from experiments in the literature. However, one should wonder the actual relevance of the parameters in the predictions; i.e. which would be the predictions if the parameters were one or two orders of magnitude different? Arguably, if the predictions are the same, the physical relevance of those parameters, and the related terms, should be minor and no relevant conclusion should be drawn. The parameters of the simulations in Fig. 4 are also those used for Fig. 2, see Table 2, but with just one parameter change in each case. The modified parameters are: Case 1 ( l O 2 = 0.2 mm ), Case 2 ( l O 2 = 0.002 mm), Case 3 ( g r = 6.0 ) and Case 4 ( g r = 0.06 ), which are variations of one order of magnitude with respect to the experimentally measured values and used in Figure 2. All those variations show a dependency on respective parameters of the global behaviour of the vasculature, not only a local fit of the result. Figure 4 shows as well the amorphous vasculatures created by the model when the parameters are far from the experimental measurements shown in the literature [14,70,72], and which are estimated as l O 2 = 0.02 mm and g r = 0.6 mm/day in the validation of Fig. 2.
In Case 1 of Fig. 4, the high l O 2 satisfy the oxygen level required by all the tumoral cells at an early stage of the vasculature growth (in blue in the figure). At that moment, the cells greatly decrease VEGF production, stopping the vessels growth. In Case 2 of Fig. 4 the length l O 2 is reduced, and the vasculature can not supply enough oxygen to the cells. The consequence is that the vessels can not significatively modify the VEGF field by supplying oxygen to the cells, and the resulting vasculature grows also slowly and is amorphous. This shows the relationship between the experimentally measured parameters and the "natural" shape of the final vasculature. In Case 3 of Fig. 4, the modified parameter is the growing rate ( g r ), which is much higher than the experimentally measured value. The observed consequence in the simulation is an abnormal growth of the vasculature and a lack of efficiency in the oxygen supply of the tumour when compared with the cases with the experimentally measured parameters. Since the model has a restriction of a minimum distance between new vessels (lateral inhibition), the final configuration is blocked and the oxygen is not efficiently distributed. Case 4 of Fig. 4 shows the opposite case, where the growing rate ( g r ) is set to a low value. It implies a very slow growth rate in which the vasculature is far from satisfying the oxygen requirements of the tumour. In summary, the predictions of the model are consistent with the physical understanding, and the parameters of the model are meaningful and have a crucial impact in the simulations.

Case 2: angiogenesis in a large and plane surface
This is the case that can be found in some artificial implants where a large vasculature is aimed to be developed [60].
Here, the parameters of the bigger capillary in Table 2 are used in the calculation. The region with hypoxic cells is a trapezoid with a height of 30 mm, one base of 10 mm and other of 60 mm. These are large dimensions for tumors, but the purpose of this case is to test the behaviour of the model at large scales as a complement of the case shown in Sect. 5.1. The setup is the one shown as Case 2 in Table 3, with 5000 groups of cells to discretize the VEGF production and the oxygen level. Figure 5 shows the temporal evolution of the TECs, anastomosis and total number of vessels, as well as the number of cells in hypoxia in the secondary axis. Figure 6 shows the key snapshots of the quantitative results of Fig. 5. Those are the initial penetration of the vessels in the hypoxic area (160 days), the increment of vessels and TECs (320 days), the initiation of anastomosis (415 days), time with more anastomosis than TECs (575 days), peak of TECs (720 days) and the final configuration (810 days), when all the area is oxygenated.

Three-dimensional analysis in a spherical tumor
The model is tested in 3D spherical tumors with three different configurations of hypoxic cells within the volume. Table 3 shows the configuration of cases 3, 4 and 5. In all the cases, the cells in hypoxia (modelled as discrete points that represent a group of cells) are in red and become transparent when oxygenated. The colors of the vessels reproduce the generation of each one, with blue being the first generation and in red the sixth generation. We consider that a vessel has one more generation than the vessel from which is bifurcated. We have again checked that given a certain minimum number of discrete points to reproduce the cells, the results are basically unchanged when that number is increased.

Case 3: concentrated hypoxic cells in the centre
This setup is the one denoted as Case 3 in Table 3 and the large capillary configuration of Table 2. The hypoxic cells are modelled with 1000 discrete points. The tumor is at a large distance from the primary vessel (Fig. 7).
In this case only a single ramification appears to supply all the cells. Figure 8 shows the snapshots when the vessel progresses to the tumor (100 days), the vessel touches the tumor (115 days), when half of the cells are oxygenated (130 days), the end of the oxygenation process (145 days) and the end of the calculation (no further changes take place). This setup is the one denoted as Case 4 in Table 3 and the large capillary configuration of Table 2. The hypoxic cells where modelled with 1000 discrete points. Figure 9 shows how the density of vessels in the tumor increases exponentially from a certain point, and how the anastomosis become important in the angiogenesis process. Figure 10 shows snapshots of the calculation when the vessel touches the tumor (100 days), the start of the ramifications (200 days), the initiation of the exponential growth of the density of vessels (300 days), the initiation of the anastomosis (350 days), and the end of the calculation when all the cells have enough oxygen supply.

Case 5: hypoxic cells distributed in the surface of the sphere
This setup is the one denoted as Case 5 in Table 3 and the Large Capillary configuration of Table 2. The hypoxic cells where modelled with 1000 discrete points. The purpose of this example is to simulate the case of a necrotic core within the tumor, where it is known that the vessels are more superficial to the tumor. In Fig. 11 we show the same measures of the angiogenesis process, which has a good resemblance with those obtained for the previous case. However, as seen in Fig. 12, the spatial distribution of the vessesls is different, closer to the surface of the tumor. Figure 12 shows the significant points of Fig. 11. Those are: the vessel touches the tumor and starts the branching (75 days), a secondary vessel touches the tumor (150 days), the exponential growth of the vessels (230 days), the growth of anastomosis (275 days) and the final configuration when any cell is in hypoxia.

Case 6: three-dimensional analysis in a cylindrical tumor
This setup is the one denoted as Case 6 in Table 3    Case study 3 (spherical tumor at a relative large distance from the primary vessel). Angiogenesis temporal evolution secondary vessels appear in multiple points and with multiple directions. Figure 14 shows the significant points of Fig. 13: the branching initiation (60 days), the exponential growth of the vessel number (125 days), the exponential growth of the anastomosis (190 days), a latter configuration (220 days) and the final configuration (280 days).

Discussion
Several models that graphically reproduce 3D tumor-induced angiogenesis processes are available in the literature, [6,18,21,63,74,85]. However, these models address the VEGF and oxygen diffusion through more complex mathematical formulations which results in models with a great number of parameters and considerable computational effort. As far as authors know, none of the models available in the literature have simulated the 3D cases shown here and none of them have reproduced the graphic evolution of the hypoxic conditions in cell locations during the angiogenesis process. Even if there had been a chance to graphically validate the model comparing our results with other microvascular networks, the qualitative comparison would have been quite difficult only considering the angiogenesis evolution or the morphology of the final configurations. This is the reason why qualitative results have been complemented by a multiparametric quantitative analysis. Then, quantitative analysis is something that many studies considered to be indispensable to provide a valid holistic representation of brain Tumors and their vasculature [26]. Use of quantitative measurements have also been proposed for experimental characterization of angiogenesis processes [28], given the observed large variability in the resulting vascular structure. The parameters used in this kind of analysis are consistent (obtained with similar figures for different similar examples) and based in density measures, instead of absolute quantities, namely, vessel density, branch point density or area/volume occupation. These density parameters allow to compare different studies [22]. Although the cases presented here have not been studied in any other research, some comparison can be performed with other models. For example, the number of vessels represented by Cai et al. [18] show the same tendency as the corresponding to this model. On the other hand, the pattern of TECs, consisting of rising and eventually decaying in favor of anastomosis, is also found in the model results of Vilanova et al. [87], see Case study 4 (spherical dense tumor at a medium relative distance). Quantitative measurements of the angiogenesis process. The number of vessels per unit volume are divided into the ones ending in a TEC and the ones ending in an anastomosis

Fig. 10
Case study 4 (spherical dense tumor at a medium relative distance). Angiogenesis temporal evolution particularly wide range of possible values for five vascular metrics, length density among them. The results presented in this paper are inside of the ranges given by these authors.
In particular, some aspects from the simulation results may be highlighted. From case 1, it is seen that quantitative aspects as number of TECs and bifurcations are similar as those found in continuum models, but with a slighly different geometric distribution. Case 2 shows that the model may be applied not only to tumors, but also to predict angiogenesis in artificial implants. The comparison of the predictions from cases 3 to 5 brings several interesting observations. First, the density of vessels and time needed to oxygenate the cells strongly depends on the distance of the tumor to the vessel and on the relative size of the tumor respect that distance. If the distance is significatively larger than the size of the tumor, only a branch will be generated, which will grow until reaching the tumor. This is expected because the maximum location of VEGF in the primary vessel is well defined. In the cases where the tumor is closer to the vessel, several primary branches may be generated. This is specially the case when there is a necrotic core (VEGF emision has a more difusse location) or the vessel surrounds the tumor (case 1); see Figs. 2 and 12. This observation also applies to case 6, Fig. 14, where the tumor is around the vessel. Second, the patterns of TECs and vessel curves, Figs. 1, 5, 9, 11 and 13, are similar. First, primal branching is developed directed by a strongly directional gradient, and when these branches reach the hypoxic zones, where the gradients are weakly directional (more constant), secondary and tertiary branches are developed. This observation is also consistent with the results in cases 1 and 3. Then, necrotic cores have also a relevant impact in the subsequent vessel structure.
Finally, we note that we have used for each problem only one class of capillary (large or small). A more realistic extension for future work is to consider simultaneously in the same problem different types of capillaries and other possible configurations. Further extensions of the model will facilitate more realistic simulations (accounting for many neglected effects) and hence, a more detailed analysis of the angiogenesis process and its possible impact in understanding GBM and the potential treatments.

Conclusions
The simulation of angiogenesis processes and of the resultant vasculature is crucial in modeling many aspects of the dynamics of tumors, as immune cells access, inflammation, drug access, cell behavior, migration and metastasis; all of them important in understanding those dynamical processes and in developing efficient therapies. The process of tumor-induced angiogenesis, including anastomosis, has been modeled through a simple 3D mathematical model amenable of very efficient computational implementation. Different cases have been considered, and one of them compared quantitatively and qualitatively with continuum-based modelling. The purpose of this model is to be combined with agent-based tumor dynamics, including proliferation and migration and FEM-based nutrient diffusion, so an important pursued advantage of the model is the very low computational cost. Although the model has been motivated in GBM, where the simulation of angiogenesis has especial relevance because of the characteristics of this tumor, it can be applied to any kind of tumor. This model is a key part of an in-silico testing platform being developed for GBM.
This approach accounts for the evolution of the hypoxic conditions in static tumor cells. To this respect, TECs play an important role at early stages whereas anastomosis is paramount in the last ones. The position of the arteriole relative to hypoxic cells, or vice versa, and the relative size to distance, is quite important in the final structure since new branches are firstly sprouted towards the hypoxic cells closer to the main arteriole dominating the gradient direction. The relative position of the initial vasculature is important if angiogenesis simulations wants to be applied to clinic cases. For that reason, an accurate mapping of the quiescent vasculature and the position of avascular tumor may be of paramount importance if a reliable reproduction of angiogenesis process is required. This aspect is important in understanding drug delivery dynamics and in addressing strategies to avoid the nutrients supply to tumor or metastatic process that could disseminate tumor cells to other part of the body.