A new model for the emergence of blood capillary networks

1.   Introduction

Networks appear ubiquitously in nature and in man-made structures. Examples of networks in living systems are the circulatory, respiratory and nervous systems of vertebrates, the branches, roots and leaves of trees, coral formations, bacterial colonies. In nature, river networks, erosion patterns, lightnings of thunder are other example of networks. The literature on network modelling is enormous and growing at an impressive speed. Yet, most works assume a given topological structure, possibly random or dynamic, but which is, for a given realization or a given time, well-defined. However, in many cases, networks are fuzzy objects which, according to the observed scale, adopt a different structure. This is particularly the case of emergent networks, where a network pattern emerges from a continuum. One can think for instance of erosion gradually sculpting an originally flat landscape into gullies, channels and ultimately valleys. Many networks in nature emerge in this way. Current network modelling methodologies are unadapted to capture this emergence phase because they need an already established network structure from the start. In this paper, we propose a novel methodology to model emergent networks and apply it to the emergence of blood capillary networks, the so-called vasculogenesis or angiogenesis phenomena. However, the approach is fairly versatile and could easily be adapted - modulo appropriate changes - to other emergent networks. In particular, we refer to earlier works developing similar ideas for ant trail formation [6] and tissue self-organization [60].

Vasculogenesis (the de novo formation of new blood vessels) or angiogenesis (the development of new blood vessels from existing ones, see e.g. [67]) play a fundamental role in living systems by influencing growth, regeneration and reparation through an increase of blood flow and tissue oxygenation and consequently by boosting the transport of nutrients as well as the disposal of waste. Angiogenesis has been the subject of intense study over the last decades for its crucial role in the development of cancerous tumours. Once a tumour has reached certain level of maturation, the surrounding tissue starts to suffer from hypoxia, or lack of oxygen, which induces the accumulation of hypoxia-inducible factors. As a consequence, hypoxic cells start to secrete tumour angiogenic growth factors (TAFs) such as the vascular endothelial growth factor (VEGF) [39,83]. The TAFs diffuse through the tissue until reaching the nearby blood vessels, triggering the formation of a new vascular network to supply the tumour with nutrients and oxygen. Angiogenesis is not unique to tumour growth but is also a key part of wound healing, the menstrual cycle, embryonic development and of several diseases such as rheumatoid arthritis [10,19,20]. New blood vessels grow through the Extracellular Matrix (ECM), which is a collection of collagen fibers, interstitial cells, proteoglycans and matrix binding molecules among many other components. The ECM provides support for cells which gather together to form organs. It also provides mechanical support for cells and offers an environment for the transmission of chemical cues.

There are many different approaches to model vasculogenesis and angiogenesis in the literature. Many of them are two-dimensional cell-based models, describing the migration of endothetial cells to form new blood vessels and including various aspects of the tissue environment [5,15,44,53,76]). Evolving network models have been developed in [70,72]. Other models use a continuum approach for the endothelial cell (EC) density [3,9,77] and the link between cell-based and continuum approaches is investigated e.g. in [63]. In most of these models, blood flow plays no or little role in determining the geometry of the network. For instance, the efficiency of blood perfusion in a predetermined network is studied in [26]. For a review of recent models, we refer to [71]. One important characteristics of our model is that blood and oxygen flows have the leading role in shaping the network structure. Another feature is that the network structure is not imposed a priori but emerges as an outcome of the evolution rules ascribed to the agents.

In the present work, we describe the emergence of vascular networks by considering agents (EC or group of such cells referred to as capillary elements) obeying elementary heuristic rules. Indeed, the biology, biophysics and biochemistry involved are incredibly complex and still incompletely known [61,74]. However, by their combined influence cells are likely to respond in a determinable (if not deterministic) way to environmental cues, such as blood flow, oxygen concentration, etc. In other words, cells behave like social agents with a behavior determined by their own state and their environment. This is the viewpoint developed here. Moreover, cell social behavior is not fully understood yet and the model aims to provide a platform to test the validity of such behavioral hypotheses. Once elucidated, cell behavior can be traced back to biochemistry through experiments whose design is facilitated by the intuition gained by the model.

In this model, unlike most models of the literature we put the flow of blood, interstitial fluid and oxygen at the heart of the capillary creation process. We assimilate the tissue and ECM as a porous medium in which blood and interstitial fluid (we make no distinction between those) flow under the influence of a pressure gradient (powered by the pressure difference between arteries and veins). We assume blood velocity follows a Darcy law, i.e. the fluid is incompressible and its velocity is proportional to the pressure gradient through a proportionality matrix called the hydraulic conductivity tensor. Oxygen is supplied by e.g. an oxygenated tissue, transported by the blood velocity and diffuses through and is consumed by the tissue. If there are no capillaries, due to low hydraulic conductivity of avascular tissue flow is very slow and oxygen is mostly transported by diffusion. Thus, oxygen concentration falls quickly off even at short distance from the supply due to quick consumption. This is where capillary creation comes at play. So far, the treatment of blood flow and oxygen transport are through a continuum model. We now introduce an agent-based description of the formation of the capillary network, making our model a hybrid model coupling continuum and agent-based descriptions of the system.

We assume that capillary elements (representing individual EC or group of such cells) can be created in response to a gradient of oxygen concentration. The creation of this capillary element could either correspond to a de novo blood vessel creation or to the recruitment of EC previously at a different location. We do not make any difference between the two processes as, for the sake of simplicity, we do not follow the motion of these EC prior to their recruitment. For the same reason, we assume that once created, a capillary element stays in place without moving, until it is eventually destroyed. A capillary element is modelled by a rod-shape particle with finite area (we assume a 2-dimensional model). Once created, it locally enhances the hydraulic conductivity in the direction of the rod, contributing to locally increase the blood velocity and the flow of oxygen in this direction. This increase is limited to the area of the particle. The process is time-dynamic: capillary elements are created following a spatio-temporal Poisson process the intensity of which is large when and where the gradient of oxygen concentration is large. Each time a capillary element is created, we assume an instantaneous adaptation of the blood velocity through solving Darcy's equation with the new value of the hydraulic conductivity. This new velocity is fed into the convection-diffusion-reaction equation satisfied by oxygen concentration and the process goes on. It leads to the appearance of more capillary elements.

Importantly, we do not assume any connectivity between the capillary elements. They could appear anywhere following the Poisson process and do not need to intersect. The influence of each of them is summed up at each update of the hydraulic conductivity tensor. In spite of the absence of connectivity constraint between capillary elements, simulations show the creation of streets of capillary elements forming what we could assimilate as a new blood capillary. Connectivity emerges as an outcome of the model but does not need to be prescribed. This is an important advantage of the model as keeping track of connectivity requires the use of adequate data structures that considerably increase the complexity of code development and its computational cost.

In our model, the creation of capillary elements in response to a gradient of oxygen concentration is a proxy for the influence of growth factors such as VEGF, which we do not include in the model. Indeed, both gradients of VEGF and of oxygen concentration carry information about the direction to be followed by new capillaries. The VEGF signaling is likely to amplify this information for improved sensing by the cells or to be a chemical intermediate that triggers a stronger reaction than mere oxygen gradient. However, both encode similar information and ignoring VEGF allows us to reduce the model complexity without distorting the phenomenology too much. We do assume several other capillary creation mechanisms. One is by reinforcement of existing branches. Another one is through wall shear stress which has proven to play a key role in the creation of new blood vessels (see for instance [22,37]). We also consider pruning mechanisms when the capillary element density is too high.

To summarize, there are two main innovations in this work. From the modelling viewpoint, we propose a new paradigm for the emergence of capillary networks where, by contrast to the literature, the flow of blood plays the central role in defining the shape of the network. From the methodological viewpoint, we present a new method to describe evolving networks which relies on disconnected discrete elements. Connectivity is recovered at the large scale when computing macroscopic quantities (such as, in this work, the hydraulic conductivity).

The outline of the article if the following. In Section 2, we provide a complete description of the model although the details of the numerical methods are deferred to the appendices. Results of the model are given in Section 3. A discussion is then developed in Section 4.

2.   Model

2.1.   Presentation

We consider a two-dimensional model. The tissue occupies a spatial domain $\Omega \subset {\mathbb R}^2$. We will specifically consider two different geometries described in Section 2.6. As outlined above, the model considers four different entities:

● The flow of blood or interstitial fluid described by its pressure $p(x, t)$ and its velocity $u(x, t)$ where $x \in \Omega$ and $t \geq 0$ and where $p \in {\mathbb R}$ and $u \in {\mathbb R}^2$. The quantities $u$ and $p$ are continuum variables which satisfy Partial Differential Equations (PDE), namely Darcy's equations, described in Section 2.2, with boundary conditions given, according to the geometrical case studied, in Section 2.6.

● The flow of oxygen, described by the oxygen concentration $\rho(x, t) \in {\mathbb R}$. The quantity $\rho$ is also a continuum variable satisfying a Partial Differential Equation (PDE), namely a convection-diffusion-reaction equation described in Section 2.3 with initial and boundary conditions given in Section 2.6.

● Capillary elements in turn are discrete quantities (particles) carrying a direction and their dynamics is given by an agent-based model described in Section 2.4.

● The tissue and the ECM mediate the interaction between the capillaries on the one hand and the blood and oxygen flows on the other hand. The tissue characteristics determine the coefficients of the PDE for the blood and the oxygen and those are dependent on the position and direction of the capillaries. The computation of these coefficients is detailed in Section 2.5.

We now describe the various elements of the model in detail. The numerical values of all the parameters are summarized in Table 1.

2.2.   Blood and interstitial fluid flow

We lump blood and interstitial fluid in one and single fluid which we will call 'blood' for simplicity. We assume that blood is an incompressible Newtonian fluid described by its pressure $p(x, t)$ and its velocity $u(x, t)$. The tissue viewed as a mixture of cells, ECM and blood vessels is assimilated to a single porous medium with hydraulic conductivity tensor $\textbf{K}(x, t)$, which is a symmetric uniformly positive-definite matrix. Because of the presence of capillaries and the possible creation/removal of capillaries, $\textbf{K}$ is dependent on both position and time. In a porous medium velocity and pressure are linked by Darcy's law

where $\nabla_x p$ denotes the gradient of $p$. This equation asserts that in spite of the time variations of the conductivity matrix $\textbf{K}$, pressure and velocity instantaneously adjust one with each other according to (1). This is consistent with the assumptions leading to Darcy's law, i.e. viscous terms dominate inertial ones, resulting in a quasi-steady flow. Darcy's law is complemented with the incompressibility condition:

where $\nabla_x \cdot \textbf{u}$ stands for the divergence of the vector field $u$. Substituting (1) into (2) we obtain $p$ as a solution of the following elliptic problem:

where $T$ is the total simulation time. The expression of $\textbf{K}$ will be given in Section 2.5. This equation is complemented by boundary conditions detailed in Section 2.6.

2.3.   Oxygen flow

We recall that the oxygen concentration at a point $x \in \Omega$ and time $t \geq 0$ is denoted by $\rho ( x, t)$. The oxygen may be either convected by the blood flow, diffused through the tissue or consumed. Thus, we assume that $\rho$ is described by the following convection-diffusion-reaction equation

where

The left-hand side of (4) together with (5) is the convection-diffusion part. The first term of (5) is the blood velocity $\textbf{u}$ and models the convection of oxygen by the blood flow. The second term is a nonlinear diffusion with diffusivity matrix $\textbf{D}(x, t) \frac{\rho(x, t)}{\rho(x, t) + \widetilde{\rho}}$. For small values of the oxygen concentration $\rho \ll \widetilde{\rho}$, this nonlinear diffusion behaves like $- \nabla_x \cdot ( \frac{\textbf{D}}{2 \widetilde{\rho}} \nabla_x \rho^2)$ and leads to the porous medium equation. In particular, it maintains $\textbf{v}$ finite (and actually equal to $0$) where $\rho = 0$. On the other hand, for high oxygen concentrations $\rho \gg \widetilde{\rho}$, it behaves like $- \nabla_x (\textbf{D} \nabla_x \rho)$, which leads to the standard linear diffusion equation. The diffusivity $\textbf{D}(x, t)$ is a symmetric nonnegative matrix which, like the hydraulic conductivity $\textbf{K}$, depends on the presence of capillaries and consequently, on space and time. It will be specified in Section 2.5. The right-hand side of (4) is a reaction term that accounts for oxygen consumption by the tissue. Formula (6) is the classical Menten-Michaelis reaction term which is widely used for biological reactions. The consumption term $\beta(\rho) \rho$ is a linear function of the oxygen concentration at low oxygen concentration and models the increase of the number of cells able to consume oxygen (through for instance, cell proliferation). The consumption term $\beta(\rho) \rho$ saturates at large oxygen concentration because the number of cells and reaction sites able to consume oxygen cannot grow indefinitely and saturate to a maximal value. The Michaelis constant $K_m$ is the oxygen concentration at which the consumption rate reaches half of its maximum value and $\beta_{\text{sat}}$ is the maximal consumption rate. Eq. (4) must be supplemented with initial and boundary conditions which will be given in Section 2.6.

