3D hybrid modeling of vascular network formation

We develop an agent-based model of vasculogenesis, the de novo formation of blood vessels. Endothelial cells in the vessel network are viewed as linearly elastic spheres and are of two types: vessel elements are contained within the network; tip cells are located at endpoints. Tip cells move in response to forces due to interactions with neighbouring vessel elements, the local tissue environment, chemotaxis and a persistence force modeling their tendency to continue moving in the same direction. Vessel elements experience similar forces but not chemotaxis. An angular persistence force representing local tissue interactions stabilises buckling instabilities due to proliferation. Vessel elements proliferate, at rates that depend on their degree of stretch: elongated elements proliferate more rapidly than compressed elements. Following division, new cells are more likely to form new sprouts if the parent vessel is highly compressed and to be incorporated into the parent vessel if it is stretched. Model simulations reproduce key features of vasculogenesis. Parameter sensitivity analyses reveal significant changes in network size and morphology on varying the chemotactic sensitivity of tip cells, and the sensitivities of the proliferation rate and sprouting probability to mechanical stretch. Varying chemotactic sensitivity also affects network directionality. Branching and network density are influenced by the sprouting probability. Glyphs depicting multiple network properties show how network quantities change over time and as model parameters vary. We also show how glyphs constructed from in vivo data could be used to discriminate between normal and tumour vasculature and, ultimately, for model validation. We conclude that our biomechanical hybrid model generates vascular networks similar to those generated from in vitro and in vivo experiments.


