## Abstract

In vertebrates, networks of capillary vessels supply tissues with nutrients. Capillary patterns are closely mimicked by endothelial cells cultured on basement membrane proteins that allow single randomly dispersed cells to self‐organize into vascular networks. Here we provide a model including chemoattraction as the fundamental mechanism for cell‐to‐cell communication in order to identify key parameters in the complexity of the formation of vascular patterns. By flanking biological experiments, theoretical insights and numerical simulations, we provide strong evidence that endothelial cell number and the range of activity of a chemoattractant factor regulate vascular network formation and size. We propose a mechanism linking the scale of formed endothelial structures to the range of cell‐to‐cell interaction mediated by the release of chemoattractants.

## Introduction

To distribute nutrients throughout the body vertebrates have evolved a hierarchical branching blood vascular system that terminates in a network of size‐invariant units, namely capillaries. Development of capillary networks characterized by typical intercapillary distances ranging from 50–300 μm is instrumental for optimal metabolic exchange (Chilibeck *et al*., 1997; Guyton and Hall, 2000; Goldman and Popel, 2001). The ability to form networking capillary tubes is a cell autonomous property of the endothelial cells (ECs), which need permissive but not instructive signals from the extracellular environment (Folkman and Haudenschild, 1980). In particular, it is well known that culturing ECs on Matrigel, a natural basal membrane matrix, markedly accelerates their morphological differentiation in geometric tubular networks, which are almost identical to vascular networks observed *in vivo* (Kubota *et al*., 1988; Grant *et al*., 1989; Risau and Flamme, 1995). In the embryo, the morphogenetic process by which free EC precursors (angioblasts) of mesodermal origin self‐assemble into a primitive vascular plexus is known as vasculogenesis (Risau and Flamme, 1995). Angioblasts either coalesce at the location where they emerge from the mesoderm (type I vasculogenesis) or they migrate through tissues and form blood vessels at distant sites (type II vasculogenesis) (Coffin and Poole, 1991). Angiogenesis is the cellular mechanism by which the primitive embryonic vasculature is remodeled into a mature vascular bed comprising arteries, capillary networks and veins (Carmeliet, 2000). Angiogenesis also includes penetration by sprouting and branching of vessels into avascular regions of the embryo, such as the central nervous system (Risau and Flamme, 1995). In addition, sprouting and branching blood vessels have been shown to receive significant contribution from free invasive angioblasts that move throughout the embryonic mesenchyme (Noden, 1989). Attempting to understand the mechanisms operating during formation of the vascular system presents a tremendous challenge. In particular, the issue of how ECs self‐organize geometrically into vascular networks is still to be elucidated.

In the study of self‐organized aggregation and pattern formation, a well‐established paradigm is a non‐linear diffusion equation known as the adhesion model, which has been utilized to describe various natural phenomena. The main feature of this equation is that it describes the dynamics of a particle fluid where initial slight inhomogeneities are strongly amplified leading to the emergence of patterns (Whitham, 1974; Kuramoto, 1984; Kardar *et al*., 1986; Zeldovich *et al*., 1990; Gurbatov *et al*., 1991; Vergassola *et al*., 1994).

In this paper we study an adhesion‐type model including chemoattraction as a fundamental mechanism for cell‐ to‐cell communication. On the basis of experimental evidence and theoretical insights, we show that non‐linear mechanics and chemotactic cellular dynamics fit into a model able to reproduce with great accuracy the formation of capillary networks *in vitro*.

## Results

### Directional EC migration in early stages of capillary network formation

Vascular network formation proceeds along three main stages: (i) migration and early network formation; (ii) network remodeling; (iii) differentiation in tubular structures (Folkman and Haudenschild, 1980; Figure 1). We performed accurate experimental and statistical analysis of *in vitro* pattern formation by plating ECs on a Matrigel film. The process was recorded by time‐lapse videomicroscopy (Figure 1A–D and Supplementary figure 1, available at *The EMBO Journal* Online) and individual cell trajectories were tracked by computer‐assisted analysis.

