Mathematical modelling of angiogenesis using continuous cell-based models

In this work, we develop a mathematical formalism based on a 3D in vitro model that is used to simulate the early stages of angiogenesis. The model treats cells as individual entities that are migrating as a result of chemotaxis and durotaxis. The phenotypes used here are endothelial cells that can be distinguished into stalk and tip (leading) cells. The model takes into account the dynamic interaction and interchange between both phenotypes. Next to the cells, the model takes into account several proteins such as vascular endothelial growth factor, delta-like ligand 4, urokinase plasminogen activator and matrix metalloproteinase, which are computed through the solution of a system of reaction–diffusion equations. The method used in the present study is classified into the hybrid approaches. The present study, implemented in three spatial dimensions, demonstrates the feasibility of the approach that is qualitatively confirmed by experimental results.


Introduction
Angiogenesis is the process through which a new blood vessel is formed from a pre-existing blood vessel network. An adequate blood vessel network is required to supply blood to the entire human or animal body. In cases of (mechanical) damage, like a wound, the small blood vessel network in the wounded area has been disrupted and needs to be restored. In other cases of damage, one can think of the re-establishment of a vascular network around cardiac (coronal) arteries that may have closed as a result of atherosclerosis. In the aforementioned processes, angiogenesis is indispensable for the survival of the organism. During the early stages of development of a tumour, growth occurs through cell division and proliferation. Subsequently, it halts as a result of lack of oxygen and even develops a necrotic core. Finally, it is able to continue growing if a vascular network around the tumour has been developed. Here, angiogenesis is responsible for turning a benign tumour into a malignant tumour, which will possibly metastasise (or spread out) to other parts of the body, often leading to morbid and mortal consequences for the patient.
In order to understand the underlying mechanisms of angiogenesis, it is important to carry out experiments both in in vitro and in in vivo settings. Since qualitative (images) and quantitative (after analysis) results are obtained from these experiments, it is important to quantify and to test the hypotheses that are formed after theoretical assessment and analysis of the results. Therefore, mathematical modelling of phenomena like wound healing, wound contraction, tumour growth, ulcer development and angiogenesis has become very important and has developed into a mature state. The maturity of the modelling can be seen from the number of approaches that are used to simulate the aforementioned processes. Some of the approaches treat cells by the use of averaged quantities where the models end up as continuum-scale partial differential equations for cell densities. Examples of such continuum-scale approaches in the context of wound healing are Britton and Chaplain (1993), Javierre et al. (2009), Prokharau et al. (2014), Valero et al. (2014) and Gaffney et al. (2002), Maggelakis (2003), Maggelakis (2004) in the context of angiogenesis. The list of examples is far from complete. Next to the continuumscale approaches, several formalisms have been developed on the smaller cell colony scale, where we start with mentioning the relevant work in Oers et al. (2014), Graner and Glazier (1992), Merks and Koolwijk (2009) on cellular Potts modelling in the context of angiogenesis. The cellular Potts models fall within the class of cellular automata models, which divide the computational domain into a discrete set of lattice points. Each lattice point is either occupied or not occupied by a cell (or by one of the subdomains) based on several biologically derived constrained optimality principles. Since cells or their boundaries move in a discrete fashion, and since intra-cellular adherence can be built in easily as a penalisation, the use of cellular Potts models has become a very natural choice for the simulation of angiogenesis where endothelial cells move and stay attached to each other. On the same cell colony scale, we mention the semi-continuous approach, where cells are treated as discrete entities, but where their migration is not restricted by any lattice points. Here, several modelling approaches have been developed in the context of wound closure, wound contraction, cell migration, and tumour growth and development. Examples of such models are the studies in Byrne and Drasdo (2009), Groh and Louis (2010), Mousavi et al. (2013), Neilson et al. (2011), Rey and Garcia-Aznar (2013), Vermolen and Gefen (2012), , . A recent review on particle methods applied to wound healing and tumour growth can be found in Vermolen (2015). Regarding cancer initiation, growth and invasion of cancer cells, and the use of cell-based modelling, we refer to the studies by Schlüter et al. (2014) and . One can distinguish between models in which the cell geometry does not change over the simulation and those models where cells actually geometrically deform. An example of 3D models where cells deform and migrate, using a probabilistic voxel finite-element method, is given in Borau et al. (2014). This voxel finite-element is stochastic, as well as discrete, whereas some cell deformation models are based on phase-field methods or on moving surface partial differential equations like in, respectively, the work by Marth and Voigt (2014) and Elliott et al. (2012).
Looking at the process of angiogenesis, one finds the classical continuum-scale models along with the cellular automata approaches; however, hybrid approaches combining cell-based approaches with finite-element simulations are very scarce in the literature. Regarding bone growth and angiogenesis, we refer to the work of our Belgian colleagues (Carlier et al. 2012) where endothelial tip cells are moving individually in a lattice-free manner and where other cells are treated in terms of cell densities. Chemotactic and haptotactic signals determine the migration of tip cells. Several chemotactic factors as growth factors are taken into account by approximating the solution of a system of diffusion-reaction equations. This interesting work also treats impaired angiogenesis in a framework with two spatial dimensions, where angiogenesis is considered in the context of bone formation. Another interesting study on modelling angiogenesis where mechanical cues were taken into account was done by Stéfanou et al. (2015). The present approach that we consider in this study deals with the modelling of an experimental in vitro setting, where a fibrin matrix is considered with a confluent monolayer of endothelial cells on top submerged in an extracellular fluid. Our approach is three dimensional where all endothelial cells are treated as soft spheres. Further, we distinguish between tip and stalk cells, and opposed to Carlier et al. (2012); stalk cells and tip cells are able to differentiate to either (sub-)phenotypes at all times, see Tammela et al. (2011), Tung et al. (2012) and Blanco and Gerhardt (2014) for an experimental justification. The concentrations of the growth factors are treated analogously to the work in Carlier et al. (2012) in terms of diffusion-reaction equations. However, we also include proteins that are secreted by the tip cells through mathematical point sources which make the stalk cells follow them based on the mechanism of chemo/haptotaxis. Next, to chemotaxis, a durotaxis term is added such that the cells preferably stay near the transition of the extracellular fluid and fibrin matrix (hence approximately at the transition from fluid to matrix, that is, approximately at the basement membrane). To the best of our knowledge, we believe that these additions are innovative and complementary to the existing literature.
The paper is organised as follows: in Sect. 2, the experimental set-up is presented, and subsequently in Sect. 3, the mathematical formalism and the numerical solution strategy are presented. This description is followed by the presentation of the simulation results in Sect. 4. Finally, in Sect. 5, the model is discussed and some conclusions are drawn.

