Self-assembly, buckling and density-invariant growth of three-dimensional vascular networks

The experimental actualization of organoids modelling organs from brains to pancreases has revealed that much of the diverse morphologies of organs are emergent properties of simple intercellular ‘rules’ and not the result of top-down orchestration. In contrast to other organs, the initial plexus of the vascular system is formed by aggregation of cells in the process known as vasculogenesis. Here we study this self-assembling process of blood vessels in three dimensions through a set of simple rules that align intercellular apical–basal and planar cell polarity. We demonstrate that a fully connected network of tubes emerges above a critical initial density of cells. Through planar cell polarity, our model demonstrates convergent extension, and this polarity furthermore allows for both morphology-maintaining growth and growth-induced buckling. We compare this buckling with the special vasculature of the islets of Langerhans in the pancreas and suggest that the mechanism behind the vascular density-maintaining growth of these islets could be the result of growth-induced buckling.


Introduction
Tubes are ubiquitous features of numerous biological systems.In humans, they form the gastrointestinal tract, the ductal network of the pancreas, the fallopian tubes, the urinary tract and so on, with the most obvious example being the entire vascular network of blood vessels.On the relevant time scales of multicellular energy consumption, diffusion is limited to delivering metabolites over length scales smaller than ∼100 μm.Instead, on larger length scales, tissue needs some form of directed transport [1].In vertebrates, this active transport is provided by the beating heart through the vascular network, which in turn has to branch into every part of the organism to nourish tissue and remove waste.
The development of the vascular network involves mainly two processes: vasculogenesis and angiogenesis [2,3].During vasculogenesis, individual endothelial cells coalesce and de novo form functional vessels [4][5][6].Studies of vasculogenesis in vitro have been mainly restricted to two dimensions [7], but recently three-dimensional vascular organoids have been produced [8].Vasculogenesis results in a randomly connected vascular plexus, which is subsequently remodelled by pruning or branching [9][10][11][12] to a mature vascular network, e.g. with a hierarchical tree-like structure.In angiogenesis, the tree-like structure is formed by branching processes involving either splitting (intussusception) or sprouting dynamics from already formed blood vessels [2,13,14].This remodelling can be guided by blood flow, pressure and vessel wall stresses [15].We shall be interested in modelling blood vessel organoids and will thus not consider this latter reorganization, which becomes relevant only in connection with certain organs (such as a pumping heart).
From a theoretical and computational viewpoint, the most intriguing feature of vasculogenesis is its three-dimensional self-assembly of tubular networks.Of equal importance is whether these self-assembled networks percolate across the tissue, i.e. whether a fully connected network of tubes is formed.What densities of endothelial cells are needed in three dimensions to ensure these criteria?We will additionally be interested in questions of growth.Once a network is formed, can this network undergo stable growth?And what are the possible mechanisms for such networks to grow while maintaining a constant space-to-vessel density?
Understanding blood vessel formation computationally has received much attention [16].Continuum models enable descriptions of density fields of chemotaxing endothelial cells during vasculogenesis [7,[17][18][19].Likewise, cellular Potts models [20] and models of individual cells [21] have been employed.These studies of vasculogenesis have focused mainly on two-dimensional systems.In this paper, we introduce a coarse-grained description of tubes in three dimensions using a formulation that resolves features of both single cells and full organs (vessel network).In particular, we are able to simulate vascular networks comprising up to hundreds of thousands of particles.
Our focus will be on emergent features of the model such as vasculogenesis and buckling during growth, but we note that our model also has the ability to describe angiogenetic sprouting [22] (budding) and intussusception, which on the cell level resembles gastrulation [23].While lumen formation in blood vessels is its own research field [5,[24][25][26], we introduce a simple mechanism for lumen formation by describing the evolution of the apical-basal polarity of cells, thus yielding a fully emergent approach to tubulogenesis.
The paper is organized as follows.In §2, we introduce the methodology and mathematics of the model and demonstrate lumen formation.Section 3 is devoted to vasculogenesis and the percolation of the vascular network.In §4, we study the growth of vascular networks for various parameters of the model and show that both morphology-maintaining and buckling growth patterns can arise.Lastly, in §5, we describe the vasculature of the islets of Langerhans in the pancreas and describe how their tortuous features could be the result of buckling during growth.In particular, we show that the vascular density during the growth of these islets could be maintained simply by a buckling mechanism without the need for angiogenetic processes.

Model, polarities and lumen formation
Contrary to two-dimensional models, in three dimensions, cell polarity is crucial to model cell sheets.Our coarse-grained model describes a collection of particles/cells each defined by their position x, their apical-basal polarity (AB) p, and their planar cell polarity (PCP) q, illustrated in figure 1.Our model is coarse grained in the sense that a collection of particles model a cell, and as such, even though each particle is a sphere, shape deformations are possible in a collection of particles.In a tube, such as a blood vessel, the AB polarity of cells will define the inside versus the outside of the vessel, while PCP defines the direction around the tube versus the direction along the tube.
To model cell behaviour, we use a slightly modified version of the model of Nissen et al. [23].In this model, particles interact pairwise only if they are line-of-sight Voronoi neighbours and their mutual potential energy is for which where We keep β = 5, which sets the inter-particle spacing to ≃2 units, and furthermore enforce λ 0 + λ 1 + λ 2 + λ 3 = 1 with λ 1 ≥ λ 3 .This strikingly simple model can describe a plethora of phenomena related to polarity-driven morphogenesis in organoids [23].Naturally, real interactions between polarized cells will be much more complex than the model used here; indeed, these depend on the precise distribution of the surface proteins that make up the polarity of the cells.The interactions used here can be thought of as the first relevant and symmetryobeying terms that give rise to polarity-aligning cells.The simplicity of the present model is thus agnostic towards the underlying microscopic details.Here we introduce a small extension to this model that permits the de novo formation of tube-like structures.
The dynamics of the model follows from taking all mobilities to be equal.Hence, with the norms of p and q kept at unity.For this study, the model was implemented using PYTORCH and run with CUDA-acceleration.
With only spherical interactions, i.e. λ 0 = 1.0, a solid tube is a meta-stable structure of this model, as shown in figure 2a.Lumen formation corresponds to the formation of AB polarity, i.e. the discrimination of the inside and outside of the tube.Various methods for lumen formation exist, e.g.
. The polarity of each particle, which stems from a distribution of proteins, is modelled as vectors.AB polarity is indicated by p and PCP by q.For a tube, as illustrated, p points away from the tube (distinguishing outside from inside) and q curls around the tube (distinguishing along versus around the tube).Interactions between cells depend on their orientation between each other, rij .In equations (2.3), S 1 favours p i and p j parallel and both orthogonal to rij , S 2 favours p's and q's orthogonal and S 3 favours q i and q j parallel and both orthogonal to rij .Note that, in reality, cells can deform based on the polarities, but we model them as point particles.Shape deformation is captured by collections of multiple particles.(Online version in colour.) extracellular cord hollowing and lumen ensheathment or intracellular vacuole fusion [27].While the specifics of these mechanisms vary, they all establish the AB polarity of the tubes.If we turn on AB polarity in the present model, that is, we let λ 0 → 0.0 and λ 1 → 1.0, the solid structure tube of figure 2a opens up into a sheet-like structure as shown in figure 2b.This behaviour occurs because of the random initialization of the AB polarity (as illustrated in figure 2a).
To allow AB polarity to form properly we introduce the potential where f(r) ≃ e Àr 2 =2' 2 , such that the total potential is This potential aligns AB polarity against local areas of high density, in correspondence with experiments suggesting cell-cell contact directs AB polarity [28].It can also be thought of as alignment along a gradient field c, where c is a molecular, diffusing field of particles nucleated at cell locations, This formulation assumes the existence of such a molecular field.Although many molecular gradients are set up during blood vessel growth, such as vascular endothelial growth factor (VEGF), which elongates and reorganizes cells [6], it is unclear if these interact with and orient the polarity of cells.It is thus easier to think of the interaction as a direct cell-cell interaction as described by equation (2.8).
With γ > 0, lumen formations occur and the solid tube becomes hollow and fully enclosed, as shown in figure 2c.While network formation and lumen formation in reality are separate processes, in this study we will consider the simplified system of them occurring simultaneously.While γ > 0 ensures enclosed structures, PCP with λ 3 > 0 is needed to control the tube thickness.That is, λ 3 > 0 creates a preference for lengthwise alignment of particles, and thus establishes convergent extension, which in turn happens through cell intercalation events.Mathematically, the only difference between AB (p) and PCP (q) is the fact that λ 1 > λ 3 .The λ 2 term keeps p and q approximately orthogonal, and the magnitude of λ 3 thus determines how much AB alignment is favoured over PCP alignment, which in turn controls the thickness of the tubes since thicker tubes will have better aligned PCP.
In equation (2.2) we include vectorial interactions of AB polarity (S 1 ), but only nematic interactions of PCP (S 2 , S 3 ), since we do not want to impose a handedness to the vascular tubes.At branch points of the vascular network, there must be defects in PCP alignments, since, in a similar fashion to how you cannot perfectly comb the hair on a sphere, you cannot have a smooth surface vector field at a tube branch point.Taking the absolute value in equation (2.2) turns vectorial − 1 charge defects into two − 1/2 defects, which establishes symmetric branch points (see electronic supplementary material).
Finally, we note that we only model the endothelial cells themselves; in the jargon of active matter research, our model is 'dry'.In organoid experiments, there will naturally also be culture medium, extracellular matrix, pericytes, etc., present, and the system will perhaps be embedded in, for example, matrigel and collagen [29].These components mitigate their own interactions between one another and with the cells and could be explicitly modelled in a similar manner to our cell-cell interactions.Such interactions would complicate our model a lot and make interpretations harder, but one should keep in mind that the parameters we use effectively include these interactions and are not solely the result of pure cell-cell interactions.For instance, the effects of the viscosity of the culture medium would effectively introduce mobilities in equation (2.5).Likewise, the effects of shaking could be modelled effectively by including external noise in equation (2.5).We have tested such effects and our results remain qualitatively unchanged.

Vasculogenesis
During vasculogenesis blood vessels form from aggregating endothelial cells [2]. Figure 3 shows the self-assembly of three-dimensional vessels in our model.From an initial random distribution (figure 3a) the cells start aggregating (figure 3b) and form a tubular network (figure 3c).In figure 3, cells are initially sampled from a uniform distribution within a sphere but any distribution works.
Naturally, a major concern in vasculogenesis is to form a network of blood vessels that is fully connected.It has previously been shown how this percolation condition, i.e. whether all particles connect to one another, depends on the density of endothelial cells [7].In our model, cells have a preferred distance to one another and can attract over long distances.At first glance, therefore, it seems that initial density might not be an important quantity.However, particles only attract if their polarizations match, and as soon as vessel structures have formed, enclosed vessels will not attract one another, since two vessels nearing each other will have opposing AB polarity on their adjacent surfaces.
Because of this polarization, the vasculogenesis process in our model is also density dependent.
The density-dependent percolation behaviour is visualized in figure 4, which shows the probability for a particle to be part of the largest cluster P as a function of the initial density ρ.This is shown for various initial radii, or, in other words, for various numbers of particles ranging from ∼250 to ∼30 000.As is clear, the vessel network percolates at around ρ c ∼ 8.2 × 10 −3 , i.e. at an initial length scale of r À1=3 c 5-the same order of magnitude as the inter-particle spacing = 2. Figure 4 also shows some finite-size effects, since for a small number of particles, even those far below the transition point, the largest cluster, albeit small, will constitute a significant fraction of the whole system.

Growth and buckling
As organisms grow, their networks of blood vessels also need to grow.The vascular system needs to grow in two distinct ways: first, blood vessels need to increase their diameter in order to deliver increased amounts of blood.However, as vessels grow, their surface area to volume fraction decreases and so does their effectiveness.Thus, they also need to grow their network structure to maintain a space-filling network with small diameter vessels, capillaries, as the 'leaves' of the network [30].This latter version of growth is called angiogenesis and, as mentioned, is not the focus of our study.In this section, we introduce the growth of the blood vessel and consider the effect of PCP strength λ 3 , which creates a preference for growth in tube length rather than in tube diameter.The next section will demonstrate a less considered alternative to space-filling growth exploiting the buckling phenomenon demonstrated in this section.
First, we demonstrate that the growth of vessel diameters follows naturally when λ 3 = 0. Cell division is implemented as a Poisson process in the sense that each cell has a constant rate of division ν.When a cell divides, a new cell is created with the same polarities p and q, but placed at a random position next to its mother cell in the plane orthogonal to p, meaning that cells divide within the cell sheet.
Figure 5 shows the results of growth dynamics with λ 3 = 0. Figure 5a is the steady-state result of a self-assembly with λ 3 = 0.08.We then let λ 3 → 0 and ν → 3 × 10 −5 .Figure 5b  shows the result of the growth from 3500 particles to 25 000 particles.As is evident, the vascular network can grow uniformly under this model.However, the structure as a whole does not grow much, and, in fact, the density (vascular volume to free space) grows as well.This happens because we only model the cell division of the vascular network.In reality, the tissue between the vessels, which are cells we are not modelling, will also be dividing and in turn grow the structure as a whole.
If we instead consider the case of λ 3 > 0 with cell division, this leads to completely different growth.Since λ 3 induces a preferred diameter, this creates a tendency to grow more in length than thickness.Figure 6a-c shows how this leads to a growth-induced buckling instability [31], i.e. growth that does not retain the structure's shape.The buckling occurs because cell division is faster than the time to relax shape perturbations on long length scales.
A simple measure of buckling of a growing ring is to compare the curve length L with the structure's effective radius r.If no buckling occurs L/r = 2π.This definition works well as long as the structure remains a simple curve.Figure 6d shows this buckling as a function of time for various division rates ν.Clearly, the buckling occurs even for exceedingly small division rates at approximately the same value of νt, i.e. approximately the same number of cells.After the transition, the buckling degree grows exponentially in time, or, in other words, linearly with the number of cells.In this regime, this can be understood as the structure growing mostly in structure length L and not in effective size r.This continues until one side of the structure meets the other.
The inset of figure 6d shows the thickness of tubes calculated as N/L, where N is the number of cells in the structure.As is clear, a high growth rate leads to thicker tubes as the time scale for growth severely outpaces relaxation.The thickness of the tubes grows until buckling starts to occur, after which the thickness decreases as there is suddenly room to grow in length instead of thickness and the thickness stabilizes.The peak in thickness occurs at slightly later νt for larger ν.

Islets of Langerhans
In the pancreas, the so-called islets of Langerhans are responsible for the production of hormones such as insulin.While these islets constitute only about 1% of the pancreatic volume they contain about 10% of the blood vasculature.This dense vessel network is needed to provide energy to  the islets; it is also needed for intercellular communication [32] and to measure and regulate the production and injection of hormones into the bloodstream.The dense vasculature of pancreatic islets is shown in figure 7. The figure also shows the growth of these islets over 44 weeks.A crucial observation is that the islets' vascular density, despite the overall growth of the islets, remains almost constant [33].Furthermore, the thickness of the blood vessels also does not change much over this course (figure 7; variations in thickness are linked to age [34]).In other words, the growth of the vessel network is unlike that of figure 5. Angiogenesis is the canonical explanation for the growth of vascular network structure and this is indeed a factor for pancreatic islets [32].
What is also clear from figure 7 is the tortuous structure of the vasculature.This is in contrast to normal vasculature, which is regular and structured, as can be seen, for example, in the surrounding vessels in figure 7. Considering only angiogenesis as a growth mechanism, the tortuous structure and the space-filling growth are independent observations.Both observations can, however, be simultaneously explained by growth-induced buckling, as we demonstrate in figure 8. Figure 8b,c shows the different structures that can be attained in our model with different growth rates.As is clear, growing vessel structures with λ 3 > 0 leads to very tortuous structures.Furthermore, as shown in figure 8d, the division rate ν also determines the vascular density.With a large division rate, the structure does not have time to relax under the division-induced stress and the vascular density increases.Conversely, the vascular density is decreased for small division rates.Hence, a simple feedback mechanism where the proliferation rate inversely depends on the density can keep vascular density constant during growth.While this effect would definitely be co-occurring with angiogenesis, it is intriguing that it simultaneously gives an explanation for the tortuousness of the vascular network.Although the AB polarity-aligning parameter γ is only used for self-assembly, during this sort of growth γ > 0 can allow for anastomosis, the fusing of separate blood vessels.

Conclusion
We have demonstrated three-dimensional vasculogenesis in a simple model that aligns cell polarities through cell-cell interactions.The initial self-assembly of enclosed structures is ensured through the alignment of AB polarity against cell density.This is the key driver of lumen formation and enables the de novo formation of tubes.This interaction also allows for fusing of vessels (anastomosis).While PCP is not needed for the formation of enclosed structures, this polarity ensures thin tubular structures and convergent extension.
The self-assembly of the vascular networks results in fully connected vessels if the particles are initialized above a critical density, in accordance with previous two-dimensional experiments [7].Enclosed, non-connected structures appear at lower densities, and these do not interact because of their opposing AB polarities.
Introducing cell proliferation in our model leads to distinct behaviour depending on the strength λ 3 of PCP, which controls the preference of tube diameter.With λ 3 = 0, we have shown that uniform growth is possible.With λ 3 > 0 the tubes buckle under growth.While blood vessel buckling is typically associated with high blood pressure [35], this shows that such behaviour can also stem from cell proliferation.
Considering this buckling mode of growth, we compared it with the vasculature of islets of Langerhans, which show a large degree of tortuousness.The vessel density of these pancreatic islets furthermore remains constant during their growth.We suggest that a simple explanation for this behaviour is growthinduced buckling in which the cell division rate is coupled to the vascular density.This rate, in turn, may be controlled by negative feedback from the blood supply to the tissue.
Tortuous blood vessels are in general also found in cancerous tumours [36].Cancer growth is often associated with angiogenesis in nearby tissue, a process that we did not explore here, but could easily be introduced via local cell shape changes [22].The abnormality of tumour vessels is thought to be, in part, due to over-expression of VEGF-A [37].We have shown how tortuousness in blood vessels can be linked to the growth rate.This could thus also play a major role in the abnormal morphology of tumour vascularization.
Data accessibility.This article does not contain any additional data.Authors' contributions.J.B.K. conceptualized the study, designed the model, wrote the model code, ran the experiments and analysis and drafted the manuscript.B.F.N. participated in model discussions, was involved in writing the model code and revised the manuscript.A.T. provided biological insight and revised the manuscript.K.S. conceptualized the study, participated in model and analysis discussions and revised the manuscript.

Figure 6 .
Figure 6.Buckling of vessels during cell growth with λ 3 > 0. (a) Initial condition of a torus of 1000 cells.(b,c) Buckling during growth.(d ) Buckling measured as the vessel contour length L over the effective radius r as a function of time re-scaled by division rate ν, which is shown by colour; its value is given in the legend.νt ∼ 0.03 corresponding to ∼2100 cells.Inset shows tube thickness (N/L) during buckling.(Online version in colour.)

Figure 7 .Figure 8 .
Figure 7. Growth of the vasculature of islets of Langerhans over 44 weeks.Images show the vasculature at weeks 2 (a), 6 (b), 35 (c) and 44 (d).The islet is growing, but its vascular density remains approximately constant, as shown in (e).The y-axis shows the percentage of the vascular volume of the total islet volume.During these 44 weeks, the islet more than triples its volume.Data and images from Berclaz et al. [33].(Online version in colour.)