In the first stage (0–3–6 h; Figure 1A–C), after being randomly seeded on Matrigel (thickness: 44 ± 8 μm, *n* = 30), human ECs (from large veins or adrenal cortex capillaries) start moving in different directions, interact and adhere with their neighbors, spread and eventually form a continuous multicellular network (Figure 1C), whose geometry is not substantially modified in the following stages. Tracking of individual trajectories (Figure 2B and C) shows a highly directional movement toward zones of higher cell concentration. In the second stage (6–9 h) the network just undergoes a slight deformation (Figure 1C and D). In the third stage (9–12 h), individual cells fold up to form capillary‐like tubes along the previously formed bidimensional network (not shown) (Kubota *et al*., 1988). Because the main geometrical features of network patterns are determined during the first stage of the process, we focused our experimental and theoretical investigations on the early migration dynamics.

First we performed a statistical analysis of EC trajectories. In Figure 2A an early snapshot of the cell deposition process is shown (see Supplementary figure 1). Trajectories longer than one cell diameter have been tracked and are shown in Figure 2B and C. Migration appears to be directed toward zones of higher cell concentration. As chemotactic movement is a key mechanism in vasculogenesis (Roman and Weinstein, 2000), we hypothesize that early EC aggregation is supported by chemotactic intercellular cross talk mediated by chemoattractants. To test this hypothesis, we computed chemoattractant concentration fields starting from the assumption that each EC secretes a chemoattractant diffusing in the medium and degrading in time (see below). Figure 2B shows a color map of the simulated concentration field produced by virtual sources placed in the same positions as ECs in Figure 2A. Cell velocity vectors have been measured over 1‐min intervals and two types of angle (φ and φ) have been computed (Figure 2D): φ is the angle between two velocities relative to the same trajectory, measured at subsequent intervals of time; φ is the angle between a velocity in a point of a given trajectory and the simulated gradient of the concentration field of chemoattractant in the same point. Histograms of φ and cos φ, which are indices of mere directional persistence of cell movement, and φ and cos φ, which measure the tendency of cells to move towards zones of higher simulated chemoattractant concentrations, are shown in Figure 3. In standard culture conditions, trajectories of ECs plated on Matrigel have a high degree of directional persistence and are directed with high probability towards zones of higher simulated chemoattractant concentration (Figure 3 and Table I).

Next, we directly interfered with the hypothesized chemoattractant mechanism and identified as candidate vascular endothelial growth factor (VEGF)‐A, which is known to induce EC growth, survival, motility and formation of capillary beds (Neufeld *et al*., 1999), also by participating in EC autocrine loops (Namiki *et al*., 1995; Seghezzi *et al*., 1998; Helmlinger *et al*., 2000). ECs both from large vessels and microvasculature plated for 3 h on Matrigel release detectable amounts of VEGF‐A (control medium: 0 ng/ml, *n* = 4; ECs from large umbilical vein: 0.24 ± 0.08 ng/ml, *n* = 7; microvascular ECs: 0.67 ± 0.05 ng/ml, *n* = 4). To inhibit the effect of this autocrine loop on vascular pattern assembly, we tested an anti‐VEGF‐A neutralizing antibody, which inhibited capillary network formation by triggering EC apoptosis (percentage of Annexin V positive cells: control, 0.61 ± 0.20; antibody anti‐VEGF‐A, 32.11 ± 2.20; irrelevant IgG, 0.75 ± 0.41; *n* = 6). To overcome the apoptotic effects elicited by direct inhibition of VEGF‐A, we extinguished VEGF‐A gradients spreading from individual ECs plated on Matrigel by adding a saturating amount of exogenous VEGF‐A_{165}. Saturation of chemoattractant gradient with VEGF‐A_{165}, but not with a heat‐inactivated molecule, resulted in a strong inhibition of network formation (Figure 4A and B). Statistical analysis of individual EC trajectories obtained by computer assisted time‐lapse videomicroscopy showed that although trajectories retained a certain degree of directional persistence, saturating VEGF‐A gradients caused ECs to be indifferent to the direction of simulated concentration gradients of chemoattractant (Figure 3 and Table I). This observation was also confirmed by checkerboard analysis of EC migration (Table II), showing that VEGF‐A_{165} triggered directed, but not random EC migration, as also reported for EC precursors (Asahara *et al*., 1999).