Introduction
In order to ensure their continued survival, many biological tissues are endowed with a network of blood vessels that delivers vital nutrients such as oxygen and glucose to cells, removes waste products and facilitates information exchange between different organs [1]. The vessel networks typically form via angiogenesis and/or vasculogenesis [2]. The process of angiogenesis has been widely studied due to its importance in wound healing and tumour growth, angiogenesis marking the transition from the relatively harmless and localised phase of avascular tumour growth to the potentially life-threatening phase of vascular growth [3]. During angiogenesis, new blood vessels emerge from pre-existing, perfused vessels, the endothelial cells that constitute the vessels being stimulated to proliferate and migrate chemotactically in response to growth factors produced by cells that lack an adequate supply of nutrients. In contrast, vasculogenesis is the de novo formation of new blood vessels from isolated endothelial cells. As such, it is not reliant upon the presence of a pre-existing vascular network and is a prominent feature of embryonic development. While in practice, both vasculogenesis and angiogenesis may simultaneously participate in vascular network formation during wound healing, tumour growth and embryonic development, their relative contributions remain keenly debated [4]. As a result, increased understanding of both processes and their interactions is urgently needed. Such understanding may also enable experimentalists and clinicians to establish how best to combine vascular targeting agents with other treatments either to stimulate healing of chronic wounds or to arrest the growth of solid tumours.
There is an extensive theoretical literature devoted to mathematical and computational modelling of vascular network formation. Models have been developed across the spectrum of physiological space and time scales, using a variety of frameworks. The most widely used methods are based on ordinary differential equations, partial differential equations and/or agent-based approaches. They differ in the geometrical resolution and detail they represent. Ordinary differential equation models describe the time evolution of global quantities, such as the number of vessels and the tumour volume [5]. Models formulated using partial differential equations typically describe the time evolution of spatially distributed, macroscale features, such as vessel volume fractions and the concentrations of oxygen and chemoattractants [6], although more recent, phase-field models provide a framework for simulating morphological features of vascular networks with continuous variables [7]. Agent-based approaches permit a more detailed study of the biological phenomena, on a scale at which the spatial and temporal evolution of individual blood vessels and/or endothelial cells may be resolved. In the paragraphs that follow, we illustrate briefly each of the aforementioned categories, focussing on agent-based approaches, which are the subject of the present work.
Hahnfeldt et al. [5] and Arakelyan et al. [8] proposed two-compartment ordinary differential equation models for the growth of a vascular tumour and its response to treatment with an anti-angiogenic chemical.
More generally, spatially-resolved, continuum models feature partial differential equations (PDEs) for the endothelial cell volume fraction that are coupled to PDEs for chemoattractants, chemorepellents and the extracellular matrix. The key advantages of such continuum models are the relatively small number of parameters that they contain and that they can be simulated efficiently. The drawbacks are that individual cell properties, geometrical details of the vascular morphology and subcellular networks cannot be included easily. Following a continuum approach, Balding and McElwain [9] developed a model including densities of blood vessels and capillary tips, that describes sprout formation and fusion from a pre-existing network in response to tumour angiogenic factors. This model applies the "snailtrail" concept under which moving capillary tips leave behind blood vessels. Flegg et al. developed a similar, continuum model with three-species (oxygen concentration, capillary tip density, blood vessel density) to study the efficacy of hyperbaric oxygen therapy for healing chronic wounds [10,11]. Several modelling papers have investigated the effects that deformation of the extracellular matrix have on network structure. For example, Edgar et al. [12] have shown how fibre orientation may guide the movement of tip cells. Other approaches based on mechano-chemical models [13], continuum [14,15], and continuum-discrete modelling [16], have tackled this issue. For example, Dyson et al. have used a multiphase approach to show how fibres embedded in the tissue matrix may bias cell movement and how cell movement may deform the fibres [17]. However, none of these models couples the mechanical 3D Vascular Network Formation stress that the endothelial cells experience to their growth, proliferation and the phenotype of their daughter cells. As an example of an agent-based off-lattice angiogenesis model we refer to Stokes and Lauffenberger [18] which includes sprout movement as a biased persistent random walk, based on a stochastic differential equation for the cell acceleration. See also Anderson and Chaplain [19] and Plank and Sleeman [20]. The persistence of motion is directly linked to chemotaxis. More information about these approaches can be found in the review articles of Mantzaris et al. [21], Ambrosi et al. [22], Merks and Koolwijk [23], and Scianna et al. [24].
We now highlight some existing agent-based models of vasculogenesis. Bentley et al. [25] resolve the shape and movement of vessel sprouts: each cell is decomposed into a number of connected agents and attention focusses on the influence of delta-notch-signalling on the initiation of vascular sprouts. Their work focusses on small segments of vascular networks. Other authors have simulated vasculogenesis and vascular network formation on larger length scales. For example, Merks et al. [26] have used the Cellular Potts Model (CPM) to investigate how contact inhibition influences the movement of capillary sprouts in response to an extracellular, diffusible chemoattractant produced by the endothelial cells. Cell movement, shape and alignment are determined by minimising an energy function that accounts for cell-cell binding, volume constraints, and chemotaxis, with chemotactic movement restricted to lattice sites that are adjacent to sites occupied by extracellular matrix (ECM). Scianna et al. [27] also use the CPM to develop a multiscale model that accounts for cell activation, migration, polarisation and adhesion, and in which each cell is decomposed into a nuclear and a cytosolic region. In the intracellular space the concentrations of arachidonic acid, nitric oxide and calcium are described by reaction-diffusion equations. The energy functional that determines cell morphology and movement accounts for the shape (area, perimeter), adhesion energy (intraellular, and nuclei and cytosol within the same cell), chemotaxis, contact inhibition, and persistence in cell movement. The extracellular concentration of vascular endothelial growth factor (VEGF) is described by a reaction diffusion equation. Building on the CPM, Szabo et al. [28][29][30] studied the roles of cell elongation and cell-matrix interactions on vascular network formation. In other work, Oers and Merks [31] coupled a CPM to a finite element model for substrate deformation. The energy function in the CPM depended on the strain and orientation-dependent stiffness of the extracellular matrix. Their model simulations yielded realistic patterns of network formation and sprouting from clusters of endothelial cells.
In contrast to the agent-based models mentioned above, where a single cell is represented by several agents, cellular automata represent one cell by one agent. For example, Stéphanou [32] developed an agent-based framework in which the movement of individual cells is influenced by chemotaxis and haptotaxis. In their model, interactions between vessels and the extracellular matrix play an important role in regulating cell movement. In addition to changes in the structure of the network, vessel radii also adapt in response to hydrodynamic, metabolic and angiogenic stimuli. Watson et al. adapted this approach to develop a hybrid two-dimensional model of retinal vascular plexus development [33]. The model accounts for individual astrocytes and endothelial tip cells that move on a lattice, and uses a continuum approach to model the distribution of biochemical signalling molecules and ECM components. Building on earlier work by Alarcòn et al. [34,35], Perfahl et al. [36], Macklin et al. [37], Shirinifard et al. [38], and Welther et al. [39,40] developed a similar cellular automaton model of vascular tumour growth.
A common feature of the cellular automaton model mentioned above is that cell-cell interactions are rule-based: mechanical forces are not included. Additionally, network nodes are restricted to lattice sites, creating angular, non-smooth networks. Off-lattice models can be used to address these limitations and also generalise naturally to three dimensions. They have been widely used to investigate cell-cell interactions in multicellular systems, including the liver [41], development [42], avascular tumour spheroids [43] and homeostasis in the intestinal crypt [44]. A mechanically-based, hybrid model for network formation was proposed by Jackson and Zheng [45]. In their approach, tip and stalk cells are modelled as interconnected, elastic agents. Chemotactic movement of the tip cells produces mechanical forces that act on the neighbouring stalk cells whose rates of proliferation depend on their mass and maturity. An alternative off-lattice agent-based model of angiogenesis and vascular tumour growth was presented in Drasdo et al. [46].
Similay to Drasdo et al. [46], and in contrast to the two-dimensional, lattice-based models described above, in this paper we develop a three-dimensional (3D), off-lattice, hybrid model of vasculogenesis, although the cells (or agents) are assumed to have fixed shape, for simplicity and to improve the computational efficiency. As such, blood flow is not considered as part of our model. Further, we do not describe the individual endothelial cells that line a small diameter blood vessel; we use spherical elements to represent small vessel sections. Our model is motivated by in vitro experiments in which isolated green fluorescent rat microvessel segments were embedded within matrigel that sits on top of a droplet containing immortalised red fluorescent MDA-231 breast cancer cells so that the tumour and endothelial cells were not in direct contact. In vitro images, acquired on days 1, 4 and 6 using a Zeiss Z1 Ob-server microscope at 5x magnification, reveal endothelial cell proliferation and migration in response to tumour-derived growth factors, and the formation of small unperfused networks (see Figure 1). Further motivation for our hybrid model comes from experiments designed to understand how endothelial cells respond to mechanical stretch [47,48]. In [47], Liu and colleagues show that endothelial cells increase their rate of proliferation under stretch and that both cell-cell adhesion and engagement of vascular endothelial cadherin are needed to transduce the mechanical stretch into proliferative signals. In separate work [48], Zheng et al. showed that endothelial cells response to static stretch by increasing their rates of cell proliferation, lumen formation and branching, and that VEGF binding is necessary to mediate these responses. Figure 1. Series of images showing how vascular networks emerge from rat microvessel segments in vitro. Isolated green fluorescent rat microvessel segments were embedded within matrigel that sits on top of a droplet containing red fluorescent breast cancer cells. The endothelial cells proliferate and migrate in response to tumour-derived growth factors, forming small networks. The images were collected on days 1, 4 and 6 (left-to-right). The green and red scale bars correspond to length scales of 500µm.
Our agent-based, off-lattice model represents the vasculature as a network of spheres whose centres are connected by linear springs and whose movement in 3D is determined by applying local force balance. We distinguish two types of endothelial cell agents: vessel elements which are contained within the network and proliferate, and tip cells which are located at the ends of vessels, do not proliferate and respond via chemotaxis to spatial gradients in angiogenic factors such as VEGF. For simplicity, we assume that tip cells and vessel elements neither change phenotype nor swap position, although there is experimental evidence for these phenomena [25]. In contrast to existing hybrid models and motivated by work of Liu et al. [47] and Zheng et al. [48], here cell proliferation and branching are assumed to depend on the degree of mechanical stretch (or compression) experienced by individual vessel elements. Thus, two key processes drive endothelial cell movement: chemotactic movement of tip cells creates a pull which acts on vessel elements contained within the network whereas cell proliferation creates a mitotic push which acts on the tip cells.
Numerical simulations reveal that our hybrid, mechano-chemical model can reproduce the key features of vasculogenesis. Extensive computations are performed in order to show how parameter changes influence network development. Since model simulations are stochastic, exact network comparisons are not possible. Instead, our parameter sensitivity analyses are based on network features extracted from multiple simulations generated with different random seeds. Metrics used to characterise the networks include the following: the total network length, the average number of branch points per unit length, the tortuosity and the distribution of vessel segment lengths. We note that these metrics can also be extracted from in vitro and in vivo experimental data and could, thereby, facilitate comparisons between our model and available data. Since variation of model parameters can affect multiple metrics, we also introduce glyphs simultaneously to visualise several network metrics. A glyph is a graphical object whose attributes are bound to data [49]. The two-dimensional glyphs that we develop enable us to clearly present multidimensional data or metrics in a single graphical entity.
The main achievements of the paper can be summarised as follows: • The development of an off-lattice hybrid model of vasculogenesis in which mechanical stretch regulates endothelial cell proliferation and capillary sprout formation; • The identification of quantitative metrics that can be used robustly to characterise and compare vessel networks and to study the impact of external perturbations on network structure; • The design and use of glyphs as an objective way of aggregating multiple network features in one diagram.
• A demonstration that mechanical stimuli alone can generate networks whose morphological features are qualitatively similar to those observed in vitro and in vivo; The remainder of the paper is structured as follows. The mathematical model is introduced in Section 2. Simulation results are presented in Section 3. The paper concludes in Section 4 where we discuss our findings and outline directions for future work.

