<figure><img src="image_1.png" /></figure>

HAL
open science

# Tumour angiogenesis as a chemo-mechanical surface instability

Chiara Giverso, Pasquale Ciarletta

► To cite this version:

Chiara Giverso, Pasquale Ciarletta. Tumour angiogenesis as a chemo-mechanical surface instability. Scientific Reports, 2016, 6, pp.22610. 10.1038/srep22610. hal-01285078

HAL Id: hal-01285078

https://hal.sorbonne-universite.fr/hal-01285078

Submitted on 8 Mar 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
SCIENTIFIC REPORTS

<figure><img src="image_3.png" /></figure>

OPEN

# Tumour angiogenesis as a chemo-mechanical surface instability

Chiara Giverso<fnref n="1" /> & Pasquale Ciarletta<fnref n="1" />,²

Received: 29 October 2015

Accepted: 17 February 2016

Published: 07 March 2016

The hypoxic conditions within avascular solid tumours may trigger the secretion of chemical factors, which diffuse to the nearby vasculature and promote the formation of new vessels eventually joining the tumour. Mathematical models of this process, known as tumour angiogenesis, have mainly investigated the formation of the new capillary networks using reaction-diffusion equations. Since angiogenesis involves the growth dynamics of the endothelial cells sprouting, we propose in this work an alternative mechanistic approach, developing a surface growth model for studying capillary formation and network dynamics. The model takes into account the proliferation of endothelial cells on the pre-existing capillary surface, coupled with the bulk diffusion of the vascular endothelial growth factor (VEGF). The thermo-dynamical consistency is imposed by means of interfacial and bulk balance laws. Finite element simulations show that both the morphology and the dynamics of the sprouting vessels are controlled by the bulk diffusion of VEGF and the chemo-mechanical and geometric properties at the capillary interface. Similarly to dendritic growth processes, we suggest that the emergence of tree-like vessel structures during tumour angiogenesis may result from the free boundary instability driven by competition between chemical and mechanical phenomena occurring at different length-scales.

Angiogenesis is the complex process by which new blood vessels develop from an existing vasculature in order to supply nutrients and/or metabolites to tissues, playing a fundamental role in many physiological and pathological conditions<fnref n="1" />⁻⁵.

For instance, angiogenesis physiologically occurs during embryogenesis, placenta formation, fetal development⁶ and during tissue-repair⁷. On the other hand, it may drive the onset and the progression of rheumatoid disease⁸, some eye diseases, duodenal ulcers and the initiation and growth of most types of solid tumors⁹,¹⁰.

Since angiogenesis involves the dynamics of the endothelial cells (ECs) forming the blood vessel wall, we investigate whether it is possible to develop a thermodynamically consistent *surface growth model* describing the onset of new vessels, unlike existing mathematical approaches requiring some empirical rules for capillary formation.

Angiogenesis invariably occurs through a well-ordered sequence of events driven by the rearrangement and migration of ECs (i.e. sprouting) forming the lining of the existing vasculature, and the subsequent EC proliferation and fusion, culminating in the formation of a network of new capillaries¹¹,¹².

In this paper, we focus on tumour-induced angiogenesis, the process in which a small avascular tumour (i.e. a colony of cancer cells that lacks its own blood supply) reaches a critical diameter of approximately 2 mm. Above this critical size, the existing vasculature can no longer sustain tumor growth only by means of nutrients and oxygen diffusion⁹. Tumour cell hypoxia is assumed to trigger the first event of tumor-induced angiogenesis, which is the secretion by tumor cells of a number of chemicals⁵, collectively called tumour angiogenic factors (TAFs)¹³. These substances diffuse through the surrounding tissue until they reach the nearby vasculature, whereupon they initiate the degradation of the basal membrane of the capillary wall⁵. Then, ECs respond to the TAF concentration gradient by proliferating near the sprout tip and by chemotactically migrating towards the TAF source, forming protrusions (sprouts)¹⁴⁻¹⁶. Following the initial small finger-like protrusion, branching later occurs (i.e. sprout branching)⁴,¹⁴,¹⁷. The growing branches move towards the tumour, following the motion of the leading EC at the sprout-tip and organize themselves into sort of a dendritic structure⁴,¹⁵⁻¹⁷ inside the extracellular matrix (ECM). Furthermore, as the sprout approaches the tumor, eventually joining the tumour mass, the branches dramatically

<footnote n="1">¹Dipartimento di Matematica - MOX, Politecnico di Milano and Fondazione CEN, Piazza Leonardo da Vinci, 32-20133 Milano, Italy. ²CNRS and Sorbonne Universités, UPMC Univ Paris 06, UMR 7190, Institut Jean le Rond d'Alembert, 4 place Jussieu case 162, 75005 Paris, France. Correspondence and requests for materials should be addressed to P.C. (email: pasquale.ciarletta@polimi.it)</footnote>
increase in number¹⁸. Once the tumour has granted access to the vasculature, it gains access to a virtually endless supply of nutrients, possibly metastasizing in distant sites¹,⁴.

Tumor angiogenesis is an active field of research not only for the biologists but also for the physical and mathematical researchers. In particular, several continuous¹⁹⁻²² and hybrid²³⁻²⁶ models have been proposed in the past 20 years²⁷⁻²⁹. These works focused mainly on the role played by ECs and by the different chemicals, including both those associated with the soluble angiogenic factors secreted by the cancer cells, and eventually those arising from insoluble molecules in the ECM². The modelling of the interactions between the ECs and both the different angiogenic factors and ECM macromolecules typically encapsulates systems of coupled nonlinear partial differential equations (PDEs) describing the migration of ECs from the parent vessel towards the solid tumour. Such PDEs continuous models, first proposed by Balding and McElwain²⁰ and later refined⁴, can describe some important features of angiogenesis occurring at the macroscopic scale, such as average sprout density, average vessel growth rates and network expansion rates¹⁷. However, they are not able to reproduce the morphology of the developing capillary network, since they use a diffuse interface approach. Therefore this kind of model can be used neither to capture the overall dendritic structure of the network nor to evaluate the inner blood flow²⁴. The finer description of these vessel networks can be reproduced with some discrete models that operate at the scale of a single EC¹⁷,³⁰,³¹. Thus some hybrid models combining a continuous representation for the chemicals and discrete elements to track the motion of individual ECs have been proposed³,⁴,²³,²⁵,²⁶,³²,³³. Although discrete/hybrid models have the advantage of describing the motion of individual ECs for simulating a realistic capillary network, empirical rules for branching should be generally defined. Moreover, a proper mechanical description of the interactions occurring between the ECs and the ECM has not been incorporated²³.

Furthermore, both the continuous and the discrete mathematical models in the state-of-the-art focused uniquely on the formation of the capillary sprout network in response to soluble and insoluble chemical stimuli (e.g. TAFs and eventually fibronectin), without giving a proper mechanical representation of the process. Notwithstanding, in addition to chemical stimuli, mechanical cues play a fundamental role in vascular sprouting and maturation, since they govern the interaction between ECs and the surrounding extracellular environment⁴. Indeed, ECs interact with the ECM components, which strongly affect the cellular migration characteristics⁴,¹⁹,²²,³⁴. An extensive description of the key chemical and mechanical processes occurring during angiogenesis can be found in recent experimental works³⁵⁻³⁷.

