The Force at the Tip - Modelling Tension and Proliferation in Sprouting Angiogenesis

Sprouting angiogenesis, where new blood vessels grow from pre-existing ones, is a complex process where biochemical and mechanical signals regulate endothelial cell proliferation and movement. Therefore, a mathematical description of sprouting angiogenesis has to take into consideration biological signals as well as relevant physical processes, in particular the mechanical interplay between adjacent endothelial cells and the extracellular microenvironment. In this work, we introduce the first phase-field continuous model of sprouting angiogenesis capable of predicting sprout morphology as a function of the elastic properties of the tissues and the traction forces exerted by the cells. The model is very compact, only consisting of three coupled partial differential equations, and has the clear advantage of a reduced number of parameters. This model allows us to describe sprout growth as a function of the cell-cell adhesion forces and the traction force exerted by the sprout tip cell. In the absence of proliferation, we observe that the sprout either achieves a maximum length or, when the traction and adhesion are very large, it breaks. Endothelial cell proliferation alters significantly sprout morphology, and we explore how different types of endothelial cell proliferation regulation are able to determine the shape of the growing sprout. The largest region in parameter space with well formed long and straight sprouts is obtained always when the proliferation is triggered by endothelial cell strain and its rate grows with angiogenic factor concentration. We conclude that in this scenario the tip cell has the role of creating a tension in the cells that follow its lead. On those first stalk cells, this tension produces strain and/or empty spaces, inevitably triggering cell proliferation. The new cells occupy the space behind the tip, the tension decreases, and the process restarts. Our results highlight the ability of mathematical models to suggest relevant hypotheses with respect to the role of forces in sprouting, hence underlining the necessary collaboration between modelling and molecular biology techniques to improve the current state-of-the-art.


Free energy
In the present model, the configuration of the vascularised tissue is represented by an order parameter φ(r) that distinguishes the two phases of the system, being equal to 1 inside the capillary and -1 outside [1]. To describe the elastic behaviour we include a displacement field u(r) that represents the local displacement with respect to the reference unstressed configuration. Akin to what happens in phase separation of elastic blends [2,3], the order parameter relaxes much slower than the mechanical environment, and therefore, we consider that all the forces are constant in time when describing the evolution of φ. Our strategy starts by defining a free-energy functional of the order parameter and the displacement field. For this specific free energy we will start by obtaining the displacement as a function of the order parameter. We can then calculate the functional derivative of the free energy with respect to the order parameter, i.e. the chemical potential of the system. The evolution of the order parameter depends of this chemical potential in the following manner [1,4]: where M is the mobility coefficient and the last term corresponds to the proliferation of the tissue cells.
The free energy of the system is then the sum of F φ , the energy associated to the system configuration, with the elastic free energy associated to the displacements, F el , and with the free energy related to the forces exerted by the cells in the system F f . Therefore: In this equation ρ φ − aφ 2 + φ 4 4 + ǫ 2 2 (∇φ) 2 is the contribution to the free energy density that is associated to the order parameter. This term is responsible for the interface between the two tissues [1], and leads to a Cahn-Hilliard term [4] in the final equation. In here, ρ φ is proportional to the free energy density of a pure phase, a is a constant and ǫ is the width of the interface, which in phase-field models is also associated to the surface tension. This term guarantees the existence of two minima in the free energy (corresponding to the two phases) and includes an energy penalty for spatial variations of the order parameter, i.e. a surface tension.
The elastic contribution for the feee-energy density is given by 1 2 σ : ε, where the strain tensor elements are We consider the tissues to be homogeneous and isotropic for simplicity, though a similar derivation of the model can be also obtained for other constitutive relations. For isotropic and homogeneous materials the stress tensor elements are given by In this equation K is the compressibility modulus, µ the rigidity modulus, d the dimensionality of the system, and δ ij the Kronecker symbol.
The free energy density associated to the force exerted by the cells, −χ∇ · u, will be discussed below.
Since the displacement field relaxes much faster than the order parameter, we can consider a situation of mechanical equilibrium, i.e. δF δui = 0. Calculating the functional derivative of F el with respect to the displacement we obtain Therefore, since the mechanical equilibrium law states that δF δui = ∂ j σ ij + f i = 0, the force density exerted by the cells has to obey −f i = δF f δui = ∂ i χ. And so, the gradient of the field χ(r) if proportional to the force density. There are two important forces in the system: the traction exerted by the tip cell, f t i = −∂ i χ t , and the force exerted when one cell adheres to its neighbour. This force of adhesion mimics the attraction between neighbouring cells that is conveyed by the adhesion molecules [5], and can be written as f a i = −∂ i χ a = α∂ i φ, where α is a constant. This force points in the direction of increasing φ and has the role of balancing the pulling exerted at the tip during angiogenesis [6]. The corresponding field χ a (r) is then given by χ a (r) = −αφ(r). Notice that while the adhesion χ a (r) is a function of the order parameter, we assume that the traction force exerted by the tip cell is not a function of φ. At the tip, the cell will contract and keep the traction force for some time while the tissue adapts to that force.
During that time, χ t (r) is a constant field.