Methods
In this section we introduce the 3D agent-based model that we have developed, using an off-lattice approach, to simulate vasculogenesis. Our goal is to establish whether physically realistic vessel networks can be generated when endothelial cell proliferation and capillary sprout formation are regulated by mechanical effects. A detailed description of the computational model is included below while information about the algorithm and parameter values used to generate numerical simulations is provided in the Supplementary Materials. In order to characterise objectively morphological changes that occur as the networks evolve and/or system parameters are varied, several metrics are introduced (e.g. total vessel length, tortuosity and number of branch points per unit vessel length) and applied to the networks.
Glyphs that simultaneously depict multiple metrics are also introduced and used to present, in a concise manner, network attributes.

The Computational Model
Our three-dimensional, off-lattice, agent-based model aims to describe the de novo formation of vascular networks that occurs during vasculogenesis. We distinguish two types of endothelial cells: tip cells and vessel elements (see Figure 2). Tip cells are located at the blunt end of each capillary sprout and all other (internal) segments are vessel elements. While recent experimental results have shown that tip cells may change position (and phenotype) with endothelial cells located in the same sprout [50], here, for simplicity, we neglect such effects and assume that the identity of the leading tip cell in a particular sprout is fixed. Tip cells and vessel elements are modelled as linearly elastic spheres and all pairs of connected cells or elements exert equal and opposite (mechanical) forces on each other.
In our model, tip cells and vessel elements differ in two important ways (see Figure 2): vessel elements can proliferate whereas tip cells cannot; tip cells are subject to a chemotactic force caused by spatial  gradients in local levels of the diffusible angiogenic factor, VEGF, whereas vessel elements are not sensitive to VEGF (cell-cell contacts have been shown to inhibit VEGF-induced signalling within a vessel [26]).
Thus, as shown in Figure 2(A), tip cells perform a persistent random walk, biased by chemotaxis in response to spatial gradients in VEGF and constrained by mechanical forces due to (linearly elastic) cell-cell interactions and drag forces due to interactions with the local tissue matrix. Vessel elements are subject to drag, mechanical and random forces and an angular persistence force, the latter two forces accounting for cell interactions with the local environment. We introduce below the functional forms used to model each force but first explain how we derive the equations of motion for each vessel element.
Suppose that at time t, the network comprises N = N (t) elements (tip cells and vessel elements) that are located within a three-dimensional, Cartesian domain of size the position of the centre of vessel segment i (i = 1, . . . , N ), and record in an adjacency matrix E the node numbers of all pairs of connected vessel elements. We use Newton's second law to derive the equations of motion. In the over-damped limit, we neglect inertial effects and obtain the following force balances Table 1. Model assumptions. Summary of the key differences between tip cells and vessel elements. Key: + and 0 indicate whether a cell type exhibits a particular property (+) or not (0).