In this work, we propose an original mechanistic approach to angiogenesis, defining a thermo-dynamically-consistent continuous model of interfacial growth, that takes into account geometrical, physical and chemo-mechanical factors.

Inspired by the striking similarity to tree-like structures found in pure liquid solidification, isothermal solidification of liquid mixtures and oil recovery by fluid injection³⁸,³⁹, we resort to the theory of interfacial growth⁴⁰ and pattern formation in crystal growth³⁹,⁴¹, proposing a continuous model with sharp interface. These theories have largely been applied to inert matter, ranging from the growth of snowflakes to the solidification of metals³⁸, and recently employed for living systems, such as bacterial colonies⁴²,⁴³, cell differentiation during morphogenesis and cancer growth⁴⁴. Compared to non-living systems, biological processes present a significantly higher complexity in the mechanisms of self-organization, even if some general principles can be still applied⁴¹,⁴². All these processes, indeed, are characterized by the occurrence of an interfacial pattern, i.e. the developing structures will form at the interface between two (meta-) stable phases/materials rather than in the bulk of the material. Furthermore, complex patterns result from out-of-equilibrium phenomena dominated by the interplay of driving forces of different physical nature³⁸,³⁹. Our understanding of this mechanism is given by the Mulins-Sekerka instability⁴¹,⁴⁵, which describes pattern formation during the solidification of a solid phase front growing into a supercooled liquid³⁹.

In the following we propose a free-boundary mathematical model of tumor-induced angiogenesis. The model encapsulates the interfacial and bulk balance laws for ECs and extracellular space, respectively, as well as the reaction-diffusion equation for the vascular endothelial growth factor (VEGF), one of the principal tumour angiogenic factors¹⁵,²⁶. Implementing a finite element code, we later present the numerical simulations of the model, together with a sensitivity analysis on the model parameters, which highlights the key factors controlling the vascular morphology. Finally, we critically discuss the numerical results of the proposed model and we add few concluding remarks.

## Mathematical Model

**Description of the biological system model.** Solid tumors can stimulate new vessel formation¹⁰ across distances of some millimeters (from 1–3 mm⁴⁶ to 5 mm⁴⁷), whilst it takes approximately 10–21 days for the growing network to link the tumour to the parent vessel¹⁴,¹⁸,⁴⁷. However, generally, mathematical models generally focus on smaller length-scales (e.g. 17–68 µm³).

The major biological components of the vessel are ECs, forming a mono-layer of flattened and extended units. The abluminal surface of the capillary is then wrapped into the basal lamina, which is a collagenous network composed by laminin, other proteins and carbohydrates. The thickness of the basal lamina is only a fraction of the endothelial layer, so that the typical dimension of the capillary wall (including both the basal lamina and the endothelial layer) is in the order of 190–270 nm, for a capillary lumen spanning from 4–5 µm up to 30–40 µm⁴⁸. Thus all the relevant biological processes occurring during angiogenesis take place within a region $\Delta V_c$ whose characteristic width is much smaller than the capillary-tumor distance and the lumen size (see Fig. 1a).

**Definition of the mathematical model.** Accordingly, let us consider a continuous system made by two different materials occupying the two adjacent regions, $V^+$, representing the extracellular space, occupied by the ECM and healthy tissue, and $V^-$ representing the capillary lumen filled with blood. In the following, we will refer to all properties related to the material in $V^+$ and $V^-$ with the superscripts + and –, respectively (see Fig. 1b).
<figure><img src="image_4.png" /><figcaption>**Figure 1.** (a) Schematic representation of the biological domain consider: the capillary lumen, the extracellular space and the tumor can be represented through the control volumes $V^-(t)$, $V^+(t)$ and $V_i(t)$, respectively, whereas the endothelial layer of the capillary is represented through the non-material interface, $\Sigma(t)$. (b) 2D domain used for the numerical simulations, representing the section of the 3D domain in (a).</figcaption></figure>

Let us consider angiogenesis as an *interfacial growth* process$^{40,44}$, i.e. as the ensemble of phase transformation phenomena occurring between the two contiguous materials in a very narrow volume $\Delta V_e$ across their interface (see Fig. 1):

$$ \Delta V_e = \bigcup (\mathbf{x} + s\mathbf{n}_{\Sigma}) \quad \forall \mathbf{x} \in \Sigma(t), \quad -\varepsilon/2 < s < \varepsilon/2 \quad (1) $$

where $\mathbf{n}_{\Sigma} = \mathbf{n}^- = -\mathbf{n}^+$ is the local outward unit normal vector of the surface, $\varepsilon$ is the small thickness of this interfacial volume and $\Sigma(t)$ is the non-material interface that separates the two regions and that represents the capillary wall.

Thus, we can define the surface fields on the interface $\Sigma(t)$ by homogenizing the volumetric physical properties defined inside the small volume $\Delta V_e^{40}$. In practice, the surface $\Sigma(t)$ behaves as a moving non-material interface, which carries the thermo-chemo-mechanical properties of $\Delta V_e$. Therefore $\Sigma(t)$ can be treated as a discontinuity moving inside the continuous biological system with parametric velocity $\bar{\mathbf{n}}_{\Sigma}$ and associated physical velocity $\mathbf{v}_{\Sigma}$. Then we can assume$^{40}$ that the projection $\mathbf{v}_{\Sigma s}$ on the surface $\Sigma$ of the physical velocity $\mathbf{v}$ is uniform inside the volume $\Delta V_e$, i.e.

$$ \mathbf{v}_{\Sigma s} = (\mathbf{I} - \mathbf{n}_{\Sigma} \otimes \mathbf{n}_{\Sigma}) \cdot \mathbf{v}(\mathbf{x}) = (\mathbf{I} - \mathbf{n}_{\Sigma} \otimes \mathbf{n}_{\Sigma}) \cdot \mathbf{v}(\mathbf{x} + s\mathbf{n}_{\Sigma}) \quad \forall \mathbf{x} \in \Sigma(t), \quad -\varepsilon/2 < s < \varepsilon/2. \quad (2) $$

The tangential component $\mathbf{v}_{\Sigma s} = \bar{\mathbf{n}}_{\Sigma s}$ is dependent on the parametrization, whilst the normal component $\bar{\mathbf{n}}_{\Sigma n}$ is not. Finally, we remark that the decomposition of the physical velocity is given by $\mathbf{v}_{\Sigma} = \mathbf{v}_{\Sigma s} + \mathbf{v}_{\Sigma n}\mathbf{n}_{\Sigma}$, so that $\mathbf{v}_{\Sigma n} = \mathbf{v}_{\Sigma} \cdot \mathbf{n}_{\Sigma}$ is unique.

