1 Introduction

The vascular system is one of the most important organs in large multicellular organisms. During embryogenesis, the need for an efficient transport of nutrients and waste products arises naturally as the increasing size of the organism makes diffusion less and less efficient. A similar situation arises in solid tumors, but instead of the de novo formation of the vasculature (vasculogenesis), tumors are able to remodel the vascular network in the surrounding tissues and promote formation of new branches from the existing ones (angiogenesis) [7]. Understanding the basic principles behind these processes could help in controlling tumor growth, as well as improving restoration of normal vasculature after trauma such as a stroke. It is tempting to speculate that the same basic principle that generally promotes network formation during vasculogenesis also drives blood vessel sprouting during angiogenesis.

To explore mechanisms behind vasculogenesis, it is beneficial to turn to computational modeling where all parameters are under control, before attempting experimental validation where unknown factors could complicate the interpretation of the results. Various hypotheses have been constructed with the aid of computational biology [3]; however, we argue that the best modeling approach to study network formation is through cell-based models [13]. Such models are well suited to describe collective behaviors emerging from cellular properties and interactions, such as collective cell motion [11], or various morphogenetic processes [14]. Vascular networks contain structures typically on the scale of a couple to tens of cells which would invalidate continuum-based approaches that are better fitted for describing cellular structures of at least hundreds of cells in size.

Previous studies using cell-based modeling have identified a handful of candidate mechanisms for sprouting and network formation [4]. By the nature of biological cells, all of these mechanisms include a form of long-range attraction that ensures cell cohesion and a short-range repulsion that is needed to represent cell incompressibility and thus prevent unrealistically high cell densities. One class of models explains the long-range attraction via mechanotaxis, whereby cells contract and deform an underlying soft substrate, pulling cells toward each other [10] or where the substrate locally stiffens in response to cellular contraction directing cells toward one another via durotaxis [16]. This mechanism suffices to explain how in some culture systems, for example bovine aortic endothelial cells on poly-acrylamide gels [16], networks form only on substrates of specific stiffness. Contractile cells on deformable substrates have been described as force dipoles revealing that such dipoles tend to form locally or globally ordered structures depending on the Poisson ratio of the substrate [1, 2]. An alternative mechanism, based on the observation that cells preferentially adhere to elongated cellular structures, explains how networks can form on bare substrates [19, 20]. Another body of work focuses on chemotaxis toward a self-secreted chemokine, in combination with either contact inhibition of chemotaxis [15] or cell elongation [12], that proves to be a robust mechanism albeit the chemokine remains to be identified. While the different mechanisms of attraction distinguish between different hypotheses, the mechanism of repulsion is less important and is either introduced explicitly (for example, in partial differential equation models [10] or particle-based models [20]) or is an implicit property of the model (as in the CPM [12, 16] for instance).

Fig. 1
figure 1

Cellular Potts simulations of chemotaxis and cell elongation models of vasculogenesis. a A simple mechanism of chemotaxis toward a self-secreted chemokine generates cell aggregates. b Chemotaxis together with cell elongation gives rise to robust networks. c Cell elongation without chemotaxis is sufficient to explain network formation. Used with permission from [17]

The model mechanism we focus on in this paper demonstrates that cell elongation and volume exclusion together with contact-dependent cell adhesion [17] or longer range cell–cell attraction [12] are sufficient for network formation (Fig. 1). Whereas slightly adherent, dispersed cells aggregate into compact clusters (Fig. 1a), elongated cells form network-like structures in the models with (Fig. 1b) or without (Fig. 1c) additional longer range chemotaxis. In these models of elongated cells, the branches are stabilized by the increasing rotational inertia of growing clusters of elongating cells: a tightly packed branch of elongated cells is harder to displace or reorganize, and therefore cells are frozen into the branched pattern in the model. A somewhat similar cellular process to elongation described by this model is the apical constriction of epithelial cells that plays a major role in, for example, gastrulation and neurulation [18]. Although apical constriction leads to cell elongation, this elongation may only be a passive result of volume conservation and cytoplasmic flow [8].