Property Tip cell Vessel element
Chemotaxis + 0 Angular persistence 0 + Directional persistence of movement + 0 Stress-dependent proliferation 0 + for tip cell i and vessel segment j respectively: In equations (1) and (2), we assume that the drag force on vessel element i is proportional to its velocity dx i /dt, the positive constant µ denoting the drag coefficient. We denote the mechanical force by F m i , the random force by F r i , the chemotactic force by F c i , and the directional and angular persistence forces by F p i and F a i respectively. These forces are prescribed as follows.
• Mechanical force, F m i (tip cells and vessel elements) The mechanical force acting on a tip cell or vessel element is the net force exerted on it by its neighbours. We assume that cells/elements i and j only interact if the distance between their centres is less than a fixed value, l c . In more detail, and following the interacting sphere approach outlined in [51][52][53][54], if |x i − x j | < l c , then the interaction force between cells/elements i and j is parallel to the vector x i −x j connecting their centres and its magnitude S is piecewise linear, taking values that range between S = S c for highly compressed cells/elements, and S = 100S c for highly stretched cells/elements. If we denote by N i those cells/elements j = i for which (i, j) ∈ E and |x i − x j | < l c , then the net mechanical force acting on cell/element i is given by: In equation (3), R c is the target cell/element radius and the spring function S(x) is a piecewise linear function so that the interaction force depends on the distance between the cells/elements and is much stronger, and attractive, for pairs of stretched cells/elements (|x i − x j | > 2R c ) than it is for pairs of compressed cells/elements. We note that l c > 2R c holds for the cut-off distance l c .