We comment on the chosen values for the parameters $\widetilde{\rho}$, $\beta_{\rm{sat}}$ and $K_m$ (see Table 1). All quantities relating to oxygen densities are expressed relative to the oxygen density in the oxygenated tissue $\rho_0$ which is a boundary condition for (4) (see Section 2.6). The particular value of $\rho_0$ is irrelevant. Indeed, since the model is two-dimensional, it can be viewed as representing a three dimensional volume of certain thickness where all quantities are independent of the third coordinate. Thus, $\rho_0$ is equal to the actual volumetric oxygen density in the oxygenated tissue times the thickness of the considered volume and that thickness is arbitrary. The chosen value of $\rho_0$ is thus purely arbitrary. The value $\widetilde{\rho}$ is chosen to be $10\%$ of $\rho_0$ meaning that we switch to nonlinear diffusion when the oxygen density is $10\%$ lower than that of the oxygenated tissue. We could not find any documentation of this effect in the literature and this value appeared as a sensible choice. As for the oxygen consumption, we refer to [75] where the rate of consumption of oxygen was measured in a rabbit eye. This rate was measured to be of the order of $10^{-5}$ ml(O$_2$) (ml(tissue) s)$^{-1}$, while the O$_2$ volume fraction was of the order of $4 \, 10^{-3}$ ml(O$_2$) (ml(tissue))$^{-1}$ which gives a consumption rate of $14\%$ min$^{-1}$. However, we have chosen the much lower value $2\%$ min$^{-1}$. If we match this value to the maximal consumption rate of formula (6), this gives us the value $\frac{\beta_{\rm{ sat}}}{K_m} = 0.02$ min$^{-1}$ which we can deduce from Table 1. The reason for choosing this low value lies in the approximations made in the model. Indeed, as already mentioned above, we use the oxygen gradient as a proxy for the VEGF gradient, which we do not include in the model. However, for a physiological value of the consumption rate, oxygen is consumed too quickly and there is not enough oxygen left to trigger the formation of capillaries (see Section 2.4). By lowering the value of oxygen consumption, we allow for enough oxygen to remain and fuel capillary growth. A physiologically realistic value of the oxygen consumption should be usable when we include VEGF in future developments of the model. Now, we have chosen $K_m$, the oxygen concentration which halves the consumption rate, to be equal to half the oxygen concentration in the oxygenated tissue. This is consistent with values given in [13] which suggest that the Michaelis constant is of the same order of magnitude as the oxygen concentration in the oxygenated tissue.

2.4.   Capillaries

2.4.1.   Capillary element

The fundamental assumption of the model is to consider capillary elements as directional information of the flow, rather than describing the exact geometry of the network. Biologically capillary elements could represent one or a group of EC depending on the size chosen for the capillary element. We represent each capillary element as a rod of fixed length $L_c$ and width $w_c$, defined by the position of its centre $\textbf{X} \in \Omega$ and direction $\boldsymbol{\omega} = ( \cos \theta, \sin \theta)$ where $\theta \in [0, 2 \pi)$ (see Fig. 1). We assume that each capillary element influences the local direction of the blood flow as well as the diffusivity properties of the tissue (see Sect. 2.5). The capillary network is formed by superposing all the capillary elements and its connectivity is not imposed at the particle level but is recovered through an averaging process at the continuum level detailed in Section 2.5. We have taken the values $L_c = 15 \, \mu$m and $w_c = 4 \, \mu$m (see Table 1) as these correspond to the typical size of human EC [23].

We assume that capillary elements are created in response to environmental cues. This process could represent both de novo formation of new capillaries as in vasculogenesis or recruitment of existing EC that are attracted to the location where this creation happens. We do not model the process of EC migration (by contrast to many other works in the literature) for the sake of modelling simplicity. Likewise, once created, capillary elements remain immobile, until they are eventually destroyed. Thus capillary element dynamics reduces to creation and destruction processes. We assume three different creation processes and one removal process according to the nature of the environmental cues considered. They are all described by spatio-temporal Poisson processes whose intensities depend on these environmental cues. These processes are:

● Creation by oxygen concentration gradient (detailed in Section 2.4.2): when the norm of the gradient of oxygen concentration at some location exceeds a certain threshold, this signals a fast drop of oxygen concentration and a potentially hypoxic region at some distance along the direction of the gradient. Then, the probability of capillary creation is increased. The newly created capillary element, directed along the direction of the gradient, will contribute to bring more blood and eventually more oxygen, feeding the hypoxic region. However, if the oxygen concentration is already high, there is no need for more blood inflow and the capillary creation is inhibited. As already noted, this mechanism is a proxy for the creation of new capillaries in response to VEGF gradients, while we ignore VEGF in this model.

● Creation by reinforcement of small capillaries (detailed in Section 2.4.3): to stabilize the newly created capillary elements, we have introduced a reinforcement mechanism. As long as the norm of the blood velocity is smaller than a threshold value and the oxygen concentration lies between two bounds, we trigger the formation of new capillary elements directed in the direction of the blood velocity. The bound on the velocity is because we want to avoid the blood flow going too fast, and 'missing' some nearby hypoxic regions. The bounds on the oxygen concentration is because there is no point in reinforcing capillaries if no or too much oxygen is already flowing in them.

● Creation by Wall Shear Stress (WSS) (detailed in Section 2.4.4): WSS is known to play a fundamental role in the development of new capillaries [36,37,38,49,73]). A fluid flowing in a duct exerts a shear stress on its walls, called WSS. There is evidence of a mechanosensing mechanism which triggers capillary sprouting when WSS exceeds a certain value. WSS is expressed as the maximal eigenvalue of the viscous stress tensor. In the model, when at a given location WSS exceeds a threshold value, the creation of a new capillary element in the normal direction to the velocity is favored. The choice of this direction is motivated by experimental observations [38] where it seems favored. But this point will be scrutinized in future work.

● Capillary pruning (detailed in Section 2.4.5): we assume that capillary elements can be destroyed if the tissue hydraulic conductivity is already large. Indeed, the creation of new capillaries increases the tissue hydraulic conductivity (see Section 2.5) but on the other hand, the conductivity should not be locally too large, to avoid channeling all the blood flow in a localized region and leaving the rest of the domain hypoxic. We use the Frobenius norm $\gamma: = \| \textbf{K}\|$ (defined in Section 2.4.5) of the hydraulic conductivity tensor $\textbf{K}$ to sense the conductivity of the tissue. We trigger the removal of capillary elements when $\gamma$ exceeds a threshold value $\gamma^*$ and the removal probability is then a quadratic function of $\gamma-\gamma^*$. The choice of the quadratic function is motivated by [31] where it is assumed that there is a cost for the maintenance of the network which is a quadratic function of a quantity which plays a similar role as the hydraulic conductivity used here. In [31], such a nonlinear cost is needed for the appearance of capillaries and we wanted to test if it is similar here.

In the following sections, these four processes are described in detail.

2.4.2.   Capillary creation triggered by oxygen gradient

We recall that $\rho$ denotes the oxygen concentration in the tissue. We assume that the creation of capillary elements is given by a spatio-temporal Poisson process with intensity

where

and $\nu_c^*$, $L_0^c$, $h_c$, $\rho_s$, $h_s$, $\rho^*$ are positive parameters. The function $g = g( \rho, \nabla_x \rho)$ defined in (9) represents a smoothed logarithmic sensing function of the oxygen gradient $\nabla_x \rho$ where the parameter $\rho^*$ is introduced in order to avoid singularities where $\rho = 0$. The function $\psi$ defined in (8) is a sigmoid. It operates as a switch which is turned on or off whether $z$ is positive or negative respectively. It smoothly transitions from 0 to 1 which introduces uncertainty in the turn on/off of the switch, with a mushy zone located between the values $z = -1$ and $z = 1$. For instance, consider the first activation function

on the right-hand-side of (7). The parameter $h_c$ prescribes the width of the fuzziness region of the activation function. We set $h_c = 0.1$ which represents a $\pm 10\%$ uncertainty around the value $1/L_0^c$ (see Fig. 2(A)). Therefore, (10) turns on the creation process when the logarithmic sensing $g( \rho, \nabla_x \rho)$ of the oxygen gradient exceeds a threshold $1/L_0^c \pm 10\%$. The second sigmoid on the right-hand-side of (7) follows the same rationale. It activates capillary creation when $\rho$ is less than the threshold value $\rho_s$ up to a fuzziness of order $\pm 10 \%$ (taking $h_s = 0.1$) (see Fig. 2(B)). This factor prevents the formation of more capillary elements if the oxygen concentration is high enough. The parameter $\nu_c^*$ is the intensity of the creation rate when both switches in (7) are totally on. If a capillary element is created at the position $\textbf{X} \in \Omega$ then it is oriented towards the oxygen gradient

which represents the optimal direction for the spreading of the oxygen.

We have already commented on the choice of $h_c = h_s = 0.1$ corresponding to $\pm 10 \%$ uncertainties on the threshold values. The constant $\nu_c^*$ is the parameter of the Poisson process where all the switches are on i.e. where the two factors involving the function $\psi$ in (7) are equal to $1$. It corresponds to the fastest rate of capillary creation where large oxygen concentration gradient and low oxygen concentration signal a strongly hypoxic region. It is difficult to obtain such data from the experimental literature and we have resorted to trial and error with the model. The choice made, $\nu_c^* = 0.05$ min$^{-1}$ $\mu$m$^{-1}$ can be best understood when relative to a capillary element area $S_c = L_c \, w_c = 60 \mu$m$^2$. Then $\nu_c^* S_c = 3$ min$^{-1}$ meaning the maximal creation rate of capillary elements due to oxygen concentration gradient is one new capillary element every $20$ seconds per surface of a capillary element. This value sounds too high and was chosen on the basis of numerical constraints. This time-scale issue will be discussed in greater details in Sect. 3. We have taken the oxygen gradient length threshold $L_0^c$ which triggers the formation of new capillaries as half the size of an EC ($L_0^c = 8 \, \mu$m) as concentration sensing should be done at the cell level. The quantity $\rho^*$ is just a regularizing constant and has been chosen as $10\%$ of the oxygen concentration in the oxygenated tissue $\rho_0$ consistently with the choice done for the other regularizing constant $\widetilde{\rho}$. The threshold oxygen concentration beyond which the capillary creation is turned off has been chosen equal to the oxygen concentration in the oxygenated tissue $\rho_0$ which means that for concentrations above $\rho_0$ the tissue is assumed fully oxygenated. The numerical values of parameters $\nu_c^*$, $L_0^c$, $h_c$, $\rho_s$, $h_s$, $\rho^*$ are summarized in Table 1.

2.4.3.   Capillary creation by reinforcement

In order to enhance the creation of small branches we introduce a spatio-temporal Poisson process with intensity function

where $\nu_f^{*}$, $\bar{u}$, $h_f$, $\underline{\rho}$ and $\bar{\rho}$ are positive parameters and $\textbf{u}$ is the blood velocity. The activation function $\psi$ is defined in (8). Following the same rationale as in Sect. 2.4.2, the intensity (11) promotes the creation of capillary elements when the magnitude of the blood velocity $|\textbf{u}|$ is below the threshold parameter $\bar{u}$ up to a fuzziness region of $\pm 10\%$ set by the parameter $h_f = 0.1$. We also impose a cut-off when the oxygen concentration $\rho$ is below $\underline{\rho}$ or above $\bar{\rho}$ up to an uncertainty of $\pm 10\%$ set by the same parameter $h_f$. If a capillary element is created at a point $\textbf{X}$ its direction is set to be the normalized velocity vector at $\textbf{X}$, i.e.