The experimental set-up
The dermatology department of the VUmc carries out several in vitro assays using primary Human Dermal Tissue Endothelial Cells (DTECs) on different substrates like fibrin. In this study, our particular focus is on the sprouting assay which uses a fibrin matrix and which is carried out in a standard 96well plate depicted in Fig. 1, see Koolwijk et al. (1996) for the development of the assay. In this assay, angiogenic responses to the angiogenic growth factor vascular endothelial growth factor (VEGF) are measured for different concentrations. On the first day, a v F = 100 µL fibrin matrix (3 mg mL fibrinogen with 0.5 µg mL thrombine IIa) is placed in a total of 39 wells on top of which a 100 µL solution is poured containing around N = 20.000 ECs. The total volume in the well then is v = 200 µL. Experimental observations show that ECs have a typical diameter of around 45 µm, and hence, a radius of R = 22.5 µm. ECs are ellipsoidal being twice as long as wide. The ECs sink and adhere to the fibrin matrix, thus forming a confluent monolayer covering the surface of the fibrin matrix as depicted in the microscopic images in Fig. 2.
On the second day, the cells are stimulated using different conditions. Three wells serve as baseline controls, where no growth factors are added. All other wells are treated with 2 All wells with a "T" have a 2 ng mL = 2 × 10 −3 µg mL (microgram per millilitre) TNF-α solution added. Wells with a "V" have VEGF added to them in the given concentrations in mg/mL. The numbers behind "VT" stand for the amount VEGF in the unit of mg/mL added to these wells ng mL TNF-α to maintain and activate the monolayer of EC. In addition, most wells are treated with additional growth factor VEGF in different concentrations. All different concentrations are replicated in threefold to compare the results, and the well numbers are used to label the microscopic images. We summarise the different concentrations in Table 1.
Depending on the donor-specific endothelial cell motility, fibrin matrices are fixated 48-72 h after stimulation. The sprouting into the fibrin matrix is observed using microscopic images like those in Fig. 3. In this figure, we see cells stimulated with VT25. The monolayer is roughly undamaged, except for a couple of circular-like structures with dark edges. These dark edges form the premises of the newly formed sprouts and are most likely the effect of the fibrous layer underneath the monolayer bending out of the focal reach of the microscope. In Fig. 4, we zoom in on one of the sprouts, where one of the sprouts has been indicated by an arrow. Inside the sprout, the fibrin matrix is degraded and this shows up slightly lighter on the microscopic image. We can see that no ECs show up in the image inside the sprout. This  Dermal ECs in a well after stimulation with 25 ng/mL VEGF and 2 ng/mL TNF-α. The circular structures form the boundaries of newly formed sprouts. This figure represents a magnification of ten times with respect to Fig. 3. One of the sprouts has been indicated by an arrow is due to the fact that the sprouts move into the matrix and get out of focus in the microscopic image. The amount of sprouting in an assay is quantified using image processing software. The darker edges of the sprouts are coloured, and the cumulative area of the coloured regions is calculated as a percentage of the total area of the image. This percentage will be called P(t), and this variable will serve as a measure of sprouting. Although it cannot be seen in these microscopic pictures, we know that sprouts are, as a rule of thumb, twice as deep as their diameter at the top of the fibrin matrix. One can conclude this by varying the focal depth of the microscope. Sprouts usually are in downwards direction, but slightly bending sprouts are also observed.

Driving forces on cells in sprouting angiogenesis
The motility of cells on the fibrin matrix is subject to many mechanical and biological factors. We identify several factors driving the movement of cells on the fibrin matrix. The mathematical formulation of these principles is covered in Section 3. For detailed cell biological descriptions of the hereafter listed phenomena, we refer the reader to the extensive work on cell movements by Bray (2001). Gamba et al. (2003) and Serini et al. (2003) describe chemotaxis as the movement of cells in response to a chemical stimulus. One speaks of positive (negative) chemotaxis if the movement is in the (opposite) direction of the gradient and the chemical is called a chemoattractant (chemorepellent). Chemoattractants can be, following the Keller-Segel model formulated by Horstmann (2003), secreted by the cells themselves, leading to the formation of isolated clusters of cells. The chemotactic process takes place thanks to pseudopodia on the cell membranes that are formed on the sides of the cell in high concentrations of the chemoattractant and "reach" towards higher concentrations, pulling the cell in the desired direction. Inflammatory mediators such as TNF-α may increase the motility of cells.

Cell-cell forces, contact mechanics
Cells can adhere to each other by physically attaching their cell membranes using surface proteins like cadherins. ECs adhere to each other using vascular endothelial cadherin (VEcadherin) bonds. VE-cadherin at the same time works as an inhibitor of haptotactic movement caused by VEGF by binding to the same receptor used in the chemotaxis signalling pathway. Merks describes this contact inhibition in his cellular Potts Model in Merks and Koolwijk (2009). ECs have a certain optimal elliptical shape induced by their cytoskeleton and will try to elastically return to this shape upon deformation. The magnitude of these forces is proportional to the elasticity of the cell and the severity of the deformation. This deformation can be caused by cells colliding into one another. We will denote this effect by contact mechanics in further chapters.

Cell-matrix forces, durotaxis
Transmembrane integrin proteins on the cell membrane adhere to fibrous scaffolds such as fibrin matrix or collagen and exert contractile forces causing cell-matrix adhesion. Since these forces are caused by physical attachment to the fibrin matrix, the net force will be in the direction of the fibrin matrix gradient. However, high-density fibrin matrix may be too stiff for the cells to move into. The same cellmatrix adhesive forces cause strain in the elastic fibrin matrix, which on its turn is sensed by other cells adhering to the matrix, and they get pulled along the stress lines. This effect is called mechanotaxis. Reinhart-King et al. (2008) conducted a series of experiments considering the interplay between cellcell adhesion and mechanotactic forces for endothelial cells. They conclude that matrix stiffness is an important factor for the cell motility and the ability to mechanically communicate through the substrate.

The mathematical model
First the model formulation is presented, and this is followed by the presentation of the numerical method.

The mathematical formulation
In this section, we present the governing equations with their boundary and initial conditions. We consider a threedimensional cylindrical domain Ω with boundary ∂Ω. Initially the domain is divided into three segments: Ω E , Ω B and Ω F , denoting the regions occupied by the, respectively, from top to bottom, extracellular fluid, basement membrane and fibrin matrix, see Fig. 5 for a sketch. The basement membrane can be considered as a somewhat stiffer top layer of the fibrin matrix. Since the tip cells will chemically create holes through the boundary membrane and fibrin matrix by degrading fibrin, the extracellular fluid will occupy the channels formed by the tip cells. To this extent, the biological problem could be considered as a moving boundary problem. This approach will not be used, and in the approach that we propose, we introduce the volume fractions of fibrin matrix, basement membrane and extracellular fluid, which are, respectively, denoted by f F (t, x), f B (t, x) and f E (t, x), where x = (x, y, z), being the coordinates of the location. Hence initially, we set Since the basement membrane and fibrin matrix are similar collagen-structured materials and since the extracellular fluid is a fluid, we introduce the variable f S (t, x) := f B (t, x) + f F (t, x) being the solid fraction. In our mixture formulation, we require that f S (t, x) + f E (t, x) = 1 at all times t and at locations in Ω.