*A percolative transition in capillary network formation*

To quantitatively characterize network structures, we then measured their metric and topological properties. Since fine‐tuning of an appropriate EC density is mandatory for development of vascular networks (Fong *et al*., 1999) we performed experiments by varying the initial cell number. We measured mean chord length, obtaining the approximately constant value ℓ̅ ∼ 200 ± 20 μm over a cell density *n̅* extending from 100–200 cells/mm^{2} (Figure 5B, C and I). By studying the network connectivity properties (Figure 5 and see Supplementary figure 2) we observed the formation of groups of disconnected structures (Figure 5A) below a critical value *n*_{c} ∼ 100 cells/mm^{2} and the formation of a single connected network for *n̅* ≥ *n*_{c} (Figure 5B and C). Increasing *n̅* over 200 cells/mm^{2}, the mean chord thickness grows to accommodate an increasing number of cells. Eventually, one observes a continuous carpet of cells interspersed with holes (‘Swiss‐cheese’ pattern; Figure 5D).

An abrupt change in the connectivity properties of a randomly formed structure through the variation of a control parameter is known in physics as percolative transition (Stauffer and Aharony, 1994). Percolative transitions have been thoroughly investigated and are classified in a small number of ‘universality classes’ (Stauffer and Aharony, 1994), corresponding to specific values of critical indices, such as the fractal dimension of the network. Such values are fingerprints of the growth process that leads to vascular structure formation and may be useful in discriminating normal and pathological vascular structures (Gazit *et al*., 1995).

*Mathematical modeling*

To encode the results of the above experiments into a mathematical model, we started from the following assumptions: (i) the cell population can be described by a continuous distribution of density *n* and velocity depending on the point on the Matrigel surface; we shall also assume the presence of a concentration *c* of chemoattractant; (ii) the cell population in the early stages of its evolution (i.e. for low densities) can be modeled as a fluid of non‐directly interacting particles; (iii) cells are accelerated by gradients of chemoattractants; and (iv) chemoattractants are released by cells, diffuse and degrade in finite time.

From these assumptions one obtains the following equations:

where *D*, α and τ are, respectively, the diffusion coefficient, the rate of release and the characteristic degradation time of chemoattractant, and β measures the strength of cell response. Equation 1 describes mass conservation, Equation 2 is a momentum balance with a chemotactic force and Equation 3 is a diffusion equation for the chemoattractant produced by ECs and degrading in time; in particular, the convective term on the right hand side of Equation 2 allows describing cell migration. For simplification, details of the model are not reported here (Gamba *et al*., 2003).

Equations 1–3 can be refined by taking into account the fact that tightly packed cells resist compression. This may be implemented by adding a density dependent pressure term −▿*p*(*n*), being zero for low densities and rapidly increasing for *n* ∼ *a*^{−2}, *a* being the average cell diameter. This term avoids cell compenetration and allows describing the crossover to the ‘Swiss‐cheese’ regime of Figure 5D.

Model (1–3) contains a natural length scale that represents the effective range of interaction mediated by chemoattractant and a natural time scale .

