Particle-based simulation of ellipse-shaped particle aggregation as a model for vascular network formation

Computational modelling is helpful for elucidating the cellular mechanisms driving biological morphogenesis. Previous simulation studies of blood vessel growth based on the Cellular Potts model (CPM) proposed that elongated, adhesive or mutually attractive endothelial cells suffice for the formation of blood vessel sprouts and vascular networks. Because each mathematical representation of a model introduces potential artifacts, it is important that model results are reproduced using alternative modelling paradigms. Here, we present a lattice-free, particle-based simulation of the cell elongation model of vasculogenesis. The new, particle-based simulations confirm the results obtained from the previous Cellular Potts simulations. Furthermore, our current findings suggest that the emergence of order is possible with the application of a high enough attractive force or, alternatively, a longer attraction radius. The methodology will be applicable to a range of problems in morphogenesis and noisy particle aggregation in which cell shape is a key determining factor.


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 tumours, but instead of the de novo formation of the vasculature (vasculogenesis), tumours 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 tumour 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 modelling where all param-eters 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 modelling approach to study network formation is through cell-based models [13]. Such models are well suited to describe collective behaviours 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 continuumbased approaches, that are better fitted for describing cellular structures of at least hundreds of cells in size.
Previous studies using cell-based modelling 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 towards each other [10] or where the substrate locally stiffens in response to cellular contraction directing cells towards one another via durotaxis [16]. This mechanism suffices to explain how in some culture systems, for example bovine aortic endothelial cells on polyacrylamide 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 towards 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).
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 Fig. 1 Cellular Potts simulations of chemotaxis and cell elongation models of vasculogenesis. a, A simple mechanism of chemotaxis towards 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] models of elongated cells, the branches are stabilised by the increasing rotational inertia of growing clusters of elongating cells: a tightly packed branch of elongated cells is harder to displace or reorganise, 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 [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 artefacts. 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 artefacts 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. 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)

Computational Model
Each cell i in the model is described as an ellipse on a 2D plane, characterized by its position x i , direction of elongation φ 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 modelled as a persistent diffusion process as in a previous study [20]. The change in velocity v i (t) for cell i at time t is described as: Here , and model parameter N v sets the amplitude of this translational noise. The last term is the sum of pair interactions between cells describing short-range core repulsion, and a long-range but finite attraction: A r is the overlap area between the two repulsive ellipses of the two cells (Fig. 2b, blue area), and A a is the overlap area of the two outer ellipses (Fig. 2b, red area), and λ r and λ a are model parameters. The overlap is calculated using a previously published method [9]. Model cells are rotated to minimize overcrowding and maximize attraction [2]. For every interval ∆t an attempt is made to change the orientation of N randomly selected cells with a random angle ξ a √ ∆t, where ξ a ∈ [−π : π] is a uniform random variable. The change is accepted with a turn-probability where F i j 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 µm × 200 µm unit area. Cell area was fixed at 250 µ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, λ a = 0.0006, λ r = 0.02, α = 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 (λ r ) is two orders of magnitude larger than the magnitude of the attraction strength (λ 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 equations 1 and 2 are integrated numerically using the forward Euler method with a fixed time-step of ∆t = 0.1 (or ∆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 on 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.

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:

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 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 multi-cellular level (S(80)), which is always higher than global order (S g ). 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.

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. 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). 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 rounded cells at all length scales (∆t = 1). Within configurations order at shorter scales is always higher than order at longer scales Increased alignment is also achieved by increasing attraction strength through parameter λ a (Fig. 8), while increased repulsion leads to dissociation and consequent loss of order (Fig. 9).

Noise is not essential for alignment
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 similar manner; high translational noise results in a decay of global order (Figs. 10a and 10c), while high angular noise allows for more energy inefficient rotations, thus decreasing global order (Figs. 10d and 10f). 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.

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 aspect ratio of s = 1.5, three scenarios with increased range of attraction (R a ), increased attraction strength (λ a ) and reduced repulsion strength (λ 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).

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]. Consistently with this previous result, we found that a longer range of attraction results in higher order and a more network-like structure (Fig. 7).
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 (α in Eq. 1) controls the movement persistence of cells. For large values of α, 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 α 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 (Figure 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, both an increase in 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 maximise 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 cellular Potts model, 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 cellular Potts model 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 to 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.
Online Resource 1 Simulation video of rounded cells for t = 2.5 × 10 4 . SImulated cells of aspect ratio s = 1.5 form aggregates. Every frame is after t = 50. Every second of video is t = 1250; see https://youtu.be/BgQwd1aGxAI Online Resource 2 Simulation video of elongated cells for t = 2.5 × 10 4 . Simulated cells of aspect ratio s = 7 form networks. Every frame is after t = 50. Every second of video is t = 1250; see https://youtu.be/9zUGEBk6pio Online Resource 3 C++ implementation of the presented particle-based system, and software for measuring the order parameter