Indeed, this choice is best to reinforce the flow in the direction of the blood velocity $\textbf{u} ( \textbf{X})$.

The parameter $\bar{u}$ was estimated by taking it slightly above the typical blood velocity in the tissue in the absence of capillaries. This velocity can be estimated through Darcy's law using the hydraulic conductivity of the tissue in the absence of capillaries $\kappa_h = 400 \, \mu$m$^2$ min$^{-1}$ mmHg$^{-1}$ (see Section 2.5 for a justification of this value). With a pressure difference of about $20$ mmHg and a typical tissue length of $1$ mm (see section 2.6), such velocity is $8 \, \mu$m/min. We have chosen $\bar u = 20 \, \mu$m/min. The values that we chose for $\underline{\rho}$ and $\bar{\rho}$ are, respectively, $10\%$ and $50\%$ of the the oxygen concentration $\rho_0$ injected to the tissue (see Sect. 2.6). Indeed below $0.1 \, \rho_0$ the oxygen density is too small and capillaries must first be created by other creation mechanisms. Above $0.5 \, \rho_0$ the oxygen density is close to that of the oxygenated region and capillaries do not need any reinforcement. The parameter of the Poisson process $\nu_f^{*}$ corresponding to the maximum intensity when all switches are on is set to be $1/5$ of the corresponding intensity for the creation process by oxygen gradient. Indeed, once capillaries are created, the threat posed by hypoxia on the tissue is reduced and the reinforcement mechanism can be slower, but still within an order of magnitude of the faster process of creation by oxygen gradient. All these parameters are summarized in Table 1.

2.4.4.   Capillary creation by Wall-Shear-Stress

For a given point $x$ we consider the deviatoric stress tensor at $x$ which, thanks to the incompressibility condition for the blood flow (2), reduces to

where $( \nabla_x \textbf{u})_{ij} = \partial \textbf{u}_j / \partial x_i$ is the tensor gradient of the vector field $\textbf{u}$, $A^T$, where $A$ is a $2 \times 2$ matrix, denotes the transpose of $A$ and $\mu$ is the dynamic viscosity of the blood. The $2 \times 2$ tensor $\sigma (x)$ is nonnegative symmetric and traceless. Hence it has eigenvalues $\pm \lambda$ with $\lambda \geq 0$. The creation by WSS is given by the spatio-temporal Poisson process with intensity function

where $\nu_w^*$, $\lambda^*$ and $h_w$ are positive parameters and the activation function $\psi$ is defined in (8). Following the same rationale as in Sect. 2.4.2, capillary elements are created if $\lambda$ is bigger than the threshold parameter $\lambda^*$ up to an uncertainty given by $h_w$ (also set to be $0.1$) and $\nu_w^*$ is the intensity when the switch is totally on. If a capillary element is created at a point $\textbf{X}$ its orientation is taken to be

where $\textbf{u}^\perp ( \textbf{X})$ denotes the rotation of $\textbf{u} ( \textbf{X})$ by $90^\circ$ in the counterclockwise direction. We chose this branching angle following the experimental data presented in [37] and [38]. The value of $\lambda^*$ was taken from [56] and the value of the dynamic viscosity of the blood $\mu$ at $37^\circ$C was taken from [21]. Following a set of numerical trials, and owing to the literature documenting the importance of this mechanism, the value of $\nu_w^*$ was chosen large and set $6$ times larger that the rate of creation by oxygen gradient. See Table 1 for a summary of the numerical values of these parameters.

2.4.5.   Capillary pruning

In order to impede excessive concentrations of capillary elements we remove them following a spatio-temporal Poisson process with intensity defined by

where $\nu_r^*$, $\gamma^*$ are positive parameters, $\gamma$ is the Frobenius norm of $\textbf{K}$ and $z_+ = \max \{ 0, z \}$ for $z$ in $\mathbb{R}$. The Frobenius norm of a $2 \times 2$ matrix with real valued entries $\{a_{ij}\}_{1 \leq i, j \leq 2}$ is given by $\| A \| = (\sum_{1 \leq i, j\leq 2} (a_{ij})^2)^{1/2}$. Formula (14) states that the removal mechanism is only turned on when $\gamma$ exceeds the threshold value $\gamma^*$. Then for $\gamma \geq \gamma^*$ the removal intensity increases quadratically. The quantity $\nu_r^*$ is the intensity of the process for $\gamma = 2 \gamma^*$. The choice of the parameter $\gamma^*$ was estimated to be $5$ times the hydraulic conductivity of individual capillary elements (estimated in Section 2.5) meaning that we do not allow more than $5$ individual capillary elements to superimpose at the same place (superposition of capillary elements modelling broader capillaries). Indeed, larger number of superposed capillaries would model larger vessels for which the Darcy approximation and consequently the whole model would lose validity. We want to promote a fast removal rate when $\gamma$ exceeds the threshold $\gamma^*$. So, the value of $\nu_r^*$ is chosen to be such that the lifetime of a capillary element when $\gamma = 2 \gamma^*$ is equal to 2 seconds which is comparable to the maximal creation rate of capillaries by WSS, equal 1 every 3 seconds per capillary element area (see Section 2.4.4). Finally, the pruning rate being a quadratic function of $\gamma - \gamma^*$ reflects the nonlinear loss rate in the continuum model of network formation of [31]. This feature will be commented in Section 4. We refer to Table 1 for a summary of the parameter values.

2.5.   Tissue

The tissue is a mixture of cells, ECM and capillary vessels. It behaves like a porous medium (see [74]) for the blood and oxygen fluids. It is characterized by its hydraulic conductivity tensor $\textbf{K}$ for the blood flow and its diffusivity tensor $\textbf{D}$ for the oxygen. We assume that in the absence of capillaries, the tissue is an isotropic homogeneous medium characterized by scalar background hydraulic conductivity $k_h$ and oxygen diffusivity $\Delta_h$. Each capillary element $( \textbf{X}_k, \boldsymbol{\omega}_k)$ induces a local increase of the tensor $\textbf{K}$ by the quantity $\kappa \, (\boldsymbol{\omega}_k \otimes \boldsymbol{\omega}_k) \, \chi_{S_k}(\textbf{X})$. The positive scalar $\kappa$ is the conductivity of a capillary element, computed below under the assumption of Poiseuille flow in the capillary. The symbol $\otimes$ denotes the tensor product of two vectors, i.e. if $A = (A_i)_{i = 1, 2}$ and $B = (B_i)_{i = 1, 2}$ are two vectors in ${\mathbb R}^2$, then $A \otimes B$ is the $2 \times 2$ matrix whose entries are $(A \otimes B)_{ij} = A_i \, B_j$. We have the formula $(\boldsymbol{\omega}_k \otimes \boldsymbol{\omega}_k) \nabla_x p = (\boldsymbol{\omega}_k \cdot \nabla_x p) \, \boldsymbol{\omega}_k$ giving that the flow velocity through a capillary element is proportional to the directional derivative of $p$ in the direction of $\boldsymbol{\omega}_k$ and points in the direction of $\boldsymbol{\omega}_k$. This is the feature of a porous medium whose pores are all oriented in the direction of $\boldsymbol{\omega}_k$ and which is impermeable in the orthogonal direction $\boldsymbol{\omega}_k^\bot$ (where $\boldsymbol{\omega}_k^\perp$ represents the counterclockwise rotation of $\boldsymbol{\omega}_k$ by $90^\circ$). The increase of hydraulic conductivity is limited to the domain occupied by the capillary element, which is reflected by the factor $\chi_{S_k}(X)$ where $S_k$ is the capillary element domain, i.e. the rectangle defined by

For a subset $S$ of ${\mathbb R}^2$, the function $\chi_S(\textbf{X})$ is the indicator function of $S$ i.e. it takes the value $1$ if $\textbf{X} \in S$ and $0$ otherwise. The contribution of each capillary element is summed up with those of the other capillary elements and with that of the bulk tissue given by $k_h \textbf{I}_2$ (where $\textbf{I}_2$ is the $2 \times 2$ identity matrix and reflects the fact that the bulk conductivity is isotropic). The same phenomenology holds for the oxygen diffusivity tensor $\textbf{D}$ with a local increase of the diffusivity due to the $k$-th capillary element by the quantity $\Delta (\boldsymbol{\omega}_k \otimes \boldsymbol{\omega}_k) \chi_{S_k}(\textbf{X})$ where $\Delta$ will be estimated below.

The resulting formula for the hydraulic conductivity and diffusivity tensors of the tissue are

The first component on the right-hand-side of (15) is the background conductivity of the tissue in the absence of capillaries and the second one combines the influence of all capillary elements whose domains contain $\textbf{X}$ (see Fig. 3). So, each capillary element $(\textbf{X}_k, \boldsymbol{\omega}_k)$ containing $\textbf{X}$ in its domain exerts a local bias by facilitating blood flow along $\boldsymbol{\omega}_k$. The same phenomenology holds for the oxygen diffusivity given by (16).

The value of $\kappa_h$ in (15) was taken from [74], whereas $\kappa$ was estimated supposing that the flow within a capillary follows Poiseuille's law. It relates the pressure drop $\Delta p$ in a duct of length $L_c$ and radius $w_c/2$ to the velocity $|u|$ of the flow by the relation $\frac{\Delta p}{L_c} = \frac{8 \mu |u|}{(w_c/2)^2}$ where $\mu$ is the dynamic viscosity of the fluid. With Darcy's law within a capillary element $|u| = \kappa \frac{\Delta p}{L_c}$, we get $\kappa = \frac{w_c^2}{32 \mu}$. With the values of $w_c$ and $\mu$ from Table 1, we get $\kappa \approx 10^6 \, \mu$m$^2$ (mmHg min)$^{-1}$. In fact this value brings too large stiffness to the elliptic problem (3) and we settle to a value about $10$ times smaller equal to $0.8 \, 10^5 \mu$m$^2$ (mmHg min)$^{-1}$. The value $\Delta_h$ deduced from from [76] is about $0.5 \, \mu$m$^2$ min$^{-1}$. We take a slightly larger value of $10 \, \mu$m$^2$ min$^{-1}$ owing to the large variability of this parameter according to the type of tissue. This larger value allows for a more uniform distribution of oxygen in the tissue. The diffusivity of oxygen in a vessel is enhanced but not as much as the hydraulic conductivity because the fluid molecules still constitute obstacles against the diffusion of the gaseous oxygen. We assume for our model that $\Delta = 20 \, \Delta_h$. Again, the entire set of parameter values is reminded in Table 1.

2.6.   Geometry and boundary conditions

2.6.1.   Geometry

For the sake of computational efficiency, we consider a 2D rectangular domain $\Omega = [0, L_x] \times [0, L_y]$. Given a capillary-free tissue at initial time, we want to explore how the network appears. We assume that there is a highly oxygenated tissue on the left border with high blood pressure bringing oxygenated blood into the tissue. We analyze two different geometrical conditions distinguished by their boundary conditions.

Case 1. In the first domain $\Omega_1$ depicted in Fig 4 (A), the blood inflow region on the left-hand boundary is surrounded by a lower blood pressure region, such as a region of uptake of blood by the venous system, modelled by the top, bottom and right-hand-side boundaries. The blood source is located in the middle of the left boundary and has small extent. The rest of the left boundary is impermeable. This case models the spread of capillaries around a blood vessel in a cross-section of the tissue normal to the blood vessel. Only half of the environment of the blood vessel is taken into account since, at least in a homogeneous tissue, the capillaries will spread in a similar fashion on the other side of the left boundary. We take domain sizes equal to $L_x = 1$ mm and $L_y = 2$ mm. This corresponds roughly to a domain which extends up to a distance of $1$ mm to the blood source, which is a typical distance covered by capillaries in wound healing for instance. The extent of the blood and oxygen sources in the middle of the left boundary is $100 \, \mu$m which corresponds to the diameter of a small blood vessel.