Fig. 5
A schematic of the chemical interaction between the tip and stalk cells in which the tip cells secrete the the DLL4 to make the stalk cells follow them. The tip cells start migrating as a result of the gradient of the VEGF. Further, the arrangement of the fibrin matrix, basement membrane and extracellular fluid is shown, as well as the degradation of the basement membrane and fibrin matrix by, respectively, the chemicals MMP and uPA. In the model itself, gravity is not dealt with; however, in the simulations, the cells are seeded on the top surface of the basement membrane by first positioning them on top of the extracellular fluid and let them "sink" (by gravity) onto the top surface of the basement membrane to get a somewhat more randomised arrangement of endothelial cells as initial configuration for the simulations. In the sprouting assay set-up used at the VUmc, gravity together with contact mechanics forms a reasonable explanation for the formation of the initial confluent mono-layer. This initial configuration is also determined by the contact forces that the cells experience when seeded on the top of the basement membrane

Cell dynamics
The endothelial cells are treated as discrete spheres with radius R. This simplification has been chosen to facilitate a non-complicated approach for the intercellular contact forces. In case of ellipsoid cells, then, contact mechanics would need the determination of the points and angles of contact to compute the resulting direction of cellular displacement. The treatment of the contact forces in the case that cells are colliding against each other has been illustrated in Fig. 6, where the cells are considered to have collided if h > 0. In the case that two cells collide, the contact force is directed in the direction of the line connecting the centres of the cells. In case of multiple cells that collide, then the net contact force is obtained from a linear combination from all contributions from the separate cells. More information regarding this topic can be found in Vermolen and Gefen (2013). We distinguish between tip (leading) and stalk endothelial cells, see Fig. 5. The only phenotype we consider here is endothelial cells. The spatial positions of the cells are denoted by x i (t), where i and t, respectively, denote the cell index and time. To distinguish between the tip cells and stalk cells, we introduce the set of all the stalk cells, Fig. 6 A schematic of the contact forces caused by partial overlapping of two spherical cells. Note that the picture displays a two-dimensional representation, whereas the implementation is in three spatial dimensions and its complement Note that the transitions of the cells between the two states "stalk" and "tip" make the sets time dependent. It is assumed that the gradient of the vascular endothelial growth factor concentration, c V (t, x), drives the chemotactic movement of the tip cells. Next to the chemotactic signal, we take durotaxis into account where we realise that the stiffness of the material is proportional to the fraction of the volume of solid, see Griffith (1921) where experimental results confirming that solid is stronger than liquids. We take durotaxis into account by considering the gradient of f S (t, x). We also postulate that the cells are not willing to move into a very dense solid. To this extent, we want the cells to move towards the centre of the solid-liquid interface, denoted by (t), which we define implicitly by the surface We interpret the cell's willingness to reside at locations near (t) as their adherence to the fluid-solid interface. To this extent, the tip cells move according to Here, x i (t) denotes the cell's centre position at time t. The first term of the right-hand side in the above equation denotes the component of migration as a result of contact forces between neighbouring cells as well as forces that cells exert on the substrate, which are sensed by the other cells if the strain energy is large enough. Here M denotes the strain energy density and z i denotes the direction of movement determined by the strain energy density. The variable α models the mobility of cell, as well as its viability and the friction forces applied onto the cell surface as it moves over the solid material inside the channel. This term has been detailed in Vermolen and Gefen (2012), but for this study we implemented the following changes with respect to Vermolen and Gefen (2012): -The formulation has been extended to three-dimensional geometry, as in Dudaie et al. (2015); -The formulation involves regions with different structures, and in the extracellular fluid, there is no long-distance communication that is associated with longdistance mechanical signals; -The cells are moving around in a non-homogeneous lattice where the elasticity modulus changes over time and location. This issue is currently dealt with by averaging the elasticity modulus used to determine how the strain energy density signal attenuates over the line between the two cells; -The friction has been adjusted to incorporate the various solid surfaces, as well as to model the difficulty cells experience in R 3 to move through cavities, the friction coefficientμ in Vermolen and Gefen (2012) has been modified to where E k stands for the elasticity modulus in each phase.
Therewith the strain energy density associated with longdistance communication is written as where in which F i represents a mechanical force exerted by viable cells, which we treat as a constant in the present study. Further, R denotes the radius of the cells. The local elasticity modulus, E S (x i ), is determined by a mixing rule, which reads as Note that as which models that if the cells are not able to adhere to a solid, then they will not generate any force, and thus, no strain energy density is generated. Further, F i denotes the cellular traction force exerted by cell i. The attenuation of the signal over the domain surrounding cell j towards cell i is modelled by where E S (x i , x j ) represents the averaged elasticity modulus between the two communicating cells, computed by Further, E C represents the elasticity modulus of the cell. Mechanical contact between cells is also adjusted to the three-dimensional case using Hertz contact mechanics for two spheres, which gives where we refer to Fig. 6 for a schematic. Finally, the overall strain energy density is computed through The determination of the translation unit vector z i is analogous to Vermolen and Gefen (2012), except for the abovementioned adaptations. The second term in Eq. (5) mimics the contribution as a result of chemotaxis in the direction of VEGF. Here, the concentration of VEGF is denoted by c V . The γ -function incorporates the mobility of a cell, as well as the resistance by the material it has to move through, and the force that the cell is able to exert on the material. To get a dimensionally consistent relation, see Bookholt (2015) for the details, we assert where Here ρ k are the densities of the separate phases fibrin matrix, basement membrane and extracellular fluid. As for the mechanical component, there is no chemotactic movement f S = 0, that is if the cell is in an environment entirely filled with extracellular fluid. The third term in Eq. (5) takes into account the migration as a result of durotaxis. The λ-function, λ : N×R + → R, is constructed such that it is zero (no durotactic movement) in the extracellular liquid (that is f S = 0), in hard solid (that is f S = 1) and on the the interface (t) (that is f S = 1 2 ), to this extent and after some normalisation and taking into account the hypotheses that migration is inversely proportional to local stiffness, cell motility and cell viability. In Vermolen and Gefen (2012), cell viability is directly coupled to the forces that they exert, and since the dimension of λ should be m 2 s −1 , we use where β i denotes the motility coefficient of cell i and wherê λ denotes the adhesive scaling factor, which is treated as a constant. The last term of the right-hand side of Eq. (5) models random walk of the cells (diffusion), where D denotes the diffusion coefficient of the endothelial cells. The vector dW(t) has three entries that are all independent normally distributed events with zero mean and a variance of dt, that is, the entries of dW(t) are Wiener Processes dW k ∼ N (0, dt), for k ∈ {1, 2, 3}. The tip cells secrete the protein delta-like ligand 4 (DLL4), which is the chemotactic signal of the stalk cells. Hence, for the stalk cells, with centre location x i , we have Here no distinction has been made whether cells are stalk or tip cells, since they are both endothelial cells, except for the chemotaxis term, which is driven by the gradient of the concentration of DLL4, denoted by c D , where DLL4 is secreted by the tip cells. Hence by the secretion of DLL4, the tip cells make the stalk cells follow them, see Fig. 5 for a sketch of the mechanism. The stalk cells can become tip cells and, vice versa, the tip cells may become stalk cells. The transitions between these states are modelled as memoryless stochastic processes, which are classified as follows: Let P(t > τ) denote the probability that the transition does not take place until time τ and let P(t > θ + τ | t > θ) denote the probability that the transition does not take place until time τ + θ , given that the observer is at time θ where the transition has not yet taken place (that is the event did not take place before time θ ), then the memoryless property is defined by The random process is memoryless in times t > 0 if and only if The probability that the stalk cells become tip cells is modelled by an exponential distribution, given by the following probability density function for t > θ: for reverse transition, we analogously have With these probability density functions, we get the following transition probabilities for t > θ In the above equation, λ Q→P is a probability rate constant.
In the experimental case with large number of cells in many experimental experimental samplings, one could measure the amounts of tip cells and stalk cells, and from these figures, one can estimate the probabilities that a cell is either in the "tip state" (that is i ∈ T (t)) or in the "stalk state" (that is i ∈ S(t)). To this extent, Bayes' theorem applied to the longtime observations gives where n and n T (t), respectively, denote the total number of endothelial cells and the number of tip cells at time t. The above relation gives an estimate of how the probability rates λ S→T and λ T →S are related. In the simulations that we will show, the probability rates depend on the chemical environment in which transitions between the two states are favoured if the VEGF concentration is high and if the DDL4 concentration is low. Some phenomenological relations have been used in this study. We finally note that the present modelling does not incorporate cell death or cell proliferation. In "Appendix", the reader will find more details regarding the input values used in this study.