In this framework, balance laws should be formulated not only for the volume fields but also for the surface fields, thus to guarantee the consistency of the thermo-chemo-mechanical properties of the whole system. Let $\rho^\alpha$, with $\alpha = \{+, -\}$, be the spatial density in the corresponding volumes, we first derive the mass balance laws. Considering that no net proliferation occurs inside the volumes $V^\alpha$, with $\alpha = \{+, -\}$, that the ECM degradation associated to matrix metalloproteinases is locally concentrated at the interface, the mass balance equation reads$^{40,44}$

$$ \frac{\partial \rho^{\alpha}}{\partial t} + \nabla \cdot (\rho^{\alpha} \mathbf{v}^{\alpha}) = 0, \quad (3) $$

where we considered that non-convective mass fluxes are zero everywhere in the bulk volumes $V^\alpha$ but on both sides of the interface, possibly having a discontinuity across $\Sigma(t)$. Since the spatial density $\rho^\alpha$ can be considered as constant, the impressibility condition in Eq. (3) is satisfied if we take $\mathbf{v}^\alpha = 0$, thus assuming that either the capillary sprout degrades and invades the extracellular space without deforming it$^{49}$, or that the growing ECs fill $V^-$ without inducing any volumetric deformation.

On the other hand, defining the homogenized surface density field $\rho_{\Sigma}$ and surface mass source $\gamma_{\Sigma}$, the mass balance of the surface density for the ECs of the capillary wall reads$^{40,44}$

$$ \frac{\delta_t \rho_{\Sigma}}{\delta_t t} + \nabla_{\Sigma} \cdot (\rho_{\Sigma} \mathbf{v}_{\Sigma s}) - K \rho_{\Sigma} \bar{\mathbf{v}}_{\Sigma n} = \rho_{\Sigma} \gamma_{\Sigma} + \llbracket \rho(\bar{\mathbf{v}}_{\Sigma n} - \mathbf{v}_n) + \mathbf{n}_{\Sigma} \cdot \mathbf{m} \rrbracket \quad (4) $$

where $\nabla_{\Sigma} \cdot (\mathbf{v}) = (\mathbf{I} - \mathbf{n}_{\Sigma} \otimes \mathbf{n}_{\Sigma}) \cdot \nabla$ is the surface divergence, $K$ is twice the local mean curvature, $\delta_t$ indicates the Thomas (convected) derivative, $\llbracket (\cdot) \rrbracket = (\cdot)^+ - (\cdot)^-$ is the jump operator and $\mathbf{m}$ is a non-convective mass flux vector. Taking $\rho_{\Sigma}$ constant in time and in space, the mass balance for the surface density of ECs (5) leads to

$$ \rho_{\Sigma} \nabla_{\Sigma} \cdot (\mathbf{v}_{\Sigma s}) - K \rho_{\Sigma} \bar{\mathbf{v}}_{\Sigma n} = \rho_{\Sigma} \gamma_{\Sigma} + \llbracket \rho \bar{\mathbf{v}}_{\Sigma n} + \mathbf{n}_{\Sigma} \cdot \mathbf{m} \rrbracket. \quad (5) $$

Inside the proliferative interface, new ECs are constantly produced by a surface mass source $\gamma_{\Sigma}$, whereas the extracellular material is degraded$^{49}$ in order to make the capillary advance and the waste material is transported
inside the capillary in order to be eliminated. Therefore, we can define the mass fluxes $\mathbf{m}^\alpha$ at the interface by setting

$$ \mathbf{m}^{\alpha} = \alpha \rho^{\alpha} (\bar{\nu}_{\Sigma n} - v_{n}^{\alpha}) = \alpha (\rho^{\alpha} \bar{\nu}_{\Sigma n}) \mathbf{n}^{\alpha} \quad \text{on } \Sigma(t). \quad (6) $$

Assuming that $\rho^+ = \rho^-$, motivated by the fact that both materials are mainly composed by water, the previous relation states that the growing material deposited in the capillary lumen is created with the same rate of the degraded ECM at the interface. Substituting Eq. (6) into Eq. (5) and, for the sake of simplicity, restricting to surface divergence-free motion of the interface, we obtain that $\gamma_{\Sigma} = -K\bar{\nu}_{\Sigma n}$. This interface balance law states that the proliferation of EC is proportional to the curvature, as observed in some biological experiments$^{50}$.

For the sake of simplicity, we here consider only the effect of the VEGF, which is the main chemical involved in the early stage of angiogenesis$^{15}$. In the proposed model, we assume that the secreted VEGF diffuses from the tumor cell located at the boundary $\partial V_t$ (see Fig. 1) through the extracellular region and the intercapillary space, with the same diffusion coefficient $D_c$, so that the balance of the VEGF concentration, $c$, reads

$$ \frac{\partial c^{\alpha}}{\partial t} - D_c \nabla^2 c^{\alpha} = - \beta^{\alpha} c^{\alpha} \quad (7) $$

being $\beta^\alpha$ the decay rate of the substance in the volume $V^\alpha$. Since the VEGF is a diffusible factor with negligible inertia, the introduction of an interfacial field is not required and Eq. (7) must be complemented with the following boundary conditions

$$ \|c\| = 0, \quad (8) $$

$$ \|D_c \nabla c\| \cdot \mathbf{n}_{\Sigma} = \gamma_c(K, c), \quad (9) $$

which represent the continuity of the chemical field and the jump of the normal gradient due to the absorption rate $\gamma_c(K, c)$ at the interface.

Let us now impose the balance of the physical linear momentum inside the biological system. Neglecting inertial terms and in absence of external forces, the balance of linear momentum inside the volume $V^\alpha$ reads

$$ \nabla \cdot \boldsymbol{\sigma}^{\alpha} = 0, \quad (10) $$

being $\boldsymbol{\sigma}^\alpha$ the Cauchy stress tensor of the material in the volume $V^\alpha$. As done for the mass balance, we have to consider also surface quantities along with volumetric ones. Therefore, defining the Cauchy stress tensor for the interface, $\boldsymbol{\sigma}_{\Sigma}$, the balance of linear momentum for the interface, under the condition (6), reads

$$ \nabla_{\Sigma} \cdot \boldsymbol{\sigma}_{\Sigma} + \mathbf{n}_{\Sigma} \cdot \|\boldsymbol{\sigma}\| = 0. \quad (11) $$

Considering that the living material in $V^+$ and $V^-$ behaves as a perfect fluid with surface tension $\omega$, being $\boldsymbol{\sigma}_{\Sigma} = \omega \mathbf{I}_{\Sigma} = \omega(\mathbf{I} - \mathbf{n}_{\Sigma} \otimes \mathbf{n}_{\Sigma})$, eq. (11) reduces to the standard Young-Laplace equation.

In order to impose the thermodynamical consistency of the model, the reduced dissipation inequalities inside the volume $V^\alpha$ and on the surface $\Sigma(t)$ should be defined (see the Supplementary Information for further details). Since the VEGF concentration is continuous at the interface, the associated volumetric chemical potentials $\mu^\alpha$ should also be continuous, i.e. $\|\mu\| = 0$. Accordingly, the quasi-static reduced dissipation inequality on $\Sigma(t)$ (see the Supplementary Information) reads