As each model implementation is a simplification with different limitations, it is worth to examine the model at hand using more than one implementation with different implicit assumptions associated with the simulation methodology. If the hypothesis is independent of the implementation, it should be able to reproduce the phenomenon in different implementations. For example, this has been performed in the case of the preferential adhesion hypothesis, which has been tested in both a lattice-free model [20] and a grid-based model [19].

The previous studies [12, 17] used the same cellular Potts model (CPM) description which could introduce model-specific artifacts. For example, highly elongated cells could breach the limitations of cell representation within the CPM, as cell elongation is achieved by maximizing the largest moments of inertia of the cell [12], which leads to unrealistic thickening of the cell body at its extremities. As artifacts may result from the grid-based nature of the model and the implementation of cell elongation, we chose to implement the cell elongation hypothesis in a lattice-free model using a generalized attraction and repulsion between cells, describing a class of models. We use this model to test if the hypothesis holds and to compare our results to previous reports.

2 Methods

2.1 Computational model

Each cell i in the model is described as an ellipse on a 2D plane, characterized by its position \(\mathbf {x}_i\), direction of elongation \(\phi _i\), area \(A_i\), and the aspect ratio of major-to-minor axis \(s_i\). This ellipse is a repulsive core representing the incompressible cell body (Fig. 2a, blue) and is surrounded by a larger concentric ellipse that is responsible for cell–cell attraction in the model (for example, via filopodial adhesion) (Fig. 2a, red). Cell motion is modeled as a persistent diffusion process as in a previous study [20]. The change in velocity \(\mathbf {v}_i(t)\) for cell i at time t is described as

$$\begin{aligned} \frac{d\mathbf {v}_i(t)}{dt}=\frac{1}{m_i}\left( -\alpha \mathbf {v}_i(t)+N_v\mathbf {\xi }_v dt^{-0.5}+\sum _{j\ne i}^{j} \frac{\mathbf {x}_i-\mathbf {x}_j}{|\mathbf {x}_i-\mathbf {x}_j|} F_{ij}\right) . \end{aligned}$$
(1)

Here \(m_i\) is the mass of cell i (\(m_i=m=1\)) and \(\alpha \) is a damping parameter. \(\mathbf {\xi }_v\) is a uniformly distributed random vector in the 2D plane (\(\mathbf {\xi }_v \in [-0.5:0.5] \times [-0.5:0.5]\)), and model parameter \(N_v\) sets the amplitude of this translational noise. The last term is the sum of pair interactions between cells describing a short-range core repulsion and a long-range but finite attraction:

$$\begin{aligned} F_{ij}=\lambda _r A_{r} - \lambda _a A_{a}. \end{aligned}$$
(2)

Here \(A_r\) is the overlap area between the two repulsive ellipses of the two cells (Fig. 2b, blue area), \(A_{a}\) is the overlap area of the two outer ellipses (Fig. 2b, red area), and \(\lambda _{r}\) and \(\lambda _{a}\) are model parameters. The overlap is calculated using a previously published method [9].

Fig. 2
figure 2

Cells and their interaction in the model. a Cells are represented as 2D ellipses consisting of a repulsive core as cell body (blue). Cell attraction occurs at a finite range around the core in an area represented as a bigger ellipse surrounding concentrically the repulsive core (red). b Cell–cell interaction is governed by the overlap area of the two inner ellipses (blue area) and the overlap of the two outer ellipses (red area). (Color figure online)

Model cells are rotated to minimize overcrowding and maximize attraction [2]. For every interval \(\Delta t\), an attempt is made to change the orientation of N randomly selected cells with a random angle \(\xi _a\sqrt{\Delta t}\), where \(\xi _a \in [-\pi :\pi ]\) is a uniform random variable. The change is accepted with a turn-probability