The proteins involved
Next we treat the concentrations of the various proteins VEGF (c V ), DLL4 (c D ), matrix metalloproteinase (MMP) (c M ) and urokinase plasminogen activator (uPA) (c U ). Note that the VEGF makes the tip cells move and further the DLL4 is secreted by the tip cells, and this chemokine makes the stalk cells follow the tip cells. In the equations, we disregard shrinkage or expansion of the total computational domain that could possibly occur due to mixing processes. All concentrations are modelled by diffusion-reaction processes, where we have for t > 0 where r V is a decay rate constant due to consumption by tip cells, D V is the diffusivity of VEGF depending on the phase (fibrin matrix, basement membrane or extracellular fluid), and δ(.) represents the Dirac delta distribution, which is defined by For all the chemokines, there is flux normal (perpendicular) to the boundary, hence for t > 0 The initial condition for the VEGF concentration is given by DLL4 is regenerated from conversion of VEGF by the tip cells, and it is consumed by the stalk cells, to this extent, we have for t > 0 Here r D and s D are regeneration and consumption rates. Initially, there is assumed to be no DLL4 in the domain of computation. The metalloprotease MMP is secreted by the tip cells by conversion from VEGF, and this chemical breaks down the basement membrane, for t > 0, and we have where the first term in the right-hand side models the decay of the MMP concentration as a result of the breakdown of the basement membrane, which is a somewhat stiffer extension of the fibrin matrix. This breakdown enables the Dermal ECs to migrate into the fibrin matrix (including the basement membrane). Further, r M and s M , respectively, are rate constants for decay and regeneration of MMPs. Initially, the concentration of MMPs is zero at all locations of the computational domain. Finally, the protein uPA breaks down the fibrin matrix. This protein is also secreted by the tip cells, and hence, we have for t > 0 also here the first term of the right-hand side models decay of uPA due to the breakdown of fibrin matrix. Furthermore, r U and s U , respectively, are decay and regeneration rate constants regarding uPA. Initially, there is no uPA in Ω. The diffusivities are modelled using a mixing rule: where D 0 p denotes the diffusivity of protein p ∈ {V, D, M, U } and D k for k ∈ {F, B, E} (fibrin matrix, basement membrane, extracellular fluid) denotes the diffusion factor corrected for the phase that is considered.
Since the proteins MMP and uPA, respectively, change the basement membrane and fibrin matrix into extracellular fluid, we have for t > 0 The above relations are consistent with the requirement that the sum over all volume fractions should be equal to one. It is noted that the above equations warrant that if there would be a sufficient number of cells that the long-time behaviour becomes f B , f F → 0, as well as f E → 1, along with c V , c D , c M , c U → 0 in Ω as t → ∞, which indicates stability of the system.