• Random force, F r i (tip cells and vessel elements)
We assume that all vessel elements experience random forces due to heterogeneity in the surrounding tissue environment and that the random force acting on vessel element i is given by: where ξ i is a unit vector (|ξ i | = 1), pointing in a direction chosen randomly from a uniform distribution, and the constant σ denotes the vessel element's sensitivity to random fluctuations.
• Chemotactic force, F c i (tip cells) Following [26,27], we suppose that cell-cell contact inhibits the chemotactic sensitivity of vessel elements so that only tip cells are chemotactic. We assume further that the chemotactic force acting on tip cell i is proportional to the gradient of the VEGF concentration c = c(t, x i ) at its centre so that where the positive constant χ denotes the chemotactic sensitivity. For simplicity, in the simulations that follow we prescribe a fixed chemoattractant gradient, with ∇c = (−c x , 0, 0) so that the constant c x denotes the VEGF gradient.
• Directional persistence force, F p i (tip cells) We suppose that tip cells have a tendency to continue moving in the same direction. Rather than retaining inertial effects in Equations (1) and (2), we model this tendency by introducing the directional persistence force. Unlike inertia, we assume that this force is induced by active cell movement (along the ECM, for example) and that it acts in the direction of cell movement over a timescale of length τ > 0. Accordingly, we denote the directional persistence force, F c i acting on tip cell i as follows: where ω p > 0 denotes the persistence coefficient. We remark that a similar approach was adopted in [18] where persistence was directly linked with chemotaxis.
• Angular persistence force, F a i (vessel elements) The angular persistence force accounts for forces that vessel elements experience due to interactions with their microenvironment. We assume that it acts to stabilise buckling instabilities induced by proliferation (see Figure 3). If vessel element i is not a branching point, and has two neighbours j and k (so that (i, j), (i, k) ∈ E), then the angular persistence force acting on vessel element i is given by In equation (7), the positive constant ω a denotes the angular spring constant, and 0 ≤ α angular ≤ π is the (smallest) angle between the vectors (x j − x i ) and (x k − x i ). If vessel element i is a branch point, connected to more then two other elements, then α angular is taken to be the smallest branching angle between the vessel elements.
We complete our model of vasculogenesis by explaining how division and sprouting are represented.
We assume that both processes are regulated by the amount of stretch that a vessel element experiences.
Thus, for example, an elongated vessel element increases its rate of progress through the cell cycle while a compressed element decreases its rate. For simplicity, and extending an approach used in Owen et al. [55], the progress of vessel segment i through the cell cycle is monitored by a phase variable φ i (t) ∈ [0, 1] whose evolution is determined by the following differential equation: where the positive constants k φ and β φ denote, respectively, the maximum rate of progress through the cell cycle and the extent to which progress is modified by mechanical compression and/or tension, while l i approximates the length of vessel element i: if vessel element i has only two neighbours, elements j Figure 3. Schematic of the angular persistence force. Suppose vessel element i is connected to vessel elements j and k and that their centres have coordinates x i , x j , and x k , respectively. Then the angle α angular is defined to be the smallest angle between the two vessel segments, as indicated. This angle and the coordinates of the three vessel elements are used in Equation (7) to determine the angular persistence force, F a . and k, say, then If element i is a branching point, with more than two neighbours, then l i is the average over all distances between i and its neighbours. Vessel element i divides if φ i = 1 and then the phases of the parent and new daughter cells are set to zero.
The daughter cell is located randomly within the ball of radius R c , centred on its parent and its fate is determined by its position and the degree of mechanical compression being experienced by its parent (see Figure 4 in main text). In more detail, we denote the position of the daughter cell relative to its parent in polar coordinates by x random = (r random , α random , γ random ) where r random ∈ (0, R c ), α random ∈ [0, π] is the smallest angle between the vessel and the random vector (Figure 2), and γ random ∈ [0, π] is the angle between the plane spanned by the vessel elements connecting cells i, j and k and the random vector, x random . Then the new cell is a tip cell which forms a new capillary sprout if where l i is defined by Equation (9) and (x) + = max (0, x). The angles α random and α that appear in Equation