Initially, cells are randomly distributed on Matrigel. We therefore impose initial conditions in the form of a set of randomly distributed bell‐shaped bumps having a finite width of the order of the average cell diameter *a*, with null velocities. Non‐zero velocities are built up by the chemoattractive term starting from inhomogeneities in the density distribution. Then, the dynamics amplifies the density variations and forms a capillary‐like network. In the fast diffusion approximation, density inhomogeneities are instantaneously translated in a landscape of concentration of chemoattractants where details of scales ≤ξ are averaged out. Particles move toward landscape crests, which are separated by valleys of width ∼ξ. The dynamics sharpens the crests and empties the valleys, eventually producing a network structure characterized by a length scale of order ξ. This way, the model provides a direct link between structure size and intercellular interaction range.

*Prediction of the characteristic size of network structures*

From available data of molecular radii (Walter *et al*., 1992; Muller *et al*., 1997) and experimental results (Pluen *et al*., 1999) we infer that the order of magnitude of the diffusion coefficient for major angiogenic growth factors is *D* ∼ 10^{−7} cm^{2}/s. In our experimental conditions the half‐life of VEGF‐A is 64 ± 7 min. This gives ξ ∼ 100 μm, which is the same order of magnitude as the observed dimensions of capillary‐like networks on Matrigel.

*Computer simulations*

We performed numerical simulations of Model (1–3), using the experimental values *D* = 10^{−7} cm^{2}/s, τ = 64 min, *a* = 30 μm and the same initial cell number as in the biological experiments. We observed a remarkable correspondence between simulations and experiments both in the dynamic evolution details (Figure 1A–C and E–G) and in the final pattern structure (Figure 5A–D and E–H). The experimentally observed percolative transition is correctly reproduced: by varying the initial cell number we switched from a phase with several disconnected structures to a phase with a single connected structure (Figure 5E–H). The introduction of the pressure term described above allows reproducing the crossover to the experimentally observed ‘Swiss‐cheese’ regime (Figure 5D and H).

We verified that parameter ξ actually controls the pattern size (Figure 6A–C). To this aim, we measured the radial part P(*r*) of density correlations on mature structures, averaging over a large set of random initial conditions. A typical graph of P(*r*) (Figure 6D) shows a marked minimum, whose position ℓ_{o} measures the characteristic scale of network structures. We plotted ℓ_{o} versus ξ, obtaining an approximate linear relation with slope of order 1 (Figure 6E), confirming previous qualitative statements.

As a further validating prediction, we verified the capability of the model to reproduce the cell organization patterns obtained experimentally (Figure 4A and B) in the presence of saturating amounts of VEGF‐A_{165}. To this aim we simulated Model (1–3) switching off the chemoattractant mechanism, i.e. imposing β = 0, and assigning random velocities with a random distribution centered around 0, with SD ≈ 5 μm/min as measured in experiments, on ECs plated on Matrigel alone and saturated by VEGF‐A_{165}. Indeed, simulations produced a clear‐cut inhibition of pattern formation (Figure 4C and D) similar to that obtained by extinguishing VEGF‐A gradients experimentally (Figure 4A and B).

## Discussion

Formation of a vascular plexus by motile‐free EC precursors is a hallmark of embryonic vasculogenesis (Risau and Flamme, 1995). In this study, we provide evidence that the motility of a strictly constrained number of ECs under the control of a chemoattractant with a defined limit of activity is the pivotal parameter for vascular network formation by free ECs. A strategy based on tightly coupling experiments, theory and computer simulation allowed us to demonstrate that: (i) ECs endowed with a specific gene profile act according to simple mathematical laws, which ultimately govern vascular network assembly and geometry; (ii) a non‐linear equation, encoding these laws and essentially depending only upon two physically measurable parameters and the initial EC number, allows realistic numerical simulations of the process of vascular formation.

The first parameter of our model is the number of ECs. Its key relevance has been inferred by the demonstration of the existence of a percolative transition depending on cell density. Actually vascular networks are obtained precisely in the proximity of the percolative transition, corresponding to a critical number of ∼100 cells/mm^{2}. Cell number over or below the range 100–200 cells/mm^{2} is not permissive for network formation. Results of numerical simulations with >200 cells/mm^{2} are very similar to observations carried out in *flt‐1 ^{−/−}* mice (Fong

*et al*., 1999), where severe vascular system disorganization is a direct consequence of an overcrowded population of ECs. Thus, our model could prove useful in predicting and characterizing pathological situations, where the altered expression of angiogenic inducers and their counter‐receptors is thought to be responsible for abnormal blood vessel formation.