The numerical method
To solve the stochastic differential equations for the spatial positions of the centres of the cells, the Euler-Maruyama method is used (Steele 2001). Further, the diffusion-reaction equations for the concentrations of the chemicals have been solved using the finite-element method in three spatial dimensions. To this extent, we give the weak (variational) formulation, where Sobolev/Bochner spaces are omitted, of the diffusion-reaction equation for VEGF as an example: Approximating where ϕ k (x) are a set of chosen basis functions (in the present study piecewise linear), and taking φ i (x) for ϕ(x), necessitates determining whether the cell centre is located in the tetrahedral element of consideration. In order to determine this, we consider a tetrahedron, e with vertices x 1 , x 2 , x 3 and x 4 . We use the barycentric coordinates of the tetrahedral element. Consider tetrahedral element e, and let ψ i (x) be the linear function that is characterised by where x j represents the vertices of e, and δ i j denotes the Kronecker delta function. Note that and only if x ∈ e and that outside the tetrahedron e the function ψ can assume values beyond the interval [0, 1]. Further, the cell centre with coordinates x j (t) is located within e if and only if 0 ≤ ψ(x j (t)) ≤ 1. An alternative treatment can be applied if the ordering of the vertices of the tetrahedral element e has been carried out such that the numbering over each face of the tetrahedron is in the positive orientation. Let Δ be the determinant given by then this determinant represents six times the volume of the tetrahedron since it is positive by the choice of the orientation. Further, we introduce the following auxiliary determinants that are constructed on the same principle, but now we replace the coordinates of vertex x k with the coordinates of the cell centre x j (t), to get Δ k . For instance, Δ 1 is given by This should be done for all vertices of the tetrahedron e. If all Δ k · Δ > 0 for all k ∈ {1, . . . , 4}, then x k (t) lies within the tetrahedron e. The finite-element implementation has been done for MATLAB, where the meshes have been constructed using the iso2mesh (Fang and Boas 2009). The meshed domain is shown in Fig. 7. The iso2mesh-package generates meshes without actually seeking for an optimal bandwidth of the discretisation matrices, and to this extent, the Cuthill-McKee algorithm has been used to optimise the bandwidth. The finite-element method uses linear elements; hence, linear tetrahedra and the mass matrix needed in the time derivative are lumped through Newton-Cotes integration to prevent spurious oscillations that could even occur when implicit methods are used. Time integration of the partial differential equations is based on a first-order IMEX scheme where the diffusion operator is evaluated at the new time step, and all the nonlinear terms at the previous time step, meaning that the right-hand side and the gradient of the concentration are all treated at the new time step and that the determination of the diffusion coefficients has been performed at the previous time step. The main advantage is that the numerical stability of the time integration is not determined by the mesh size and despite this feature, one does not have to solve a complete nonlinear problem at each time step, and hence, no inner iterations are needed, and therewith, the time integration is relatively cheap. The choice of basis simple tetrahedral elements is justified because no advection terms are to be discretised, and hence, no SUPG discretisation or flux corrections are needed to suppress spurious oscillations. Higher-order finite-element methods will not improve the accuracy because of the used point sources in the partial differential equations for the concentrations (since c p / ∈ H 2 (Ω) for a fixed t > 0, p ∈ {V, U, D, M}). If one aims at improving the efficiency and accuracy, then adaptive finite-element mesh strategies could be helpful, where the iso-concentration surface, implicitly defined by f S = 0.5, represents the moving interface between the fluid and solid phases.
Finally, it is noted that the numerical solution of the volume fractions is performed using Euler backward method with the concentrations that were computed earlier. One  Fig. 7 The three-dimensional finite-element mesh used in the current study. The mesh was constructed using the isomesh-mesh generator (Fang and Boas 2009) could also interchange the order: first compute the volume fractions with the use of concentrations at the previous time step and subsequently use a complete implicit Euler time integration for the numerical time integration of the concentrations. This variant has not been studied in the present study since our approach gave satisfactory results. Another alternative time-integration method is the fully coupled implicit approach, which needs an inner iteration loop within each time step. This lastmentioned approach is thought to be more expensive, and therefore, it has not been applied either.

Simulation results
First we show the visualisation of the simulation in terms of field plots and cell plots. Subsequently, we show results in terms of quantitative measures as well as a comparison with experimental outcomes. We further carry out a sensitivity analysis on the simulations. The default input data have been listed in "Appendix".

Visualisation of the simulations
Initially, the cells are located on the basement membrane, and the initial volume fractions of the three phases, fibrin matrix, basement membrane and extracellular fluid are shown in Fig. 8, where the three-dimensional nature of the problem necessitates to represent the phases in a slice plot. We used slices that are perpendicular to the coordinate axes. In Fig. 8, left, the initial fibrin fraction is shown. It can be seen that the initial fibrin is distributed on the bottom. In the middle of Fig. 8, the initial fraction of basement membrane is plotted on the same slices in the domain of computation. It can be seen that the initial basement membrane is localised on a horizontal layer adjacent to the extracellular fluid and fibrin. The initial fraction of extracellular fluid is localised on positions above the membrane, as can be seen in Fig. 8 on the right. On the bottom of Fig. 8, we show the histograms of the volume fractions that are experienced by the cells. It can be seen that initially these histograms give a polarised behaviour reflecting that cells are located either in the fibrin or in the extracellular fluid to a lesser extent. The cells are modelled to migrate through chemotaxis, durotaxis (adhesion), contact mechanics and random walk. The contact mechanics prevents cells from coinciding with one another. Further, as an example, we show the positions of the cells after 3354 s in Fig. 9, where the green and red cells are the stalk and tip cells, respectively. On the left, the three-dimensional representation is shown, whereas on the right some projections are shown so that it is clearer to see how the cells fit in the channels through the fibrin matrix. It can be seen in this figure that the tip cells are localised on positions closest to the On the bottom, histograms are given of the volume fractions of fibrin, basement membrane and extracellular fluid as experienced by the cells bottom since they indeed take the lead in chemically creating holes in the basement membrane and fibrin. The tip cells are also shown in the projections on the right, where they are represented by the red crosses. In Fig. 10, on the left side, the level surface of f S = 0.5 is shown, which is the surface on which we let the cells adhere to. This figure shows how the channels have been formed by the cells through the release of MMP and uPA that convert the basement membrane and fibrin matrix into extracellular liquid. Hence, this iso-surface of f S = 0.5 shows a time instant in the time evolution of the region that is occupied by the solid phases, and it can be used to visualise the evolution of angiogenesis if different time frames are shown after one other. Next to this figure, we show the positions of the cells that are residing on the top surface of the basement membrane, that is in the vicinity of the initial interface between the extracellular liquid and the basement membrane in Fig. 11. Like in the experimental setting, we show cells that are within a layer of three cell diameters around the initial position of the interface between the basement membrane and the extra-cellular fluid. The top, middle and lower layer cells have been plotted in red, blue and green, respectively, in Fig. 11. On the projections, the tip cells are represented by the red crosses. In Fig. 11, several gaps arise (see the white regions surrounded by the black lines and where we indicated two of them by arrows, see also Fig. 4 in the in vitro experiments). These gaps coincide with sprouts formed by the tip cells. We compute their areas and compare these areas that have been computed for the same setting in the in vitro experiments. The concentrations of all the chemicals are obtained by the finite-element approximation of the solution to the three-dimensional diffusion-reaction equations in which the cells either consume or regenerate the chemicals at their spatial positions through point sources or sinks. Sometimes, one observes some cells that are located within the circumferences of sprouts in Fig. 11, such as in the white patch on the top left. The probable reason for this observation is that these cells just detach from their neighbours that are still on the top surface and start migrating downwards into the sprout that is being formed. Note that here this sprout is very premature The "camera" is in an angle slightly lower than the x, y plane We also see tip cells (denoted by red x-marking) that have not formed a sprout. The iso-lines are calculated in the surface at z = 1039.4 directly beneath the initial placement of the cells. The total number of tip cells is n t = 9. Two sprouts have been indicated by arrows and that here the F S = 1/2 level curve has not formed yet on the top layer of the solid, which means that the summed solid fractions are still above 0.5 there. Since the sprout is occupied with extracellular fluid, the migration of the cell does not proceed instantaneously, and hence, they remain visible on the top (though they are not located on the interface between the solid and fluid phases) for a while. The concentrations of the VEGF, DLL4, uPA and MMP after 3766 s are shown in Fig. 12 in terms of slice plots on the planes perpendicular to the coordinate axes, where it can be seen that the gradients are largest in the vicinity of f S = 1/2 where the tip cells chemically create holes in the solid phases. On the bottom of the figure, histograms of the concentration that all the cells experience have been shown. It can be seen that at the top of the fibrin matrix, the VEGF concentration has decreased a bit as a result of consumption by the tip cells, whereas further away to the bottom, the concentration has not changed considerably. The concentration of DLL4 increases slightly in the vicinity of the top surface of the fibrin matrix since this ligand is secreted by the tip cells that are predominantly at the interface between the extracellular fluid and the solid substances. This ligand makes the stalk cells follow the tip cells. A similar behaviour is observed for the uPA and MMP concentrations, which, respectively, are responsible for the degradation of the fibrin matrix and basement membrane. Next to the concentration plots, we plot histograms of the values of all the concentrations that are experienced by the cells, which determine to what extent the stalk cells will follow the tip cells. Many cells have not moved yet and at their positions, the concentration of all proteins is between zero and very low values. The other cells that did move (tip cells and stalk cells with positions initially adjacent to the tip cells) are entering the region where VEGF has higher values. Therewith there is a considerable portion of cells that experience high values of VEGF. Furthermore, in the course of time VEGF stimulates secretion of the other proteins by tip cells that make the stalk cells follow them and that facilitate the degradation of the fibrin matrix and basement membrane.