Analysis and Visualisation of Simulation Results
Since our model is stochastic, for each choice of parameter values, multiple simulations should be performed, using different random seeds, and suitable statistics extracted from the simulations to determine how the network structure evolves over time and how it is affected by changes in the parameter values.
In order to characterise the vessel networks, and to facilitate future comparisons with experimental data, the following metrics are calculated for each simulation:  x random = (r random , α random , γ random ) within the three-dimensional ball of radius R c centred on its parent cell. The fate of the daughter cell depends on the degree of stretch experienced by the parent cell, more compressed cells having greater probability of producing tip cells than mechanically stretched ones. The volume of the cells shaded red and green, respectively, indicate the probability that, when a cell divides, its offspring will be incorporated within the same vessel (red region) or form a new capillary sprout (green region), larger red regions corresponding to cells that are more stretched. evolve over time (or as particular model parameters vary). We remark that a variety of metrics could be used to characterise the synthetic networks that our model generates. For example, in a series of papers Konerding and coworkers [56][57][58] have measured inter-vessel distance, inter-branch distance, mean vessel diameter, vessel diameter and branching angle in three-dimensional corrosion casts of tumour networks.
The five metrics we use are chosen for simplicity: we postpone consideration of alternatives for future work.
We use two complementary approaches to visualise the results of our statistical analyses. When focussing on a single feature (e.g. total network length), we include plots showing how the mean and variance of that metric change over time (see Figure 6) or as system parameters vary (see Figure 7 and 8) [49]. Since variation of model parameters may affect multiple metrics, we have designed a glyph so that we may visualise simultaneously, in a simple and comprehensive manner, several network metrics ( Figure 5). In order to facilitate comparison of glyphs, they contain information on mean values only (where appropriate, information about standard deviations is presented elsewhere).  this work. We perform sensitivity analysis to determine the influence of individual parameters on the network morphology. The impact of chemotaxis is investigated by varying the assumed constant chemotactic gradient, c x (for simplicity, we assume that the chemoattractant varies only in one direction, and that its gradient in this direction is constant; see Supplementary Material, Equation (5)). The way in which mechanical stimuli may regulate endothelial cell proliferation and, thereby, vessel morphology, is studied by varying the mechanical sensitivity parameter β φ that relates the rate of cell cycle progress to the degree of stretch that a cell is experiencing, larger values of β φ indicating greater mechano-sensitivity (see Supplementary Material, Equation (8)). In what follows, unless otherwise stated, all simulations are generated using the default parameter values detailed in Table S.1 and the initial network comprises two, straight and parallel vessels, each with three elements (tip-vessel-tip).
Before presenting our parameter sensitivity analysis, we pause to consider a typical dynamic simulation in which chemotaxis is neglected. In Figure 6 we present results generated from a single realisation  Having introduced typical network simulations and shown how glyphs can be used to characterise network features, we proceed to investigate how network morphology changes when model parameters vary. In the absence of detailed data with which to estimate system parameters, parameter ranges were chosen individually in order to obtain results that appear physically realistic (a comprehensive parameter sensitivity analysis is beyond the scope of the current work and, thus, we cannot guarantee that our choice of parameter values is the only feasible combination: we postpone such considerations to future work). For example, large increases in the chemotactic sensitivity (χ in Equation (5)), or the sensitivity to random fluctuations (σ in Equation (4) The results presented in Figures 7, 8 and 9 illustrate the variability in size and morphology of the simulated networks for given sets of parameter values and as particular parameter values vary. In Figure 7 we vary the chemotactic sensitivity, χ, and the mechanical sensitivity, β φ , while in Figure 8 we vary the sensitivity to random fluctuations σ, and the sprouting probability k spr . For each set of parameter values, we present the structure of 12 networks generated from different random seeds at t = 900 timesteps (ts) to illustrate the variability in the networks that can be generated from the same parameter values.
As expected, the networks increase their directionality as the chemotactic sensitivity, χ, increases, while varying the mechanical sensitivity, β φ , has a significant influence on the overall length of the networks and the number of branches per unit length of network. Equally, increasing σ, the parameter measuring the cells' sensitivity to random fluctuations, produces networks that are more tortuous while varying the branching probability, k spr , influences both the number of branches per unit length of the network and its total length. (These results are reinforced and made more precise in the glyphs presented in Figures 9 and S.3.) When we vary both χ and β φ we find that higher mechanical sensitivities lead to longer and more branched networks that cover a larger area while increasing the chemotactic strength decreases the area over which the network spreads, and increases the number of branches per unit length.
The chemotactic strength strongly influences the centre of mass displacement (results not shown). The above parameter sensitivity analyses enable us to determine the influence of particular parameters on the structure of the emerging vessel network (see Table 2). Not surprisingly, changes in the chemotactic sensitivity are predicted to shift the centre of mass of the vessel network towards the source of a chemoattractant. The mechanical sensitivity influences the temporal dynamics of network formation and the sprouting probability has a strong effect on the area spanned by the network area, the total network length and the number of branches per unit vessel length. Additionally, these results yield testable  The influence of varying the chemotactic sensitivity, χ, on the network morphology at t = 900 ts is depicted in Figure 10 we choose this timepoint to ensure that any transients associated with the initial conditions have disappeared and that metrics associated with the networks' internal structure (e.g. the average number of branch points per unit length or tortuosity) have stabilised at constant values. While the distribution of segment lengths, the total network length and the number of branches per unit length do not change significantly, as expected, the centre of mass shifts significantly, in the direction of the chemoattractant. Figure 11 shows, in more detail than either Table 2 or Figure 9, that increasing the mechanical sensitivity, β φ , has a negligible effect on most network metrics, except for the total network length which is significantly larger. In addition to random fluctuations, drag and spring forces, vessel cells are also subject to an angular persistence, which accounts for the characteristic persistence time during which cells maintain memory of their former direction and therefore no significant changes of direction are observed [59]. In this case, persistence is a consequence of the interaction with the surrounding   Figure 11. Statistical analysis showing how network metrics depend on the mechanical sensitivity, β φ . As β φ varies, the total network length at t = 900 ts changes significantly; changes in the other metrics are less pronounced, although the tortuosity increases slightly as β φ increases.  (10) and (11)) has a similar effect on network morphology (see Figure S.7). Increasing k spr shifts the distribution of segment lengths, leading, as we would expect, to more shorter vessel segments. The number of branches per unit length increases slightly and has a marked effect on the overall length of the network. In the absence of chemotaxis, the networks are isotropic and the tortuosity does not change as k spr varies.

Discussion
Understanding the mechanisms by which new blood vessels form has long fascinated and challenged experimental researchers, their interest stimulated in large part by the important roles played by angiogenesis and vasculogenesis during development, wound healing and their contributory roles in diseases such as rheumatoid arthritis, retinopathy of prematurity and cancer. Improvements in imaging techniques mean that it is now possible to collect high-resolution data showing how the spatial extent and morphology of vessel networks change over time and how the behaviour of the endothelial cells that form the networks is influenced by biophysical stimuli.
In general terms, mathematical and computational modelling represent a natural framework for integrating diverse sources of experimental information about a biological system, for investigating the ways in which different processes may interact and generate observed behaviours, for testing hypotheses and generating new predictions about the system's response to perturbations. In this paper we have developed a new, off-lattice model of vascular network formation. Where most existing hybrid models assume that chemotaxis is the dominant process driving endothelial cell proliferation and capillary sprout formation, our hybrid model allows us to investigate the possible roles played by mechanical stimuli. In particular, motivated by in vitro experiments performed by Liu et al. [47] and by Zheng et al. [48], we assume that elongated cells increase their proliferation rate and that compressed cells are more likely to form new capillary sprouts when they divide. Using these rules, our hybrid models generates vascular networks whose morphologies are qualitatively similar to those of real vascular networks cultivated in vitro (compare Figures 1 and 6).
We have performed extensive 3D numerical simulations to establish how different model parameters influence the size and morphology of the vascular networks. Since model simulations are stochastic, our parameter sensitivity analyses were based on network features extracted from multiple simulations that were generated using different random seeds. The quantitative measures used to characterise the networks included the distribution of vessel lengths, the total network length, the average number of branches per unit length, tortuosity and the displacement of the network's centre of mass from its initial position. For example, increasing the chemotactic sensitivity of the tip cells does not appear significantly to affect the distribution of vessel lengths, the total network length or the average number of branches per unit length (see Figures 9(A) and 10), although it shifts the centre of mass in the direction of the chemoattractant.
Equally, increasing the sensitivity of proliferating cells to mechanical effects leads to markedly larger networks (see Figures 9(B) and 11). Increasing the sprouting probability has a significant effect on the area spanned by the network, total network length and the average number of branches per unit length An alternative to the use of angular persistence to model the effect of the extracellular medium on network structure is explicitly to account for such interactions [60][61][62][63]. In Figure 12 we present simulation results obtained by embedding the developing vessel network in a 3D gel which comprises non-proliferating, linearly elastic particles that possess the same mechanical properties as the endothelial cells. Force balances are applied to the vessel and gel elements, as before, except that forces due to cell-gel interactions supercede the angular persistence forces. The simulation results presented in Figure 12 reveal that this model extension yields realistic network structures and could be used to investigate how the mechanical properties of the gel or tissue in which the networks develop affect their size and morphology [5].
We have shown that mechano-transduction, whose effects on solid tumour growth have been studied by Byrne and Preziosi [64] and Chaplain et al. [65], is a key factor in our model, preventing compressed vessels from proliferating and stimulating sprout formation. Another important feature of our model is the active movement of tip cells that stabilises the network and establishes network morphology. Without tip cells, and/or a balance between tip-cell movement, vessel cell proliferation and cell-cell adhesion (as modelled here via a spring constant), our simulated endothelial cells do not form a spatially distributed structure. Indeed, if the chemotactic sensitivity of the tip cells is too high, then they break away from the vessel cells which then form dense aggregates. Experimental studies have shown that tip-and vessel cells can change their phenotype [66], that such changes may be reversible, and may be regulated by subcellular signalling pathways [67]. In future work, by building upon previous work on multiscale modelling [34,35,63,68,69], we aim to increase the biological realism of our model by embedding, within the vessel cells, mathematical models of subcellular signalling pathways (e.g. delta-notch signalling [25], vascular endothelial growth factor signalling [48,70] and Rac1 [47]) and using them to determine whether a particular cell adopts a tip-or vessel phenotype and/or whether a new daughter cell forms a new capillary sprout.
In this paper, attention has focussed on characterising geometric features of the simulated networks: functional properties, such as average blood flow, wall shear stress and oxygen delivery rates, and their influence on network evolution have been neglected. Equally, the metabolic demands of the cells contained within the tissue, and their cross-talk with the evolving vessel networks, have been ignored. Thus, in practice, our model focuses on the early stages of network development (e.g during vasculogenesis), before flow has been established, so that only stresses exerted on the vessels by the surrounding microenvironment need be considered. The impact of wall shear-stress and the metabolic demands of the perfused tissue can be incorporated in future model extensions that account for intravascular flow [71,72], oxygen delivery and consumption [32,[34][35][36]. We can then extract from these simulations similar statistics and construct glyphs which are similar to those generated in this paper to investigate how inclusion of these processes affects the morphology of the resulting networks. By further extending the model to account for the action of anti-angiogenic compounds such as endostatin or combretastin we could also identify how such drugs should act in order to transiently normalise the tumour vasculature and, thereby, increase the tumour's sensitivity to radiotherapy and standard chemotherapeutic agents [73,74] In our analysis of vascular networks, glyphs have served as objective tools with which to compare averaged properties of synthetic networks. They reveal how particular physical quantities change as the networks evolve over time, and show how model parameters and combinations thereof influence the network morphologies. The advantage of these glyphs lies in their high information content: they combine several network characteristics in a single image in an objective and quantitative manner. In addition to characterising synthetically generated data sets, the glyphs could also be used to characterise in vivo images of vascular networks from animals and patients and, in so doing, provide a first step towards validation and parameterisation of angiogenesis and/or vasculogenesis. As an illustrative example, in Figure 13 we present a high-resolution image from [75] which shows two different tissue areas, one associated with reference, healthy tissue and the other associated with a tumour. Vascular graphs extracted from the two regions of interest are used to calculate the total length of the networks and the average number of branches per unit length. The corresponding glyphs are normalised so that the properties associated with the reference tissue take half of their maximum values (and so that undervascularised necrotic areas and denser vascularised tissues can be easily identified). Our analysis reveals that, compared to the reference tissue, the tumour vasculature is longer, with thicker vessels that are more densely distributed across the tissue segment, and with more branch points per unit length of vessel. Figure 13. Generating glyphs from experimental data. These glyphs were generated using imaging data (A) from cancerous and reference tissue areas [75]. The vascular graphs (B) were extracted and used to construct glyphs. In contrast to the glyphs introduced before, the shading of the central circle indicates the vascular volume fraction in the region of interest. All glyph metrics were normalised so that the values for the reference (healthy) tissue are 0.5.
strates that mechanical regulation of endothelial cell proliferation and capillary sprout formation can generate realistic 3D vascular networks, even when chemotaxis is not active. Additionally, by introducing a range of statistical measures to characterise the networks and designing glyphs with which to visualise multiple network metrics, we are able to predict how network properties evolve over time and as key system parameters vary. Additionally, comparison of glyphs generated from synthetic data with those generated from in vitro and in vivo data reveal how, in the longer term and with suitable dynamic data, it may be possible to validate and parameterise our model and other hybrid models of vasculogenesis, angiogenesis and vascular remodelling.

Supplementary Material
The Supplementary Material comprises three sections in which we describe the computational algorithm and parameter values that were used to generate the numerical results (S.1 and S.2, respectively), before presenting additional simulation results (S.3).

S.1 The Computational Algorithm
We now outline the algorithm that is used to generate numerical simulations of our hybrid, agent-based model (see, also, Figure S.1). We consider a three-dimensional Cartesian domain, of size W X × W Y × W Z , discretised uniformly, with spacings ∆X, ∆Y and ∆Z in the x, y and z directions, respectively. While an off-lattice approach is used to model the evolution of the vascular network, the domain discretisation is used to sort the cells (to increase computational efficacy in finding interacting neighbours [76]) and, thereby, to increase the speed of the numerical simulations. As stated above, the total number of segments (both tip and vessel segments) at time t is denoted by N = N (t). The vascular network structure is stored in an undirected graph in which cell centre coordinates N and cell-cell connections E are recorded.
The adjacency matrix E contains 2-tuples of the node numbers of all connections in the vessel network.

Initialise the simulation (t = 0)
The initial network is defined by specifying a fixed number of vessels, the number of segments within each vessel, their positions and connectivity. Unless otherwise stated, the initial network comprises two parallel vessels, each formed of three elements (tip-vessel-tip). The initial cell cycle phase of each vessel element and the chemoattractant gradient are also prescribed (∇c = (−c x , 0, 0)).    (1) and (2)). We assume that these forces act on two different time-scales: a slow time-scale is associated with chemotactic and random movement, and a faster one with mechanically induced cell movement. We exploit this separation of time-scales by using a two-step method to integrate the equations of motion and update x i , the position of element i.
• Step 1. Calculate the random force F r i using Equation (4). If cell i is a tip-cell, then calculate the chemotactic and persistence forces F c i and F p i using Equations (5) and (6).
Update x i , the position of cell i, using an explicit Eulerian method to integrate Equation (1) for tip cells or (2) for vessel cells, and neglecting the mechanical and angular persistence forces.
• Step 2. Calculate the mechanical force F m i using Equation (3) and if cell i is a vessel cell then calculate the angular persistence force F a i using Equation (7 (3)).
3. If t ≥ T then end the simulation; otherwise, return to step 2.
A flowchart of the computational algorithm is given in Figure S

S.2 The Default Parameter Values
In the