Evolution of the order parameter
Considering a small difference between the rigidity moduli of the two phases, such that µ = µ 0 − µ 1 φ, with µ 1 ≪ µ 0 , we will expand the derivative of the free energy with respect to the displacement in powers of µ 1 . In the same way, the displacement field u(r) and the free energy function can also be split into components of order zero and components of order one: u(r) = u 0 (r) + u 1 (r) and F [φ(r), w(r)] = F 0 [φ(r), w(r)] + F 1 [φ(r), w(r)].
Collecting the terms of first order, we obtain: where L 0 = K − 2µ0 d + 2µ 0 is a constant and we define a field w such that u 0 (r) = ∇w. The derivative of the whole free energy yields δF δui = 0, giving −L 0 ∂ ijj w + ∂ i χ = 0. This equation is satisfied if Using equation (2), and the definition of w, we can re-write F 0 (1) to get: Repeating the process for the first order term in the free energy we observe that all the terms in this order that are proportional to u 1 (r) cancel out, giving simply Putting all together, we can finally calculate the partial derivative of the free energy with respect to the order parameter: where [∇] −2 is the inverse laplacian operator, and the constant a can be set for convenience such that a − α 2 L0ρ φ = 1, such that the two minima of the order parameter fall at φ = 1 and φ = −1.
Finally we can obtain the equation for the evolution of φ by ∂φ ∂t = M ∇ 2 δF δφ (to which we will add the proliferation term below), yielding