Quantitative measures from the simulations
In order to quantify angiogenesis, several measures have been introduced. The first measure is the total area of the sprouts on the initial top layer of the basement membrane divided by the total basement membrane area. This measure is denoted by A(t). The second measure is called the sprouted perimeter, which is obtained by computing the total perimenter of the sprouts on the top of the basement membrane multiplied by the cell diameter and subsequently divide this result by the total area of the basement membrane. We denote this measure by P(t). The third measure is the total volume fraction of the sprouts computed by the integral over f E over the initial fibrin matrix domain. The lastmentioned measure is denoted by V (t). Since the model contains a stochastic nature, through migration and differentiation, it is important to determine the amount of uncertainty for the set of parame-  Fig. 13, where the mean curves for the respective quantities have been plotted over time as well as the 95 % interval of confidence (exceeding probability of p = 0.05) for all these quantities using 12 runs with identical choices for all parameters. It is clear that all the measures go up as the vascularisation process continues. Further, there is an incubation time, which is a result of the following sub-processes:

-Endothelial cells become tip cells by a random selection
in the model and only after a short, nonzero length period, the first tip cells appear; -The VEGF concentration has to reach the tip cells in order to be able to degrade the basement membrane and fibrin matrix; -The (tip) cells can only migrate quickly to the bottom provided the basement membrane and fibrin matrix have decayed.
Besides the incubation time, a small jump (at t ≈ 14000) in the plot for the percentage of the sprouted area is observed. After having examined the video for this simulation, it turned out that two sprouts merged and that the routine to compute the total area of the sprouts counted the merged sprouts twice. This small jump is not visible in the computation of the volume fraction of the vessels. Simulations have been done with different adhesion valuesλ in the durotaxis term, and the dependence did not seem to be significant since the behaviour was not monotonic and the variations were not larger than the the variations in different runs with identical input parameters. The input values forλ ranged between 6 and 100. Therefore, these results are not shown here. Possibly using lower values could show more dependency but this is not certain since the adhesion is not the main trigger for the vascularisation process. Lowering theλ-value only allows cells to move towards the f S = 0.5 in a slower manner and the time-integration method would allow for less overshoot. Furthermore, the initial VEGF concentration has been varied and higher initial VEGF concentrations predict a higher volume of vascularisation, see the bottom of Fig. 14.
Note that we only plot the mean of all the 12 runs that were carried out. Furthermore, the other measures were computed over time for the various initial VEGF concentrations, and the monotonic behaviour was also observed except for the Fig. 13 The estimators for the mean area A(t) and perimeter P(t) of all sprouts at h F as a percentage of the total area and the percentage of the total area and the percentage degraded substrate V (t) and 95 % confidence intervals based on 12 runs with identical parameters highest VEGF concentrations. This deviation is attributed to the fact that higher initial concentrations give a larger concentration gradient. Since the total movement of the cells is determined by contact forces, protein signals, durotaxis (for the adherence) and random walk, the chemotactic and haptotactic components to the overall movement of the cells increases as the initial concentration increases. Thereby the path that the cells follow towards the bottom will be more according to a straight line, and hence, the cross sectional will be lower than in the case that the movement of the cell is distorted more from all the other signals such as contact forces and random walk. In Fig. 15, we show the influence Fig. 14 The area A(t) and perimeter P(t) of all sprouts at h F as a percentage of the total area and the percentage degraded substrate V (t) for varying initial concentrations VEGF of the variation of the probability that stalk cells become tip cells. It can be seen that the amount of vascularisation in all the measures increases monotonically with the probability of stalk cells to become tip cells. Finally, we show the influence of the regeneration of uPA and MMP's in the evolution of the vascularisation over time in Fig. 16. Since both concentrations act in a similar way, we took the regeneration constants equal. It can also be seen that the vascularisation rate monotonically increases with the regeneration constants. The decrease in speed at the latest times is due to flattening out of the VEGF signal due to diffusion, which is the main trigger for the further growth of the sprouts. This is also Fig. 15 The area A(t) and perimeter P(t) of all sprouts at h F as a percentage of the total area and the percentage degraded substrate V (t) for varying maximal sprout to tip probability λ PT confirmed by the absence of a dependence of the sprouting dynamics upon changing theλ-parameter.

Comparison with experiments
Finally, we compare the simulation results to the outcomes obtained from the experiments in terms of the metric P(t), which accounts for the perimeter over time. The results can be seen in Fig. 17, where human dermal tissue (DTECs) endothelial cells are considered. In both the in vitro and simulation experiments, the concentration of VEGF was varied and the sprout perimeter P(t) was computed after the third day in the in vitro experiments and after 14,400 s (4 h) in the Fig. 16 The area A(t) and perimeter P(t) of all sprouts at h F as a percentage of the total area and the percentage degraded substrate V (t) for varying MMP and uPA sourcing rates simulations. This discrepancy in times was caused because the right parameter values were not yet available. Firstly, it is noted that the trends of increasing vascular perimeter are observed in both simulations for increasing values of VEGF concentration. However, in the simulations it can be seen that for the largest value of initial VEGF concentration, the amount of sprouting seems to drop. This drop is attributed to the migration mechanism of the endothelial cells: for larger values of the initial concentration in the fibrin matrix region, the magnitude of the VEGF concentration is larger. Thereby the chemotaxis movement becomes larger. Since for all the cases the other mechanisms that contribute to migration (cell-cell contact, random walk and durotaxis) are approximately the same, the relative portion of chemotaxis is larger if the initial VEGF concentration in the fibrin matrix is larger. Therefore, the paths and holes that the cells will make through the solid will be more straight with fewer migrational components perpendicular to normal of the initial interface between solid and fluid. Hence, for larger initial concentrations, the chemotaxis component dominates and the cells will migrate straight towards the bottom, by which the small vessels will have a smaller diameter and a smaller perimeter. If the concentration of VEGF would be much larger in the experimental case, then one possibly observes the same behaviour for the relation between P and the initial stimulation with the VEGF concentration, since results of in vitro experiments often reveal a bell-shaped curve in the response of cells to increasing concentrations of a stimulus. Probably it is a matter of adjusting the parameters to more appropriate values to have a drop in the sprouting percentage at a higher concentration. This behaviour cannot be attributed to a chemical saturation effect. Furthermore, adjusting the parameter values will also lead to a decrease of the discrepancy between the times at which we determined the sprouting percentages in Fig. 17. A more efficient implementation of the present model will be needed if one wants to match the in vitro experiments to the simulations through inverse modelling. This inverse modelling will make a better fit between simulations and the in vitro experiments possible. Since the aim of the current paper is to introduce the mathematical model which is original in its kind, being a cell-based model for angiogenesis, this inverse modelling is omitted here.