$$\begin{aligned} P_i= \min \left[ 1,\exp \left\{ \frac{1}{N_a}\left( \sum _{j\ne i}F_{ij} - \sum _{j\ne i}F_{ij}'\right) \right\} \right] , \end{aligned}$$
(3)

where \(F_{ij}'\) is the pair interaction in a configuration after the proposed rotation and \(N_a\) is the angular noise parameter.

Simulations were initiated with \(N=1000\) overlapping cells distributed at random within a \(200\;\mu \mathrm {m}\times 200\;\mu \mathrm {m}\) unit area. Cell area was fixed at \(250\;\mu \mathrm {m}^{2}\) throughout the study. The size of the outer “adhesive” ellipse is set by the attraction radius parameter defined as \(R_a = a_{a} / a\), where \(a_a\) is the semi-major axis of the adhesive (outer) ellipse and a is the semi-major axis of the repulsive (inner) ellipse, and the aspect ratio of the outer and inner ellipses is the same. Baseline parameters used are \(R_a=1.5\), \(\lambda _a=0.0006, \lambda _r=0.02, \alpha =0.4, N_v=0.4\), \(N_a=1\), and \(s=7\). These parameters were estimated by macroscopic inspection of the model to yield “realistic” pattern formation as a measure of validation. To prevent excessive cell overlap, the magnitude of the repulsion strength (\(\lambda _r\)) is two orders of magnitude larger than the magnitude of the attraction strength (\(\lambda _a\)). The model’s sensitivity to the parameters was tested by altering the parameter values until unrealistic results, either concerning cell movement or cell interactions, were produced. Model Eqs. 1 and 2 are integrated numerically using the forward Euler method with a fixed time step of \(\Delta t=0.1\) (or \(\Delta t=1\) for the longer simulations of up to \(t=100,000\) shown in Fig. 5), and orientations are updated after each iteration synchronously using Eq. 3. Integration was stopped at \(t=25,000\) (or at \(t=100,000\) for simulations shown in Fig. 5). In each case, the system reached a quasi-stationary state at the end of the simulations. An implementation of this particle-based system is provided in Online Resource 3.

2.2 Order parameter

Local alignment of elongated cells has been shown to play an important role in the formation of networks previously [17]. To measure alignment, a local orientational order parameter is used:

$$\begin{aligned} S(r)=\langle \cos \left( 2\theta _i(r)\right) \rangle _i. \end{aligned}$$
(4)

Here \(\theta _i(r)\) is the polar angle between \(\phi _i\), the orientation of cell i, and \(\eta _i(r)=\langle \phi _j \rangle _{\left\{ j: |\mathbf {x}_j - \mathbf {x}_i| \le r\right\} }\), the average cell orientation within r distance of cell i. S(r) is isotopic and takes the value of 0 for randomly oriented cells and 1 for perfectly aligned cells. Smaller values of the order radius r describe the alignment in the close vicinity of the cell (\(r=20\)), or in the local structures (\(r=80\), Fig. 3), while any global order is captured by the global order parameter, \(S_g=S(r \rightarrow \infty )\), where all the cells are considered for the calculation of \(\eta _i\). The software used for measuring the order parameter is provided in Online Resource 3.

Fig. 3
figure 3

Ranges of local order. Two values of r for the order parameter calculation shown for the center cell within a typical cell arrangement. While \(r=20\) describes order at the local cell neighborhood, \(r=80\) is an indication of order in multicellular structures in the vicinity

3 Results

3.1 Elongated cell shape induces network formation

Rounded cells (\(s=1.5\)) aggregate into less-ordered clusters (Fig. 4a and the animation in Online Resource 1), while more elongated cells with an aspect ratio \(s=7\) in the simulation interconnect to form elongated branches and networks (Fig. 4b and the animation in Online Resource 2). The time evolution of the order parameters calculated for two simulations of round and elongated cells shown in Fig. 4 demonstrates that the initially disordered cells align within \(t=10^4\) and keep a similar level of order for at least 10 times longer at finite length scales (S(20) and S(80) on Fig. 5a). Positive value of the global order parameter (\(S_g\)) in the elongated cell configurations results from finite system size and indicates the emergence of a system-wide alignment (Fig. 4c). Elongated cells (\(s=7\)) produce a significantly higher order than more round cells (\(s=1.5\)) at all length scales (Fig. 5b). Order at the cell-to-cell range (S(20)) always supersedes order at the multicellular level (S(80)), which is always higher than global order (\(S_g\)).

Fig. 4
figure 4

Cell elongation induces network formation. a Short cells with an aspect ratio \(s=1.5\) form compact aggregates. b Elongated cells with an aspect ratio of \(s=7\) form networks. c Temporal evolution of network formation for cells with an aspect ratio of \(s=7\). Global order is established at about \(t=10^4\). The system reaches a fluctuating steady state: the morphological details change over time, but the morphological features remain statistically constant

Fig. 5
figure 5

Order parameters of elongated and rounded cell simulations. a Temporal evolution of the order parameter, S(r), for \(r=20\) (green curves), \(r=80\) (red curves), and \(r\rightarrow \infty \) (blue curves) for two systems of more rounded cells (\(s=1.5\), dotted lines) and elongated cells (\(s=7\), solid lines; \(\Delta t=1\)). Ribbons show standard deviations over 10 runs. Finite scale order reaches a plateau at \(t=10^4\) and remains approximately stable up to at least \(t=10^5\). b, Order parameters from 10 independent simulations show that alignment at \(t=10^5\) is higher for elongated cells than for rounded cells at all length scales (\(\Delta t=1\)). Within configurations, order at shorter scales is always higher than order at longer scales. (Color figure online)

Simulations with cells of increasing aspect ratios and constant areas reveal that order is gradually increased in simulations with more elongated cells (Fig. 6). When simulations with increasing s are compared, order emerges first at short length scale (S(20)) followed by longer range order. Global order appears in simulations with \(s\ge 4\). The high variation of the global order among the different simulations is consistent with the fluctuation shown in Fig. 5a.

Fig. 6
figure 6

Emergent order in the model is an increasing function of cell elongation. a Simulations with uniform cells of various aspect ratios show that slightly elongated cells align locally, while more elongated cells give rise to order at longer length scales. Ribbons show standard deviations over 10 runs after \(t=2.5\times 10^4\). b–d Configurations at \(t=2.5\times 10^4\) for the aspect ratios \(s=1.5\), \(s=4.5\), and \(s=7\)

3.2 Increased range of cell attraction enhances cell alignment

Previous studies showed that chemotaxis could play an important role in vascular network formation [12, 15]. Since chemotaxis may be interpreted in our model as an adhesive interaction over a longer range, we investigated how the range of attraction affects the observed cell alignment.

Fig. 7
figure 7

Increasing range of attraction increases order. a Local, intermediate, and global orders are higher with larger attraction ranges for cells with an aspect ratio of \(s=7\). Ribbons show standard deviations over 10 runs after \(t=2.5\times 10^4\). Configurations at \(t=2.5\times 10^4\) for b short, c medium, and d long attraction ranges \(R_a=1.1\), \(R_a=1.6\) and \(R_a=2\)

As expected, a short-range attraction (\(R_a=1.1\)) results in a dissociated cell configuration and low order (Fig. 7a and Fig. 7b, \(R_a=1.1\)), compared to the control parameters (Fig. 4b). A long attraction range, however, results in a marked increase in ordering, consistent with the longer range attraction of the chemotaxis studies. Interestingly, global order increases markedly at an interaction range of less than twice the cell size, where long, aligned strands of cells are formed (Fig. 7d, \(R_a=2\)).

Increased alignment is also achieved by increasing attraction strength through parameter \(\lambda _a\) (Fig. 8), while increased repulsion leads to dissociation and consequent loss of order (Fig. 9).

Fig. 8
figure 8

Stronger attraction leads to increased order. a Order parameters from simulations with increasing attraction parameters for cells with an aspect ratio of \(s=7\). b Simulation configurations with low attraction strength \(\lambda _a=4\times 10^{-4}\). c Simulation configurations with high attraction strength \(\lambda _a=2\times 10^{-3}\)

Fig. 9
figure 9

Repulsion acts against order. a Order parameters from simulations with increasing repulsion parameters for cells with an aspect ratio of \(s=7\). b Simulation configurations with a low repulsion strength of \(\lambda _r=10^{-2}\). c Simulation configurations with a high repulsion strength of \(\lambda _r=9\times 10^{-2}\)

3.3 Noise is not essential for alignment

Fig. 10
figure 10

Effect of noise on alignment. a Order parameters as a function of vectorial noise for cells with an aspect ratio of \(s=7\). b Simulation configurations in the absence of translational noise. c Simulation configurations with increased translational noise (\(N_v=3\)). d Order parameters as a function of angular noise for cells with an aspect ratio of \(s=7\). e Simulation configurations with low angular noise (\(N_a=0.1\)). f Simulation configurations with increased angular noise (\(N_a=2\)). Order parameters and configurations are shown at \(t=2.5\times 10^4\) simulation time

We assessed the importance of noise for cell alignment in our model (Fig. 10). In the absence of translational noise (Fig. 10b) or low rotational noise (Fig. 10e), cells align and form networks. At high noise levels, the two types of noise act in a similar manner; high translational noise results in a decay of global order (Fig. 10a, c), while high angular noise allows for more energy inefficient rotations, thus decreasing global order (Fig. 10d, f). At low noise levels, the two types of noise act differently; while low translational noise has no or very little effect on the orientational order (Fig. 10a), low angular noise results in higher local order and lower global order (Fig. 10e), suggesting the formation of isotropic branches that point in different directions.

3.4 Cell elongation is essential for alignment

Finally, we tested whether cell elongation is required for alignment, by influencing the parameters that had positive impact on cell order in the previous simulations. Using rounded cells with an aspect ratio of \(s=1.5\), three scenarios with increased range of attraction (\(R_a\)), increased attraction strength (\(\lambda _a\)), and reduced repulsion strength (\(\lambda _r\)) were tested. Order parameters from these simulations together with order parameters from simulations of round and elongated cells (Fig. 4a, b) are summarized in Table 1. These indicate that neither increased adhesion strength or range, nor decreased repulsion, is able to order the cells to an extent comparable with elongated cell simulations. Without elongation, cell clusters are unable to deviate from their compact aggregates (Fig. 11).

4 Discussion

A simple particle-based model is introduced to demonstrate that elongation can indeed aid cells to aggregate into a network. The mechanism for network formation is simple: Cell elongation, together with attraction and core repulsion, leads to local cell alignment and increase in orientational order. The fact that this prediction holds in two completely different model implementations, our particle-based model and the CPM [12, 17], gives confidence that the predictions are due to the explicit model assumptions, not due to unintended properties of the simulation methodology. Time evolution of order and dependence on elongation (Fig. 5b) are similar in both implementations. Increasing the attraction range of cells in the current implementation allows the study of the transition between the pure elongation hypothesis [17] to the chemotaxis and elongation hypothesis [12]. Extended attraction range, in the form of chemotaxis, has been shown to help in the formation of regular networks with a more defined pattern size across the whole system [12]. Consistent with this previous result, we found that a longer range of attraction results in higher order and a more network-like structure (Fig. 7).

Table 1 Order parameters for round (\(s=1.5\)) and elongated (\(s=7\)) cells, compared with cells with an aspect ratio of \(s=1.5\) at \(t=2.5\times 10^4\) tested with parameters that favor cell ordering
Fig. 11
figure 11

Effect of parameters on alignment without cell elongation. a Simulation configurations with increased radius of attraction \(R_a=2\). b Simulation configurations with increased attraction strength \(\lambda _a=2\times 10^{-3}\). c Simulation configurations with low repulsion strength \(\lambda _r=10^{-2}\). Configurations shown at \(t=2.5\times 10^4\). All cells have an aspect ratio of \(s=1.5\)

Cell movement dynamics is a marked difference between the CPM and our model. In the CPM, cell movement emerges from the displacement of the cell boundaries, yielding a more amoeboid cell movement. By contrast, in particle-based models model cells translate as a single unit, including translocation and rotation in our case. Cell movement in the model is described by an overdamped dynamics, where the damping factor (\(\alpha \) in Eq. 1) controls the movement persistence of cells. For large values of \(\alpha \), the left-hand side of Eq. 1 vanishes and Eq. 1 can be approximated by an algebraic equation giving the cell velocities as a function of the external forces acting on the cells. We have chosen here to work with a differential equation description, thus keeping the damping factor in our model as an explicit parameter. Although we have here only studied overdamped kinetics, smaller values of \(\alpha \) could mimic persistent cell motility: it takes time for a cell to change its direction if, e.g., the chemoattractant gradients change, the actin cytoskeleton needs to reorganize which takes time. Such persistent motion was also explicitly included in early chemotaxis-based partial differential equation models of vascular network formation [5]. Although persistence of motility was later shown to be unnecessary for network formation, it was argued that in the CPM the cell shape may cause some directional persistence (Fig. 11 of Ref. [12]). In the present work directional persistence is severely reduced and is largely independent of the cell shape, suggesting that persistence is not required for network formation.

Stochasticity in our model is introduced through the angular noise similar to the noise of the Vicsek model [21] and the translational noise of the Grégoire–Chaté model [6]. Note that these noise factors describe fluctuations in the system, and therefore the considered noise factors are independent from one another and from the noise in the past. Our results from the study of the angular noise show that noise hinders the global order, in good agreement with previous particle-based simulations using only angular noise and force dipoles to describe multicellular structure formation [2].

Interestingly, an increase in both attraction range (Fig. 7) and attraction strength (Fig. 8) leads to a higher tendency to global ordering as packing of the cells becomes tighter. This is similar to the previous observation in the CPM models [12, 17], where the more adhesive cells become more compact and highly ordered in domains, but without the emergence of global order. In the CPM, cell shape is an emergent property of the cells and therefore cells are able to deform in order to maximize the contact surface and compactness (see Fig. 4a of [17]). In contrast, cells in our implementation are unable to deform such that they might be less able to accommodate the boundaries of such locally ordered domains. This may force them to align globally through the compacting force of the strong attraction. Secondly, when attraction range is increased to a distance of two cell diameters in the current implementation, alignment order spans through domains of at least 1000 cells, with strands of cells aligned in parallel even without tight packing (Fig. 7d, \(R_a=2\)). An effect contributing to this global ordering may be the implementation of noise in the current particle-based implementation. In the CPM, at high densities, cells hinder each other’s translation and rotation, resulting in an increasingly slow development of the pattern as the branches grow [17]. By contrast, in the particle-based model translation is not hindered by adjacent cells and occurs for cells as a whole (Eq. 1), while cells continue to rotate independently as a whole as long as conflicts with adjacent cells persist (Eq. 3).

In conclusion, here we introduced a particle-based model to re-examine a hypothetical mechanism for the formation of microvascular networks, i.e., that elongated vascular cells tend to aggregate into branches of a network structure. Our previous work [12, 17] simulated this potential mechanism using the CPM, demonstrating that indeed elongated cell shape, in combination with mutual attraction or adhesion, suffices for the formation of network-like patterns. Here we have shown that this phenomenon also occurs in a lattice-free particle-based model, adding confidence that the effect is not caused by artifacts of the CPM or of the present particle-based model. Thus our models suggest that network formation is a natural emergent property of elongated, adhesive objects in a stochastic system. Because a range of alternative mechanisms for network formation and angiogenic sprouting have been suggested (reviewed in [14] and [3], see also [16]), to what extent the present mechanism contributes to angiogenesis in vivo or in vitro at this point must remain an open question. These and similar studies, however, help generalize and categorize the main requirements for network formation and angiogenic sprouting, thus contributing to ongoing efforts to identify the controlling factors of angiogenesis from a biophysical point of view.