Numerical Implementation
In a square domain of 120 × 120 µm, with a lattice unit of 1 µm, we start with a single vessel of width 10 µm, where φ = 1 inside that vessel and φ = −1 outside. The VEGF initial concentration V (r, 0) is set to have a constant gradient, with V = 1 at the system's right boundary, and V = 0 at its left boundary.
The VEGF concentration in the system evolves according to the following diffusion equation: where D V is the diffusion constant of VEGF, α V is the VEGF consumption rate by the endothelial cells, and Θ(φ) is the Heaviside function. The penetration width of VEGF in the endothelial tissue is given by . This equation is solved by finite differences and with fixed values of V at the left and right boundary conditions.
We start by running the simulation for a sufficient number of time steps until the VEGF field becomes stationary.
Then, at the vessel's mid-height, on its right surface, we will include the traction force field originated by the tip cell action.
The evolution of the order parameter, according to equation (3), depends on the zeroth order displacements through the field w. This field is a function of φ and of the tip cell traction (according to equation (2)), and so we start by obtaining χ t (r) numerically. The tip cell exerts a compressible force [7,8], and, hence, we consider a putative compressible traction force with the direction of the gradient of VEGF applied at the vessel tip, f t 1 . We set the magnitude of this component to decrease exponentially (as a gaussian function) far away from the application point.
In the next figure we plot the magnitude of f t 1 as a function of the position. In the way that the model is designed, the traction is proportional to the gradient of a localised field χ t . That is not the case for the putative traction f t 1 , which only has a component in the direction of the VEGF gradient, implying that the field χ t would be constant and different from zero (and therefore, not localised) along the direction perpendicular to the gradient of VEGF. We correct the traction field f t 1 by adding a small non-constant component in the direction perpendicular to the gradient of VEGF. This small component is calculated indirectly by first obtaining a field χ t (r) from f t 1 by solving the following equation: This equation is solved using finite differences with the boundary condition of χ t = 0 at the boundaries. In this way, the χ t obtained is localised, and the traction given by f t = −∇χ t that we are in fact simulating, is also localised, compressible, and has a small component perpendicular to its main direction of action (see Figure 2A, main article).
For other more complex force distributions (e.g. the traction fields measured experimentally [7,8]), the χ t (r) can be also obtained using equation (4), by replacing f t 1 for the traction field of interest. The field χ t (r) is kept constant for a certain number of time steps t cell while the order parameter evolves. After that time, the capillary wall has moved and the field χ t (r) is again calculated using equation (4), but centred at the new position of the capillary wall.
To simulate the evolution of the order parameter, we solve equation (3) with a term describing proliferation. In this equation w is obtained by solving equation (2) numerically at every time step. We use periodic boundary conditions and set the spatial average of w to be equal to zero (notice that the value of the average of w does not affect the displacement field).
As explained in the main text, the proliferation α p [V, φ, w] at a point is the average of the proliferation in an area the size of a cell, being given by where Ω is the area of a cell (we consider a circle of radius 5 µm), Θ(φ) is the Heaviside function, and p(V, w) is the proliferation function used. In the main article, we consider that the proliferation may depend either on V or on the strain, i.e. on ∇ · u 0 + α = ∇ 2 w + α = α(1 − φ) + χ t . For the first case we characterise the proliferation, i.e. p(V, w) = p(V ), by two constants: the limit VEGF concentration, L V and the maximum proliferation, M P .
Therefore, p(V ) is given by In the definition of p(V ), the maximum cell proliferation M P is multiplied by a factor 2 because the order parameter has to be increased by two units (from φ = −1 to φ = 1) for a new capillary region to form. Equation (5) is the dependence for the proliferation used in Figure 6 of the main text. For Figure 7, the proliferation of equation (5) is multiplied by the Heaviside function Θ(∇ 2 w(r) + α − S m ), where S m is the minimum strain that is required to exist proliferation.
When the proliferation depends on ∇ 2 w ( Figure 4) the proliferation function used is similar to equation (5) but with the VEGF concentration replaced by the strain, and the L V replaced by the limit strain L S . For Figure 6 of the main text we require that the proliferation is only different from zero if there is a minimum concentration of VEGF In Figures 5 and 7 the values used for V m and S m follow the reasoning that a small VEGF concentration (in Figure   5) or a small strain (in Figure 7) is sufficient to trigger endothelial cell proliferation. In this way, the values used for Having w and α p [V, φ, w], equation (3) is finally integrated using a finite different scheme with periodic boundary conditions.