Discussion and conclusions
In this section, we discuss the model and give various recommendations for further study, and the final conclusions are drawn.

Discussion
We based most parameters and scaling factors in the derivation of the cell movement model on physical or biological principles. The only parameter that forms the exception is the dimensionless scaling factorλ. This parameter mediates the contiguity of the monolayer of cells; however, in our sensitivity analysis we have seen that this parameter does not play an important role in angiogenesis. The addition of natural protein decay over time would be an interesting property to add to the system of PDEs. We have performed a sensitivity analysis on five parameters, and we did twelve model runs for the estimator of the mean of all sprouting metrics. Extra computational power or parallelisation could reduce the cost of the simulations so that these numbers can be improved. Such a parallelisation has recently been carried by Woods et al. (2014). Another interesting approach is to apply a full continuum model in terms of a system of partial differential equations like by Maggelakis (2003), Maggelakis (2004) or by Gaffney et al. (2002). The latter approach allows to consider angiogenesis on a larger scale such as on a tissue scale. It would be of great interest to apply some of the homogenisation techniques that, for instance, are currently applied in porous media applications to link the microscopic, cellular, scale to the macro, tissue scale. One could also use an optimisation in the cheaper continuum modelling approach as a start for the cell-based model. The optimisation could be refined in a cell-based formalism, which in total gives a hybrid (multiscale) regression technique.
Obtaining appropriate values for the physical parameters is often very hard and in many cases impossible in complicated models. An example is the diffusive speed D V . For most of the parameters, we used values that were reported in the literature, and if no values in literature were reported, then we used educated guesses so that the model predictions are reasonable. The difference of the VEGF concentration between the in vitro study and the initial concentrations used in the simulations was compensated for by scaling down the chemotactic response accordingly. A modification of these values will not change the model outcomes significantly. As mentioned earlier, one could use sound inverse modelling techniques to get access to modified parameter values which could reduce the gap between experimental and simulation outcomes. At this stage, we are satisfied with having developed a new sound hybrid cell-based model for angiogenesis.
In the current model, the edges of the domain do not exert any forces on the cells, making it possible for cells to move out of the computational domain. Cells that have moved out of the domain are not within an element of the FEM mesh and therefore cannot sense or source any proteins or react to the substrate properties. This fact on itself forms no problem for the rest of the computational model, but does require much useless computational effort, predominantly in trying to find the (non-existing) element a cell is located in. An improvement would be to remove these cells from the computations. Another approach could be to give the boundaries contact mechanical properties or to lay a monolayer of ghost cells on the boundaries that provide the contact mechanical forces to keep the cells in the problem domain.
In the formulation of the rates λ S→T and λ S→T of tip cell selection, we normalise the VEGF and DLL4 concentrations with the initial concentration VEGF c 0 v . In hindsight, this is no reasonable assumption since tip cell selection does not depend on the absolute VEGF concentration, but rather on a saturation with respect to the initial condition. It would be an improvement to remove this normalising factor and reconsider these rates. It would be very interesting to see whether we can verify more of our simulation results with laboratory measurements. This incorporates measurement of other metrics from the in vitro experimental results than only P(t) and maybe 3D visualisation of the sprouting assay using multi-focal plane microscopy or other techniques.
In the future, we want to combine the current angiogenesis model with modelling of cancer development, such as in Vermolen (2015), where the necrotic core of the tumour releases growth factors (Tumour Necrosis Growth factor) that trigger the angiogenesis response of the endothelial cells.