$$ \mathbf{n}_{\Sigma} \cdot [D_c \nabla c] \mu - \mathbf{n}_{\Sigma} \cdot [\boldsymbol{\sigma} \cdot \mathbf{v}_{\Sigma}] \geq 0. \quad (12) $$

Under the assumption $v_{\Sigma n} = 0$, the substitution of (9) and (11) into (12) leads to $\gamma_c \mu + K \omega v_{\Sigma n} \geq 0$. In particular in the following we will consider the thermodynamically admissible condition

$$ \gamma_c \mu + K \omega v_{\Sigma n} = 0. \quad (13) $$

Furthermore, we will assume that the VEGF activates ECs in a region of a nearby capillary where the concentrations of the tumor angiogenic growth factor reaches a given threshold$^{4,35}$, namely $c_0$. Thus, the consumption term $\gamma_c$ takes the following form

$$ \gamma_c = - \gamma K H (-K) \mu (c - c_0)_+ \quad (14) $$

where $H$ is the Heaviside function (i.e. $H(-K) = 1$ if $K \leq 0$, else $H(-K) = 0$), $\langle \cdot \rangle_+$ is the positive part of its argument and $\gamma$ is the uptake of chemical energy at the interface, that in this work we assume to be a constant. The previous constitutive assumption allows rewriting eq. (13) as

$$ v_{\Sigma n} = \frac{\mu^2 \gamma}{\omega} (c - c_0)_+ \quad \text{if } K \leq 0, \quad (15) $$

which allows deriving the normal velocity of the capillary wall. We remark that for $K=0$ the relation (13) is automatically satisfied, so that $v_{\Sigma n}$ has been extended for continuity also for $K=0$. The normal velocity depends on the curvature of the interface, since the latest influences the chemical potential at the interface, being

$$ \mu = \mu_0(1 - \eta K), \quad (16) $$
where $\mu_0$ is the nominal VEGF chemical potential and $\eta$ is a characteristic microscopic length, in accordance with the Gibbs-Thompson condition$^{51}$.

In summary, the mathematical model is given by eq. (7), with boundary conditions (8–9) on the moving discontinuity complemented by the absorption term (14), a fixed concentration $c_t$ on the boundary of the tumor and no-flux conditions on the remaining border of the domain, i.e.

$$c = c_t \quad \text{on} \quad \partial V_t, \qquad (17)$$

$$(D_c \nabla c) \cdot \mathbf{n} = 0 \quad \text{on} \quad \partial V_{ext}^{\alpha}, \qquad (18)$$

whereas the interface $\Sigma(t)$ moves with normal velocity

$$v_{\Sigma n} = \frac{\gamma \mu_0^2}{\omega} (1 - \eta K)^2 (c - c_0)_+ \quad \text{for} \quad K \le 0. \qquad (19)$$

In order to perform some numerical simulations, it is useful to derive the non-dimensional system of governing equations and to obtain the fundamental parameters of the model. Considering the following characteristic time $t_c$, length $l_c$, velocity $v_c$ and chemical concentration $c_c$:

$$t_c = \beta^{+1}, \quad l_c = \sqrt{D_c \beta^{+1}}, \quad v_c = \sqrt{D_c \beta^{+1}}, \quad c_c = c_t$$

and using the barred notation to denote dimensionless quantities, the system to be solved reads

$$\frac{\partial \bar{c}^{+}}{\partial t} = \nabla^2 \bar{c}^{+} - \bar{c}^{+} \quad \text{in } V^{+}, \qquad (20)$$

$$\frac{\partial \bar{c}^{-}}{\partial t} = \nabla^2 \bar{c}^{-} - \delta \bar{c}^{-} \quad \text{in } V^{-}, \qquad (21)$$

$$[\bar{c}] = 0 \quad \text{on} \quad \Sigma(t), \qquad (22)$$

$$[\nabla \bar{c}] \cdot \bar{n}_{\Sigma} = -\lambda(1 - \bar{\eta}K)KH(-K)(\bar{c} - \bar{c}_0)_+ \quad \text{on } \Sigma(t), \qquad (23)$$

$$\bar{v}_{\Sigma n} = \xi \lambda (1 - \bar{\eta} K)^2 (\bar{c} - \bar{c}_0)_+ \quad \text{if } K \le 0 \quad \text{on } \Sigma(t), \qquad (24)$$

$$\bar{c} = 1 \quad \text{on} \quad \partial V_t, \qquad (25)$$

$$(\nabla \bar{c}) \cdot \bar{n} = 0 \quad \text{on} \quad \partial V_{ext}^{\alpha}, \qquad (26)$$

where $\delta = \beta^{-}/\beta^{+}$, $\bar{c}_{0} = c_{0}/c_{t}$, $\lambda = \gamma\mu_{0}/D_{c}$, $\bar{\eta} = \eta/l_{c}$, $\xi = \mu_{0}c_{t}\sqrt{D_{c}}/(\omega\sqrt{\beta^{+}})$.

The nondimensionalization procedure leads to the definition of five dimensionless parameters, $\delta, \bar{c}_0, \lambda, \bar{\eta}, \xi$. In particular:

* $\delta$ represents the ratio between the decay of VEGF inside and outside the vessel.

* $\lambda$ is related to the ratio between the rate of absorption and the rate of diffusion of VEGF.

* $\xi$ represents the ratio between chemical and mechanical energies driving capillary growth.

* $\bar{c}_0$ is merely the ratio between the threshold required for EC activation and the VEGF concentration at the tumor border.

* $\bar{\eta}$ is the ratio between the microscopic length regulating chemical absorption at the interface and the VEGF diffusive length.

Finally the size characterizing the domain of the simulations are made dimensionless with respect to the characteristic length $l_c$.

For sake of simplicity, in the following we will omit the barred notation to denote dimensionless quantities.

## Results

**Finite Element Simulations.** The system of eqs (20–26) has been numerically implemented developing a finite element code with the open-source program FreeFem++$^{52}$. The jump condition (23) at the moving boundary has been introduced in the variational formulation of the problem. The equations for the chemical species (20–21) are solved on a triangular grid, fitting at every iteration the moving interface. Given the VEGF concentration $c_t$ at time $t$, the curvature of the boundary $K_t$ and the value of the chemical potential $\mu_t$, the normal velocity of the boundary $v_{\Sigma n,t}$ is computed from eq. (24). The nodes belonging to the interface are then explicitly moved accordingly to the computed velocity in order to obtain the new position of the discontinuity at time $t_{i+1}$. Once the new position is known, the curvature $K_{i+1}$, and the chemical potential $\mu_{i+1}$ are updated in order to compute
<figure><img src="image_5.png" /><figcaption>**Figure 2.** (a)Angiogenesis dynamics in numerical simulations: the initial sprout evolves and soon splits in second-generation vessels, which in turn split until a complex tree-like network is formed. The simulations have been obtained setting $ξ = 10$, $λ = 1$, $δ = 500$, $η = 0.0316$, $c_0 = 0.143$ in a domain with dimensionless length and height equal to 0.316 and the interface initially placed at $x_S = 0.079$. The total dimensionless time of this simulation is $T = 1.558$. (b) Morphological diagram of the simulated capillary morphology for different values of the parameter $ξ$. The simulations were obtained setting $λ = 1$, $δ = 100$, $η = 0.0316$, $c_0 = 0.3432$. We report at the bottom the dimensionless time at which each snapshot has been taken. The simulations stopped because of either the capillaries reached the tumour (left snapshot) or they formed anastomosis (central and right snapshots).</figcaption></figure>