Case 2. In the second one denoted by $\Omega_2$ and shown in Fig. 4 (B), a highly oxygenated high blood pressure region stands on the left-hand boundary. It faces a low blood pressure region corresponding to the right boundary. The top and bottom boundaries are assumed periodic: they represent a tissue that extends in the vertical direction on large distances in a somewhat similar fashion and for which only a representative fraction is simulated. The creation of a new capillary is triggered by an inflow of oxygen on a small portion of the left boundary. This case depicts the spread of a blood capillary longitudinally, i.e. when a cross-section of the tissue is made in a plane containing the blood vessel making the blood source. This situation also mimics the geometry of a corneal micropocket angiogenesis assay [28]. We take domain sizes equal to $L_x = 2$ mm and $L_y = 1$ mm, which is the order of magnitude encountered in such experiments. The oxygen sources in the middle of the left boundary extends along $100 \, \mu$m similar to the previous case.

In both cases, the difference of pressures between the high and low blood pressure regions produces a flow of the interstitial fluid which triggers the formation of a new vascular network. The bulk of the domains are constituted of a tissue assimilated to a porous medium as described in Sect. 2.5.

2.6.2.   Boundary conditions on the pressure and boundary/initial conditions on the oxygen concentration

Since we initially suppose a capillary-free tissue, we impose for consistency that the oxygen concentration is initially zero, i.e. for both geometries $\Omega_1$ and $\Omega_2$ we impose $\rho ( x, 0) = 0$ for all $x$. For the boundary conditions, we distinguish between the two cases depicted in the previous section.

Case 1. The labeling of the different parts of the boundary of $\Omega_1$ is given in Fig. 5. We impose the following boundary conditions for the pressure:

and for the oxygen concentration:

where $\widehat{n}$ is the outward unit normal vector to the boundary. Boundary $\Gamma_{1, D}$ has length $100 \, \mu$m and is located at the centre of the left boundary of the rectangle $\Omega_1$ i.e. it starts at the height $L_{min} = \, 950 \, \mu$m and ends at $L_{max} = \, 1050 \, \mu$m. We set $p_0 > p_1$. Condition (17) expresses that $\Gamma_{1, D}$ is the high blood pressure region with pressure $p_0$ (for instance a blood vessel). $\Gamma_{1, D}$ is also the source of oxygen and $\rho_0$ stands for the oxygen concentration in this highly oxygenated region, hence (20). The subset $\Gamma_{2, D} \cup \Gamma_{3, D} \cup \Gamma_{4, D}$ of the boundary is the low blood pressure region, hence the Dirichlet condition (18) with pressure $p_1$. It surrounds the high pressure source to mimic the conditions depicted in Fig 4 (A). The uptake of blood by the venous system takes place in this region and the oxygen flows out of the domain. Along this boundary, we assume a zero normal gradient of the oxygen concentration (21), meaning (with (5)) that the normal oxygen flow velocity $v \cdot \widehat{n}$ is equal to the normal blood flow velocity $u \cdot \widehat{n}$. We expect the latter to be large and positive (i.e. leaving the domain). Thus, oxygen will exit the domain freely without any possible inflow. So, (21) expresses that the incoming flux of oxygen is zero. The Neumann boundary condition (19) for the pressure along $\Gamma_{1, N} \cup \Gamma_{2, N}$ is a rigid wall condition stating that the normal blood velocity to the wall is zero. We also assume that there is no oxygen, leading to (22). Because the problem is two-dimensional, the quantity $\rho_0$ is a surface density. As explained in Section 2.3, its numerical value is irrelevant because it depends on the height of the tissue in the third dimension, which is arbitrary. The boundary conditions on the pressure are those that prevail between the inflow and outflow of the capillary system as measured in [81]. The parameters are given in Table 1.

Case 2. We refer to Fig. 6 for the labeling of the various parts of the boundary of $\Omega_2$ and consider the following boundary conditions for the pressure:

where $e_2$ is the unit vector in the vertical direction i.e. $e_2 = (0, 1)$. For the oxygen concentration, the boundary conditions are set to:

Again, we impose $p_0 > p_1$ and the Dirichlet condition (23) sets the pressure along the left boundary $\Gamma_{1, D}^{\mbox{ top}} \cup \Gamma_{1, D}^{\mbox{ mid}} \cup \Gamma_{1, D}^{\mbox{ bot}}$ to the high pressure value $p_0$, while the Dirichlet condition (24) sets it to the low value one $p_1$ on the right boundary $\Gamma_{2, D}$. Condition (25) simply states that the values of the pressure on $\Gamma_{1, \mbox{ per}}$ must be equal to those of $\Gamma_{2, \mbox{ per}}$ expressing the periodicity of the solution in the vertical direction. To trigger the formation of a capillary we assume a large oxygen concentration along $\Gamma_{1, D}^{\mbox{ mid}}$, hence the Dirichlet condition (26) setting the oxygen concentration to $\rho_0$ on this boundary. Along the remaining part of the left boundary $\Gamma_{1, D}^{\mbox{ top}} \cup \Gamma_{1, D}^{\mbox{ bot}}$, we assume no oxygen is present, hence the Dirichlet condition (27). We assume that $\Gamma_{1, D}^{\mbox{ mid}}$ has length $100 \, \mu$m and is located at the centre of the left boundary i.e. it starts at the height $L_{min} = \, 450 \, \mu$m and ends at $L_{max} = \, 550 \, \mu$m. Like for the pressure, condition (29) along $\Gamma_{1, \mbox{ per}}$ expresses the periodicity constraint on $\rho$ in the vertical direction. Finally, (28) is an outflow condition for the oxygen concentration along the right boundary $\Gamma_{2, D}$. Its interpretation is the same as for (21) and we refer to the previous case for a detailed explanation.

2.7.   Numerical method

In this section, we provide a brief summary of the numerical methods used and refer to the Appendices A, B, C for a detailed description. A synopsis of the numerical treatment of the problem is given in Algorithm 1 below and the choices for the numerical parameters are summarized in Table 2.

The elliptic equation describing blood flow (1)-(2) is approximated using a $\mathcal{Q}_1$ finite element method on a rectangular spatial grid [7] of uniform mesh sizes $\Delta x$ and $\Delta y$ in the horizontal and vertical directions respectively. Finite-element methods are regarded as methods of choice for the resolution of elliptic problems with complex boundary conditions. A $\mathcal{Q}_1$ finite element method provides an approximation of the solution by polynomials linear in each variable on quadrangular numerical cells. We favored $\mathcal{Q}_1$ elements over more classical $\mathcal{P}_1$ elements (which provide approximations by globally linear polynomials on triangular numerical cells) to avoid too much numerical directional bias. For all numerical simulations shown in this paper, otherwise stated, we set $\Delta x = \Delta y = 1.25 \, \mu$m, which resolves the scale of capillary elements. If we take finer mesh sizes we do not note any significant changes in the network formation.

The diffusion-advection-reaction equation (4) describing the evolution of the oxygen concentration $\rho$ is discretized using the Smoothed Particle Hydrodynamics (SPH) methodology [47]. The choice of a SPH method is motivated by the initially zero oxygen concentration. Dealing with regions of zero concentration with classical grid methods, such as finite volume methods, is delicate due to the high risk of appearance of negative values, which often leads to a breakdown of the simulation. The SPH method is not subject to this risk and is the method of choice for problems involving jets or injections of gases in vacuum. SPH is a classical method whose accuracy has been practically assessed, and thoroughly studied theoretically (see e.g. [16,42,66,78]).

It consists of approximating the function $\rho(x, t)$ by a sum of $N$ Dirac deltas located at positions $\{\textbf{Y}_\ell(t)\}_{\ell = 1, \ldots, N}$ and of equal and time-constant masses $m$. Here, since the oxygen concentration scale is irrelevant (see Sect. 2.3), we take $m = 1$. These particles move in space with a velocity $\textbf{V}_\ell(t)$ given by formula (5) evaluated at $\textbf{Y}_\ell(t)$. However since (5) involves $\nabla_x \rho$, it does not make sense with $\rho$ equal to a sum of Dirac deltas. To make sense of it, we approximate all Dirac deltas by smooth regularizations. Their role is to spread the Dirac deltas over regions of finite extent $\eta$. We take $\eta = 5 \, \mu$m which spreads the Dirac deltas over $10 \, \mu$m, i.e. approximately the scale of a capillary. The process of representing diffusion of a cloud of particles by their convection along their regularized concentration gradient is known as the diffusion velocity method, introduced in [17]. It stands as an alternative to stochastic methods for the treatment of diffusion which introduces lower noise uncertainty. The regularization kernel is taken to be the poly6 kernel (see Appendix) which has well documented reliability [48]. To account for the consumption term (right-hand side of (4)), we use a splitting method. We first solve for the left-hand-side of (4) by the SPH method and then we use a stochastic death process to discretize the Michaelis-Menten reaction term.

The simulations of the spatio-temporal Poisson point processes (7), (11), (13) and (14) are performed using a classical acceptance-rejection method (see for instance [27]). For each creation process and each discrete time $t$, we pick $N_c$ points $x$ in the domain independently and uniformly and compute the value of the corresponding Poisson parameter $\nu(x, t)$. We can then estimate the probability of creation of a capillary element by this process. The value $N_c = 10^5$ was selected by trial and error. Larger values of $N_c$ do not improve the quality of the solution but require longer computer time. Lower values of $N_c$ deteriorate the quality of the results. For the removal process, we just loop over all the existing capillaries and compute the value of the removal Poisson rate at their location. Again, this allows us to estimate the probability of removal of the considered capillary. We could have accelerated the treatment of these processes by using Metropolis-Hastings or Monte-carlo Markov Chain methods but this step was so quick compared with the resolution of the elliptic equation that it did not seem to be worth the effort.

We implemented the code in Fortran 95 and used the Mersenne twister generator [43] for the generation of random numbers. The plots were generated in Gnuplot. The simulations were conducted on two workstations, one with two Intel Xeon E5-2637v3 quad core processors clocked at 3.5 GHz with 12 MB of on-CPU cache and the other one with two Intel Xeon Gold 5122 quad core processors clocked at 3.6 GHz with 16 MB of on-CPU cache. The time taken for each simulation is in average less than 72 hours for Geometry 1 and for Geometry 2 it takes in average 6 days depending on the creation/deletion mechanisms involved. We summarize the choice of numerical parameters in Table 2.

3.   Results

3.1.   Simulations for geometry $\Omega_1$

In this section we consider Geometry $\Omega_1$ which mimics the branching of a new capillary network from a blood vessel in the plane perpendicular to this blood vessel (see Fig. 4). We explore the influences of the different creation/deletion mechanisms of capillary elements introduced in Sect. 2.4 on the network structure. All the parameters used are those given in Tables 1 and 2.

3.1.1.   Including all capillary creation/deletion mechanisms

We first turn 'On' all the mechanisms described in Sect. 2.4 i.e. creation of capillary elements by oxygen gradient, reinforcement by blood flow and by WSS as well as capillary pruning. Since all these processes are of stochastic nature, different realizations of the same model will not give the same result. However, they are qualitatively similar. The result of one given realization is displayed in Fig. 7. On this figure, the rectangle represents the domain $\Omega_1$ i.e. corresponds to $1$ mm in length in the horizontal direction and $2$ mm in the vertical direction. The positions of the oxygen particles in this domain are represented by red spots and those of the capillary elements by tiny blue rods. As red spots overlay the blue rods, capillary elements lying below the red oxygen particles are present although not seen. For the same realization, the isolines and heatmap of the pressure $p$ in $\Omega_1$ are shown in Fig. 8 and a heat map of the Frobenius norm $\gamma$ of the conductivity matrix $\textbf{K}$ is displayed in Fig. 9 (see Section 2.4.5 for the definition of the Frobenius norm). In both figures, the colour code stands from blue (low value of $p$ or $\gamma$) to red (high values) through green and yellow (intermediate values). The six pictures (A) to (F) in these figures show various stages of the evolution of the network, $2$ min (A), $4$ min (B), $6$ min (C), $8$ min (D), $10$ min (E), $12$ min (F) after initialization. The initial condition, not depicted, corresponds to no oxygen and no capillary element at all, i.e. an empty domain for Fig. 7, regularly spaced concentric pressure isolines emanating from the middle of the left boundary in Fig. 8 and a constant value equal to $k_h$ for Fig 9.