Conclusions
Our model is qualitatively successful in describing the in vitro angiogenesis sprouting assay as performed by the VUmc dermatology department. We modelled the degrading of the substrate by proteases secreted by ECs as a continuous process dependent on the properties of the substrate itself. Cell motility is modelled using a cell-based formalism based on mechano-biological principles that are well established in cell biology. A probabilistic model based on local chemical conditions is proposed to model the differentiation of ECs into tip cells and stalk cells.
The proposed metrics of the amount of sprouting seem to align with the in vitro results on a qualitative level. Quantitative comparison is hard due to many uncertainties, both in the proposed computational model and in the measurement techniques used for the in vitro experiment. The morphology of the sprouts is similar to the experimental setting.
The metrics over time produced by the model respond to variation in parameters as we would expect from biological reasoning. Only the variation of the VEGF concentration is performed in the laboratory setting, and the results are comparable. The area of sprouted perimeter ranges from 1 to 5 % in the in vitro experiments and ranges from 0.5 to 6 % in the simulations for varying concentrations VEGF, which are of the same order of magnitude.
VEGF concentrations, protease secretion rates and the probabilistic model for tip cell selection are important factors in sprout formation process. We postulate that it is the tip cell's ability to degrade the substrate in its surroundings that drives the success of producing a viable sprout as well as give rise to the sprouts proliferation speed and its final depth. This factor is at least as important as the chemotactic response to a higher concentration VEGF or the adhesive properties. The success rate of sprout formation for a tip cell is between 50 and 60 %, independent of the number of tip cells present. We postulate that also in vitro the number of tip cells is larger than the number of sprouts.
The model was constructed to simulate angiogenesis. Vasculogenesis is another process witnessed in studies containing ECs as described by Nany et al. (2004) and Merks et al. (2004). Since our formalism describes EC behaviour in a general sense, we also witness vasculogenesis-like structure formation for varying values of the substrate elasticity. Since modelling vasculogenesis was not the scope of this study, we leave further investigation of this phenomenon for future research.
programme Coupling Geometric Particle Differential Equations at the Isaak Newton Institute.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecomm ons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Appendix: Input data
ECs can change from stalk cells into tip cells (and reversibly) based on the local protein balance. Higher concentrations of VEGF will increase the chance of a stalk cell to become a tip cell. However, a nearby tip cell that sources DLL4 inhibits a stalk cells from becoming tip cells.
We model the changing time of a stalk cell to become a tip cell using an exponential distribution with a rateλ ST . We also model the changing time of a tip cell to change back into a stalk cell using an exponential distribution with rateλ T S . The probability density function (PDF) of the exponential distribution is given by The changing times are then distributed according to exponential distributions. The ratesλ ST andλ T S are dependent on cell-specific factors and local chemical conditions. Since the concentration VEGF c V only decreases over time, we know that we have c max V = c 0 V and we normalise using c V c max V . The concentration DLL4 c D can take arbitrary positive values but in practice almost never exceeds the concentrations VEGF so again we "normalise" by setting c D c max V . For ideal circumstances, i.e. c V = c max V and c D = 0, we define the maximal rate of changing from stalk to tip per second to be p ST . The rateλ ST increases for increasing c V and decreases for increasing c D . We model The strength of the DLL4 induced inhibition is governed by the parameter p i (i for inhibition) and higher values of p i constitute stronger inhibition. The strength of the stimulation by VEGF is governed by the parameter p s (s for stimulation), and lower values of p s constitute stronger stimulation even for low values of c V . Tip cells can change back into stalk cells at a rate governed by the decreasing of a VEGF gradient. We model and set p T S = 1 × 10 −6 . The rates in Eqs. 37 and 38 are per second. During a time step of length t, we see that we have the following probability of a cell type change, where we approximate the exponential with its first-order Taylor This approximation is suitable for small enough time steps t, and we have to specifically ensure thatλ T S t ≤ p T S t ≤ 1. Bayes theorem now dictates that P(S|T ) P(T |S) = P(S) and we end up with an equilibrium ratio between tip cells and stalk cells given the local chemical conditions. The values for p m , p s and p i should be found based on experimental data. Microscopic images show us that for large concentrations of VEGF (25 ng mL ), we see around one sprout for every fifty ECs after two days. We assume that every sprout is led by one tip cell. Furthermore, we assume that all tip cells are selected at 12 h after stimulation with VEGF and that no additional tip cells are formed due to inhibition by the existing ones afterwards. This means that a cell has a probability of 0.02 per 12 h of becoming a tip cell not changing back again. This is equivalent to 0.02 12×3600 = 4.6 × 10 −7 per second. The VEGF concentration is not maximal over the entire domain, so we choose p m = 4 × 10 −6 . We know that the inhibiting effect of DLL4 is stronger than the stimulating effect by VEGF, so we model p i = 10 and p s = 4 and conclude that these values are suitable for our modelling purposes. Better estimations for these values might be obtained by fitting the model results to the experimental observations, but this is cumbersome due to the long execution time of the model. An investigation of the magnitude of these parameters reasoning from a more biochemical point of view is a useful recommendation.

Parameters and Domain
Our model uses a number of parameters, most of a physical or chemical nature. We try to find accurate values for these parameters in the literature as much as possible although some values are hard to estimate. We list all our parameters with description, symbol, programming code name, value, dimension and where possible a source in Tables 2, 3, 4 and 5. The elastic modulus of endothelial cells E c is approximately 10 kPa according to Kuznetsova et al. (2007). They found this value in a study using atomic force microscopy probing.
According to the work of Ganz et al. (2006), traction forces F i are also higher than the 1 nN used by Vermolen and Gefen (2012). Reinhart-King et al. (2003) conducted in vitro studies measuring the forces exerted by ECs on    Streeter et al. (1998) polyacrylamide substrates. On page 1578, we see see the relation between cell area and exerted force of ECs. Since we model our cell radius as R = 22.5 µm, we have an area of 22.5 2 π = 1590 µm 2 = 1.59×10 −5 cm 2 . The graph gives a force of 0.1dyne = 10 −6 N = 1000 nN, and we choose to use this value in our model. Elasticity of fibrin matrix E F is approximately 10 kPa according to Rowe et al. (2007). Zhu et al. (2011) propose values for the elastic modulus of collagen-chitosan scaffolds in the order of 10 kPa. Since BM consist of collagen and more stiff components, we use an estimate E B = 20 kPa. Plank et al. (2002) set the diffusion coefficient for VEGF D V = 3.6 × 10 −3 mm 2 h = 1 µm 2 s (in matrigel). We have found no references on the diffusion coefficients of the concentrations DLL4, uPA and MMP so we estimate them using their molecular weights. We assume that substances diffuse more slowly for larger molecular weights. All four substances are sold commercially, and molecular weights m are specified very accurately. We see that m V ≈ 38.2 kDa, m D ≈ 75 kDa, m U ≈ 31 kDa and m M ≈ 72 kDa. We estimate  (2009) measured in an in vitro experiments that the VEGF diffusion constant in matrigel (similar to fibrin matrix) is equal to 278 µm 2 s . The coefficient used by Plank et al. (2002) seems to make the most physical sense and produces reasonable results in our computational model. Further research into the difference between these values might be useful. Plank et al. (2002) furthermore set the VEGF uptake rate r V within cells to 8.66 × 10 −5 mm 2 h = 0.024 µm 2 s . Note that Plank et al. work in a 2D setting and that in the 3D setting, the dimension would be mm 3 s . We set the uptake rates r D , r U and r M equal to this value since we could not find any references concerning these quantities. The degrading rate r F of the fibrin matrix fraction by uPA is, according to Lutolf et al. (2003), of magnitude 1.21 s −1 , and we use this value in our model. We estimate that the degrading rate r B of the BM by MMP is of the same magnitude.
Viscosity of water is 0.6531 × 10 −3 Ns m 2 at 40 degrees centigrade and 0.7978 × 10 −3 Ns m 2 at 30 degrees centigrade according to Streeter et al. (1998). Our experiment is conducted in an incubator at 37 degrees, and we therefore take μ E = 0.6965 × 10 −3 Ns m 2 . Viscosity fibrin matrix is approximately 7.000 × 10 −3 Ns m 2 according to Ehrlich et al. (1952). According to Urbanchek et al. (2001), the density of myocytic cells is 1.060×10 −3 kg cm 3 . Water at 37 degrees centigrade has a density of 0.9933×10 −3 g cm 3 . We model one well in a 96-well plate. We take the same dimensions of diameter 7 mm (i.e. d top = 7000µm) and a height determined by filling the well with v F = 100µL fibrin matrix and 100 µL fluid to a total of v = 200 µL with N = 20.000 suspended EC's. This mimics the laboratory setting one-to-one. To decrease computational time and load on memory, we introduce a linear scaling factor ξ , and we have scaled domain parameters ξ 2 andÑ = N ξ 2 to keep an equal aspect ratio of the domain and a monolayer of cells that is constant in density (cells per µm 2 ).