the chemical field at time $t_{i+1}$, using an implicit-Euler scheme. Then, the mesh is adaptively refined, generating an increasing number of grid points in those region close to the higher curvature of the moving interface.

The spatio-temporal dynamics for the vessel formation resulting from numerical simulations is reported in Fig. 2a. Initially, a sprout arises from the parent vessel and it soon splits in two or more branches. Some of this second-generation branches do not grow sufficiently to give rise to third-generation branches (because the concentration of VEGF at their tips is not high enough), whereas some others soon split. The process continues until a dendritic structure is formed and the solid tumor is reached. Accordingly to biological observations$^{14,18,47}$, the branching velocity increases as the vessels approach the tumor, due to the increase in the VEGF concentration close to the tumour region and the formation of an increased number of branches with higher absolute curvatures. The visible higher frequency of branching at the edge of the network as the capillary sprouts approach the tumor is biologically known as the “brush border” effect$^{14,18}$.

**Sensitivity analysis.** The growth of the network evolves both temporally and spatially in response to the combined effects of angiogenic factors, migratory cues via the ECM, mechanical factors acting at the vessel interface (i.e. the surface tension) and geometric factors, that are summarized in the model dimensionless parameters. Therefore, we conducted a sensitivity analysis changing one dimensionless parameter at a time, while keeping the other parameters fixed. In particular, Fig. 2b reports the final simulated morphologies for different values of the parameter $ξ$, that drives, along with $λ$, the velocity of the network. We find that the vessel sprouting is strongly favored by smaller values of the parameter $ξ$, i.e. by low velocities, which is a common feature of growing biological systems in a diffusion-limited regime$^{41}$. Since $ξ$ is proportional to the ratio between the chemical energy of VEGF diffusion and the mechanical energy of the capillary tip ($ω$), small values of $ξ$ correspond to situations in which the capillary surface tension is dominant on the diffusion of chemicals. In contrast to the patterns emerging in other non-living and living systems$^{39,42,43}$, here the surface tension does not act as a stabilizing effect at small wavelengths, but it rather promotes the formation of lateral ramified branches. Indeed, the instability is governed in this problem by the curvature of the interface, similarly to dendritic growth problems$^{41}$, and the only short-scale cutoff is the size of the ECs in the capillary wall.

The parameter $η$ is also very influential for the formation of the capillary network, since it weights the increase of VEGF absorption, thus regulating the velocity of the front. The dimensionless $η$ represents the ratio between the microscopic length related to the chemical absorption at the interface and the characteristic diffusive length inside the volumes. Figure 3a reports the different morphologies obtained varying the parameter $η$ for different values of $ξ$, showing that more branched patterns with sharper tips arise for higher values of $η$. Indeed, the microscopic length related to the chemical dynamics at the interface is known to give the dimensional information to set the characteristic microscopic scale of a pattern in pure solidification$^{41}$. Actually, smaller values of $η$ significantly reduces the contribution of the curvature to the velocity, whereas higher values of this parameter amplify even small variations in the curvature, setting the formation of sharper tips.

Sharper and thinner vessel tips can be obtained also increasing the parameter $λ$ (see Fig. 3b), which weights the absorption of VEGF at the vessel interface: since the absorption is proportional to the curvature, higher values of $λ$ generates a sort of “tip-effect”, sequestering the VEGF on the tips of the vessel and thus hindering the formation of long later protrusion behind the front-head tip. For small values of $λ$, the finger-like capillary sprouts tend to bend toward each other, as observed in experiments$^{53}$. This process will eventually lead to numerous tip-to-tip and tip-to-sprout fusions known as anastomoses$^{15}$. In the present model, we do not implement the formation of anastomoses, for increased technical complications, and the simulations stops as soon as two branches touch.
<figure><img src="image_6.png" /><figcaption>**Figure 3.** (a) Morphological diagram of the simulated capillary morphology for different values of the parameter η and ξ, setting λ = 1, δ = 100, c₀ = 0.33838 in a domain with dimensionless length and height equal to 0.32 and the interface initially placed in x<sub>Σ</sub> = 0.08. (b) Morphological diagram of the simulated capillary morphology for different values of the parameter λ and ξ, setting δ = 100, η = 0.05, c₀ = 0.33838 in a domain with dimensionless length and height equal to 0.32 and the interface initially placed in x<sub>Σ</sub> = 0.08. At the bottom of each snapshot we reported the dimensionless time at which it was taken.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 4.** (a) Morphological diagram of the simulated capillary network for different values of the parameter δ and ξ, setting λ = 1, η = 0.032 in a domain with dimensionless length and height equal to 0.32 and the interface initially placed in x<sub>Σ</sub> = 0.08. The value of c₀ changes in order to keep constant the initial sprout size. (b) Morphological diagram of the simulated capillary morphology for different dimension of the domain of simulation (*H*), setting λ = 1, η = 0.032 and ξ = 10. The interface is placed at distance *H*/4 from the left-side of the domain and c₀ is chosen in order to keep constant the ratio between the initial sprout size and the domain dimension *H*. At the bottom of each snapshot we reported the dimensionless time at which it was taken.</figcaption></figure>

The effect of varying the parameters $\beta$, which is the ratio between the VEGF decay inside the vessel and inside the extracellular environment, strongly influence the evolution of the chemical field and the initial condition on the VEGF (obtained solving the stationary problem, without absorption on the vessel wall). Therefore in order to obtain the same initial thickness of the initial sprout, we have to change the initial value of $c_0$: in particular, as $\delta$ increases we have to set a smaller threshold $c_0$ in order to observe the formation of branches. In all cases, the time required to form the branched vessel network decreases with increasing $\delta$.

Finally, we consider the effect of varying the size of the domain, keeping constant the ratio between the size of the capillary, *d*, and the total size of the square domain, *H* (see Fig. 1b). Figure 4 shows that relatively thicker branches (with respect to the normalized initial sprout thickness) can be obtained as *H* increases, highlighting the existence of a size effect in the pattern selection.