Parameters
The goal of the present article is to introduce a compact and simple to implement phase-field model for sprouting angiogenesis able to describe sprout morphology as a function of the mechanical characteristics of the tissue and cells.
The model is also used to provide insight to the possible factors regulating endothelial cell proliferation in the sprouts.
In the present article we are not attempting to model specific experiments, and therefore we will use rough estimates for the different parameters of the model.
We first present the table of parameters used (denoting the extra-cellular matrix by ECM and the endothelial cells by EC), and below we discuss the choices made when setting these parameter values.  The constants reflecting the mechanical characteristics of the tissues, i.e. the Young's and rigidity moduli, are used to calculate the parameters µ 0 (average of both rigidity moduli), µ 1 (difference between the rigidity moduli), K (average compressibility modulus) and L 0 (given by L 0 = K + 2µ 0 − 2µ0 d , where d = 2 is the dimensionality of the system). In the particular choice of parameters used in this example, µ 1 µ 0 ≪ L 0 . Though the terms of higher order that were disregarded in the derivation of equation (3)  The role of the first term in the free energy, F φ [φ(r)], is the formation of the interface between the capillary and the extracellular matrix. Small variations in the profile of the order parameter across the interface with respect to equilibrium are rapidly corrected by this term. The parameter ρ φ determines how fast the equilibrium shape of the interface is achieved. However, if ρ φ is very large, the dynamics of the system is governed by the Cahn-Hilliard part of equation (3) and not by the elastic component or the endothelial cell traction. Therefore, ρ φ should be large enough for the interface to be formed rapidly, but small enough for sprout extension being determined by the mechanical properties of the tissues and not by the capillary wall surface tension. This balance occurs when the order of magnitude of F φ [φ(r)] is the same as the order of magnitude of F el [u(r)] [2,3], and therefore, we set ρ φ = L 0 .
The value of the mobility M is chosen to match approximately the timescale of sprout growth ex vivo. In Figure   3 of the main text we observe that individual cells are able to migrate alone in the matrix. The maximum distance these cells are able to cover in 24 hours is approximately independent of the initial amount of VEGF in the matrix (provided it is not too large, so that it saturates the cell receptors). All these cells travel at approximately 9.5 µm/hr.
We simulated a single sprouting event, with no proliferation, setting the adhesion to α = 4.7 × 10 2 Pa and maximum traction force 3.0 × 10 3 Pa. With these parameters, and according to Figure 2C from the main document, the tip cell separates from the main vessel and migrates in the direction of increasing VEGF levels. We observed the cells migrate at constant velocity and we altered the time scale of our simulation (given by the mobility M ) so that the cell migration velocity matches the value observed experimentally.
Finally, the VEGF diffusion occurs in a faster timescale than endothelial cell rearrangement and migration. The large value for the VEGF diffusion constant [1,11] guarantees that the VEGF profile achieves the equilibrium for a particular sprout morphology. The value of VEGF consumption α V is chosen so that the VEGF penetration width D V /α V is on the order of the cell diameter [1].

Experimental Model: Matriplug Assay
Male Wistar rats (6 months old) were obtained from local breeding colonies maintained at the animal facility of the Faculty of Medicine -University of Coimbra. The animals were subjected to a constant daily cycle of 12 hours of light and 12 hours dark with constant temperature (22 -24 o C) and humidity (40 -70%). Rats were given free access to water and to standard commercial pellet chow diet (AO4; Panlab, Barcelona, Spain).
To model VEGF-driven in vivo angiogenesis, a subcutaneous matrigel plug assay was performed, as previously Ex vivo study of angiogenesis: live confocal imaging. Live imaging of growing sprouts was performed using zeiss confocal microscopy (633 nm) with 10x of magnification and at stable temperature. The acquisitions were made at different time points of the experience and always in the same position to allow sprout follow.
Statistical analysis. Results are presented as mean ±SEM. One way ANOVA-test was used to determine statistical differences. Post-hoc Tukey test was performed to assess the differences between groups. P < 0.05 was considered significant.

Examples of Single Runs
In this section we include the some graphics of single runs for different parameters.
First we analyse the growth of a single sprout for different values of the proliferation rate. In our model there is a balance between the force of the tip cell and the surface tension term (the first term of equation (1)  In this figure we observe that when there is no proliferation the vessel grows very slowly (proportional to ∼ t 1/2 ).
After 12 hours the tip cell is already moving at a very slow speed: only 10 µm in 8 hours, and it is still slowing down.
Though the vessel does not effectively stop, the velocity decreases so much that in an experiment the cell would end up moving less than its length per day. Compare with what occurs when there is proliferation, where the vessel is able to maintain a constant growing velocity. The proliferation is always able to increase the velocity of the growing sprout.
Then we observe the sprout length for different values of ρ φ . As we increase ρ φ the surface tension increases and the vessel's final length will be shorter. Lower values of ρ φ would make easier for the tip cell to leave the growing sprout. Finally we show the length of the sprout for different values of the ECM's Young's modulus. The ability of endothelial cells to migrate and orient depends on the value of the Young's modulus in a non-monotonous manner [6] (i.e. there is an optimal value for the Young's modulus for network formation). In the figure below we present the growth in time of the position of the front of a growing sprout with no proliferation for four different values of the Young's modulus of the ECM: E = 0.2E 0 , E = 0.6E 0 , E = E 0 and E = 1.5E 0 , with E 0 = 3.0 KPa, the value used in the article. The values for the adhesion and traction are the same used in the previous two figures. High values and low values of the Young's modulus lead to shorter vessels, as expected.