As the pressure drops (see isolines in Fig. 8) between the portion $\Gamma_{1, D}$ of the left boundary (see Fig. 5 for the nomenclature) representing a blood vessel, and the top, right and bottom boundaries $\Gamma_{2, D} \cup \Gamma_{3, D} \cup \Gamma_{4, D}$ where uptake of blood happens, there is a flow of blood from the former to the latter. Oxygen, which can enter the domain through $\Gamma_{1, D}$, is transported by blood flow. Blood flow and oxygen transport trigger the formation of capillaries through the various creation/deletion mechanisms described in Sect. 2.4. This results in the initiation of a capillary network emanating from $\Gamma_{1, D}$ which gradually develops radially towards the surrounding boundaries $\Gamma_{2, D} \cup \Gamma_{3, D} \cup \Gamma_{4, D}$ as seen in Fig. 7. Capillaries are made by the aggregation of many capillary elements along distinctive branches and contribute to a dramatic increase of the hydraulic conductivity along these branches as seen in Fig. 9. In Fig. 7 and 9 we observe spontaneous branch sprouting, the merging of existing branches (forming so called 'anastomoses') and the spontaneous emergence of capillary tortuosity, which are characteristic features of actual vascular networks (see for instance [24] or [49]). Fig. 9 which displays the norm of the conductivity matrix $\textbf{K}$ informs us on the density of capillary elements. In particular, the big red spot of oxygen particles in the trunk of the network seen in Fig. 7 overlays a uniform bed of capillary elements (hidden by the red spot) as indicated by the uniform green color in the same region in Fig. 9. The creation of a capillary network instead of the expansion of a homogeneous cloud of capillary elements is due to an instability akin to the viscous fingering (or Saffman-Taylor) instability in fluids [69]. Such an instability arises when a viscous incompressible fluids with large mobility penetrates a porous medium occupied by another fluid of low mobility. If the two fluids are incompressible, they are separated by a sharp interface which gradually forms indentations or fingers. This instability has been mathematically studied in e.g. [52] and is recognized as a powerful morphogenetic mechanism in biology [8]. Here, the creation of capillaries modifies the mobility of the medium. Hence, blood flowing in a region filled with capillaries behaves like the higher mobility fluid compared with blood flowing in the background medium. Within this analogy, the creation of the blood capillary network can be viewed as a fingering instability.

The pressure seems to remain constant and equal to its boundary value on $\Gamma_{1, D}$ throughout the region covered by the capillary network and which roughly corresponds to its convex hull (see Fig. 8). Indeed, due to the very large hydraulic conductivity along the branches of the network, there is almost no pressure drop between $\Gamma_{1, D}$ and the tip of the branches. Even if the hydraulic conductivity between the branches is much lower, there is no reason for the pressure to drop significantly between these branches. So, all the pressure drop between $\Gamma_{1, D}$ and the outer domain boundary occurs between the boundary of the convex hull of the network (which we will call the 'envelope' of the network) and the outer boundary. As the envelope of the network moves towards the outer boundary, the pressure isolines become closer and closer to each other, indicating a dramatically increasing pressure gradient (see e.g. 8 (F)). When a network branch touches the outer boundary (not depicted here), there is a sudden 'short circuit' of the pressure difference and the flow becomes virtually infinite, which does not make biological sense. So, the model currently cannot describe the last phase of the network expansion, when it connects to the outer boundary. Various mechanisms which may take place when the network transitions from an expansion phase to an established phase are being investigated.

The results show that the network evolves over a time scale of the order of $10$ minutes. This is too fast and experimental results from e.g. [38] suggest that the right time-scale would rather be $10$ hours. So, the network dynamics is between $1$ and $2$ orders of magnitude faster in the model than in reality. The issue there is a numerical issue. There is a time-step limitation for the stability of the SPH resolution of the oxygen concentration equation. It requires that a particle does not move more than a smoothing length $\eta$ during one time step. The particle velocity is roughly that of the blood velocity. A slower capillary creation/deletion dynamics would not dramatically affect this velocity, and thus, the stability constraint. So, the time-step would be roughly similar to that used in the current simulations. But the simulation would have to run for longer times, with 10 to 100 times more time-steps than in the current simulations. This would lead to unaffordable simulation times. Resorting to implicit particle methods would be a possible cure. Some implicit particle methods have been proposed in the context of plasma physics [40] but they are cumbersome and would require to be adapted and validated for the present case. Alternately, one could use the stationary form of the oxygen concentration equation. Indeed, the time-scale for equilibration of the oxygen concentration is very fast compared to the speed of network evolution. Therefore, we could safely assume an adiabatic evolution of the oxygen concentration, where it would instantaneously adjust to the new network structure exactly like the blood flow does in the current version of the model. The issue here is the resolution of a degenerate stationary convection-diffusion-reaction equation where the concentration takes zero values in large parts of the simulation domain. There is no established method for this kind of problem and there are indeed potential issues such as non-uniqueness of the solutions. This is why, in the present state, we have used the SPH method and 'accelerated' the evolution of the network to make solutions computable.

3.1.2.   Influence of individual capillary creation/deletion mechanisms

In this section, we attempt to analyze the roles of individual creation/deletion mechanisms. In Fig. 10 we display a set of simulations where the various creation mechanisms are turned 'On' or 'Off'. When a mechanism is turned 'Off', this simply means that the corresponding Poisson parameter $\nu^*$ is set to zero, while when it is turned 'On', $\nu^*$ is set to the value given in Table 1. In Fig. 10 we have kept the deletion (capillary pruning) mechanism (see Sect. 2.4.5) 'On'. If we turn this mechanism 'Off' (not shown in the text but a video of this situation is provided in the supplement), the magnitude of (the Frobenius norm of) the conductivity tensor $\mathbf{K}$ reaches higher values in the middle of the branches, without this having a perceivable influence on the network morphology itself, except for the dynamics being slightly quicker (by about $10 \, \%$) than if the pruning mechanism is turned 'On'.

Fig. 10 takes the form of two binary trees placed back-to-back, the upper one (Ⅰ) corresponding to the reinforcement mechanism turned 'On' and the lower one (Ⅱ) to that mechanism turned 'Off'. Within these trees, each branching corresponds to the choice of one of the two other mechanisms (creation by oxygen gradient or by WSS) being 'On' or 'Off'. The first branching (the topmost one for (Ⅰ) and bottommost one for (Ⅱ)) corresponds to choosing whether creation by WSS is 'On' (this is depicted by a green edge in the tree) or 'Off' (red edge). Then the second branching (the bottommost one for (Ⅰ) and the topmost one for (Ⅱ)) chooses whether creation by oxygen gradient (referred to as '$O2\nabla$' on the picture) is 'On' (greed edge) or 'Off' (red edge). Each leaf of the tree is a typical result of the model in the situation corresponding to the various mechanisms being turned 'On' or 'Off' according to the path followed in the tree. Like in Fig. 7, the picture shows the positions of the oxygen particles (red spots) and those of the capillary elements (tiny blue rods). We refer to the previous section for a detailed description.

First we examine what happens when only one of the creation mechanisms is turned 'On', the other ones being turned 'Off'.

● Creation of capillary elements by oxygen gradient 'On', the other mechanisms 'Off': this corresponds to Fig. 10 (G). We notice a fully developed network. So, creation by Oxygen gradient alone provides enough positive feedback to trigger an instability. Let us describe a possible mechanism for this. First, due to the large hydraulic conductivity in the branches, the blood flow velocity there is large. Neglecting oxygen diffusivity, consumption and time variation of the oxygen concentration which are believed to play minor roles here, Eq. (4) roughly reduces to $u \cdot \nabla_x \rho = 0$, showing that oxygen concentration is nearly constant in the branch and equal to that in the inflow boundary $\Gamma_{1, D}$, i.e. $\rho_0$. Ahead of the branch tip, the hydraulic conductivity drops, and so does the blood velocity and consequently the oxygen concentration, which cannot be transported by the blood flow any more. So, the creation mechanisms detects a large oxygen gradient at this place, oriented in the direction along the branch. Therefore, it triggers the creation of a new capillary element ahead of the branch tip in the direction of the branch, thereby increasing its length. However, there are also oxygen gradients across the boundaries of the branches. This results in capillary creation across the branches and contributes to branch widening and the initiation of many new small branches across the main branch. Indeed, we can see from Fig. 10 (G) that the trunk and main branches of the network are thick and that branches terminate in a sponge-like structure. Several anastomoses can also be seen. It takes about $20$ min to observe a fully developed network.

● Creation of capillary elements by WSS 'On', the other mechanisms 'Off': this corresponds to Fig. 10 (F). Again, the WSS mechanism alone is able to generate a full network. Here, the mechanism of network formation is a bit more mysterious as WSS creates capillaries normal to the main blood flow. So, a new capillary created at a branch tip will pursue the branch in the direction normal to the existing branch. However, the next capillary will be created in the direction normal to the direction of the previous capillary and thus again in the direction of the main branch. In this way the branch advances a bit like a sailing boat forced to tack in headwind. Indeed, the observation of the videos shows that the motion of the branch tip is undulatory while it is straighter in the previous case. Behind the tip, the branch consolidates into larger and straighter channels than in the previous case. The branch width is very large at the trunk and decreases quickly in the secondary branches. Obviously WSS also favors sprouting but most of these branches are not pursued because they are not used by blood flow which finds an easier way towards the outer boundary through the main branches. The formation of the network is faster than in the previous case, with a fully developed network after only $12$ min.

● Creation of capillary elements by reinforcement 'On', the other mechanisms 'Off': this corresponds to Fig. 10 (D). Here, we cannot see much of network formation. Rather, we observe the expansion of oxygen and capillary elements in a diffusive way from the inflow boundary $\Gamma_{1, D}$. The reinforcement mechanism does not carry enough positive feedback to trigger a significant instability at least within the $30$ min time of the simulation (compared with the $12$ min which are sufficient for a fully developed network with WSS).

From this analysis, we can infer that the mechanisms that influence the network structure the most are the creation by oxygen gradient and by WSS. The creation by reinforcement seems to play a minor role. Let us examine this question more closely by comparing what happens when we turn the reinforcement 'On' and 'Off' with all other parameters unchanged. This corresponds to comparing homologous pictures in the upper (I) and lower (II) binary trees of Fig. 10, i.e. (A) with (E), (B) with (F) and (C) with (G) (we leave (D) aside because in this case, if we turn reinforcement 'Off', there is no capillary creation mechanism at all and nothing happens).

● Creation of capillary elements by WSS and oxygen gradients both 'On'; comparison between reinforcement mechanism 'On' and 'Off': this corresponds to comparing Figs. 10 (A) and (E). Fig. 10 (A) is the same as Fig. 7 (F) but is reproduced again for the sake of comparison. The network structures in the two plots are fairly different. While in the presence of reinforcement (Fig. 10 (A)) there are thick branches with significant sprouting near their tip, in the absence of it (Fig. 10 (E)) branches are thinner and there is much less sprouting at the tip. In the latter case, there are more thin branches directly branching off the trunk of the network. Also, the network seems less developed in the absence of reinforcement, suggesting that network formation is slightly slower. The reinforcement mechanism consolidates small branches thereby giving them more chances to grow. When reinforcement is absent, branching looks more difficult.

● Creation of capillary elements by WSS 'On' and oxygen gradient 'Off'; comparison between reinforcement mechanism 'On' and 'Off': this corresponds to comparing Figs. 10 (B) and (F). Roughly speaking the same comparisons as in the previous case can be made here, although the difference is less striking.

● Creation of capillary elements by oxygen gradient 'On' and WSS 'Off'; comparison between reinforcement mechanism 'On' and 'Off': this corresponds to comparing Figs. 10 (C) and (G). Here the two network structures are qualitatively almost indistinguishable, except for the extension of the network, which is larger when reinforcement is on, than when it is off.

So reinforcement seems to have little yet perceivable influence on the network structure. In general it gives slightly more regular shape of branches so, we are going to keep it in the next round of comparisons, thereby discarding the lower tree (II) of Figs. 10. Assuming that we include reinforcement and looking only at the upper tree (I), we are going to compare the cases where both WSS and oxygen gradient creation are 'On' to the cases of either one or the other is 'On', i.e. compare Figs. 10 (A), (B) and (C) with one another