**Physical interpretation of the model results.** The tree-like network structure emerging in the numerical simulations is somewhat similar to the morphology of crystal fronts in solidification problems, where a stable thermodynamic phase propagates into a metastable one. In both cases, the competition is between the diffusion of a given field (e.g. thermal or chemical) on one hand, and the microscopic dynamics occurring at the interface (such as the surface tension and the chemical kinetics of absorption at the interface) on the other, modulated by the ratio between the typical distances involved in the process and the diffusive length. The result of this competition is the onset of branches and the formation at the macroscopic scale of intricate tree-like structures, which closely resemble to the ones experimentally observed⁵⁴,⁵⁵. Thus, this work shows that the onset of vessel sprouting from the existing vasculature can be reproduced without defining ad-hoc empirical rules for branching, as in
previous mathematical approaches, but it simply results from the interactions between chemical and mechanical phenomena, generating an unstable growth process at the vessel interface.

Unfortunately, a direct quantitative comparison between the numerical simulations and the biological experiments is not straightforward, since not all the data required by the mathematical model, even though measurable in principle, are reported in literature.

However, it is useful to compare our results with the available experimental data with respect to the size of the domain and the time required for forming the vascular network. The diffusion coefficient inferred from biological data, for the vast majority of angiogenic growth factors, is in the order of $10^{-600 \mu m^2/s}$$^{4,56,57}$, while the decay rate of the VEGF in the surrounding tissue is in the order of $0.456-0.65 \text{ h}^{-1}$$^{56,58}$. Accordingly, taking $D_c = 10 \mu m^2/s$ would result in a computational domain having a width of about 90µm, which is comparable to the value used in previous mathematical models$^3$; whereas $D_c = 600 \mu m^2/s$ would represent a width of about 0.7 mm, that is closer to the tumor-capillary distance reported in literature$^{5,58}$. For what concern the dynamics, experiment on the avascular cornea have shown that angiogenetic processes take 2–5 days to form an initial sprout$^{59}$. Once the sprouts are developed, the growth rate of new capillaries is much faster, with an initial growth rate of about 0.22 mm/day which increases up to 0.61 mm/day when the network becomes highly branched$^{59}$. An increasing velocity over time has also been reported in previous mathematical models$^3$, precisely from ~0.20 mm/day to ~1 mm/day, highlighting the long time required for the initial spouting. These effects are consistent with our numerical results, where the initial sprout formation requires a long time, whereas branching accelerates as the new capillaries approach the tumor (see the Supplementary Information for further details). Even if our results at this stage do not have any intent of being quantitative, the proposed model should be regarded as a proof-of-concept to outline a new mechanism of angiogenesis formation.

## Discussion

In this work we have presented a theoretical and numerical analysis for studying the well-orchestrated sequence of biological events occurring during tumour angiogenesis. In particular, we have proposed a thermodynamically-consistent growth model describing the morphology of the new developing vasculature, depending on the different geometrical and chemo-mechanical factors involved in the initial stages of network formation. The proposed continuous approach couples the diffusion of angiogenic molecules (VEGF) with the chemically-driven migration and proliferation of ECs. Unlike existing modelling approaches$^{19–22}$ we explicitly reproduce both the temporal and spatial network evolution taking into account for the chemo-mechanical cues (i.e. the surface tension, the variation of the chemical potential across a curved interface) and geometric factors (i.e. the capillary to tumour distance). In addition, this application of a mechanistic interface growth model to angiogenesis is completely new. Our approach brings novel insights on the role played by physical forces along with chemical factors in directing angiogenesis. In particular, five dimensionless parameters, encapsulating the geometric, chemical and mechanical cues, are found to characterize the model dynamics: $\delta_t, \bar{\epsilon}_0, \lambda, \eta, \xi$. The effects of the different parameters on the predicted morphology are studied in the numerical simulations. Angiogenesis initiates with a first vessel tip sprouting when the local VEGF concentration reaches a given threshold. Such a protrusion undergoes an unstable growth process, with the curvature having a destabilizing effect on the capillary motion. Indeed the simulations demonstrate that both the microscopic chemical kinetics at the interface of the capillary, following the Gibbs-Thompson law, and the size effects of the biological domains compete with the reaction-diffusion of VEGF in order to determine the occurrence of complex branched patterns. The predicted morphologies resulting from the numerical simulations strikingly resemble the tree-like vascular structures experimentally observed *in vivo* and *in vitro*$^{54,55}$. Nonetheless, this work should be regarded as a proof-of-concept of the fundamental role of physical forces during angiogenesis, since it neglects some molecular processes acting at the cellular level. Thus, future studies should focus on the incorporation of the molecular and other cellular mechanisms observed in biological experiments (e.g. the Notch signaling$^{16,30,60}$, the interaction between stalk cell and tip cell$^{16,60}$). Moreover, we remark that in this simulations the concentration of VEGF is kept constant at the tumor boundary, whereas time-dependent conditions at the tumor interface can be useful to consider the physiological feedback between the tumor oxygenation due to the onset of the new vasculature and the secretion of VEGF in response to hypoxic condition. Taking this feedback mechanism into account might affect the dynamics of vessel formation as the branches approaches the tumor. It is also worth noticing that the proposed numerical implementations should be refined in order to describe the formation of anastomosis and the role played by the environmental stress in influencing capillary branching. In fact, despite the mathematical model has been derived considering the role of the stress field in the surrounding environment, the simulations have been performed considering an external inviscid fluid and a surface tension acting at the interface. By formulating proper hypotheses on the mechanical stress exerted by the surrounding tissue, future works should focus on the role played by external mechanical cues on the overall process. Indeed, there is a growing recognition that the balance between internally generated and externally applied forces, along with ECM remodelling and mechanical factors connected with blood flow or extravascular mechanical stress are key determinants of a cell's fate and function and are important regulators in postnatal physiological angiogenesis$^{61–65}$. Nevertheless, mature ECs can develop capillary-like networks in cell culture even in the absence of flow or any other externally applied stresses$^61$, demonstrating that extrinsic stresses are not strictly necessary for triggering angiogenesis. Furthermore, even if the majority of the existing mathematical models of angiogenesis are 2D$^{1,3,4,25}$, angiogenesis is typically a 3D process, with tips sprouting in directions other than that of the propagating vascular front. Thus, future works will certainly focus on the 3D implementation of the proposed model. In spite of these limitations, the model provides original insights about the influence of the physical and chemical effects on the pattern dynamics during angiogenesis. In conclusion, the proposed mechanistic approach, possibly combined with biologically more detailed diffusion-based mathematical models, has the potential to foster our understanding on the process of vessel formation. A deeper comprehension of the key factors directing angiogenesis is fundamental for many clinical
applications, since vessel tortuosity is known to strongly affect the anti-tumor treatment response⁶⁶,⁶⁷. Finally, even though we focus exclusively on tumour-induced angiogenesis, the proposed model can be useful to model other biological processes, such as wound healing⁶⁸ and tissue optimization in engineering scaffolds⁶⁹.

## References

1. Alarcon, T., Byrne, H., Maini, P. & Panovska, J. Studies in Multidisciplinarity, vol. 3, chap. *Mathematical modelling of angiogenesis and vascular adaptation*, 369–387 (Elsevier, 2006).

