Percolation, Morphogenesis, and Burgers Dynamics in Blood Vessels Formation

Experiments of in vitro formation of blood vessels show that cells randomly spread on a gel matrix autonomously organize to form a connected vascular network. We propose a simple model which reproduces many features of the biological system. We show that both the model and the real system exhibit a fractal behavior at small scales, due to the process of migration and dynamical aggregation, followed at large scale by a random percolation behavior due to the coalescence of aggregates. The results are in good agreement with the analysis performed on the experimental data.

The problem of morphogenesis and self-organization in biological systems is both an actual and a long standing one [1,2]. In living beings, complex structures arise from the aggregation of separate constituents. Such structures, although irregular, often show some sort of hidden order, like self-similarity and scaling laws [3,4]. Scaling laws are known to emerge not only in the physics of phase transition [5] and in percolation [6], but also from several kinds of aggregation dynamics, and are often considered as fingerprints of the process which led to the formation of the structure itself [7,8]. A typical example of natural structure characterized by non trivial scaling laws is the vascular network [9]. In recent years many experimental investigations have been performed on the mechanism of blood vessel formation [10]. Cells are cultured on a gel matrix and their migration and aggregation observed through videomicroscopy. Tracking of individual trajectories shows marked persistence in the direction, with a small random component superimposed [11]. The motion is directed towards zone of higher concentration of cells, suggesting that chemotactic factors play a role. Cells migrate over distances which are an order of magnitude larger than their radius and aggregate when they get in touch with one of their neighbors. In a time of the order of 10 hours they form a continuous multicellular network which can be described as a collection of nodes connected by chords. The mean chord length is approximately independent on the initial cell densityn, with an average valuel = 200 ± 20 µm for n between 100 to 200 cells/mm 2 . Observing the formed structures, one finds that, by varying the initial cell density, there is a percolative transition. Below a critical value n c ∼ 100 cells/mm 2 groups of disconnected structures are formed (Fig. 1a). Forn ≥ n c a single connected network is instead visible (Fig. 1b,c). For higher values ofn one observes a sort of "swiss cheese" pattern ( Fig.  1d). In the biological system the percolating property is of physiological relevance, since it is directly linked with the functionality of blood vessels.
In this paper we propose a theoretical model which turns to be in good agreement with experimental observations. It allows to well reproduce both the observed percolative transition and the typical scale of observed vascular networks. With respect to standard percolative models, here an important role is played by migration and dynamical aggregation of particles. The model appears thus as a possibly new and physically interesting representative in the class of percolation models. We characterize the phase transition from the point of view of scaling laws, both in the model and in the real system, computing critical indices and the fractal dimension of the percolating cluster through extensive numerical simulations and analysis of the experimental data. The model describes the cell population as a continuous density field n and velocity v; it also assumes the presence of a concentration field c of soluble factor. Cells are modeled as a fluid accelerated by gradients of the soluble factor. The latter is supposed to be released by cells, diffuse, and degrade in finite time, in agreement with experimental observations. These assumptions give rise to the following equations: where D, α, τ , are respectively the diffusion coefficient, the rate of release and the characteristic degradation time of soluble mediators, and µ measures the strength of cell response. Initial conditions are given in the form of a set of randomly distributed bell-shaped bumps in the density field, representing single cells initially spread in the system, with zero velocities and zero concentration of the soluble factor. In Eqs. (1) a multidimensional Riemann (1b) equation is coupled to the diffusion equation (1c), while (1a) is just a continuity equation expressing conservation of the total cell number. Multivaluedness of solutions to (1b) can be prevented by adding a viscous term ∇ 2 v in the r.h.s. One than gets the analog of Burgers' equation, which has been used to describe growth of molecular interfaces [12] and the emergence of network-like patterns [13] in the large-scale distribution of masses of the Universe [14,15]. Instead of viscosity, we introduce a phenomenological, density dependent pressure term −∇φ(n), in the r.h.s. of (1b), with φ(n) zero for low densities, and rapidly increasing for densities above a threshold σ −2 , where σ = 30 µm is the mean cell radius, in order to describe the fact that cells do not interpenetrate.
The model provides a natural lenght scale r 0 = √ Dτ , which is the effective range of the interaction mediated by the soluble factor. This suggests that, starting from the above mentioned initial conditions, Eqs. (1) should develop network patterns characterized by the same scale. Direct measurements of D and τ give D ∼ 10 −7 cm 2 /s, τ ∼ 1h, from which one gets r 0 ≃ 200 µm, a value which is in good agreement with the mean chord lengthl experimentally observed.
We performed numerical simulations of model (1) on square boxes of size L = 1, 2, 4, 8 mm, with periodic boundary conditions, using a finite volume method [16]. We used the experimental values for D and τ , while the unknown parameters µ, α were set to one (which amounts to an appropriate rescaling of time and the concentration field c). The initial conditions were given by throwing the same number of cells as in the biological experiments in random positions inside the box, with zero velocities and zero concentration of the soluble factor, with a single cell given initially by a Gaussian bump of width σ and unitary weight in the integrated cell density field n. Starting from these condition, Eqs. (1) were numerically integrated. The simulation was stopped when the vascular network was formed, or a stationary state was reached. We observed a remarkable correspondence between simulations and experiments both in the the dynamic evolution details (not shown) and in the reproduction of the percolating patterns (Fig. 2).
In order to study the percolative transition we computed between 100 to 200 realizations of the system for each box size and cell density, with different initial positions of the cells. The simulation gives as a result a continuous density field n. In order to study percolation properties we partitioned the box in little square boxes of size 2 −6 mm, forming the sites of a square lattice, and set each site as filled if the mean cell density inside the corresponding little box was larger than the threshold σ −2 , empty otherwise. We then identified clusters of nearest neighbor filled sites. In Fig. 2  represent the underlying substrate. As in the case of experimental data, a percolative transition is observed at some critical cell density n c . Furthermore the model reproduces quite well the typical structure of the vascular network, with chords of lengthl ≃ 200 µm for a wide range of cell densities, in good agreement with the expected value r 0 . For each set of realizations, corresponding to a value of the box size L and of the mean cell densityn, we identified clusters of nearest neighbor filled sites, and measured the following quantities: 1) the percolation probability Π, defined as the fraction of realizations in which a percolating cluster appears; 2) the infinite cluster density P , defined as the density of sites belonging to the infinite (greatest) cluster; 3) the mean cluster size S = (1/N ) s n s s 2 , where N is the total number of sites and n s is the number of clusters of size s, excluding from the sum the infinite cluster. In percolation models, in presence of a second order transition at some mean cell density n c , these quantities show a characteristic singular behavior near the percolation point. Namely, the percolation probability is zero below the transition and one above it; the infinite cluster density, which plays the role of order parameter, vanishes as |n − n c | β from above; and the mean cluster size has a singularity |n − n c | −γ on both sides of the transition. For finite size systems, the same quantities obey near the transition to the following relations: Using Eq. (2a), one can find the critical density n c as the point where the curves Π(n, L) for different sizes L cross. The estimated value is n c = 94 ± 1 cells/mm 2 . Then, the critical index ν is found by plotting Π(n, L) as a function of (n − n c )L 1/ν , where n c is held fixed to the value previously computed, and looking for the value of ν that gives the best data collapse. The error on the computed value of the index is estimated looking for the range of values that gives an acceptable data collapse.
Holding the value of n c and of ν fixed, the other two relations Eq. (2b) and (2c) are used to estimate the value of the indices β and γ, by looking for the data collapse of L β/ν P (n, L) and L −γ/ν S(n, L) as a function of (n − n c )L 1/ν . We observed that while n c is sensitive to the choice of the threshold in a neighborhood of the natural value σ −2 , the critical exponents do not depend on it. In Fig. 3 the data collapses for Π and S are shown. The estimated values of the critical indices are reported in the second column of Tab. I. Within the errors, the values obtained are the same that one expects for the simple random percolation model, that is when the sites of the lattice are uncorrelated, and each one occupied with the same probability. As critical indices are linked to the large scale structure of the percolating cluster, this means that, on such a large scale, the structure of the vascular network is mainly determined by the initial random positioning of the cells,  and is not altered by the dynamical process of migration and aggregation. On the other hand, on smaller distances clusters are quite different from the ones of random percolation. For example the mean fraction of occupied sites, which is a quantity sensible to the small scale structure of clusters, is about 0.2 at the percolation point, much lower than the corresponding figure for random percolation, that is 0.59. A quantity that can give us informations about the structure of the percolating cluster at different scales is the density ρ(r) as a function of the radius. This is defined as the mean density of sites belonging to the percolating cluster, inclosed in a circle of radius r centered at one site belonging to the cluster, and averaged over different centering sites and different realizations of the system. For a fractal object with fractal dimension D, this should scale as ρ(r) ∼ r D−d . For the percolating cluster of random percolation at the critical point, one expects a fractal dimension D = 1.896. We computed ρ(r) for our model at the percolation point n = n c (see Fig. 4), and found two different regimes depending on the scale. For r > r c , with r c = 0.77 ± 0.08 mm, a behavior compatible with random percolation is found, with fractal dimension D = 1.87 ± 0.03. On the other hand, for r < r c a different exponent is found, D = 1.50 ± 0.02. This exponent may be the signature of the dynamic growth process driven for r < r c by the rapidly oscillating components of the concentration field. We have repeated the measures described above on a set of experimental data, consisting of 28 digital photographs (with a resolution of 1024 2 pixels) of mature structures obtained on a matrix of size L = 2 mm, starting from initial densities of 50, 75, 100, . . ., 200 cells/mm 2 . The photographs are the results of two experimental sessions, each performed in duplicate, so that 4 photographs are available for each density value. The data show the presence of a percolative transition at n c ∼ 125 cells/mm 2 , and are compatible with the values γ/ν = 1.83 obtained from numerical simulations. The available data are not enough to obtain ν and β/ν with a reasonable precision. Estimated values of the critical indices are reported in the third column of Tab. I. Finally, we measured the density of the percolating cluster as a function of the radius at the percolation point, and found a behavior surprisingly similar to the one found in the numerical model (see Fig. 4): for r > r c the fractal dimension of the percolating cluster is D = 1.85 ± 0.10 which is compatible with D = 1.87±0.03 from the numerical simulations and with the value D = 1.896 of random percolation. For r < r c instead D = 1.48 ± 0.05.
In conclusion, we have introduced a model which describes the formation of the vascular network with a typical scale lenght r 0 , and is in good agreement with experimental data. The percolative phase transition appears to be second order, at least to the accuracy allowed by the limitations in size which are intrinsic to our numerical experiments. Comparison of critical indices suggests that the phase transition falls in the universality class of random percolation, even in presence of migration and dynamical aggregation. This is confirmed by the fractal dimension of the percolating cluster, on scales larger than r c ≃ 0.8 mm, that is in agreement with the one expected for random percolation. On the other hand we found that, if observed on scales smaller than r c , the vascular network shows a different scaling behavior, with a fractal dimension D ≃ 1.5, that may be the signature of the dynamical process that led to the formation of the network. Thus, the model appears to be quite successful in describing in vitro experiments, where all the parameters are under control and one can easily tune the cell density. It would also be interesting to understand the relevance of this picture to vascular networks in living beings. It is known [10] that during embriogenesis blood vessels form through local (in situ) differentiation of a fraction of endothelial cells, that afterwards assemble in a vascular labyrinth, a process similar to the one described in this paper. However, in this case it is more difficult to have complete control over all the relevant parameters.
In particular, one would expect that the cell density is fixed through some feedback mechanism, in order to maximize efficiency [9]. Indeed, some experiments [17] on in vivo capillaries show compact space-filling structures, with fractal dimension D = 2, which would indicate that cell density is tuned by Nature to a value larger than the percolation threshold. More experiments, focusing on the relevance of the percolation mechanism, and the nature of cell density self-regulation, are needed to shed light on these issues.