● Comparison between creation of capillary elements by WSS and oxygen gradient both 'On' to cases where either one is 'On' and the other is 'Off' (creation by reinforcement 'On'): this corresponds to comparing Figs. 10 (A), (B) and (C). We can readily single out Fig. 10 (C) (where WSS is absent) for the sponge-like structure of the network already noticed above. The WSS mechanism present in Figs. 10 (A) and (B) provides a neater network with a well defined branching structure. This may be related to the way the branch progresses by creating a wider capillary path due to the undulatory motion of its tip as discussed above. On the other hand, in terms of homogeneity of the oxygen perfusion, the structure shown in Fig. 10 (C) seems to be more satisfactory as a wider surface of the tissue seems to have access to oxygen. On the other hand, Figs. 10 (A) and (B) look quite similar. A closer inspection shows that less branches are generated and the generated ones are thicker in the absence of the oxygen gradient mechanism (Fig. 10 (B)) than when it is present (Fig. 10 (A)). Overall, when all mechanisms are present, the network seems to be better balanced, with a gradual decrease of the thickness of the branches when going from the trunk to the periphery. For this reason, we conclude that all mechanisms seem to concur to the formation of a harmonious network.

3.1.3.   Robustness of the simulations with respect to mesh size

We have realized that too coarse a mesh-size for the finite-element method may impair the quality of the results. For all the simulations shown so far, we took $\Delta x = \Delta y = 1.25 \, \mu$m as mesh-size (see Table 2). This mesh-size resolves the scale of each capillary (tubes of length $15 \, \mu$m and width $4 \, \mu$m). By taking smaller mesh sizes we did not notice any significant difference. Coarser mesh sizes on the other hand gave significant differences which we attributed to a bad spatial resolution. In Fig. 11 we display the results of a realization taking mesh size $\Delta x = \Delta y = 0.625$, i.e. half the size of the previous simulations, with all other parameters in Tables 1 and 2 being the same as before. Figs. 11 (D) is to be compared with Fig. 7 (E) roughly corresponding to the same time. We do observe some differences due to stochastic nature of the model, which never provides exactly the same results twice. But qualitatively, the features of the network are comparable, with a gradual decrease of the thickness of the branches from trunk to tip. Hence, we concluded that the mesh size $\Delta x = \Delta y = 1.25 \, \mu$m is a good compromise to resolve the small scales without excessive computer time.

3.1.4.   Influence of the pressure gradient and capillary element size

We first investigate the influence of the pressure gradient. Fig. 12 displays simulation results for three different boundary values for the pressure. Fig. 12 (B) is the same as Fig. 7 (E) while Fig. 12 (A) and (C) are for pressure gradients reduced by $10 \, \%$ and increased by $10 \, \%$ respectively. All the creation/deletion mechanisms of capillary elements are turned 'On'. All the parameter used are those given in Tables 1 (for geometry 1) and 2 except for Figs (A) and (C) where the boundary conditions for the pressure are modified as follows:

● Fig. 12 (A): $p_0$ is replaced by $p_0' = p_1 + 0.9*(p_0-p_1) = 35.4$ mmHg, and $p_1 = 14.6$ mmHg is unchanged;

● Fig. 12 (C): $p_0$ is replaced by $p_0" = p_1 + 1.1*(p_0-p_1) = 40.0$ mmHg, and $p_1 = 14.6$ mmHg is unchanged.

As seen in Fig. 12, the obtained capillary networks are qualitatively similar (but for differences due to the stochastic nature of the model) but their spatial extension increases with increasing pressure gradient. These results reflect that the blood velocity is proportional to the pressure gradient. So, a smaller pressure gradient results in a smaller blood velocity and in a slower expansion of the network (Fig. 12 (A)) compared with the control case Fig. 12 (B). By contrast, an increase of the pressure gradient results in a faster expansion of the network (Fig. 12 (C)). These results also document the stability of the model with respect to the physical parameters. A small variation of the pressure gradient gives rise to a small variation of the network speed of expansion, without affecting the other qualitative features of the network. It gives confidence that the model is well-posed in the Hadamard sense, i.e. it responds continuously (in a probabilistic sense) to a small variation of its parameters.

We next investigate the influence of capillary element size. Fig. 13 displays simulation results for three different values of the capillary length $L_c$. Fig. 13 (B) is the same as Fig. 7 (D) while Fig. 13 (A) and (C) are for capillary lengths divided by $2$ and multiplied by $2$ respectively. All the creation/deletion mechanisms of capillary elements are turned 'On'. All the parameter used are those given in Tables 1 (for geometry 1) and 2 except for Figs (A) and (C) where the capillary length $L_c$ is modified as follows:

● Fig. 13 (A): $L_c$ is replaced by $L_c' = L_c/2 = 7.5$ mm;

● Fig. 13 (C): $L_c$ is replaced by $L_c" = 2 \, L_c = 30$ mm.

We can see in Fig. 13 that the expansion speed of the network increases with capillary size. Fig. 13 (A) obtained with half-sized capillary elements compared to the control case (Fig. 13 (B)) shows a network extension roughly the same as that of Fig. 7 (C) which corresponds to the control case at an earlier time ($t = 6$ min) than the one corresponding to Fig. 13 ($t = 8$ min). Qualitatively the trunk and branches look thicker with a lower number of bifurcations in the small capillary length case than in the control case (compare Fig. 13 (B) with Fig. 7 (C)). Symmetrically, in Fig. 13 (C) which corresponds to double-sized capillary elements compared to the control case (Fig. 13 (B)), the network extension is comparable to an intermediate between Figs. 7 (E) and (F) which corresponds to the control case at a later time ($t \approx 11$ min). However, in Fig. 13 (C) on sees a branch that extends up to the boundary and signals the end of the simulation, which corresponds to time $t = 12$ min in the control case (Fig. 7 (F)). The morphology of the network is also quite different with thinner trunk and branches. It has also a more well-defined network morphology with less "foamy" structures at the tip of the branches. Thus, acting on the capillary length allows us to tune the network morphology from moss-like to tree-llke structures. It demonstrates the ability of the model to generate a wide range of different network-like patterns. A more systematic parametric exploration is left to future work.

3.2.   Simulations for geometry $\Omega_2$

We now consider Geometry $\Omega_2$ which represents the branching of a new capillary network from a blood vessel in the plane parallel to itself (see Fig. 4). It is also relevant for experimental situations such as [28,51]. We refer to Sect. 2.6 for the description of the geometrical setting and the boundary conditions.

A realization of the model with all capillary element creation/deletion mechanisms turned 'On' and with parameters given in Tables 1 and 2 is shown in Fig. 14. The graphics are the same as in Fig. 7, i.e.: (ⅰ) the rectangle represents the domain $\Omega_2$ which is $2$ mm long in the horizontal direction and $1$ mm long in the vertical direction and (ⅱ) the red spots are the positions of the oxygen particles and the tiny blue rods depict the capillary elements. Again, the red spots overlay the blue rods, so there are capillary elements below the red oxygen dots even if they are not apparent. We remind that the boundary conditions along the horizontal boundaries are periodic. Fig. 14 shows six different snapshots of the network evolution, $6.8$ min (A), $13.6$ min (B), $20.4$ min (C), $27.2$ min (D), $34$ min (E), $46.8$ min (F) after the initial time. The same instability as in geometry $\Omega_1$ takes place and contributes to initiate the capillary network from the injection boundary $\Gamma_{1, D}^{\mbox{ mid}}$. Because of the presence of the WSS capillary creation mechanisms, the network branches are very distinct and of fairly constant width between two bifurcations. There is a gradual decrease of the branch width from trunk to branch tip as already observed in Fig. 7. There are several anastomoses. The periodic boundary conditions along the horizontal boundaries prevent the generation of a pressure gradient in the vertical direction. The blood flow is predominantly horizontal, from left to right and does not provide any drive for the development of network branches in the vertical direction. Therefore, most of the network branches develop horizontally. We also observe in Fig. 14 (F) that the early made branches (those bifurcating away from the trunk close to the injection region) do not seem to be used by oxygen particles. This is because most of the blood flows along the central trunk due to its high hydraulic conductivity until it reaches the tip of the trunk. There, the pressure gradient is large because of the short gap between the tip and the right-hand outflow boundary. It outcompetes the low hydraulic conductivity in the gap and produces a large blood flow resulting in a large perfusion in the whole trunk. This large blood flow drags the oxygen particles away from the side branches which, after some time, become depopulated. In actual tissues, unused capillaries are pruned after some time. Therefore, the long term structure of the network would gradually become reduced to the mere trunk itself and would be close to the structures shown in [28] or in [51,55,57]. We note that the time scale is of the order of one hour to see a fully developed network, and is too fast by at least one order of magnitude compared to actual observations. We refer to Sect. 3.1.1 for the rationale behind this discrepancy.

This geometry allows for comparisons with results of the literature using other models. Ref. [15] considers a cell-based model for EC (based on a cellular Potts methodology) moving along the gradient of the VEGF concentration. Motion of EC is also influenced by the ECM (a phenomenon known as haptotaxis) but ECM can be degraded by matrix degrading molecules (MDM) secreted by the EC. VEGF, ECM and MDM concentrations are described by continuum equations. Fig. 3 of [15] shows that motion of EC under VEGF alone (with no consideration of ECM or MDM) generates a single stem. However, when both the ECM and MDM fields are added, the network has a very similar tree-like structure as that shown by our model in Fig. 14. In [77] similar mechanisms as in [15] are considered but the model uses a continuum model for the EC except for the tip cells which are kept discrete (in a similar fashion as in [44], the difference being that [77] uses a phase-field approach). Again, the tree-like structures shown in [77] (such as e.g. Fig. 2 C and D) are quite similar to those of our model in Fig. 14. Ref. [70] once more uses similar mechanisms as [15,77] but relies on a three-dimensional evolving network description of the vascularization. Here also, the results display a tree-like structure with initial stems gradually dividing into smaller branches (see e.g. Fig. 6 of [70]). The three-dimensional nature of [70] does not allow further comparisons with our two-dimensional results. By contrast, [44] produces more complex structures, where initial stems divide in smaller branches but branches initiated from different stems may merge when they become close (see e.g. Fig. 3 of [44]). It seems that [44] favors more mergers than [15,70,77] and than the present model. In summary, different network formation mechanisms, powered by the flow of blood in our model and by chemo- and haptotactism in [15,70,77] produce qualitatively similar network structures. It would be interesting to understand if common underpinning principles are at play which could explain this similarity. Also, it will be difficult to decide which of the network formation mechanism is the biologically relevant one on the basis of sole comparisons with experiments. On the other hand, although similar mechanisms are at play in [44], the qualitative features of the generated networks vary significantly from those obtained in [15,70,77]. This highlight the fact that these models are very sensitive to small details of the modelling principles and the values of the parameters. We should keep this in mind and avoid overstating the abilities (or inabilities) of a model to generate specific structure types.

4.   Discussion

In this model, the emergence of blood capillary networks relies on the positive feedback between blood circulation and oxygen transport on the one hand and capillary formation on the other hand. It is the first of this kind which gives a central role to the circulation of blood. As new capillary branches are formed, they elongate or are reinforced as a consequence of this positive feedback. Both the network topology (i.e. its connectivity) and its geometry (e.g. the branch widths, the branching angles, etc.) are emerging properties, i.e. they are not directly encoded in the capillary creation rules but rather emerge from the interactions between the various entities of the model. This is made possible by the main originality of the model: the description of the network through elementary entities, the capillary elements, which are not required to connect to one another. Rather, the connectivity is recovered when the characteristics of the tissue, its hydraulic conductivity, is computed through the summation of the individual contributions of each capillary element. In this way, the capillary elements are free to appear in the adequate locations and directions in response to environmental cues, such as the local oxygen gradient or shear stress tensor. This flexible approach, which does not require to keep track of the network connectivity has first been proposed in the context of ant trail formation in [6].