2. McDougall, S. R., Anderson, A. R. A. & Chaplain, M. A. J. Mathematical modelling of dynamic adaptive tumour-induced angiogenesis: Clinical implications and therapeutic targeting strategies. *J. Theor. Biol.* **241**, 564–589 (2006).

3. Levine, H. A., Pamuk, S., Sleeman, B. & Nilsen-Hamilton, M. Mathematical modeling of capillary formation and development in tumor angiogenesis: Penetration into the stroma. *Bulletin of Mathematical Biology* **63**, 801–863 (2001).

4. Anderson, A. R. A. & Chaplain, M. A. J. Continuous and discrete mathematical models of tumor-induced angiogenesis. *Bull. Math. Biol.* **60**, 857–900 (1998).

5. Hillen, F. & Griffioen, A. W. Tumour vascularization: sprouting angiogenesis and beyond. *Cancer Metastasis Rev.* **26**, 489–502 (2007).

6. Graham, C. H. & Lala, P. K. Mechanisms of placental invasion of the uterus and their control. *Biochem. Cell Biol.* **70**, 867–874 (1992).

7. Arnold, F. & West, D. C. Angiogenesis in wound healing. *Pharmacol. Ther.* **52**, 407–422 (1991).

8. Walsh, D. Angiogenesis and arthritis. *Rheumatology* **38**, 103–112 (1999).

9. Folkman, J. Tumor angiogenesis. *Adv. Cancer Res.* **43**, 175–203 (1985).

10. Folkman, J. Angiogenesis in cancer, vascular, rheumatoid and other disease. *Nat. Med.* **1**, 21–31 (1995).

11. Risau, W. Mechanisms of angiogenesis. *Nature* **386**, 671–674 (1997).

12. Madri, J. A. & Pratt, B. M. Endothelial cell-matrix interactions: *in vitro* models of angiogenesis. *J. Histochem. Cytochem.* **34**, 85–91 (1986).

13. Folkman, J. & Klagsbrun, M. Angiogenic factors. *Science* **235**, 442–447 (1987).

14. Ausprunk, D. & Folkman, J. Migration and proliferation of endothelial cells in preformed and newly formed blood vessels during tumour angiogenesis. *Microvasc. Res.* **14**, 53–65 (1977).

15. Eilken, H. M. & Adams, R. H. Dynamics of endothelial cell behavior in sprouting angiogenesis. *Curr. Opin. Cell Biol.* **22**, 617–625 (2010).

16. Gerhardt, H. Simulation of vessel morphogenesis using cellular automata. *Organogenesis* **4**, 241–246 (2008).

17. McDougall, S. R., Anderson, A. R. A., Chaplain, M. A. J. & Sherratt, J. A. Mathematical modelling of flow through vascular networks: Implications for tumour-induced angiogenesis and chemotherapy strategies. *Bulletin of Mathematical Biology* **64**, 673–702 (2002).

18. Muthukkaruppan, V. R., Kubai, L. & Auerbach, R. Tumor-induced neovascularization in the mouse eye. *J. Natl. Cancer Inst.* **69**, 699–705 (1982).

19. Liotta, L. A., Saidel, G. M. & Kleinerman, J. Diffusion model of tumor vascularization. *Bull. Math. Biol.* **39**, 117–128 (1977).

20. Balding, D. & McElwain, D. L. S. A mathematical model of tumour-induced capillary growth. *J. Theor. Biol.* **114**, 53–73 (1985).

21. Byrne, H. M. & Chaplain, M. A. J. Mathematical models for tumour angiogenesis: numerical simulations and nonlinear wave solutions. *Bull. Math. Biol.* **57**, 461–486 (1995).

22. Orme, M. & Chaplain, M. A mathematical model of the first steps of tumour related angiogenesis: Capillary sprout formation and secondary branching. *IMA J. Math. Appl. Med. Biol.* **13**, 73–98 (1996).

23. Stokes, C. L. & Lauffenburger, D. A. Analysis of the roles of microvessel endothelial cell random motility and chemotaxis in angiogenesis. *J. Theor. Biol.* **152**, 377–403 (1991).

24. Alarcon, T., Byrne, H. & Maini, P. A cellular automaton model for tumour growth in inhomogeneous environment. *J. Theor. Biol.* **225**, 257–274 (2003).

25. Vilanova, G., Colominas, I. & Gomez, H. Capillary networks in tumor angiogenesis: From discrete endothelial cells to phase-field averaged descriptions via isogeometric analysis. *Int. J. Numer. Meth. Biomed. Eng.* **29**, 1015–1037 (2013).

26. Milde, F., Bergdorf, M. & Koumoutsakos, P. A hybrid model for three-dimensional simulations of sprouting angiogenesis. *Biophysical Journal* **95**, 3146–3160 (2008).

27. Mantzaris, N., Webb, S. & Othmer, H. Mathematical modeling of tumor-induced angiogenesis. *J. Math. Biol.* **49**, 111–187 (2004).

28. Preziosi, L. *Cancer modelling and simulation* (Chapman and Hall/CRC, 2003).

29. Scianna, M., Bell, C. & Preziosi, L. A review of mathematical models for the formation of vascular networks. *J. Theor. Biol.* **333**, 174–209 (2013).

30. Bentley, K., Mariggi, G., Gerhardt, H. & Bates, P. A. Tipping the balance: Robustness of tip cell selection, migration and fusion in angiogenesis. *PLoS Comput. Biol.* **5**, e1000549 (2009).

31. Markus, M., Böhm, D. & Schmick, M. Simulation of vessel morphogenesis using cellular automata. *Mathematical Biosciences* **156**, 191–206 (1999).

32. Nekka, F., Kyriacos, S., Kerrigan, C. & Cartilier, L. A model of growing vascular structures. *Bull. Math. Biol.* **58**, 409–424 (1996).

33. Olsen, L., Sherratt, J., Maini, P. & Arnold, F. A mathematical model for the capillary endothelial cell-extracellular matrix interactions in wound-healing angiogenesis. *IMA J. Math. Appl. Med. Biol.* **14**, 261–281 (1997).

34. Holmes, M. J. & Sleeman, B. D. A mathematical model of tumour angiogenesis incorporating cellular traction and viscoelastic effects. *J. Theor. Biol.* **202**, 95–112 (2000).

35. Paweletz, N. & Knierim, M. Tumor-related angiogenesis. *Crit. Rev. Oncol. Hematol.* **9**, 197–242 (1989).

36. Carmeliet, P. Angiogenesis in health and disease. *Nature Med.* **9**, 653–660 (2003).

37. Jain, R. Molecular regulation of vessel maturation. *Nature Med.* **9**, 685–693 (2003).

38. Ben-Jacob, E. From snowflakes formation to growth of bacterial colonies. part i: Diffusive patterning in azoic systems. *Contemp. Phys.* **34**, 247–273 (1993).

39. Cross, M. & Hohenberg, P. Pattern formation outside of equilibrium. *Rev. Mod. Phys.* **65**, 851–1112 (1993).

40. Romano, A. Thermomechanics of Phase Transitions in Classical Field Theory. Series on advances in mathematics for applied sciences (World Scientific, 1993).