The second crucial element is directed cell migration governed by the presence of a chemoattractant, which mainly acts in a paracrine manner. This is proved by: (i) the production of VEGF‐A by ECs during the first hour of formation of capillary‐like structures; (ii) biological experiments showing that abrogation of a VEGF‐A‐mediated chemotactic field is detrimental to capillary‐like formation; (iii) time‐lapse videomicroscopy and statistical analysis of individual cell trajectories during the early phase of vascular pattern formation, which show how ECs are directed with high probability towards zones of higher simulated chemoattractant concentration in control, but not in VEGF‐A_{165}‐saturating conditions. This is consistent with the role of EC movement regulated by VEGF‐A gradients in the establishment of vascular patterns observed in several animal species (Cleaver and Krieg, 1998; Poole *et al*., 2001; Ylikorkala *et al*., 2001). Moreover, activation of the angiomotin–angiostatin system interferes with murine EC motility, inhibiting the formation on Matrigel of early EC aggregates and the subsequent formation of capillary‐like structures (Troyanovsky *et al*., 2001). Interestingly, the control of EC self‐assembly by chemotaxis is not an exception in nature. In fact, in many living organisms chemical field‐based movement is a key mechanism in the generation of self‐organized complex patterns and shapes, such as during bacterial colony formation (Ben‐Jacob *et al*., 1998), *Dictyostelium* morphogenesis (Dormann *et al*., 2000) and collective phenomena observable in social insects (O'Toole *et al*., 1999; Beekman *et al*., 2001).

Our model predicts that the diffusivity and the half‐life time of the chemoattractant determine the size of the vascular structure. Intriguingly, this seems to be in agreement with the observation that mice lacking heparin‐binding isoforms of VEGF‐A, characterized by smaller diffusivity, form vascular networks with a larger mesh (Ruhrberg *et al*., 2002). Although the vascular patterns observed *in vivo* by Ruhrberg *et al*. (2002) are mainly thought to be generated by angiogenesis rather than by vasculogenesis, one cannot exclude the possibility of a simultaneous occurrence of both phenomena, as previously reported (Cleaver and Krieg, 1999). Moreover, it has been shown that invasive angioblasts are also capable of freely moving throughout the embryonic mesenchyme and significantly contribute to sprouting and branching of angiogenic blood vessels (Noden, 1989). An additional issue is represented by the fact that *in vivo*, both in physiological and pathological conditions, several non‐endothelial cell types produce long‐range paracrine chemoattractant gradients that participate in the vascular networks' remodeling (Bussolino *et al*., 1997). To more properly describe this complex *in vivo* scenario, further technical refinement of our model and new experimental data are needed.

A theoretical model for *in vitro* angiogenesis has been previously proposed (Manoussaki *et al*., 1996), based on experiments performed with bovine aortic ECs (BAECs), focusing on cellular network reorganization driven by cellular traction (Vernon *et al*., 1992). In these experiments, where BAECs form *in vitro* angiogenesis after >20 h, a free‐migrating regime is not observed and late cellular motility is only along strain‐produced tracks of aligned matrix. Our observations and Vernon's clearly correspond to different regimes and to different, possibly complementary, biological features, including those of arterial, capillary and venous ECs (Garlanda and Dejana, 1997). Our observations focus on the early stages of vascular network differentiation and to exclude the role of mechanical forces we used a thinner Matrigel surface (44 ± 8 μm, see Materials and methods) which does not allow BAEC forming *in vitro* vascular networks (Vernon *et al*., 1992). In our experimental conditions, both human vein and capillary ECs rapidly formed a defined vascular network, whose geometry was largely established in the early, free migration stage of the process. The viscoelastic model of Manoussaki *et al*. (1996) is not adequate to describe the early stages of human vein or capillary EC network formation, where directed migration is the predominant issue. Moreover, the model presents several shortcomings: (i) simulations are performed starting from a nearly uniform initial cell density, with a very small perturbation around the uniform steady state, physically corresponding to a continuous initial cell monolayer. These initial conditions are very far from the more realistic ones we considered in our simulations; (ii) since cell traction is assumed to be the main motor of network formation, we suspect that with more realistic initial conditions the model would not lead to network formation at all; (iii) the main theoretical results are strongly affected by the specific form of the pressure dependence on EC density, which is not determined experimentally. Manoussaki's model, or some modification of it, is probably complementary to ours and useful to describe the stage of network remodelling where mechanical forces are clearly essential.

We conclude that vascular network assembly by ECs follows a mathematical law, which ultimately governs its size and geometry. We discovered the existence of a percolative transition depending on the cell density. Anatomical and physiological evidence shows that the small capillary wall is constituted by a single folded EC; thus vascular networks having minimal chord thickness are optimal (Folkman and Haudenschild, 1980; Fawcett and Raviola, 1994; Guyton and Hall, 2000). Such networks are obtained precisely in the proximity of the percolative transition. This is a really intriguing fact, since it connects a fundamental characteristic of a vital process with a very peculiar physical condition. One could wonder whether this is an isolated fact or if other similar examples can be found.

## Materials and methods

*In vitro experiments*

Chemotaxis, chemokinesis and apoptosis were performed on human ECs from cord veins and adrenal cortex (BD Biosciences Clontech, Palo Alto, CA) as described (Bussolino *et al*., 1992; Strasly *et al*., 2001). To analyze *in vitro* angiogenesis, glass coverslips (2 cm^{2} growth area) were coated with 0.12 ml of growth factor‐reduced Matrigel (8.8 mg/ml; BD Biosciences), that was immediately aspirated and the remaining film was allowed to solidify. The thickness was measured by confocal laser‐scan microscopy (Radiance 2100; Bio‐Rad Labs, Hercules, CA) on a Matrigel film made fluorescent by adding fluorescinated‐dextran (200 μg/ml, mol. wt 464 kDa, Sigma‐Aldrich, Milan, Italy) and was of 44 ± 8 μm (*n* = 30). ECs at the appropriate cell dilution in M199 containing 0.5% fetal calf serum (Invitrogen‐Gibco, Milan, Italy) were incubated at 37°C in a 5% CO_{2} humidified atmosphere for 12–14 h on glass coverslips and observed with an inverted photomicroscope (model DM IRB HC; Leica Microsystems, Bensheim, Germany). Phase contrast snap photographs and 12–14 h long movies (one frame every 5 min) were taken with a cooled digital CCD Hamamatsu ORCA camera (Hamamatsu Photonics Italia, Rome, Italy), recorded and analyzed with ImageProPlus 4.0 imaging software (Media Cybernetics, Leiden, The Netherlands). Cell paths were generated from centroid positions and migration parameters were computed with DIAS software (Solltech Inc., Oakdale, IA). For blocking experiments, an anti‐VEGF neutralizing antibody (20 μg/ml; Pepro Tech Inc., Rocky Hill, NJ) or 3 U/ml heparin (Sigma‐Aldrich) plus 20 nM VEGF‐A165 (R&D System, Abingdon, UK), which saturates both VEGF receptors ‐1 and ‐2 (Neufeld *et al*., 1999), were added to the medium where cells were resuspended before plating them on Matrigel. Irrelevant immunoglobulins (Sigma‐Aldrich) of the same species and heat‐inactivated VEGF‐A165 (5 min at 95°C) were employed for control purposes. To measure the half‐time of VEGF‐A, 5 ng of VEGF‐A_{165} and 10 nCi of [^{125}I]VEGF (Amersham‐Pharmacia‐Biotech; Milan, Italy; 1200 Ci/mmol) were incubated for different intervals of time with EC plated on Matrigel. VEGF‐A was immunoprecipitated from the medium pre‐cleared with an anti‐goat Ig (Sigma) with a polyclonal anti‐VEGF‐A antibody (Santa‐Cruz Biotech., Santa Cruz, CA). Radioactivity corresponding to the VEGF‐A band in SDS–PAGE (12%) was counted and used to calculate the half‐time by EnzFitter software (Biosoft, Cambridge, UK). The half‐time calculated was not significantly modified by the cell number, being 64 ± 7 (*n* = 6) and 60 ± 3 min (*n* = 3), respectively, with 125 and 400 cells/mm^{2}.

To measure VEGF‐A production ECs (2.5 × 10^{6}/150 mm diameter Petri dish) were plated on Matrigel‐coated surface as described above. As control M199 supplemented with 0.5% fetal calf serum was treated in the same conditions. After 3 h, supernatants were recovered, centrifuged (13 000 *g*, 5 min at 4°C) and stored at −80°C. Pool of three dishes were concentrated by Centriprep YM‐10 (Millipore S.p.A., Milan, Italy) and the VEGF‐A was quantified by QuantiGlo Chemiluminescent Sandwich ELISA (R&D Systems, sensitivity 1.7 pg/ml).

*Measurement of capillary network dimension*

Chord lengths were measured directly from node center to node center on experimental photographs and the resulting values averaged. Alternatively, from digital photographs a density function *n _{i,j}* was obtained, where

*i,j*are pixel coordinates and

*n*is 1 above a suitably chosen luminosity threshold, 0 below. The radial part of the correlator

_{i,j}*P*

_{kl}= Σ

_{ij}n

_{ij}n

_{i+lj+k}has been computed and the position of its absolute minimum ℓ

_{o}taken as a characteristic measure of the structure size. We checked the approximate coincidence ℓ̅ ≈ ℓ

_{o}suggesting that for vascular‐type structures the mean chord length can be extracted directly from

*l*

_{o}. A similar procedure has been followed in the analysis of numerical experiments where it was used to obtain results averaged over a great number of random realizations.

*Numerical simulations*

Numerical simulations were performed discretizing the system of Equations 1–3 on squares of 128^{2}, 256^{2} and 512^{2} equispaced nodes with periodic boundary conditions. A finite volume numerical scheme of Godunov type (LeVeque, 1990) has been adopted for the hyperbolic part of the equation (i.e. transport and pressure) and a simple centered discretization used for the chemotactic term. Equation 3 has been discretized using a spectral method and solved by Fast Fourier Transform.

*Statistical analysis*

Trajectory parameters were fitted by the maximum likelihood method to the Beta distribution *p*(*x*) = *C*(*x* + 1)^{a − 1} (*x* −1)^{b − 1} adapted to the interval −1< *x* <1, where *C* is the corresponding normalization factor (Grant and Ewens, 2001). Differences were weighted by χ^{2} goodness‐of‐fit test. Where indicated, significance was estimated by *t*‐test.

*Supplementary data*

Supplementary data are available at *The EMBO Journal* Online.

## Supplementary Information

Supplementary data [emboj7595085-sup-0001.pdf]

Supplementary movie [emboj7595085-sup-0002.mov]

## Acknowledgements

The authors thank N.Bellomo, M.Gasparini (Polytechnic of Torino), I.Kolokolov (Landau Institute, Moscow), P.A.Netti (University of Naples) and M.Vergassola (CNRS, Nice) for discussions. The work was supported by A.I.R.C., C.N.R., I.S.S. (Programma Nazionale di Ricerca sull'AIDS: grants 40B.19 and 30B.9; Programma Terapia dei Tumori), M.I.U.R. and EU.

## References

- Copyright © 2003 European Molecular Biology Organization