There is a vast literature about angiogenesis and it is not feasible to compare our model with all previous approaches. However, most of them focus on the migration of EC as the main driving mechanism of capillary network formation (see e.g. [3,5,9,15,41,44,53,59,63,76,77,79,80,82] and the review [71]). More importantly, most of them let the tip cell play a special role. The mechanism is that an EC, in response to a chemical signal emitted by a tumour or a hypoxic region, starts to move upwards the gradient of this chemical. Doing so, it lays down a chemical trail like a slime trail left behind by a snail (hence the nickname 'snail-trail model' for this kind of models). This chemical trail is then followed by other EC that form the stem of the blood vessel. An early model of this type is [3], further elaborated in [9,63]. In the multidimensional cell-based versions of this model the rules for branching and creation of a new stem are often postulated with characteristics borrowed from the observations. Hence, the resulting patterns cannot be qualified as emergent in that they are directly resulting from the rules imposed to the agents. Importantly, our model does not give any special role to the tip capillary element. All capillary elements obey the same rules and this is enough to produce a blood vessel. To this extent our model is more parsimonious than the snail trail model. Likewise, there are no special instructions to tell the branches to bifurcate and how they should bifurcate. The resulting patterns emerge in a non-direct way from the individual rules. Finally while in most previously cited models, the driving force of network development is cell motion upwards a chemotactic gradient, in ours, the leading role is taken by the blood flow. Likely, both phenomena are important and the two types of mechanisms should be combined to accurately describe the biological phenomenology. Indeed, there is evidence of both EC motion during angiogenesis and sprouting of new blood vessels through WSS generated by blood flow.

Another viewpoint developed in the literature is that of the optimality of the network with respect to some cost functional. This approach dates back to [50]. Recently, inspired by this approach and the discrete network models of [35], a series of papers derive and analyse continuum models of network formation [1,30,31]. These are PDE models for a vector quantity akin to a magnetization in ferromagnetism and which encompasses both the density and orientation of the network. This quantity co-evolves with the hydraulic conductivity of the porous medium generating a positive feedback loop in a very similar fashion as in our model. There are however two important features that are needed to generate a network in these models: network diffusion on the one hand, and more importantly a nonlinear loss term interpreted as the metabolic cost of maintenance of the network on the other hand. In particular it is shown in [31] that there is no network generation unless a specific power law range is used for this nonlinear loss term. Our model does not involve any network diffusion (however, one could view the noise carried by the stochastic processes as a proxy for network diffusion). As reported at the beginning of Sect. 3.1.2, it is insensitive to the presence of a nonlinear capillary pruning process, itself being analog to the metabolic cost term of [31]. It is not surprising that continuum models do not generate instabilities in situations where microscopic models do (see another example in [6]). Noise involved in particle models is a powerful instability trigger which is absent from continuum models. This also brings some questions about the validity of the continuum models, although some rigorous results have been proved [29].

An analogy worth being made is with Diffusion-Limited Aggregation (DLA) models [2,32,33] which have been used to model electrical breakdown [62]. In the DLA process, a network is triggered by the large gradient of a potential function satisfying Laplace's equation and feedbacks on this potential by changing the boundary conditions. It could be seen as related to our model through a singular limit in which the capillary element conductivity would be infinite. Our model bears also analogies with canalization models in leaf venation processes (see [45,46] and the review [68]). These models describe how auxin transport (a hormone essential for plant development) shapes the vein network of leaves. Similarly to our model, they are based on a diffusion equation in which the diffusion coefficient evolves self-consistently with the flux of auxin. This generates a positive feedback triggering an instability. In another domain, geosciences, models of landscape evolution are based on a similar feedback between water flux and soil elevation [11,12] although in this case, soil elevation feedbacks into the advection part of the water height equation rather than in its diffusion part.

Obviously there are many directions by which the model could be improved to better account for biological observations. We have already commented on the need to include growth factors such as VEGF and other signaling chemical species as well as to allow for capillary element mobility and chemiotactic sensing ability. This would require the addition of another convection-diffusion-reaction equations describing the transport of these chemicals and its coupling with the capillary creation processes. Another aspect is active migration of EC from existing vessel walls to form new blood vessels. Currently capillary elements appear from nowhere and are immobile. In reality, EC are recruited among those lining the existing blood vessels and move in procession (the snail trail) to form new blood vessels. This active migration could be easily added to the model. In this motion EC many be impeached or guided by the structure of the ECM, a process named haptotaxis. In doing so, EC may also contribute to remodel the ECM (see [34,44,54] for models of haptotaxis). In [60], the self-organization of ECM and functional cells is modeled using a strategy similar to that presented here and the two models could be coupled to offer a comprehensive description of capillary formation in a developing tissue. The restriction to two-dimensions should be removed as biological reality is fully three-dimensional. Here, the real issue is computational efficiency and fast resolution methods should be sought. The time-scale problem already pointed out in Sect. 3 must also be addressed, by e.g. directly looking for a stationary solution of the oxygen transport equation rather than using its time-dependent version as in the current model. This requires the development of robust stationary solvers able to deal with regions of zero oxygen density. One can also dispute the validity of Darcy's law in the tissue (see e.g. [64,65]). One improvement would be to consider Brinkman's equation, i.e. adding the influence of blood viscosity. This would allow for a larger span of possible boundary conditions for the blood velocity and the possible emergence of more complex flow patterns. An extra step in increasing complexity would be distinguishing between blood in the capillaries and the interstitial fluid outside. This is needed if one wants to take specific account of oxygen transport across the capillary walls [58]. More generally, tissues have poroelastic character (i.e. there are interactions between tissue deformations and blood flow) which might be interesting to take into account.

Moving towards identifying which structures are biologically relevant, there is a need for network shape quantification. One possibility would be to skeletonize the capillary network and extract quantifiers such as statistics of branch lengths, branch widths and branching angles at junctions. Other shape quantifiers could involve the surface ratio between the network and its convex envelope for instance. Once a set of relevant quantifiers are selected, a parametric exploration of the model could be attempted. However this would first require significant algorithmic improvements as discussed below. Due to the large dimension of parameter space, sensitivity analysis would be needed. After these prerequisites, it would be possible to compare the model with biological images on which the same quantifiers would have been calculated. In doing this model calibration, the use of recent machine-learning techniques may prove necessary. Of course, these comparisons will be limited by the two-dimension character of the model, which makes the development of a three-dimensional model all the more necessary. Among algorithmic improvements needed for the model to be of practical use, the first one is the development of a faster elliptic solver. Due to the simple domain shape, a Fast Fourier Transform (FFT) solver could be an option, provided attention is paid to the aliasing problem in view of the large spatial inhomogeneity of the problem [4]. A possible remedy to that large inhomogenity could be multi-scale elliptic methods [18]. The acceptance/rejection for the capillary creation/deletion process could be improved by the use of Markov Chain Monte Carlo methods such as the Metropolis-Hastings algorithm. As pointed out above, the resolution of stationary oxygen transport would be necessary to recover time-scales compatible with observations. Such algorithmic improvements will not resolve the computational complexity bottleneck facing simulations of large tissues or organs. Here a paradigmatic shift is necessary and consists of resorting to a continuum model for the capillary network as well. A rigorous passage from the discrete model to the continuum one is however quite involved (see e.g. [29]) and faces conceptual problems, such as the greater difficulty to generate networks with a continuum model as discussed above. Without solving these difficulties, any phenomenological continuum model of capillary network would be subject to caution.

As summarized here this new network formation model opens many different exciting research avenues. It offers a new paradigm for capillary network creation by placing the flow of blood at the central place in the process. This paper provides a proof of concept of this approach and elaborates a road map by which the model can be gradually improved towards a fully fledged simulator of blood capillary network formation. Such simulator would have huge potential for biological or clinical applications in cancer, wound healing, tissue engineering and regeneration. Besides biological or clinical sciences applications the approach could also be adapted to plant biology (for leaf venation or root formation), physics (lightnings of thunder) or engineering (dielectric breakdown).

Appendices

A.   Finite Element Method for blood flow

A.1.   Weak formulation and approximation

Let us first consider the elliptic problem (3) posed on $\Omega_1$ with boundary conditions (17), (18), (19). Let $\mathcal{T}^h$ be a decomposition of the domain $\Omega$ into non-overlapping identical rectangle elements with mesh size $\Delta x$ and $\Delta y$ in the $x$ and $y$ directions and let $h = \max\{\Delta x, \Delta y\}$. Let $\mathcal{V}_h$ be the conforming finite element space associated with the partition $\mathcal{T}^h$, i.e. $\mathcal{V}^h : = \{ v \in \mathcal{C} ( \bar \Omega) : v|_{\mathcal R} \in \mathcal{Q}_1 \text{ for all } \mathcal{R} \in \mathcal{T}^h \}$, where $\mathcal{Q}_1$ denotes the space of polynomials of degree $1$ in each direction in ${\mathcal R}$ and ${\mathcal C} ( \bar \Omega)$ is the space of continuous functions defined on the closure $\bar \Omega$ of $\Omega$. Introduce the decomposition of the boundary $\partial \Omega = \Gamma_{N} \cup \Gamma_{D}$ with $\Gamma_{D} = \Gamma_{1, D} \cup \Gamma_{2, D} \cup \Gamma_{3, D} \cup \Gamma_{4, D}$ and $\Gamma_{N} = \Gamma_{1, N} \cup \Gamma_{2, N}$. The set of test functions is given by $\mathcal{V}^h_0 : = \{ u \in \mathcal{V}^h \, : \, u|_{\Gamma_D} = 0\}$. Multiplying (3) by a test function $v \in \mathcal{V}^h_0$, integrating over $\Omega$ and using Green's formula we obtain

Thanks to the assumptions on the test function and the boundary condition (19), the second term of the left-hand-side of (30) is zero. Let $\{ \varphi_1, \varphi_2, \ldots, \varphi_N \}$ be a basis of $\mathcal{V}^h_0$ and introduce $p_{\rm{ lift}}$ the function of $\mathcal{V}^h$ which interpolates the Dirichlet boundary values (17), (18) of $p$ at the nodes on $\Gamma_D$. We introduce $\tilde p = p - p_{\rm{ lift}}$ and denote by $\tilde p^h$ an approximation of $\tilde p$ in $\mathcal{V}^h_0$ satisfying

Introducing $A = ( A_{ij})_{i, j = 1, \ldots, N}$ the $N \times N$ symmetric positive-definite matrix whose entries are $A_{ij} = \int_\Omega ( \textbf{K} \nabla_x \varphi_i) \cdot \nabla_x \varphi_j \text{d}x$, equation (31) gives rise to the linear system

where $b = (b_i)_{i = 1, \ldots, N}$ is the vector with entries $b_i = - \int_\Omega ( \textbf{K} \nabla_x \varphi_i) \cdot \nabla_x p_{\rm{ lift}} \text{d}x$ and where $\textbf{p} = (p_i)_{i = 1, \ldots, N}$ is the vector of unknown values $p_i$ such that $\tilde p^h = \sum_{i = 1}^N p_i \, \varphi_i$. To solve (32) we use a preconditioned conjugate gradient method using the inverse of the diagonal of $A$ as a preconditioner (see [14]). In order to compute an approximation for each of the entries of the matrix $A$ we approximate $\textbf{K}(x)$ in each element $\mathcal{R}$ of the triangulation $\mathcal{T}^h$ by a constant function $\widetilde{\textbf{K}}$ equal to the average of the values of $\textbf{K}$ at the vertices of $\mathcal{R}$. With this approximation we compute the local stiffness matrix $A^{\mathcal{R}}$ in the element $\mathcal{R}$ given as (where the nodes are ordered in counterclockwise order starting from the lower left corner):

For Geometry 2, we proceed following a similar procedure as for Geometry 1 taking into account the periodic boundary conditions. For more details on the theory and implementation of the finite-element method we refer to [7,25].

A.2.   Velocity computation

We compute the gradient of the pressure using a second order finite difference method and for the sake of completeness we state the formulas below. For the partial derivatives with respect to $x$ and assuming that $( x_i, y_j)$ is a node not lying on the boundary of $\Omega$ (where $\Omega$ is $\Omega_1$ or $\Omega_2$) we use the centered difference scheme

For nodes on the left boundary of $\Omega$ we use a forward finite difference scheme

and a symmetric formula for nodes on the right boundary. The partial derivatives with respect to $y$ are computed in a similar fashion.