41. Langer, J. Instabilities and pattern formation in crystal growth. *Rev. Mod. Phys.* **52**, 1–30 (1980).

42. Ben-Jacob, E. From snowflakes formation to growth of bacterial colonies ii: Cooperative formation of complex colonial patterns. *Contemp. Phys.* **38**, 205–241 (1997).

43. Giverso, C., Verani, M. & Ciarletta, P. Branching instability in expanding bacterial colonies. *Journal of The Royal Society Interface* **12**, doi: 10.1098/rsif.2014.1290 (2015).

44. Ciarletta, P., L.Preziosi & Maugin, G. A. Mechanobiology of interfacial growth. *Journal of the Mechanics and Physics of Solids* **61**, 852–872 (2013).

45. Mullins, W. W. & Sekerka, R. F. Morphological stability of a particle growing by diffusion or heat flow. *J. Appl. Phys.* **34**, 323–329 (1963).

46. Cavallo, T., Sade, R., Folkman, J. & Cotran, R. S. Tumor angiogenesis. rapid induction of endothelial mitoses demonstrated by autoradiography. *J. Cell. Biol.* **54**, 408–20 (1972).
47. Gimbrone, M. A., Leapman, S. B., Cotran, R. S. & Folkman, J. Tumor angiogenesis: iris neovascularization at a distance from experimental intraocular tumors. *J Natl Cancer Inst.* **50**(1), 219–228 (1973).

48. Lash, J., Sherman, W. & Hamlin, R. Capillary basement membrane thickness and capillary density in sedentary and trained obese zucker rats. *Diabetes* **38**, 854–860 (1989).

49. Kabelic, T., Ganbisa, S., Glaser, B. & Liotta, L. Basement membrane collagen: degradation by migrating endothelial cells. *Science* **221**, 281–283 (1983).

50. Rumpler, M., Woesz, A., Dunlop, J. W., van Dongen, J. T. & Fratzl, P. The effect of geometry on three-dimensional tissue growth. *Journal of The Royal Society Interface* **5**, 1173–1180 (2008).

51. Cermelli, P. & Jabbour, M. E. Possible mechanism for the onset of step-bunching instabilities during the epitaxy of single-species crystalline films. *Phys. Rev. B* **75**, 165409 (2007).

52. Hecht, F. New development in freefem++.*J. Numer. Math.* **20**, 251–265 (2012).

53. Paweletz, N. & Knierim, M. Tumor-related angiogenesis. *Crit. Rev. Oncol. Hematol.* **9**, 197–242 (1989).

54. Auerbach, R., Lewis, R., Shinners, B., Kubai, L. & Akhtar, N. Angiogenesis assays: A critical overview. *Clinical Chemistry* **49**, 32–40 (2003).

55. Nguyen, D.-H. *T. et al.* Biomimetic model to reconstitute angiogenic sprouting morphogenesis *in vitro*. *PNAS* **110**, 6712–6717 (2013).

56. Serini, G. *et al.* Modeling the early stages of vascular network assembly. *EMBO J.* **22**(8), 771–1779 (2003).

57. Miura, T. & Tanaka, R. *In vitro* Vasculogenesis Models Revisited-Measurement of VEGF Diffusion in Matrigel. *Math. Model. Nat. Phenom.* **4**(4), 118–130 (2009).

58. Plank, M. J., Sleeman, B. D. & Jones, P. F. A Mathematical Model of an *In Vitro* Experiment to Investigate Endothelial Cell Migration. *Journal of Theoretical Medicine* **4**(4), 251–270 (2002).

59. Brem, H. & Folkman, J. Inhibition of tumor angiogenesis mediated by cartilage. *J Exp Med.* **141**(2), 427–439 (1975).

60. Siekmann, A. F., Affolter, M. & Belting, H. The tip cell concept 10 years after: New players tune in for a common theme. *Exp. Cell. Res.* **319**, 1255–1263 (2013).

61. Shiu, Y. T. *et al.* The role of mechanical stresses in angiogenesis. *Critical reviews in biomedical engineering* **33**, 431–510 (2005).

62. Ichioka, S. *et al.* Effects of shear stress on wound-healing angiogenesis in the rabbit ear chamber. *Surg. Res.* **72**, 29–35 (1997).

63. Zhou, A. L. *et al.* Capillary growth in overloaded, hypertrophic adult rat skeletal muscle: an ultrastructural study. *Anat. Rec.* **252**, 49–63 (1998).

64. Egginton, S. *et al.* Capillary growth in relation to blood flow and performance in overloaded rat skeletal muscle. *J. Appl. Physiol.* **85**, 2025–2032 (1998).

65. Korff, T. & Augustin, H. G. Tensional forces in fibrillar extracellular matrices control directional capillary sprouting. *J. Cell Sci.* **112**, 3249–3258 (1999).

66. Bullitt, E., Wolthusen, P. A., Brubaker, L., Lin, W., Zeng, D. & Van Dyke, T. Malignancy-associated vessel tortuosity: a computer-assisted, MR angiographic study of choroid plexus carcinoma in genetically engineered mice. *American journal of neuroradiology* **27**(3), 612–619 (2006).

67. Bullitt, E. *et al.* Blood vessel morphological changes as visualized by mra during treatment of brain metastases: A feasibility study. *Radiology* **245**, 824–830 (2007).

68. Vermolen, F. J., Gefen, A. & Dunlop, J. W. C. *In Vitro* "Wound" Healing: Experimentally Based Phenomenological Modeling. *Advanced Engineering Materials* **14**, B76–B88 (2012).

69. Bidan, C. *et al.* Geometry as a factor for tissue growth: Towards shape optimization of tissue engineering scaffolds. *Advanced Healthcare Materials* **2**, 186–194 (2013).

## Acknowledgements

This work is funded by the AIRC MFAG grant 17412, partially supported by the “Start-up Packages and PhD Program” project, co-funded by Regione Lombardia through the “Fondo per lo sviluppo e la coesione 2007–2013 -formerly FAS”, by the “Progetto Giovani GNFM 2015”, funded by the National Group of Mathematical Physics (GNFM-INdAM) and by the European Commission FP7 grant 260566. We are grateful to Luigi Preziosi and Marco Verani for helpful discussions.

## Author Contributions

C.G. and P.C. developed the mathematical model, C.G. performed the numerical simulations, C.G. and P.C. analyzed the data, C.G. and P.C. wrote the paper.

## Additional Information

Supplementary information accompanies this paper at http://www.nature.com/srep

Competing financial interests: The authors declare no competing financial interests.

**How to cite this article:** Giverso, C. and Ciarletta, P. Tumour angiogenesis as a chemo-mechanical surface instability. *Sci. Rep.* **6**, 22610; doi: 10.1038/srep22610 (2016).

<figure><img src="image_8.png" /></figure>

This work is licensed under a Creative Commons Attribution 4.0 International License. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in the credit line; if the material is not included under the Creative Commons license, users will need to obtain permission from the license holder to reproduce the material. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/