3D hybrid modeling of vascular network formation
Abstract
We develop an offlattice, agentbased model to describe vasculogenesis, the de novo formation of blood vessels from endothelial progenitor cells during development. The endothelial cells that comprise our vessel network are viewed as linearly elastic spheres that move in response to the forces they experience. We distinguish two types of endothelial cells: vessel elements are contained within the network and tip cells are located at the ends of vessels. Tip cells move in response to mechanical forces caused by interactions with neighbouring vessel elements and the local tissue environment, chemotactic forces and a persistence force which accounts for their tendency to continue moving in the same direction. Vessel elements are subject to similar mechanical forces but are insensitive to chemotaxis. An angular persistence force representing interactions with the local tissue is introduced to stabilise buckling instabilities caused by cell proliferation. Only vessel elements proliferate, at rates which depend on their degree of stretch: elongated elements have increased rates of proliferation, and compressed elements have reduced rates. Following division, the fate of the new cell depends on the local mechanical environment: the probability of forming a new sprout is increased if the parent vessel is highly compressed and the probability of being incorporated into the parent vessel increased if the parent is stretched.
Simulation results reveal that our hybrid model can reproduce the key qualitative features of vasculogenesis. Extensive parameter sensitivity analyses show that significant changes in network size and morphology are induced by varying the chemotactic sensitivity of tip cells, and the sensitivities of the proliferation rate and the sprouting probability to mechanical stretch. Varying the chemotactic sensitivity directly influences the directionality of the networks. The degree of branching and thereby the density of the networks is influenced by the sprouting probability. Glyphs that simultaneously depict several network properties are introduced to show how these and other network quantities change over time and also as model parameters vary. We also show how equivalent glyphs constructed from in vivo data could be used to discriminate between normal and tumour vasculature and, in the longer term, for model validation. We conclude that our biomechanical hybrid model can generate vascular networks that are qualitatively similar to those generated from in vitro and in vivo experiments.
D Hybrid Modelling of Vascular Network Formation
Holger Perfahl, Barry D. Hughes, Tomás Alarcón, Philip K. Maini, Mark C. Lloyd, Matthias Reuss, Helen M. Byrne
1 Center Systems Biology, University of Stuttgart, Stuttgart, Germany
2 School of Mathematics and Statistics, University of Melbourne, Australia
3 ICREA, Pg. Lluís Companys 23, 08010 Barcelona, Spain
4 Centre de Recerca Matemàtica, Campus de Bellaterra, Barcelona, Spain
5 Departament de Matemàtiques, Universitat Autònoma de Barcelona, Barcelona, Spain
6 Wolfson Centre for Mathematical Biology, Mathematical Institute, University of Oxford, Oxford, UK
7 H. Lee Moffitt Cancer Center and Research Institute, Tampa, USA
Email: holger.perfahl@ibvt.unistuttgart.de
1 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 lifethreatening phase of vascular growth [3]. During angiogenesis, new blood vessels emerge from preexisting, 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 preexisting 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 agentbased 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, phasefield models provide a framework for simulating morphological features of vascular networks with continuous variables [7]. Agentbased 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 agentbased approaches, which are the subject of the present work.
Hahnfeldt et al. [5] and Arakelyan et al. [8] proposed twocompartment ordinary differential equation models for the growth of a vascular tumour and its response to treatment with an antiangiogenic chemical. More generally, spatiallyresolved, 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 preexisting 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 threespecies (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 mechanochemical models [13], continuum [14, 15], and continuumdiscrete 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 stress that the endothelial cells experience to their growth, proliferation and the phenotype of their daughter cells. As an example of an agentbased offlattice 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 agentbased 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 deltanotchsignalling 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 cellcell 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 reactiondiffusion 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 cellmatrix 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 orientationdependent 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 agentbased 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 agentbased 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 twodimensional 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 cellcell interactions are rulebased: mechanical forces are not included. Additionally, network nodes are restricted to lattice sites, creating angular, nonsmooth networks. Offlattice models can be used to address these limitations and also generalise naturally to three dimensions. They have been widely used to investigate cellcell interactions in multicellular systems, including the liver [41], development [42], avascular tumour spheroids [43] and homeostasis in the intestinal crypt [44]. A mechanicallybased, 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 offlattice agentbased 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 twodimensional, latticebased models described above, in this paper we develop a threedimensional (3D), offlattice, 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 MDA231 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 Observer microscope at 5x magnification, reveal endothelial cell proliferation and migration in response to tumourderived 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 cellcell 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.
Our agentbased, offlattice 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, mechanochemical 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 twodimensional 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 offlattice 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;
2 Methods
In this section we introduce the 3D agentbased model that we have developed, using an offlattice 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.
2.1 The Computational Model
Our threedimensional, offlattice, agentbased 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 (cellcell contacts have been shown to inhibit VEGFinduced 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) cellcell 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.
Property  Tip cell  Vessel element 
Chemotaxis  +  0 
Angular persistence  0  + 
Directional persistence of movement  +  0 
Stressdependent proliferation  0  + 
Suppose that at time , the network comprises elements (tip cells and vessel elements) that are located within a threedimensional, Cartesian domain of size . We denote by the position of the centre of vessel segment , and record in an adjacency matrix the node numbers of all pairs of connected vessel elements. We use Newton’s second law to derive the equations of motion. In the overdamped limit, we neglect inertial effects and obtain the following force balances for tip cell and vessel segment respectively:
(1) 
(2) 
In equations (1) and (2), we assume that the drag force on vessel element is proportional to its velocity , the positive constant denoting the drag coefficient. We denote the mechanical force by , the random force by , the chemotactic force by , and the directional and angular persistence forces by and respectively. These forces are prescribed as follows.

Mechanical force, (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 and only interact if the distance between their centres is less than a fixed value, . In more detail, and following the interacting sphere approach outlined in [51, 52, 53, 54], if , then the interaction force between cells/elements and is parallel to the vector connecting their centres and its magnitude is piecewise linear, taking values that range between for highly compressed cells/elements, and for highly stretched cells/elements.If we denote by those cells/elements for which and , then the net mechanical force acting on cell/element is given by:
(3) In equation (3), is the target cell/element radius and the spring function 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 () than it is for pairs of compressed cells/elements. We note that holds for the cutoff distance .

Random force, (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 is given by:(4) where is a unit vector (), pointing in a direction chosen randomly from a uniform distribution, and the constant denotes the vessel element’s sensitivity to random fluctuations.

Chemotactic force, (tip cells)
Following [26, 27], we suppose that cellcell 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 is proportional to the gradient of the VEGF concentration at its centre so that(5) where the positive constant denotes the chemotactic sensitivity. For simplicity, in the simulations that follow we prescribe a fixed chemoattractant gradient, with so that the constant denotes the VEGF gradient.

Directional persistence force, (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 . Accordingly, we denote the directional persistence force, acting on tip cell as follows:(6) where denotes the persistence coefficient. We remark that a similar approach was adopted in [18] where persistence was directly linked with chemotaxis.

Angular persistence force, (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 is not a branching point, and has two neighbours and (so that ), then the angular persistence force acting on vessel element is given by(7) In equation (7), the positive constant denotes the angular spring constant, and is the (smallest) angle between the vectors and . If vessel element is a branch point, connected to more then two other elements, then 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 through the cell cycle is monitored by a phase variable whose evolution is determined by the following differential equation:
(8) 
where the positive constants 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 approximates the length of vessel element : if vessel element has only two neighbours, elements and , say, then
(9) 
If element is a branching point, with more than two neighbours, then is the average over all distances between and its neighbours. Vessel element divides if 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 , 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 where , is the smallest angle between the vessel and the random vector (Figure 2), and is the angle between the plane spanned by the vessel elements connecting cells and and the random vector, . Then the new cell is a tip cell which forms a new capillary sprout if
(10) 
where
(11) 
is defined by Equation (9) and max (). The angles and that appear in Equation (10), and the probabilities with which sprouting or incorporation into the parent vessel happen, are illustrated in Figure 4. Thus, incorporates the influence that mechanical stretch of the has on the fate of its offspring, with less stretched vessel elements having smaller values of , larger values of and, therefore, being more likely to form a new tip cell. In particular, if does not satisfy Equation (10) then the new cell becomes a vessel element and contributes to elongation of the parent vessel. In order to prevent rapid and nonsmooth cell movement following proliferation, spring constants connecting each new cell/element to its neighbours are initially fixed at and then incremented by over a short period (10 timesteps), until they reach .
For details regarding the computational implementation of our model and parameter values, we refer the reader to the Supplementary Materials, Sections S.1 and S.2.
2.2 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: (i) histograms showing the distribution of vessel element lengths, (ii) the total network length, (iii) the number of branches per unit length, (iv) the area covered by the network, (v) the displacement of the initial centre of mass and (vi) the tortuosity of the network (vessel tortuosity is determined by decomposing the network into its component vessels and calculating the ratio of the sum of the lengths of all vessels to the sum of the distances between the endpoints of each vessel element, so that tortuosity ). Metrics (i)–(v) are calculated at fixed time points and the results aggregated to produce summary statistics showing how their mean and variance 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 intervessel distance, interbranch distance, mean vessel diameter, vessel diameter and branching angle in threedimensional 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).
Figure 5 defines the glyphs that we use to visualise network properties. Each glyph has a bounding ellipse which is proportional to the mean bounding ellipse that circumscribes the vascular networks. The main axes of the ellipse are determined by the difference between maximum and minimum values in and direction. One could also argue to take a bounding rectangle, but since the networks develops in a round way, an ellipse is more appropriate. Due to the isotropy of the growth process we only consider the bounding ellipse as a rough measure for the spatial extension. The red pie chart that surrounds the ellipse is proportional to the total length of the network. The red star at the centre of the ellipse indicates the degree of branching of the network while the black hairpin indicates the initial (or default) position of the centre of mass of the network and indicates the extent to which the centre of mass of the network has been displaced during its evolution. The size of the black triangle on the left hand side indicates the chemotactic sensitivity of the network. We remark that chemotactic sensitivity and the displacement of the initial centre of mass are related: stronger chemotaxictic sensitivity will lead to larger displacement of the network’s centre of mass. As stated above, the main motivation for introducing glyphs is to facilitate comparison of networks generated using different parameter values. To effect these comparisons, the metrics that appear in the glyphs are normalised against values from a suitable series of reference simulations.
3 Results
In practice, vessel networks develop in heterogeneous environments and in response to multiple biophysical stimuli that include biochemical signals from angiogenic factors, such as VEGF, and mechanical forces from cellcell and cellsubstrate interactions. Integration of these signals by each cell regulates its phenotype (here, whether it is a tip cell or a vessel element), its rate of progress through the cell cycle, its movement and, thereby, the evolution of the entire network. In this paper we show how mechanical and chemical factors may influence the evolving morphology of vascular networks. We focus here on understanding the behaviour of the mathematical model and how simulation results involving multiple metrics can be presented in a clear and concise manner. Although our work is motivated by results obtained in in vitro experiments, a detailed comparison with experimental data is beyond the scope of 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, (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 mechanosensitivity (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 (tipvesseltip).
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 in which all parameters, except the chemotaxis coefficient, are fixed at their default values (see Supplementary Material, Table S.1). Over time, the network increases in size and develops an intricate pattern of interconnected vessel segments. The glyphs that accompany these images were obtained by averaging over 20 simulations. They reveal how the total vessel length and area spanned by the network increase over time and also how the average number of branches per unit length increases. For these simulations, since chemotaxis is inactive (), the mean and variance of the network’s centre of mass does not change significantly over time and nor do the mean and variance of the tortuosity. By contrast, the mean and variance in the total network length increase as the network expands, as does the mean number of branch points, although the variance decreases over time (see Figure S.2 in Supplementary Material for information). We remark further that the simulation results presented in Figure 6 suggest that realistic vascular networks can be generated in the absence of chemotaxis when mechanical stimuli alone drive cell proliferation and movement.
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)), can lead to ”tip cell tear off.” Here abnormally strong chemotactic forces cause tip cells to detach from the networks which then form degenerate clusters of vessel cells, because the directionality provided by the tipcell is absent (results not shown). These findings highlight the important role that active tipcell movement plays in network formation. We remark further that assigning a tipcell fate to the end cell in a vessel segment following tearoff does not resolve this issue, leading instead to total disintegration of the network as new tip cells rapidly detach from the branch (results not shown). This (unphysical) behaviour is due to a mismatch between tip cell movement, vessel cell proliferation and the strength of cellcell interactions, as represented by the spring constant (see Equation (3)).
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 . For each set of parameter values, we present the structure of networks generated from different random seeds at timesteps () 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, , 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).
length distr.  total length  branches per unit length  disp. of centre of mass  
chemo. sens.  0  0  0  + 
mech. sens.  0  +  0  0 
sprouting prob.  +  +  +  0 
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 predictions about the impact on network morphology of varying particular system parameters.
Motivated by the qualitative comparisons that the glyphs provide and to facilitate future comparisons with experimental data, we now focus on more objective ways of quantifying the impact of parameter variations on network size and morphology. As mentioned in Section 2, the following metrics are calculated for each simulation: (i) histograms showing the distribution of vessel segment lengths, (ii) the total network length, (iii) the number of branches per unit length, (iv) the area spanned by the network, (v) the displacement of the initial centre of mass of the network, and (vi) the tortuosity of the network. For (B)–(E) the mean values are plotted as dots and the standard deviation as bars.
The influence of varying the chemotactic sensitivity, , on the network morphology at 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 tissue matrix and suppresses small scale buckling instabilities (see Figures S.4 and S.5 in Supplementary Material). Changing , the strength of the angular persistence force, has a significant effect on vessel morphology, the stronger angular persistence producing networks which are less dense, have few branch points and are less tortuous (results not shown). The statistics presented in Supplementary Figure S.6 show clearly that increasing skews the distribution of vessel lengths towards longer, straighter vessel segments, that it reduces tortuosity and the overall length of the network since the cells are more compressed in stiffer networks. Variation of the sprouting parameter (see Equations (10) and (11)) has a similar effect on network morphology (see Figure S.7). Increasing 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 varies.
4 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 highresolution 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, offlattice 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 (see Figures S.3(B) and S.7). Finally, increasing the random sensitivity or decreasing the strength of the angular persistence force (i.e., the stiffness of the gel in which the vessels are embedded) produces more tortuous networks (see Figures S.3(A) and S.6).
There are many ways in which the work presented in this paper could be extended. For example, by applying the statistical measures used to characterise the networks generated from our hybrid model to networks generated using other approaches (e.g. cellular automata and Cellular Potts frameworks), it should be possible to identify those network features which are robust to the choice of cellbased framework and those which are not.
A shortcoming of our hybrid model is the introduction of an angular persistence force that mimics the forces that endothelial cells experience due to interactions with their microenvironment and that dampen buckling instabilities induced by cell proliferation (see Figures 3, S.4 and S.5 in Supplementary Material). 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 nonproliferating, 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 cellgel 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 mechanotransduction, 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 tipcell movement, vessel cell proliferation and cellcell 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, 68, 35, 69, 63], we aim to increase the biological realism of our model by embedding, within the vessel cells, mathematical models of subcellular signalling pathways (e.g. deltanotch signalling [25], vascular endothelial growth factor signalling [70, 48] 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 crosstalk 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 shearstress 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 [34, 32, 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 antiangiogenic 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 highresolution 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.
While there are many ways in which our hybrid model could be improved and extended, it demonstrates 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.
Acknowledgements:
This publication was based on work in part sponsored by the “Federal Ministry of Education and Research” Germany, funding initiative “e:Med” Systems Medicine (MultiscaleHCC) and funding initiative “e:bio” Systems Biology (PREDICT). TA gratefully acknowledges the Spanish Ministry of Economy (MINECO) for funding under grant MTM201571509C21R and Generalitat de Catalunya for funding under grant 2009SGR345. PKM was partially supported by the National Cancer Institute, National Institutes of Health grant U54CA143970. MCL acknowledges funding from Moffitt PSOC NIH/NCI U54CA143970. BDH acknowedges funding from the Australian Research Council (DP110100795).
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).
Appendix S.1 The Computational Algorithm
We now outline the algorithm that is used to generate numerical simulations of our hybrid, agentbased model (see, also, Figure S.1). We consider a threedimensional Cartesian domain, of size , discretised uniformly, with spacings , and in the , and directions, respectively. While an offlattice 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 is denoted by . The vascular network structure is stored in an undirected graph in which cell centre coordinates and cellcell connections are recorded. The adjacency matrix contains 2tuples of the node numbers of all connections in the vessel network.

Initialise the simulation ()
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 (tipvesseltip). The initial cell cycle phase of each vessel element and the chemoattractant gradient are also prescribed (). 
Increment time () and loop over all cells ()
The following steps are performed for each element in the network.
Test for cell proliferation (vessel elements only)
If element is a vessel element and then it proliferates and its daughter cell is placed at a random position within a sphere of radius centred on the parent cell. The daughter cell forms a new capillary sprout if its location is such that inequalities (10) are satisfied; otherwise, it becomes a vessel cell and contributes to elongation of the parent vessel. 
Update the position of element
Different forces act on vessel and tip cells (see Equations (1) and (2)). We assume that these forces act on two different timescales: a slow timescale is associated with chemotactic and random movement, and a faster one with mechanically induced cell movement. We exploit this separation of timescales by using a twostep method to integrate the equations of motion and update , the position of element .
Step 1. Calculate the random force using Equation (4). If cell is a tipcell, then calculate the chemotactic and persistence forces and using Equations (5) and (6). Update , the position of cell , 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 using Equation (3) and if cell is a vessel cell then calculate the angular persistence force using Equation (7). If element is a new cell, then its spring constant is increased in increments of on consecutive time steps until it reaches the value . Then , the position of cell , is updated using an explicit Eulerian method to integrate Equation (1) for tip cells or (2) for vessel cells. This step is repeated until the displacement on consecutive steps is less than a threshold or the number of iterations has reached the maximum value, .


Update the cell cycle phase of cell using Equation (8) (vessel elements only)

Update vessel connectivity
New connections are added to the vascular graph, and old ones removed when tip cells anastamose with other vessel segments or detach from their branch (detachment occurs if a tipcell moves beyond the cutoff distance , see Equation (3)).


If then end the simulation; otherwise, return to step 2.
A flowchart of the computational algorithm is given in Figure S.1. The algorithm illustrates how other species (e.g. normal cells, tumour cells, immune cells and tissue matrix) could be incorporated into the computational framework. We remark further that the order in which the algorithmic steps are carried out is arbitrary; reordering generates identical, mean results (results not included).
Appendix S.2 The Default Parameter Values
In the table that follows we state the default parameter values that were used to generate the numerical simulations. We stress that these parameter values were chosen to generate physically realistic results and not by fitting the model to experimental data (an investigation is beyond the scope of the current work).
Parameter  Physical meaning  Default value 
target cell radius  
cell interaction radius  
drag coefficient  
sensitivity to random fluctuations  
spring constant for compressed cells  
spring constant for stretched cells  
chemotactic sensitivity  
chemoattractant gradient  
directional persistence coefficient  
directional persistence time  
angular persistence spring constant  
progression rate through the cell cycle  
cell cycle progression sensitivity to cell elongation and compression  
probability of sprouting  
parameter for sprouting probability  
domain length (direction)  
domain length (direction)  
domain length (direction)  
grid spacing (direction)  
grid spacing (direction)  
grid spacing (direction)  
time step  
max. number of iter. for calculation of tissue deformation per time step 
Appendix S.3 Additional Simulation Results
References
 Folkman J (1995) Angiogenesis in cancer, vascular, rheumatoid and other disease. Nature Med 1: 27–30.
 Risau W (1997) Mechanisms of angiogenesis. Nature 386: 671–674.
 Carmeliet P (2005) Angiogenesis in life, disease and medicine. Nature 438: 932–936.
 Drake C (2003) Embryonic and adult vasculogenesis. Birth Defects Research Part C: Embryo Today: Reviews 69: 73–82.
 Hahnfeldt P, Panigrahy D, Folkman J, Hlatky L (1999) Tumor development under angiogenic signaling: a dynamical theory of tumor growth, treatment response, and postvascular dormancy. Cancer Research 59: 4770–4775.
 Hubbard M, Byrne H (2013) Multiphase modelling of vascular tumour growth in two spatial dimensions. J Theor Biol 29: 1015–1037.
 Vilanova G, Colominas I, Gomez H (2013) Capillary networks in tumour angiogenesis: from discrete endothelial cells to phasefield averaged descriptions via isogeometric analysis. Int Jl Num Meth Biomed Eng 316: 70–89.
 Arakelyan L, Merbl Y, Daugulis P, Ginosar Y, Vainstein V, et al. Multiscale analysis of angiogenic dynamics and therapy. In: Preziosi L, editor, Cancer Modelling and Simulation, CRC Press, LLC, (UK).
 Balding D, McElwain D (1985) A mathematical model of tumourinduced capillary growth. J Theor Biol 114: 53–73.
 Flegg J, McElwain D, Byrne H, Turner I (2009) A three species model to simulate application of hyperbaric oxygen therapy to chronic wounds. PLoS Comp Biol 5: e1000451.
 Machado M, Watson M, Devlin A, Chaplain M, McDougall S, et al. (2011) Dynamics of angiogenesis during wound healing: a coupled in vivo and in silico study. Microcirculation 18: 183–197.
 Edgar L, Sibole S, Underwood C, Guilkey J, Weiss J (2013) A computational model of in vitro angiogenesis based on extracellular matrix fibre orientation. Comp Meth Biomech Biomed Eng 16: 790–801.
 Manoussaki D, Lubkin S, Vernon R, Murray J (1996) A mechanical model for the formation of vascular networks in vitro. Acta Biother 44: 271–282.
 Namy P, Ohayon J, Tracqui P Critical conditions for pattern formation and in vitro tubulogenesis driven by cellular traction fields. J theor Biol 227: 103–120.
 Tosin A, Ambrosi D, Preziosi L Mechanics and chemotaxis in the morphogenesis of vascular networks. Bull Math Biol 68: 1819–1836.
 Stephanou A, Le Floch S, Chauviere A A hybrid model to test the importance of mechanical cues driving cell migration in angiogenesis. Math Model Nat Phenom 10: 142–166.
 Dyson R, Green J, Whiteley J, Byrne H (2015) An investigation of the influence of extracellular matrix anisotropy and cellmatrix interactions on tissue architecture. J Math Biol 23: 1–35.
 Stokes C, Lauffenburger D (1991) Analysis of the roles of microvessel endothelial cell random motility and chemotaxis in angiogenesis. J Theor Biol 152: 377–403.
 Anderson A, Chaplain M (1998) Continuous and discrete mathematical models of tumourinduced angiogenesis. Bull Math Biol 60: 857–899.
 Plank MJ, Sleeman BD (2004) Lattice and nonlattice models of tumour angiogenesis. Bull Math Biol 66: 1785–1819.
 Mantzaris N, Webb S, Othmer H (2004) Mathematical modeling of tumorinduced angiogenesis. J Math Biol 49: 111–187.
 Ambrosi D, Bussolino F, Preziosi L (2005) A review of vasculogenesis models. J Theor Med 6: 1–19.
 Merks R, Koolwijk P (2009) Modeling morphogenesis in silico and in vitro: towards quantitative, predictive, cellbased modeling. Math Model Nat Phenom 4: 149–171.
 Scianna M, Bell C, Preziosi L (2013) A review of mathematical models for the formation of vascular networks. J Theor Biol 333: 174–209.
 Blanco R, Gerhardt H (2013) Vegf and notch in tip and stalk cell selection. Cold Spring Harb Perspect Med 3: a006569.
 Merks R, Perryn E, Shirinifard A, Glazier J (2008) Contactinhibited chemotaxis in de novo and sprouting bloodvessel growth. PLoS Comp Biol 4: e1000163.
 Scianna M, Munaron L, Preziosi L (2011) A multiscale hybrid approach for vasculogenesis and related potential blocking therapies. Prog Biophys Mol Biol 106(2): 450–62.
 Szabo A, Perryn E, Czirok A (2007) Network formation of tissue cells via preferential attraction to elongated structures. Phys Rev Lett 98: 38102.
 Szabo A, Mehes E, Kosa E, Czirok A (2008) Multicellular sprouting in vitro. Biophys J 95: 2702–2710.
 Szabó A, Varga K, Garay T, Hegedűs B, Czirók A (2012) Invasion from a cell aggregateâthe roles of active cell motion and mechanical equilibrium. Phys Biol 9: 016010.
 van Oers R, Rens E, LaValley D, ReinhartKing C, Merks R (2014) Mechanical cellmatrix feedback explains pairwise and collective endothelial cell behavior in vitro. PLOS Comp Biol 10: e1003774.
 Stephanou A, McDougall S, Anderson A, Chaplain M (2005) Mathematical modelling of flow in 2d and 3d vascular networks: applications to antiangiogenic and chemotherapeutic drug strategies. Math Comp Model 41: 1137–1156.
 Watson M, McDougall S, Chaplain M, Devlin A, Mitchell C (2012) Dynamics of angiogenesis during murine retinal development: a coupled in vivo and in silico study. J Roy Soc Interface 13.
 Alarcon T, Byrne H, Maini P (2005) A multiple scale model for tumour growth. Multiscale Model Simul 3: 440–475.
 Owen M, Alarcón T, Maini P, Byrne H (2009) Angiogenesis and vascular remodelling in normal and cancerous tissues. J Math Biol 58: 689–722.
 Perfahl H, Byrne H, Chen T, Estrella V, Alarcón T, et al. (2011) Multiscale modelling of vascular tumour growth in 3d: the roles of domain size and boundary conditions. PLoS One 6: e14790.
 Macklin P, McDougall S, Anderson ARA, Chaplain MAJ, Cristini V, et al. (2009) Multiscale modelling and nonlinear simulation of vascular tumour growth. J Math Biol 58: 765–798.
 Shirinifard A, Gens JS, Zaiden BL, Poplawski NJ, Swat M, et al. (2009) 3d multicell simulation of tumour growth and angiogenesis. PLoS One 4: e7190.
 Welter M, Bartha K, Rieger H (2008) Emergent vascular network inhomogeneities and resulting blood flow patterns in a growing tumor. J Theor Biol 250: 257–280.
 Welter M, Bartha K, Rieger H (2009) Vascular remodelling of an arteriovenous blood vessel network during solid tumour growth. J Theor Biol 259: 405422.
 Hoehme S, Brulport M, Bauer A, Bedawy E, Schormann W, et al. (2010) Prediction and validation of cell alignment along microvessels as order principle to restore tissue architecture in liver regeneration. Proc Natl Acad Sci U S A 107: 10371–10376.
 McLennan R, Dyson L, Prather K, Morrison J, Baker R, et al. (2012) Multiscale mechanisms of cell migration during development: theory and experiment. Development 139: 2935–2944.
 Schaller G, MeyerHermann M (2005) Multicellular tumor spheroid in an offlattice VoronoiDelaunay cell model. Phys Rev E 71: 051910.
 PittFrancis J, Pathmanathan P, Bernabeu M, Bordas R, Cooper J, et al. (2009) Chaste: a testdriven approach to software development for biological modelling. Comp Phys Commun 180: 2452–2471.
 Jackson T, Zheng X (2010) A cellbased model of endothelial cell migration, proliferation and maturation during corneal angiogenesis. Bull Math Biol 72: 830–868.
 Drasdo D, Jagiella N, I RamisConde I, VignonClementel I, Weens W (2010) Modeling steps from a benign tumor to an invasive cancer: examples of intrinsically multiscale problems. In: Chauviere A, Preziosi L, Verdier C, editors, Cell Mechanics: From Single ScaleBased Models to Multiscale Modeling, Chapman & Hall/CRC. pp. 379–417.
 Liu W, Nelson C, Tan J, Chen C (2007) 3d microvascular architecture of precancerous lesions and invasive carcinomas of the colon. Circ Res 101: e44e52.
 Zheng W, Christensen L, Tomanek R (2008) Differential effects of cyclic and static stretch on coronary microvascular endothelial cell receptors and vasculogenic/angiogenic responses. Am J Physiol Heart Circ Physiol 205: H794H800.
 Maguire E, RoccaSerra P, Sansone S, Davies J, Chen M (2012) Taxonomybased glyph design–with a case study on visualizing workflows of biological experiments. IEEE T Vis Comp Gr 18: 2603–2612.
 Hellström M, Phng LK, Hofmann J, Wallgard E, Coultas L, et al. (2007) Dll4 signalling through notch1 regulates formation of tip cells during angiogenesis. Nature 445: 776–780.
 Drasdo D, Höhme S (2005) A singlecellbased model of tumor growth in vitro: monolayers and spheroids. Phys Biol 2: 133.
 Drasdo D, Loeffler M (2001) Individualbased models to growth and folding in onelayered tissues: intestinal crypts and early development. Nonlinear Analysis – Theory Methods and Applications 47: 245–256.
 Meineke F, Potten C, Loeffler M (2001) Cell migration and organization in the intestinal crypt using a latticefree model. Cell Prolif 34: 253–266.
 Walker D, Southgate J, Hill G, Holcombe M, Hose D, et al. (2004) The epitheliome: agentbased modelling of the social behaviour of cells. Biosystems 76: 89–100.
 Owen M, Stamper I, Muthana M, Richardson G, Dobson J, et al. (2011) Mathematical modeling predicts synergistic antitumor effects of combining a macrophagebased, hypoxiatargeted gene therapy with chemotherapy. Cancer Res 71: 2826.
 Konerding M, Malkush W, Klapthor B, van Ackem C, Fait E, et al. (1999) Evidence for characteristic vascular patterns in solid tumours: quantitative studies using corrosion casts. Brit J Cancer 80: 724–732.
 Konerding M, Fait E, Gaumann A (2001) 3d microvascular architecture of precancerous lesions and invasive carcinomas of the colon. Brit J Cancer 84: 1354.
 Folarin A, Konerding M, Timonen J, Nagl S, edley (2010) Evidence for characteristic vascular patterns in solid tumours: quantitative studies using corrosion casts. Microvasc Res 80: 89–98.
 Selmeczi D, Mosler S, Hagedorn PH, Larsen NB, Flyvbjerg H (2005) Cell motility as persistent random motion: Theories from experiments. Biophys J 89: 912–931.
 Schlueter D, RamisCondeI, Chaplain M (2012) Computational modeling of singlecell migration: The leading role of extracellular matrix fibers. Biophys J 103: 1141–1151.
 Lee B, Zhou X, Riching K, Eliceiri K, Keely P, et al. (2014) A threedimensional computational model of collagen network mechanics. PLoS One 9: e111896.
 Riching K, Cox B, Salick M, Pehlke C, Riching A, et al. (2014) 3d collagen alignment limits protrusions to enhance breast cancer cell persistence. Biophys J 107: 2546–2558.
 Schlueter D, RamisConde I, Chaplain M (2015) Multiscale modelling of the dynamics of cell colonies: insights into cell adhesion forces and cancer invasion from in silico simulations. J R Soc Interface 12: 20141080.
 Byrne H, Preziosi L (2003) Modelling solid tumour growth using the theory of mixtures. Math Med Biol 20: 341–366.
 Chaplain M, Graziano L, Preziosi L (2006) Mathematical modelling of the loss of tissue compression responsiveness and its role in solid tumour development. Math Med Biol 23: 197–229.
 Guarani V, Deflorian G, Franco C, Krüger M, Phng L, et al. (2011) Acetylationdependent regulation of endothelial notch signalling by the sirt1 deacetylase. Nature 473: 234–238.
 Bentley K, Mariggi G, Gerhardt H, Bates P (2009) Tipping the balance: robustness of tip cell selection, migration and fusion in angiogenesis. PLoS Comp Biol 5: e1000549.
 RamisConde I, Drasdo D, Chaplain M, Anderson A (2008) Modelling the influence of the ecadherin  betacatenin pathway in cancer cell invasion and tissue architecture: A multiscale approach. Biophys J 95: 155–165.
 Andasari V, Roper R, Swat M, Chaplain M (2012) Integrating intracellular dynamics using compucell3d and bionetsolver: Applications to multiscale modelling of cancer cell growth and invasion. PLoS One 7: e33726.
 Stefanini M, Qutub A, MacGabhann F, Popel A (2012) Computational models of vegfassociated angiogenic processes in cancer. Math Med Biol 29: 85–94.
 Pries A, Secomb T, Gaehtgens P (1998) Structural adaptation and stability of microvascular networks: theory and simulations. Amer J Physiol – Heart and Circulatory Physiol 275: H349–H360.
 Alarcon T, Byrne H, Maini P (2005) A design principle for vascular beds: the effects of complex blood rheology. Microvasc Res 69: 156–172.
 Jain R (2005) Normalisation of tumour vasculature: an emerging concept in antiangiogenic therapy. Science 307: 58–62.
 Carmeliet P, Jain R (2011) Principles and mechanisms of vessel normalisation for cancer and other angiogenic diseases. Nature Reviews Drug Discovery 10: 417–427.
 Department of Radiation Oncology MGH (2013) E.L. Steel Laporatory Research Report 2013. Harvard Medical School.
 Green S (2012). Particle simulation using CUDA. Website, CUDA documentation. http://docs.nvidia.com/cuda/samples/5_Simulations/particles/doc/particles.pdf.