B.   Particle approximation of the oxygen flow

We use a splitting method to solve (4). In the first splitting step, we use the Smoothed Particle Hydrodynamics (SPH) method to solve for the advection and diffusion terms, i.e. the left-hand side of (4). In the second step, we solve for the reaction term i.e. the right-hand side of (4).

B.1.   SPH approximation of the convection-diffusion step

Here, we give the details for the first step, i.e. we assume $\beta = 0$ throughout this section. Supposing for a while that the vector field $\textbf{v}$ in (5) is given and smooth, then the density

where $\delta$ is the Dirac delta at $0$ and $m>0$ is any positive constant, is a solution of Eq. (4) if and only if $\textbf{Y}_\ell(t)$ satisfies the following ODE:

A measure satisfying (35), (36) is called a particle solution, $\textbf{Y}_\ell(t)$ is the position of the $\ell$-th particle at time $t$ and $m$ is the particle mass (here chosen identical for all particles). There is flexibility in choosing $m$ and the initial particle positions $\textbf{Y}_\ell(0)$ which is used to best approximate initial and boundary conditions. However, as such, the formula does not make sense. Indeed, $\rho$ being a sum of Dirac deltas, the right-hand side of (36) is not defined. The SPH methodology consists of introducing a mollifier kernel $W$: $x \in \Omega \mapsto W(x, \eta) \in [0, \infty)$ where $\eta>0$ satisfying $\int_{\Omega} W(x, \eta) \, dx = 1$ and $W(\cdot, \eta) \to \delta$ as $\eta \to 0$ in the distributional sense. Then $\rho$ can be approximated by

which is now a smooth function. Thus, it can be composed with a nonlinear function $U$: $r \in [0, \infty) \mapsto U(r) \in {\mathbb R}$ giving

and it can be differentiated:

This strategy is at the core of the SPH method [47]. Here, we use this strategy to replace (36) by

which is now completely well-defined. In the present context of diffusion equations, this method is called the diffusion velocity method and was first introduced in [17]. Several choices for the mollifier kernel $W$ are possible depending on the context. We have chosen the poly6 kernel for its simplicity and its reliability in the numerical solution of the Navier-Stokes equations (see [48]). This kernel is defined in ${\mathbb R}^2$ by

We solve (40) using the forward Euler scheme

where $\textbf{Y}_\ell^n$; $\textbf{V}^n_\ell$ are approximations of $\textbf{Y}_\ell (t^n)$, $\textbf{V}_\ell (t^n)$ and where $t^{n}$ is the $n$-th time discretization point. The time-step $\Delta t^n = t^{n+1} - t^n$ is chosen respecting the CFL condition

to preserve the stability of the scheme, with $C \leq 1/2$ (see Table 1). In addition, in order to have a good acceptance-rejection sampling for the creation of capillaries, we also impose the constraint $\Delta t^n \leq 0.01$. The blood velocity $\textbf{u}({\textbf{Y}}_\ell(t), t)$ is obtained through a linear interpolation from the values of $\textbf{u}$ given by the finite element calculation at the nodes of the finite element mesh. The diffusivity matrix $\textbf{D}$ is updated at the end of the capillary creation-removal process (described below) through (16). For better efficiency, it is first computed at the nodes of the finite-element mesh. Then, the diffusion coefficient $\textbf{D}({\textbf{Y}}_\ell(t), t)$ is recovered by linear interpolation exactly like the blood velocity $\textbf{u}$.

B.2.   Death process for the reaction step

Now, we consider the right-hand side of (4), i.e. we assume $\textbf{v} = 0$ in (4). We approximate this equation by a simple death process on the particle approximation (35). This step is described in Algorithm 2 below.

B.3.   Boundary conditions

B.3.1.   Periodic boundary conditions

These are boundary conditions (29) in Case 2 of Sect. 2.6.2. If the $y$-coordinate of a particle exits the range $[0, L_y)$, it is changed to $y + k L_y$ where $k$ is the unique integer in ${\mathbb Z}$ such that $y + k L_y \in [0, L_y)$. The $x$-coordinate is unchanged.

B.3.2.   Outflow boundary conditions

As described in Sect. 2.6.2, the boundary conditions (21) and (28) guarantee that the oxygen flow is outgoing across the corresponding boundaries. They are taken into account by removing the particles that exit the domain through these boundaries.

B.3.3.   Blood vessel boundary condition

These are conditions (20) and (26). We define a ghost domain $\Omega(L_{\hbox{ ghost}}) = [-L_{\hbox{ ghost}}, 0] \times [L_{\rm{ min}}, L_{\rm{ max}}]$, where $L_{\hbox{ ghost}}$ is changed dynamically as explained below. The ghost domain intersects $\Omega$ along the boundaries where the Dirichlet conditions (20) or (26) need to be enforced, i.e. $\Gamma_{1, D}$ (Case 1) or $\Gamma_{1, D}^{\mbox{ mid}}$ (Case 2), which here, for the ease of notation will be collectively denoted by $\Gamma_{1, D}$. Supposing that in $\Omega(L_{\hbox{ ghost}})$ there is a density $\rho_0$ of oxygen particles and that their velocity in the $x$-direction is $\bar{v}$, their flux $j$ across $\Gamma_{1, D}$ is then $j = \rho_0 \bar v$. A sensible estimate of $\bar{v}$ is given by the blood velocity $\textbf{u}$ along $\Gamma_{1, D}$. Thus, we take $\bar{v}$ equal to the average of $\textbf{u}$ along $\Gamma_{1, D}$ (in practice, we take the average of $\textbf{u}$ at the finite-element nodes along $\Gamma_{1, D}$). But this flux $j$ is also given by $j = N / ( L \Delta t)$ where $N$ is the number of oxygen particles (of mass $m = 1$) crossing $\Gamma_{1, D}$ during a small time interval $\Delta t$ and $L = L_{\rm{ max}} - L_{\rm{ min}}$ is the length of $\Gamma_{1, D}$. We can then estimate the number $N^n$ of particles to be created during time step $\Delta t^n$ as $N^n = \rho_0 \, L \, \Delta t^n \, \bar v^n$ (where $\bar v^n$ is the estimate of $\bar v$ at time $t^n$). We take $L_{\hbox{ ghost}}^n$ (the width of the ghost domain at time $t^n$) to be $L_{\hbox{ ghost}}^n = \Delta t^n \, \bar v^n$ (we recall that $\Delta t^{(n)}$ is determined by (43)) and put $N^n$ particles randomly uniformly in the corresponding ghost domain $\Omega(L_{\hbox{ ghost}}^n)$. During the time step $\Delta t^n$, they are moved with velocity $\bar v^n$ and end up inside the actual domain $\Omega$. In this way, we ensure that the flux of particles entering the domain corresponds to the density $\rho_0$. As $N^n$ is in general not an integer, we insert $\lfloor N^n \rfloor$ particles if $N^n \geq 1$ (where $\lfloor N \rfloor$ denotes the greatest integer smaller than or equal to $N$). If $N^n < 1$, we pick $r$ in $[ 0, 1]$ with uniform probability and insert one particle in $\Omega(L_{\hbox{ ghost}}^n)$ randomly uniformly if $r \leq N^n$ and create no particle otherwise. In the unlikely event where there are no particles in the computational domain $\Omega$, (43) cannot be use. Instead, we set $\Delta t = C h / \bar v$, where $C$ is the same CFL number as that used in (43).

B.3.4.   No oxygen boundary condition

This is condition 22 along $\Gamma_{1, N} \cup \Gamma_{2, N}$ in Case 1 of Sect. 2.6.2. We use the same methodology as in Sect. B.3.3 but with $\rho_0 = 0$, simply meaning that we inject no particle across this boundary.

C.   Acceptance-Rejection sampling method for the capillary creation

C.1.   Creation of new capillaries and modification of the tissue

We use an acceptance-rejection method in order to sample the spatio-temporal Poisson processes involved in the creation of capillary elements. Set time equal to $t^n$ and $\Delta t^n$ given by (43). We first choose the total number of sampling points to be $N_c$ (see Table 1 for the value used in our simulations and Sect. 2.7 for a discussion of this choice). We set $S = \frac{|\Omega|}{N_c}$ where $|\Omega|$ is the area of the full domain $\Omega$ and we define $\textbf{u}^\perp ( \textbf{X})$ to be the rotation by $90^\circ$ in counterclockwise direction of $\textbf{u} (\textbf{X})$. The algorithms for capillary element creation by oxygen gradient, by reinforcement or by WSS, (see Sect. 2.4) are respectively given at Algorithms 3, 4, 5.

 

 

C.2.   Removal of unused capillaries

The capillary pruning algorithm is a Poisson time process only, not a spatio-temporal Poisson process. So we do not need the area discretization parameter $S$. The capillary pruning algorithm is given in Algorithm 6.

D.   Description of the videos

Videos can be found on https://doi.org/10.6084/m9.figshare.c.4287575.v1

The videos have been obtained by running the model with the parameters set to the values given in Tables 1 and 2 and with the following choices for the geometry and the capillary element creation/pruning mechanisms:

$1)$ Movies (a)-(d): Geometry $\Omega_1$. All the capillary element creation/pruning mechanisms "ON", corresponding to Figs. 7-9 and 10(A).

$\underline{\rm Movie\; (a)}$ Positions of oxygen particles (red spots) and capillary elements (blue rods), corresponding to Fig. 7.

$\underline{\rm Movie\; (b)}$ Isolines and heatmap of the pressure $p$, corresponding to Fig. 8.

$\underline{\rm Movie\; (c)}$ Heatmap of the Frobenius norm of the hydraulic conductivity tensor, corresponding to Fig. 9.

$\underline{\rm Movie\; (d)}$ Same as movie (a) but with mesh size divided by $2$, corresponding to Fig. 11.

$2)$ Movies (e)-(j): Geometry $\Omega_1$. Positions of oxygen (red spots) and capillary elements (blue rods) with some capillary creation/pruning mechanisms turned off (corresponding to Figs. 10(B) to (G)).

$\underline{\rm Movie\; (e)}$ (corresponds to Fig. 10(B))

WSS "on", oxygen gradient "off", reinforcement "on" and pruning "on".

$\underline{\rm Movie\; (f)}$ (corresponds to Fig. 10(C))

WSS "off", oxygen gradient "on", reinforcement "on" and pruning "on".

$\underline{\rm Movie\; (g)}$ (corresponds to Fig. 10(D))

WSS "off", oxygen gradient "off", reinforcement "on" and pruning "on".

$\underline{\rm Movie\; (h)}$ (corresponds to Fig. 10(E))

WSS "on", oxygen gradient "on", reinforcement "off" and pruning "on".

$\underline{\rm Movie\; (i)}$ (corresponds to Fig. 10(F))

WSS "on", oxygen gradient "off", reinforcement "off" and pruning "on".

$\underline{\rm Movie\; (j)}$ (corresponds to Fig. 10(G))

WSS "off", oxygen gradient "on", reinforcement "off" and pruning"on".

$3)$ Movie (k): Geometry $\Omega_1$. Positions of oxygen (red spots) and capillary elements (blue rods) with all the capillary creation mechanisms "on" and pruning turned "off" (Corresponding figure not shown in the text).

$4)$ Movie (l): Geometry $\Omega_2$. Positions of oxygen (red spots) and capillary elements (blue rods) with all the capillary creation/pruning mechanisms turned "on" (corresponding to Fig. 14).

Acknowledgments

This project was supported by the Région Midi-Pyrénées, France, under contract ref. APRTCN 14050455. PD acknowledges support by the Engineering and Physical Sciences Research Council (EPSRC) under grants no. EP/M006883/1 and EP/N014529/1, by the Royal Society and the Wolfson Foundation through a Royal Society Wolfson Research Merit Award no. WM130048 and by the National Science Foundation (NSF) under grant no. RNMS11-07444 (KI-Net). PD is on leave from CNRS, Institut de Mathématiques de Toulouse, France. BA gratefully acknowledges the hospitality of the Department of Mathematics, Imperial College London, where part of this research was conducted. PAS was supported by EPSRC under grant no. EP/M006883/1.

Data statement: no new data were collected in the course of this research.

Link to videos: https://doi.org/10.6084/m9.figshare.c.4287575.v1