PLOS COMPUTATIONAL BIOLOGY RESEARCH ARTICLE A multiscale model of complex endothelial cell dynamics in early angiogenesis Daria Stepanova ID1,2*, Helen M. Byrne ID3, Philip K. Maini ID3, Tomás Alarcón ID1,2,4,5 1 Centre de Recerca Matemàtica, Bellaterra (Barcelona), Spain, 2 Departament de Matemàtiques, Universitat Autònoma de Barcelona, Bellaterra (Barcelona), Spain, 3 Wolfson Centre for Mathematical Biology, Mathematical Institute, University of Oxford, Oxford, UK, 4 Institució Catalana de Recerca i Estudis Avançats (ICREA), Barcelona, Spain, 5 Barcelona Graduate School of Mathematics (BGSMath), Barcelona, Spain a1111111111 * dstepanova@crm.cat a1111111111 a1111111111 a1111111111 a1111111111 Abstract We introduce a hybrid two-dimensional multiscale model of angiogenesis, the process by which endothelial cells (ECs) migrate from a pre-existing vascular bed in response to local environmental cues and cell-cell interactions, to create a new vascular network. Recent OPEN ACCESS experimental studies have highlighted a central role of cell rearrangements in the formation Citation: Stepanova D, Byrne HM, Maini PK, of angiogenic networks. Our model accounts for this phenomenon via the heterogeneous Alarcón T (2021) A multiscale model of complex endothelial cell dynamics in early angiogenesis. response of ECs to their microenvironment. These cell rearrangements, in turn, dynamically PLoS Comput Biol 17(1): e1008055. https://doi. remodel the local environment. The model reproduces characteristic features of angiogenic org/10.1371/journal.pcbi.1008055 sprouting that include branching, chemotactic sensitivity, the brush border effect, and cell Editor: Roeland M.H. Merks, Universiteit Leiden, mixing. These properties, rather than being hardwired into the model, emerge naturally from NETHERLANDS the gene expression patterns of individual cells. After calibrating and validating our model Received: June 10, 2020 against experimental data, we use it to predict how the structure of the vascular network Accepted: November 19, 2020 changes as the baseline gene expression levels of the VEGF-Delta-Notch pathway, and the composition of the extracellular environment, vary. In order to investigate the impact of cell Published: January 7, 2021 rearrangements on the vascular network structure, we introduce the mixing measure, a sca- Peer Review History: PLOS recognizes the lar metric that quantifies cell mixing as the vascular network grows. We calculate the mixing benefits of transparency in the peer review process; therefore, we enable the publication of measure for the simulated vascular networks generated by ECs of different lineages (wild all of the content of peer review and author type cells and mutant cells with impaired expression of a specific receptor). Our results responses alongside final, published articles. The show that the time evolution of the mixing measure is directly correlated to the generic fea- editorial history of this article is available here: tures of the vascular branching pattern, thus, supporting the hypothesis that cell rearrange- https://doi.org/10.1371/journal.pcbi.1008055 ments play an essential role in sprouting angiogenesis. Furthermore, we predict that lower Copyright: © 2021 Stepanova et al. This is an open cell rearrangement leads to an imbalance between branching and sprout elongation. Since access article distributed under the terms of the Creative Commons Attribution License, which the computation of this statistic requires only individual cell trajectories, it can be computed permits unrestricted use, distribution, and for networks generated in biological experiments, making it a potential biomarker for patho- reproduction in any medium, provided the original logical angiogenesis. author and source are credited. Data Availability Statement: All relevant data are within the manuscript and its Supporting information files. Funding: This work is supported by a grant of the Obra Social La Caixa Foundation on Collaborative PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1008055 January 7, 2021 1 / 44 PLOS COMPUTATIONAL BIOLOGY A multiscale model of complex endothelial cell dynamics in early angiogenesis Mathematics awarded to the Centre de Recerca Matemàtica through a scholarship awarded to D.S. Author summary D.S. and T.A. have been partially funded by the CERCA Programme of the Generalitat de Catalunya. Angiogenesis, the process by which new blood vessels are formed by sprouting from the They also acknowledge MINECO (https://www. pre-existing vascular bed, plays a key role in both physiological and pathological pro- ciencia.gob.es/) for funding under grants cesses, including tumour growth. The structure of a growing vascular network is deter- MTM2015-71509-C2-1-R and RTI2018-098322-B- mined by the coordinated behaviour of endothelial cells in response to various signalling I00. D.S. and T.A. participate in project cues. Recent experimental studies have highlighted the importance of cell rearrangements 2017SGR01735 which was awarded by AGAUR (https://agaur.gencat.cat/en/inici/index.html) but as a driver for sprout elongation. However, the functional role of this phenomenon with no actual funding. The funders had no role in remains unclear. We formulate a new multiscale model of angiogenesis which, by study design, data collection and analysis, decision accounting explicitly for the complex dynamics of endothelial cells within growing angio- to publish, or preparation of the manuscript. genic sprouts, is able to reproduce generic features of angiogenic structures (branching, Competing interests: The authors have declared chemotactic sensitivity, cell mixing, etc.) as emergent properties of its dynamics. We vali- that no competing interests exist. date our model against experimental data and then use it to quantify the phenomenon of cell mixing in vascular networks generated by endothelial cells of different lineages. Our results show that there is a direct correlation between the time evolution of cell mixing in a growing vascular network and its branching structure, thus paving the way for under- standing the functional role of cell rearrangements in angiogenesis. Introduction Angiogenesis, the process by which new blood vessels are generated from a pre-existing vascu- lar network, has been extensively investigated in recent years [1–10]. However, understanding of the complex behaviour of endothelial cells (ECs) during angiogenesis remains incomplete. ECs adjust their gene expression profiles in response to a variety of extracellular cues, such as the structure and composition of the extracellular matrix (ECM), angiogenic factors such as vascular endothelial growth factor (VEGF), and changes in cell-cell interactions. Broadly, the extracellular stimuli act upon the VEGF-Delta-Notch signalling pathway, which drives ECs to acquire either a tip or stalk cell phenotype [11, 12]. Tip cells, with elevated expression of Delta ligand, VEGF receptor 2 (VEGFR2) and decreased expression of Notch receptor, are more active than their stalk cell counterparts [13]; they extend filopodia, release matrix metallopro- teinases (MMPs) that degrade the ECM and, together with the pericytes that they recruit, secrete basal lamina components that stabilise growing vessels. By contrast, stalk cells, charac- terised by low expression of Delta and VEGFR2, and high Notch, migrate along the paths explored by tip cells, and proliferate to replenish the elongating sprout. EC phenotypes are not static and strongly depend on the local microenvironment (signalling cues, interactions with other ECs, among others) [2, 7]. Furthermore, ECs within sprouts exchange their relative posi- tions via a phenomenon called cell mixing [2, 7]. Every time a cell rearrangement takes place, the cells’ microenvironment changes. This, in turn, leads to recurrent (re-)establishment of phenotypes and variations in gene expression within the phenotype. As a consequence, a dynamic coupling between cell phenotypes and rearrangements is established. The functional role of cell mixing, and how it is affected by variations in the gene expression patterns of ECs, are unclear, although it is acknowledged that cell rearrangements greatly influence the pattern of the vascular network and its functionality [1, 7, 14, 15]. Angiogenic sprouting has been extensively studied from a theoretical modelling perspective in numerous physiological and pathological contexts, including tumour growth [10, 16–33]. Many models [18–20, 22, 24–28, 32] fall within the so-called snail-trail paradigm [31]. In this approach, EC phenotype acquisition is assumed to be irreversible; a subset of cells, associated with the tip phenotype and situated at the leading edge of growing sprouts, responds PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1008055 January 7, 2021 2 / 44 PLOS COMPUTATIONAL BIOLOGY A multiscale model of complex endothelial cell dynamics in early angiogenesis chemotactically to gradients of signalling cues such as VEGF. The rest of the cells are assumed to be stalk cells which passively follow cells at the leading edge of the sprouts. Cell rearrange- ments are neglected and properties such as branching are hardwired via ad hoc rules, where the probability of branching depends on environmental factors [18–20, 24, 25, 28]. Further- more, in the majority of models, chemotactic behaviour is coded into the model via an explicit bias of the migration probabilities towards regions of higher VEGF concentration [18–20, 24– 28]. As a result, these models are unable to explain how variations in subcellular signalling, which determine cell phenotype, modify the structure of the growing network, or how specific mutations might lead to pathological networks. A more extensive review on the mathematical and computational models of angiogenesis can be found in [34–36]. Several recent computational models have started to address these issues [3, 37]. For exam- ple, Bentley et al. [3] studied the mechanisms underlying cell rearrangement behaviour during angiogenic sprouting. Using computational modelling and experimentation, they identified the Notch/VEGFR-regulated dynamics of EC adhesion as a key driver of EC rearrangements. Furthermore, simulation and quantitative image analysis indicated that abnormal phenotype synchronisation exists under pathologically high VEGF conditions, in retinopathy and tumour vessels [3]. However, due to the computational intensity of the model, only small scale simula- tions (a single vessel with 10 ECs) were performed. Boas and Merks [37] have also studied EC rearrangement within small vascular networks. They investigated how a cell at the leading position in a sprout can be overtaken by another cell, a phenomenon called cell overtaking. Their simulation results, obtained using the cellular Potts framework, suggest that cell overtak- ing occurs in an unregulated fashion, due to the stochastic motion of ECs. Moreover, their findings suggest that the role of the VEGF-Delta-Notch pathway in cell overtaking is not to select a cell with the highest expression of Delta ligand and VEGFR2, corresponding to the tip phenotype, that will shuffle up towards the sprout tip. Rather, cell overtaking is proposed to ensure that a cell that ends up occupying the leading position in a vessel, due to its spontaneous migration backward and forward within the sprout, acquires the tip phenotype. In their model, time-dependent ordinary differential equations (ODEs) are used to simulate the dynamics of subcellular signalling pathways within each cell. Cells are assigned a discrete phe- notype, tip or stalk, and variations in their gene expression levels that might influence the behaviour of the cell were neglected. Cell rearrangements have been proposed to be the main driver for early sprouting angio- genesis [1–3, 7]. In particular, reduced cell mixing leads to formation of pathological networks characterised by superimposed aberrant layers of vessels [1]. In order to characterise this phe- nomenon in a quantitative manner and determine its influence on the growth of vascular net- works, we introduce a multiscale mathematical model of early angiogenesis (i.e. on a time scale of hours where cell birth/death are negligible). The model accounts for individual cell gene expression patterns associated with the VEGF-Delta-Notch signalling pathway that define two distinct cell phenotypes. Gene expression patterns are dynamically updated upon cell migration and in response to fluctuations in the cell microenvironment. Since cell behaviour has been shown to be heterogeneous and vary significantly depending on phenotype, cell- ECM interactions (ECM proteolysis, ECM remodelling, alignment of ECM fibrils) are assumed to be phenotype-dependent. In turn, the ECM configuration directly influences EC polarity and migration, thus changing the local environment of cells and providing coupling between the subcellular, cellular and tissue scales of the model. One of our main results is that the model is capable of reproducing the typical behaviour of ECs; branching, chemotactic sensitivity and cell mixing are emergent properties of the model (instead of being hardwired into it) that arise as a result of cell-ECM interactions involving cells with dynamic subcellular gene expression. By calibrating and validating our model PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1008055 January 7, 2021 3 / 44 PLOS COMPUTATIONAL BIOLOGY A multiscale model of complex endothelial cell dynamics in early angiogenesis against experimental data, we can investigate the role that abnormal subcellular signalling has on cell-matrix interactions, cell rearrangements and the general structure of the growing vas- cular network. Notably, we show how externally imposed changes in the extracellular environ- ment (ECM density, VEGF concentration, VEGF gradient) and mutations in gene expression of ECs can alter the branching pattern of growing sprouts. To quantify cell mixing within the vascular network, we introduce a mixing measure. This scalar metric is calculated for cells that are immediate neighbours at a given time moment, characterising how far they spread within the vascular network migrating during a fixed time interval. Our results suggest that in a formed network, temporal evolution of the mixing mea- sure reaches a steady state that depends on the relative proportion of EC tip and stalk pheno- types. However, during the early stages of network formation, the time evolution of the mixing measure varies as the extracellular environment and cells’ expression levels of genes such as VEGFR1 and VEGFR2 change. In particular, we show that networks created by mutant cells with impaired gene expression of VEGF receptors exhibit a delayed mixing compared to the networks formed by wild type cells with unaltered gene expression. The networks formed by these mutants demonstrate an imbalance between effective sprout elongation and branching. Therefore, our results confirm the crucial role of cell rearrangements in the formation of func- tional vascular networks in sprouting angiogenesis [1, 7]. The remainder of this paper is organised as follows. In the Experimental motivation sec- tion, we summarise the setups and results of several experimental studies which motivated the formulation of our multiscale model. The Materials and methods section contains a summary of the relevant biological information, a description of our multiscale model and metrics used to analyse vascular network evolution. In the Results section we compare our simulation results with data from Arima et al. [2] and Shamloo & Heilshorn [38]. These data were extracted from in vitro experiments, and enable us to define a set of basal parameter values for our model. Further model validation is performed against experimental results involving mutant cells, with modified expression of VEGF receptors, carried out by Jakobsson et al. [7]. Finally, we present results on EC mixing quantification and show how it relates to the different branching patterns of the growing vascular networks. We conclude by drawing together our findings and outlining possible avenues for future work in the Discussion section. Experimental motivation The model we develop is motivated by in vitro experiments in which an aortic ring assay was embedded into a collagen matrix with a uniform VEGF concentration (0, 5 or 50 ng/ml) [2, 8]. Computational analysis of dynamic images, collected using time-lapse microscopy, revealed complex dynamical cell rearrangements within growing sprouts, a phenomenon termed cell mixing. The authors concluded that over short periods of time (e.g. 22.4 h averaged over all experiments in [2]) cell rearrangements are the main driver of sprout elongation. Interestingly, successful sprout growth was seen in a uniform distribution of VEGF across the substrate— sprout elongation velocity was observed to vary as the concentration of the external VEGF was changed, but no VEGF gradient was necessary for coordinated migration of ECs. Dynamic cell rearrangements within elongating sprouts are a direct consequence of cells continuously updating their phenotype (i.e. adapting their gene expression pattern, depending on their environment [1, 3, 7]). Jakobsson et al. [7] identified the VEGF-Delta-Notch signal- ling pathway as the key pathway controlling this phenomenon. Using mutant cells heterozy- gous for VEGFR1 (VEGFR1+/-) and VEGFR2 (VEGFR2+/-) with halved (compared to wild- type (WT) cells) gene expression of the corresponding VEGF receptor, they investigated how differential levels of VEGF receptors affect the probability that a WT or mutant EC will occupy PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1008055 January 7, 2021 4 / 44 PLOS COMPUTATIONAL BIOLOGY A multiscale model of complex endothelial cell dynamics in early angiogenesis the leading position in a growing sprout. Embryoid bodies (three-dimensional spherical aggre- gates of cells) derived from WT cells mixed with one of the populations of mutant cells (50% and 50%) were indistinguishable from those formed by WT cells only, however, the contribu- tion of each cell line to the leading position differed. In particular, VEGFR1+/- (VEGFR2+/-) cells demonstrated enhanced (reduced) competition for the leading cell position. The role of the VEGF-Delta-Notch signalling pathway in establishing competitive advantage was rein- forced by experiments with the DAPT inhibitor, which abolishes Notch signalling in all EC lineages. Treatment with DAPT, although leading to hyper-sprouting, restored balance in competition for the leading cell position. Motivated by these results and our interest in study- ing cell rearrangements, we account for the VEGF-Delta-Notch signalling pathway of individ- ual cells in our model. The coordinated migration of ECs is a result of cell-ECM interactions [39–43]. Specifically, cell migration depends on the EC ability to degrade ECM proteins via proteolysis in order to form ECM-free vascular guidance tunnels for effective sprout elongation. This was confirmed by experimental results performed by Shamloo et al. [38]. Therein, ECs were cultured within a microfluidic device with a maintained gradient of VEGF concentration and the response of ECs to variations in ECM components, specifically collagen, was considered. Their results showed that there is a “sweet spot” of collagen density for the formation of angiogenic sprouts. In low collagen densities, ECs migrate freely into the ECM without forming sprouts; at inter- mediate concentrations, structures resembling angiogenic sprouts form; at high collagen den- sities, ECs are unable to migrate into the matrix significantly and form thick short protrusions. These experimental results motivated us to include cell-ECM interactions in our model. Most in vitro experiments are performed on flat substrates in which the depth of the sub- strate can be considered negligible compared to its length and width [2, 8, 38, 40]. Thus, we formulate our model in a two-dimensional framework. Materials and methods Biological summary Angiogenic sprouting is initiated when an active cytokine, such as VEGF, reaches the existing vasculature and activates quiescent ECs. Active ECs express proteases that degrade the base- ment membrane (BM) so that EC extravasation and migration can occur [44, 45]. In vivo and in vitro, coordinated migration is guided by local chemical and mechanical cues [42, 46, 47]. While VEGF-activated ECs at the front of the vascular network explore new space and anasto- mose, to form branching patterns, previously created vascular sprouts mature, form lumens and, with the help of murine cells (pericytes and smooth muscle cells), stabilise [48]. In our mathematical model, we focus on EC migration, and its regulation by local environmental (mechanical and chemical) cues. Due to the short time scales considered in the model, pro- cesses such as proliferation, vessel maturation and lumen formation are neglected. VEGF-induced actin polymerisation and focal adhesion assembly result in substantial cyto- skeletal remodelling as required for cell migration [49]. This includes actin remodelling to form filopodia and lamellipodia, stress fibre formation, and focal adhesion turnover [41]. A crucial consequence of such remodelling is that ECs acquire front-to-rear polarity which coin- cides with the direction of their migration [50]. Membrane protrusions, such as filopodia and lamellipodia, increase the cell surface area at the leading edge. More VEGF receptors and integrins, which bind to ECM components, become activated at the cell front to further rein- force the cytoskeleton remodelling and polarisation [51]. This mechanism is known as taxis (e.g. chemotaxis, haptotaxis) and allows ECs to sense extracellular cues and migrate towards them. The cells’ chemotactic sensitivity leads to the so-called brush-border effect which is PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1008055 January 7, 2021 5 / 44 PLOS COMPUTATIONAL BIOLOGY A multiscale model of complex endothelial cell dynamics in early angiogenesis characterised by increased rates of branching and higher EC densities closer to the source of the VEGF stimulus which can be sensed by cells and biases their migration towards higher VEGF concentrations. Whilst interactions with ECM components play a central role in EC polarisation and move- ment, EC migration, in turn, reorganises and remodels the ECM [45]. Prior to assembly of the basement membrane of the newly formed vessels, the ECM microenvironment consists mostly of collagen I and elastin fibers. Activated ECs secrete matrix type 1 metalloproteinases (MT1-MMPs) that degrade the ECM [44, 46, 52]. This process generates ECM-free tunnels into which the sprouts can elongate [52]. As sprouts grow, they assemble a basement mem- brane which contains, among other things, fibrous components (collagen IV, fibronectin and various laminins [44]) that are secreted by the ECs and whose function is to promote cell-cell and cell-ECM contact and to limit EC migration [46]. Further matrix reorganisation occurs in response to the mechanical forces generated by migrating cells as they realign the collagen fibrils of the ECM. Several experiments have shown that cells with extended filopodia and lamellipodia form focal adhesions with the ECM compo- nents and realign them by pulling in the direction parallel to their motion [39, 40, 51, 53]. They also move the fibrils from the local neighbourhood closer to their surfaces. As a result, collagen fibrils accumulate and align along the direction of sprout elongation. Since cell-fol- lowers form focal adhesions with these aligned fibrils they automatically polarize and migrate in the direction of sprout elongation [42, 46, 47]. In this way sprout integrity is maintained and the coordinated motion of ECs emerges. The default phenotype of an EC that has been activated by VEGF is a “tip” cell [12, 13]. Tip cells are characterised by exploratory behaviour and low proliferative activity. Thus, although the ECM regulates EC behaviour to a great extent [45], if all cells were to move at the same time in an exploratory fashion, sprout integrity would be lost. In non-pathological angiogene- sis, typical vascular network morphology is obtained when EC proliferation and migration (i.e. the exploration of space by active ECs in response to signalling cues) are properly balanced. This balance is mediated by the Delta-Notch signalling pathway, which provides ECs with a contact-dependent cross-talk mechanism by which they can acquire a “stalk” cell phenotype [12, 13]; stalk cells are characterised by reduced migratory and higher proliferative activities. In this way, cell-cell communication allows cells to adjust their gene expression patterns (and, hence, phenotype) and to coordinate their behaviour [54]. The processes influenced by cell phenotype in our model are illustrated in Fig 1A. The ratio between the number of tip and stalk cells during sprouting plays a key role in the integrity of the developing vascular network and its functionality. This was confirmed by experiments in which Notch activity was inhibited or completely blocked. This caused all cells to adopt a tip cell phenotype and led to pathological network formation [7]. Mathematical modelling of random cell migration Stochastic approaches have been successfully used to mathematically model systems with ran- dom, heterogeneous behaviour, including spatially extended systems [27, 55–64]. We model cell migration at the cellular scale stochastically to account for random cell motility, cell mix- ing, and branching dynamics. This approach allows us to compare our model with in vitro experiments which, as any biological system, are noisy. Within the context of spatially extended systems, a stochastic model can be formulated as an individual-based, or a compartment-based, model. In the former approach, the Brownian dynamics of each individual agent (cells, molecules, etc.) are simulated, which becomes com- putationally expensive as the number of agents increases. The latter is more numerically PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1008055 January 7, 2021 6 / 44 PLOS COMPUTATIONAL BIOLOGY A multiscale model of complex endothelial cell dynamics in early angiogenesis Fig 1. Schematic summary of our model. (A) The gene expression profile of each cell (and, hence, its phenotype) is accounted for at the subcellular scale of our model. It influences such processes as I. cell polarity and branching, II. overtaking, III. extracellular matrix (ECM) proteolysis, IV. basement membrane (BM) assembly and V. ECM realignment. The cartoons on the left illustrate each process, whereas the text- boxes on the right provide brief descriptions of how they are influenced by tip and stalk cells. (B) The structure of our multiscale model. The diagram illustrates the processes that act at each spatial scale. Arrows illustrate coupling between the scales. https://doi.org/10.1371/journal.pcbi.1008055.g001 PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1008055 January 7, 2021 7 / 44 PLOS COMPUTATIONAL BIOLOGY A multiscale model of complex endothelial cell dynamics in early angiogenesis efficient for the purpose of simulating small vascular networks containing at most hundreds of cells. In this method the domain is partitioned into non-overlapping compartments, or voxels so that the position of each agent is known up to the voxel size. Agents within the same voxel are considered indistinguishable and reactions between them occur independently of those in other compartments. Movement of agents between the voxels is modelled as a continuous time random walk (RW). This is also known as the Reaction Diffusion Master Equation (RDME) approach. A weakness of the RDME approach is that it does not converge to its individual-based Brownian dynamics in dimensions greater than one for compartments with size less than a specific lower bound (of the order of the reaction radius in the Smoluchowski interpretation) [65]. The RDME breaks down as the voxel size tends to zero because the waiting time for multi-molecular reactions becomes infinite. One strategy to address this problem, the so-called convergent RDME (cRDME), is an approximation of Doi’s model for the binary reactions of the form A + B ! C proposed in [66]. In the cRDME, agents from different voxels may interact via multi-molecular reactions provided they lie within a predefined interaction radius. In this way the artefact of vanishing reaction rates as voxel size tends to zero is avoided. We formulate the cellular scale of our model using the cRDME approach. We introduce an interaction radius, Rc, within which cells are assumed to interact. This approach also fits natu- rally with the VEGF-Delta-Notch driven cell cross-talk that we incorporate at the subcellular scale of our multi-scale model (see below). For the numerical implementation of the cRDME we use the Next Subvolume (NSV) method [67], which is a computationally efficient implementation of the standard Stochastic Simulation Algorithm [68] to simulate RDME/cRDME. Summary of the multiscale model The model we develop is a two-dimensional stochastic multiscale model of migration-driven sprouting angiogenesis. Its structure is shown in Fig 1B. Briefly, the model operates on three distinct spatial scales: • The subcellular scale defines the gene expression pattern of individual cells. Since ECs con- tain a finite number of proteins, some level of noise is always present in the system. Thus, we implement a stochastic model of the VEGF-Delta-Notch signalling pathway to describe the temporal dynamics of the number of ligands/receptors for each cell. This pathway is known to produce bistable behaviour. Cells exhibit either high Delta and VEGFR2, and low Notch levels (the tip phenotype) or low Delta and VEGFR2, and high Notch levels (the stalk pheno- type). Stochasticity allows random transitions between these phenotypes in regions of bist- ability, behaviour which cannot be achieved in a deterministic model. • The cellular scale accounts for cell migration. It is formulated as a variant of an on-lattice persistent random walk (PRW) of ECs. • The tissue scale keeps track of the local ECM environment of the cells. Local ODEs track the evolution of the concentrations of the existing ECM and BM components, whereas ECM fibril alignment driven by EC movement is updated using a phenomenological model. An illustration of the model geometry can be found in Fig 2. In general, an EC has an arbitrary shape which depends on its cell-cell and cell-matrix focal adhesions. In our model we do not keep track of the exact cell shape and assume that cell posi- tion is known up to the position of its nucleus. Thus, when referring to a cell position, we refer to the position of its nucleus and assume that the cell has some arbitrary shape centred on the PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1008055 January 7, 2021 8 / 44 PLOS COMPUTATIONAL BIOLOGY A multiscale model of complex endothelial cell dynamics in early angiogenesis Fig 2. Model geometry. (A) To account for cell migration in the framework of a PRW, we cover the domain, D � R2 , by a uniform SNI pffiffi lattice L ¼ fvk : k¼1 vk � Dg of non-overlapping hexagonal voxels, vk, of width h (or hexagon edge, a, h ¼ 3a). We denote by T T I ¼ fk ¼ ðkx ; ky Þ : vk 2 Lg the set of 2D indices, with cardinal, jIj ¼ NI . The coordinates of the centre of voxel vk are denoted by qk ¼ ðqxk ; qyk Þ . We assume that there is a constant supply of ECs to the domain, coming from the vascular plexus. These ECs enter the domain at fixed boundary voxels, defined as a set, I VP (coloured in green). (B) A detailed illustration of an individual voxel, vi. s = h−1(qj − qi) denotes a normalised vector of the migration direction between neighbouring voxels, vi and vj. There are at most 6 possible migration directions for each hexagonal voxel. Each migration direction, s, can be characterised by an equivalent angle interval ½�smin ; �smax �. https://doi.org/10.1371/journal.pcbi.1008055.g002 voxel containing its nucleus (see Fig 3A). A consequence of this approach is that, since cells can extend membrane protrusions and interact with distinct cells beyond their first neigh- bours, interactions between cells in our model are non-local. We introduce two interaction radii, Rs and Rc, for the subcellular and cellular scales, respectively, and assume that a cell can interact with any neighbouring cell partially overlapped by a circular neighbourhood of these interaction radii (see Fig 3A). In particular, at the subcellular scale trans-binding between a Delta ligand on one cell and a Notch receptor on another can occur if the distance between their cell centres is less than the interaction radius, Rs. Thus, the total amount of ligand/recep- tor (belonging to a neighbour/neighbours) to which a cell is exposed is proportional to the sur- face area of the overlap region between the circular neighbourhood and the neighbouring voxel/voxels (not necessarily first neighbours). A similar modelling technique is employed at the cellular scale with the interaction radius, Rc, which we use to account for cell-cell adhesion. These radii are not necessarily equal (although they are of the same order of magnitude) since different types of interactions are considered at each scale. We perform our numerical simulations over time periods that are commensurate with the duration of in vitro experiments (hours) and for which proliferation is negligible [2, 7, 54]. For these reasons the model focuses only on the coordinated migration of ECs and assumes that proliferation occurs only at the vascular plexus [54], the initial vascular bed. This effect is implemented by introducing ECs into the domain at a specific set of boundary voxels, I VP (see Fig 2A). PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1008055 January 7, 2021 9 / 44 PLOS COMPUTATIONAL BIOLOGY A multiscale model of complex endothelial cell dynamics in early angiogenesis Fig 3. Illustrations for the subcellular scale of the model. (A) Since we track the position of cells only up to the positions of their nuclei, interactions are assumed to be non-local within some interaction radius. The subcellular and cellular scale interaction radii are denoted by Rs and Rc, respectively. (B) Illustration of values of the weight coefficients, αij (see Eq (5)). (C) Schematic illustration of the VEGF-Delta-Notch signalling pathway. NICD stands for Notch intracellular domain, VEGFR2 for VEGF Receptor 2 and VEGFR2� for activated, i.e. bound to VEGF, VEGFR2. In this case, the local environment of Cell 2 is Cell 1, and Dext and Next correspond, respectively, to the Delta and Notch concentrations in Cell 1 (and vice versa). Circled numbers correspond to the kinetic reactions listed in (D). (D) Kinetic reactions used for the VEGF-Delta-Notch pathway. See Table 1 for description of the model variables. HS(var) indicates that the transition rate of gene expression of a protein is transcriptionally regulated by the signalling variable, var. Here, HS(�) is the shifted Hill function (see caption of Table 2). Simple arrows indicate reactions with constant rates. (E) Bifurcation diagram of Notch concentration, N, as a function of external Delta ligand, Dext, PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1008055 January 7, 2021 10 / 44 PLOS COMPUTATIONAL BIOLOGY A multiscale model of complex endothelial cell dynamics in early angiogenesis corresponding to the system of equations Eq (1). Full lines denote stable steady states; dashed lines—unstable steady state; yellow filled dots— saddle-node bifurcation points. (F) Phenotype diagram as a function of external Delta, Dext, and external VEGF, V, corresponding to the system of equations Eq (1). (G) In simulations, the local environment of a cell (Delta and Notch levels in a neighbourhood of the cell, Dext and Next (Eq (6))) dynamically changes with time due to cell migration. This leads to phenotype switches. https://doi.org/10.1371/journal.pcbi.1008055.g003 Since the experimental data we use for model calibration and validation were extracted from experiments carried out with constant VEGF concentration supplied externally, we assume that the distribution of VEGF is maintained at a (prescribed) constant value at all times. Interested readers will find more detail on model formulation later in this section and in S1 Text. We list the variables of our multi-scale model in Table 1. Table 1. The description of the model variables. The variables in the table are organised by spatial scales. Subcellular variables, used to simulate the VEGF-Delta-Notch signalling pathway, define the gene expression pattern of individual cells (their phenotypes). Cellular scale variables define the occupancy of the lattice by ECs. Tissue scale variables define the composition and structure of the ECM. Bold letters denote vector variables specifying variable configuration for the whole lattice; normal font letters correspond to the variables associated with a particular voxel, the index of which is specified by a subscript. Variable Description Subcellular N ¼ ðN1 ; . . . ; NNI Þ The distribution of Notch receptor among voxels. If there is no cell scale nucleus in a voxel vi, i.e. Ei = 0, then Ni = 0. D ¼ ðD1 ; . . . ; DNI Þ The distribution of Delta ligand among voxels. If there is no cell nucleus in a voxel vi, i.e. Ei = 0, then Di = 0. I ¼ ðI1 ; . . . ; INI Þ The distribution of Notch intracellular domain (NICD) among voxels. If there is no cell nucleus in a voxel vi, i.e. Ei = 0, then Ii = 0. R2 ¼ ðR21 ; . . . ; R2NI Þ The distribution of VEGFR2 receptor among voxels. If there is no cell nucleus in a voxel vi, i.e. Ei = 0, then R2i = 0. R2� ¼ ðR2�1 ; . . . ; R2�NI Þ The distribution of activated VEGFR2 receptor (VEGFR2 bound to extracellular VEGF) among voxels. If there is no cell nucleus in a voxel vi, i.e. Ei = 0, then R2�i ¼ 0. Cellular scale E ¼ ðE1 ; . . . ; ENI Þ The distribution of EC nuclei among voxels. Ei = 1 if a cell nucleus is present in the voxel vi, Ei = 0, otherwise. At most one cell nucleus is allowed per voxel. EN ¼ ðEN1 ; . . . ; ENNI Þ The neighbourhood nucleus distribution. ThisP variable is completely �ij Ej jvj \BRc ðiÞj defined by the configuration of E. Each ENi ¼ Pj6¼i � , where �ij ¼ jvj j , j6¼i ij BRc ðiÞ is a circular neighbourhood of interaction radius, Rc, centred in the voxel vi, and j�j denotes area. Tissue scale c ¼ ðc1 ; . . . ; cNI Þ The ECM density, consisting mostly of collagen I and elastin fibers. It is degraded by cells via a process termed ECM proteolysis (Fig 1A III.). We assume 0 � ci � cmax for all i 2 I, where cmax > 0 is a parameter characterising the maximum density of the ECM. m ¼ ðm1 ; . . . ; mNI Þ The concentration of basal lamina components (collagen IV, fibronectin and various laminins) newly deposited by cells that are used for BM assembly (Fig 1A IV.). We assume 0 � mi � 1 for all i 2 I; mi = 0 if no BM components have been deposited yet, whereas mi = 1 if a BM has been assembled around the sprout segment situated in voxel vi. l ¼ ðl1 ; . . . ; lNI Þ The orientation landscape (OL) variable representing the alignment of ECM fibrils (Fig 1A V.) within the voxel. For a hexagonal lattice li ¼ flis gs2S for all possible jumping directions s 2 S ¼ fr; ur; ul; l; dl; drg (r— right, ur—upward-right, ul—upward-left, l—left, dl—downward-left, dr— downward-right), s 2 R2 (see Fig 2). An example of possible orientation landscape configurations is shown in Fig 5A. https://doi.org/10.1371/journal.pcbi.1008055.t001 PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1008055 January 7, 2021 11 / 44 PLOS COMPUTATIONAL BIOLOGY A multiscale model of complex endothelial cell dynamics in early angiogenesis Subcellular scale. VEGF-Delta-Notch model of phenotype selection. At the subcellular scale we account for the VEGF-Delta-Notch signalling pathway which determines the gene expression pattern (phenotype) of each EC (see Fig 3C). This pathway mediates inter-cellular cross-talk and typically produces alternating patterns of tip and stalk cell phenotypes within growing sprouts [11]. In our model, the pathway is simulated via a bistable stochastic system which accounts for intrinsic noise and, in particular, noise-induced random transitions between the stable steady states of the system (phenotypes) [69–72]. We posit that such pheno- typic switches are essential for understanding the complex dynamics of ECs within growing sprouts [2]. Our subcellular model is based on previous work [73–76]. Following [75, 76], we combine the lateral inhibition model of the Delta-Notch signalling pathway introduced in [73, 74] with the VEGF signalling pathway. The Delta-Notch model accounts for cis-inhibition when a Delta ligand and Notch receptor from the same cell inhibit each others’ activity. We include this interaction since cis-inhibition has been shown to substantially speed up pheno- type specification [77]. Individual cell system. The kinetic reactions acting on individual cells system are illus- trated in Fig 3D. We account for trans-activation of Notch receptor (production of an NICD) when it trans-binds to a Delta ligand belonging to a neighbouring cell, Dext (reaction 1a ). If a Delta ligand trans-binds to a Notch receptor on a neighbouring cell, Next, it is either endocyto- tically recycled or degraded (reaction 1b ). In this reaction, we assume that the active Notch signal is produced in the neighbouring cell, the dynamics of which are irrelevant for the cell of interest. Once cleaved from the Notch receptor, active Notch signal, NICD, is translocated to the cell nucleus where it down-regulates gene expression of VEGFR2 (reaction ②) and up-reg- ulates gene expression of the Notch receptor (reaction ③). Cis-inhibition is accounted for in reaction ④ in which mutual inhibition is assumed for a Delta ligand and a Notch receptor interacting within the same cell. External VEGF, V, can bind to and activate a VEGFR2 (reac- tion ⑤). This leads to up-regulation of Delta production (reaction ⑥). Reaction ⑦ corre- sponds to degradation of NICD. An essential feature of our subcellular model is that it exhibits bistability. To demonstrate this, we derived the mean-field limit equations associated with the kinetic reactions shown in Fig 3D (see Eq (1) in Table 2) and performed a numerical bifurcation analysis (see S1 Table for a list of the parameter values). The results presented in Fig 3E show how the steady state value of the Notch concentration, N, changes as the concentration of external Delta ligand, Dext, var- ies. For small (large) values of Dext, the system supports a unique steady-state corresponding to the tip (stalk) phenotype. For intermediate values of Dext, the system is bistable: both pheno- types coexist. The combined effect of the external VEGF, V, and Dext on the system is shown in Fig 3F. We see that varying V does not alter the qualitative behaviour shown in Fig 3E, although the size of the bistable region decreases as V decreases. The stable states of the subcellular system are characterised by distinct gene expression pat- terns (for example, high (low) Delta level, D, corresponds to tip (stalk) phenotype). Thus, one variable suffices in order to effectively define the phenotypes. We use Delta level, D, as a proxy variable in the following way tip phenotype : D � bD ; ð3Þ stalk phenotype : otherwise; where bD is baseline gene expression of Delta ligand in ECs (see S1 Table). Multicellular system. In order to account for cell-cell cross-talk via the VEGF-Delta- Notch pathway, we extend the individual cell system (see Fig 3D) to a multicellular environment PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1008055 January 7, 2021 12 / 44 Table 2. Mean-field equations associated with the stochastic system of the VEGF-Delta-Notch signalling pathway for the individual cell system (left column) and the multicellular system (right � 1þlX;Y ðX =X0 ÞnY column). Here HS(�) is the so-called shifted Hill function [78]. Its functional form is given by H S ðXÞ ¼ H S X; X0 ; lX;Y ; nY ¼ 1þð X = ÞnY X , where X0, λX,Y and nY are positive parameters (see S1 Text 0 for more details). Description and values of parameters can be found in S1 Table. Individual cell system (1) Multicellular system (2) dN dNi ¼ bN H S ðI; I0 ; lI;N ; nN Þ gN kt Dext N kc ND; ¼ bN H S ðIi ; I0 ; lI;N ; nN Þ gNi kt Ni D�i kc Ni Di ; dt dt dD dDi ¼ bD H S ðR2� ; R2�0 ; lR2� ;D ; nD Þ gD Zkt Next D kc ND; ¼ bD H S ðR2�i ; R2�0 ; lR2� ;D ; nD Þ gDi Zkt Di N�i kc Ni Di ; dt dt dI dIi ¼ kt Dext N ge I; ¼ kt Ni D�i ge Ii ; dt dt dR2 dR2i ¼ bR2 H S ðI; I0 ; lI;R2 ; nR2 Þ gR2 kv VR2; ¼ bR2 H S ðIi ; I0 ; lI;R2 ; nR2 Þ gR2i kv R2i V; dt dt dR2� dR2�i ¼ kv VR2 ge R2� ; ¼ kv R2i V ge R2�i ; dt dt PLOS COMPUTATIONAL BIOLOGY Dext ; Next are constant input parameters: Dext ¼ D�i ; Next ¼ N�i are given by ð6Þ; i 2 I: https://doi.org/10.1371/journal.pcbi.1008055.t002 PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1008055 January 7, 2021 13 / 44 A multiscale model of complex endothelial cell dynamics in early angiogenesis PLOS COMPUTATIONAL BIOLOGY A multiscale model of complex endothelial cell dynamics in early angiogenesis by specifying for each cell the external (i.e. belonging to neighbouring cells) amount of Delta and Notch to which it is exposed. As mentioned above, since cell positions in our model are only known up to the position of their nuclei, we assume non-local interactions between cells within a reaction radius, Rs. Thus, we define the local environment of a cell whose nucleus is sit- uated in voxel vi as the set of voxels with a non-zero overlap region with a circular neighbour- hood of radius Rs centred at voxel vi, BRs ðiÞ, HðiÞ ≔ fvj : vj \ BRs ðiÞ 6¼ ;; j 6¼ i; j 2 Ig: ð4Þ The weights, αij, assigned to each voxel vj 2 H(i), (see Fig 3B) are defined as follows jvj \ BRs ðiÞj aij ¼ ; i; j 2 I; ð5Þ jvj j where j�j denotes 2D area. The external Delta (Notch) concentration, Dext (Next), for a cell situated in a voxel vi is defined as follows: X aij Dj vj 2 HðiÞ Dext ¼ Di ¼ X ; aij vj 2 HðiÞ X ð6Þ aij Nj vj 2 HðiÞ Next ¼ Ni ¼ X ; aij vj 2 HðiÞ where Dj (Nj) denotes the Delta (Notch) concentration in voxel vj. Therefore, our multicellular stochastic system at the subcellular scale consists of the same kinetic reactions as in Fig 3D formulated for each voxel vi, i 2 I with Dext, Next given by (6). It is important to note that, in simulations of angiogenic sprouting, the quantities Dext and Next dynamically change due to cell migration, thus leading to phenotype re-establishment (see Fig 3G). When simulated within a two-dimensional domain, our multicellular system produces an alternating pattern of tip/stalk phenotypes. As the reaction radius, Rs, changes the system dynamics do not change but the proportion of tip and stalk cells does. To illustrate this, we ran simulations of the stochastic multicellular system for a 10 × 12 regular monolayer of cells for different values of the interaction radius, Rs (see S1 Fig). These simulations revealed that the distance between tip cells increases (i.e. the proportion of tip cells decreases) as Rs increases. We also investigated how phenotype patterning changes within a monolayer as cis-inhibition intensity varies (see S1 Text). For low values of the cis-inhibition parameter, κc, typical chess- board tip-stalk pattern is produced; as κc increases, ECs with tip phenotype can become adjacent to each other, thus increasing the time to patterning, since the lateral inhibition is weakened (see S1 Text). The mean-field equations associated with the multicellular stochastic system are given by Eq (2) in Table 2. More detail on the derivation and analysis of the subcellular model can be found in S1 Text. Cellular scale. Persistent random walk and cell overtaking. At the cellular scale we account for EC migration and overtaking. In more detail, we consider a compartment-based PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1008055 January 7, 2021 14 / 44 PLOS COMPUTATIONAL BIOLOGY A multiscale model of complex endothelial cell dynamics in early angiogenesis Fig 4. A cartoon illustrating cell migration. Migration transition rate (Eq (7)) is illustrated for a cell coloured in green. Its nucleus is marked with a green star. The motility of this cell depends on the ECM density, ci, in the voxel where its nucleus is situated (accounted for by the S(�) function (Eq (8)). The green cell forms cell-cell adhesions with other cells coloured in red. In this work, we assume a circular neighbourhood for cell-cell interactions within the so-called interaction radius, Rc, for the cellular scale (here it is drawn to have the same value as in our simulations, Rc = 1.5h, where h is the voxel width). This allows cells to interact beyond their immediate neighbours as, for example, the cell with the nucleus marked by a red star and the focal green cell. Other geometries (e.g. elliptical neighbourhood aligned with cell polarity vector) are unlikely to alter the model behaviour significantly since, in sprouting structures, lateral regions of the circular neighbourhood, which would be ignored by an elliptic neighbourhood, are typically empty. Cell-cell adhesion is accounted for by the neighbourhood function, F(�), Eq (9). The individual cell polarity, ϕ, is sampled from the von Mises distribution (Eq (12)) with the mean value given by the mean polarization direction, ~ p (calculated as a function of local ECM fibril alignment, li, Eq (13)). The distribution spread, κ, is assumed to depend on the focal cell phenotype and the concentration of the BM components, Eq (14). https://doi.org/10.1371/journal.pcbi.1008055.g004 model of a persistent random walk (PRW) [79], in which transition rates depend on the phe- notypic state of individual cells and their interactions with ECM components. We denote by s ¼ h 1 ðqj qi Þ 2 S the migration direction, as a unit vector pointing from qi towards qj, the centres of neighbouring voxels vi and vj, respectively (see Fig 2B). In order to formulate the PRW of ECs, we introduce, ω(i ! j), to denote the probability that a cell moves from voxel vi to a neighbouring voxel vj (i.e. along the direction s) where cell cell-cell occupancy adhesion cell polarity ∙Z 2¼ # ¢ !(i ! j) = f! Ei S (ci ) F E N ½ij (Di ) D s W (Á)fvM (Ájμ; ∙ )dÁ : ð7Þ i 0 diffusion ECM overtaking coefficient density probability We explain below each of the terms that appears in Eq (7) (see also Fig 4 for a graphical illustration and Table 1 for a description of the model variables). Diffusion coefficient. Transition rates in the framework of a PRW must be appropriately scaled with the size of a voxel (depending on the particular lattice used to discretise the domain). For a uniform hexagonal lattice, L, the diffusion coefficient of the transition ω(i ! j) is scaled as Dfo ¼ Do =h2 [80]. Here h[μm] is the width of a hexagonal voxel (see Fig 2) and h 2i Do mmmin is the macroscale diffusion coefficient for the ECs. Cell occupancy. If a cell nucleus is present in voxel vi, then Ei = 1 (at most one cell nucleus is allowed per voxel). If the voxel vi is empty, then Ei = 0 and ω(i ! j) = 0. ECM density. The function S(ci) accounts for the effect of the local ECM density, ci, on cell motility. In general, S(�) is a decreasing function of its argument. We assume further that ECs cannot move if the ECM concentration exceeds a threshold value, cmax (S(ci) = 0 for PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1008055 January 7, 2021 15 / 44 PLOS COMPUTATIONAL BIOLOGY A multiscale model of complex endothelial cell dynamics in early angiogenesis ci � cmax). In our simulations we fix ( ci 1 cmax if 0 � ci < cmax Sðci Þ ¼ ð8Þ 0; otherwise: Cell-cell adhesion. The effect of cell-cell adhesion on cell migration is incorporated via the so-called neighbourhood function, FðENi Þ. Its argument, ENi , represents the number of ECs in a cell’s local neighbourhood (red-coloured cells in Fig 4). The functional form of the neigh- bourhood function was chosen in order to phenomenologically capture the way in which EC behaviour depends on the cell-cell contacts (see Fig 5E). In biological experiments, it was shown that when a cell loses contact with its neighbours (laser ablation experiments in [8]) it halts until the following cells reach it. This is captured by the increasing part of FðENi Þ; when the number of neighbours around a cell is below the first threshold, EF1, the probability of cell movement decreases rapidly to zero, thus the migration transition (Eq (7)) goes to 0 as well. Similarly, when there are many ECs in a cell’s neighbourhood, its movement slows down. This is accounted for by the decreasing part of FðENi Þ when the number of neighbouring cells exceeds the second threshold value, EF2, cell movement is slowed down and eventually halts in regions of high cell density. � �þ N � 1 1 F Ei ¼ þ 1 ; ð9Þ 1 þ exp ð sF1 ðENi EF1 ÞÞ 1 þ exp ðsF2 ðENi EF2 ÞÞ where (x)+ = max(0, x), and the parameters EF1, EF2, sF1 and sF2 characterise the shape of the curve. Overtaking probability. This term accounts for cell overtaking and excluded volumes. A jump occurs with overtaking probability, ρij(Di) = 1, if the target voxel, vj, is empty (Ej = 0). Otherwise, the cells in voxels vi and vj switch their positions with probability ρij(Di) = pswitch(Di). We consider this probability to be phenotype dependent. In particular, we assume that tip cells (see Eq (3)) are more motile because their filopodia are stronger (see Fig 1A II.). Thus, the switching probability, pswitch(Di), is assumed to be an increasing function of the Delta level, Di, of the migrating cell ( 1; if Ej ¼ 0; rij ðDi Þ ¼ ð10Þ pswitch ðDi Þ; otherwise: where pmax pswitch ðDi Þ ¼ ; ð11Þ 1 þ exp ð sp ðDi Dp ÞÞ the parameters sp and Dp characterise, respectively, the slope and position of the sigmoid, and pmax denotes its maximum value (see Fig 5F). Cell polarity. Prior to migration, cells develop a polarity which depends on their local environment. Following [81, 82], we consider the local cell polarity, ϕ, (see Fig 4) to be a ran- dom quantity sampled from the von Mises distribution. The probability density function (pdf) of the von Mises distribution reads exp ðk cos ð� mÞÞ fvM ð�jm; kÞ ¼ : ð12Þ 2pI0 ðkÞ PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1008055 January 7, 2021 16 / 44 PLOS COMPUTATIONAL BIOLOGY A multiscale model of complex endothelial cell dynamics in early angiogenesis Fig 5. Series of sketches illustrating the components of the cellular scale persistent random walk. (A) An example of the orientation landscape, l, configurations for a hexagonal lattice. 1. All fibrils are aligned in the upward-right direction; this would be an example of a strongly aligned part of the ECM. 2. Half of the fibrils are aligned to the right, half in the upward-right direction. This would correspond to a branching point. 3. As in 2. but with some additional fibrils aligned in the left direction. (B) The window function, Ws(ϕ), (Eq (15)) has been defined as an indicator function over an angle interval corresponding to each possible migration direction s 2 S (lattice-dependent). The diagram illustrates these intervals for a hexagonal lattice. (C) Illustration of the probability distribution function of the von Mises distribution, fvM(x|μ, κ), centred at μ = 0 for different fixed κ (Eq (12)). (D) Illustration of the κ function as a function of local Delta ligand level, Di, and concentration of BM components, mi, (Eq (14) with K = 13.6, km = 2.2 and kD = 0.0002). (E) An example of the neighbourhood function, FðENi Þ (Eq (9), with EF1 = 0.15, EF2 = 0.6, sF1 = 30, sF2 = 10). (F) A sketch showing how the switching probability, pswitch(Di), changes with the level of Delta in voxel vi, Di (Eq (10) with pmax = 0.26, sp = 0.0015, Dp = 1500). https://doi.org/10.1371/journal.pcbi.1008055.g005 PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1008055 January 7, 2021 17 / 44 PLOS COMPUTATIONAL BIOLOGY A multiscale model of complex endothelial cell dynamics in early angiogenesis Its shape is characterised by two parameters: the mean value, μ, and the distribution spread, κ (see Fig 5C for a sketch of this pdf for different values of κ). I0 is the modified Bessel function of the first kind of order 0. In the context of EC migration, we view μ as the mean polarisation angle, and κ as cell exploratoriness. We integrate ECM structure and composition and cell phenotype into our cell migration model by assuming that μ and κ depend on these quantities. In particular, we assume that the mean polarisation angle, μ, depends on the ECM fibril alignment which is rep- resented in our model by the orientation landscape variable, li (see Table 1). From a biological point of view, this is substantiated by experimental observations of cells forming focal adhe- sions with ECM fibrils and, consequently, aligning along them [43]. We introduce the mean polarisation direction vector, ~p 2 R2 , and compute μ as its principle argument (see Fig 4) !T X X dir dir ~ p¼ Ha;n ðl Þdirx ; i Ha;n ðl Þdiry i ; dir2S dir2S ð13Þ m ¼ Argð~ p Þ: In Eq (13) the summation is taken over all possible directions for movement, dir ¼ T ðdirx ; diry Þ 2 S (in a hexagonal lattice there are at most 6 possible directions, see Fig 5A). The Hill function, Ha,n(�), is used to reflect the natural saturation limit to alignment and deforma- tion of the ECM fibrils, with a and n being fixed positive parameters. Details about how fibril orientation (orientation landscape, li) is calculated are given below. Cell exploratoriness, κ � 0, is directly related to the EC phenotype (see Fig 1A I.). If κ � 0, then the effect of polarity on migration is weak, and cells can explore many directions (this behaviour is typical of exploratory tip cells, see Eq (3)). By contrast, when κ � 1, the von Mises pdf is concentrated around μ (this behaviour is characteristic of stalk cells, see Eq (3)). To account for such phenotype-dependent behaviour, we assume κ to be a monotonic decreas- ing function of Di, the Delta level of the migrating cell. Similarly, increased concentration of the BM components deposited by ECs (see Fig 1A IV.) reduces the exploratory capability of both tip and stalk cells. We therefore propose κ to be an increasing function of the local con- centration of the BM components, mi. Combining these effects, we arrive at the following functional form for κ = κ(Di, mi): k ¼ kðDi ; mi Þ ¼ K exp ðkm mi kD Di Þ: ð14Þ Here K, km, kD are positive parameters (see Fig 5D for a sketch of κ(Di, mi)). Since our model of cell migration is formulated on a lattice, transition rates of type ω(i ! j) (jumps from vi into vj, in the direction s 2 S) are associated with an angle interval ½�smin ; �smax � (see Fig 2B). This corresponds to the angle between the vectors connecting the centre of the voxel vi with the endpoints of the voxel edge shared by vi and vj. For example, in a hexagonal lattice, the right direction, r 2 S, is associated with an interval ½�rmin ; �rmax � ¼ ½ p =6 ; p =6 �. Therefore, we restrict the von Mises pdf, fvM(�), in Eq (7) to this interval by multiplying it by a corresponding indicator function, Ws(ϕ), ( 1; if � þ 2pk 2 ½�smin ; �smax �; k 2 Z; s W ð�Þ ¼ ð15Þ 0; otherwise: This function is 2π-periodic due to the fact that its argument is an angle. We refer to Ws(ϕ) as the window function (see Fig 5B). PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1008055 January 7, 2021 18 / 44 PLOS COMPUTATIONAL BIOLOGY A multiscale model of complex endothelial cell dynamics in early angiogenesis Tissue scale. Modelling cell environment: ECM structure and composition. We account for the alignment of ECM fibrils, l, density of the ECM, c, and concentration of the basal lamina components, m (see Table 1) in the following way. Local alignment of the ECM fibrils, l, serves as a scaffold for the orientated migration of ECs [43]. Thus, we refer to l as the orientation landscape. Traction forces exerted by migrating ECs realign the ECM fibrils, so that they move closer to the growing sprouts. Furthermore, since tip cells have more filopodia than stalk cells, they exert a greater influence on the orienta- tion landscape [39] (see Fig 1A V.). We account for phenotype-dependent ECM realignment by assuming active stretching and accumulation of the fibrils upon cell movement between the voxels i ! j (in the direction s 2 S), i.e. when a transition of type ω(i ! j) occurs, lis ¼ lis þ Dl Di ; ð16Þ ljs ¼ ljs þ Dl Di : Here, the parameter Δl > 0, which quantifies the linear response of ECM fibrils to cell migration, depends on the substrate stiffness. Besides active stretching induced by cell locomotion, we also consider passive relaxation of the orientation landscape. We assume that relaxation follows a simple elastic model so that the orientation landscape decays exponentially at a constant rate, ηl, lis ðt þ tÞ ¼ lis ðtÞ exp ð Zl tÞ; ð17Þ where τ is the waiting time of the occurred migration transition and the update is done for all voxels i 2 I and all directions s 2 S. The time evolution of the ECM density, c, and BM components, m, is modelled via local ODEs. We assume phenotype-dependent ECM proteolysis induced by ECs (see Fig 1A III.). Since tip cells exhibit higher proteolytic activity than stalk cells [83, 84], we assume that the ECM at voxel vi is degraded at rate ηc(Di), which is an increasing function of its argument ( dci Zc ðDi Þ; if ci > 0; ¼ ð18Þ dt 0; otherwise; Zmax Zc ðDi Þ ¼ : ð19Þ 1 þ exp ð sc ðDi Dc ÞÞ Here, in order to account for the natural saturation in EC proteolytic ability, we assume a sigmoidal functional form for ηc(Di) with positive parameters Dc and sc (which correspond to the threshold level of Delta for initiation of ECM proteolysis and sharpness of EC response, respectively) and maximum value ηmax. Similarly, BM assembly (i.e. deposition of basal lamina components), is assumed to be phe- notype-dependent (see Fig 1A IV.). Tip cells are known to secrete BM components and to recruit and activate pericytes which secrete basal lamina components around the sprout [83, 84]. Thus we assume that the rate of secretion of BM components, γm(Di), is an increasing function of Di ( dmi gm ðDi Þ; if mi < 1:0; ¼ ð20Þ dt 0; otherwise: PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1008055 January 7, 2021 19 / 44 PLOS COMPUTATIONAL BIOLOGY A multiscale model of complex endothelial cell dynamics in early angiogenesis gmax gm ðDi Þ ¼ ; ð21Þ 1 þ exp ð sm ðDi Dm ÞÞ where, again, a sigmoidal functional form is assumed for γm(Di) with positive parameters Dm and sm (which correspond to the threshold level of Delta for initiation of BM assembly and sharpness of EC response, respectively) and maximum value γmax. Here we assume the decay of BM components is negligible on the timescale of our simulations. Description of quantitative metrics In order to calibrate and compare our simulation results to experimental data (see [2, 8]), we computed a number of metrics defined below. Displacement. The displacement statistic is associated with the average distance travelled by a cell in 15 minutes [8]. Orientation. The orientation statistic measures the average persistence of cells as they move. It is computed as the average (over ECs in all performed realisations) quantity of ratios of the length of a smoothed trajectory to the actual trajectory travelled by a single cell during simulation [2]. Directionality. This metric measures the average proportions of cells moving in the direc- tion of sprout elongation (anterograde), cells moving in the direction opposite to sprout elon- gation (retrograde), and cells that do not move during 20 minutes (still) [2]. Tip cell proportion. This metric is computed as the ratio of cells characterised by tip phe- notype (see Eq (3)) to the total number of cells in the system at a given time point. Mixing measure. This metric is motivated by the experimental observation that the trajec- tories of individual ECs, which initially form clusters with their immediate neighbours (see Fig 6 for an illustration), at later times diverge so that cells can find themselves at distant regions of the angiogenic vascular network [2, 8]. This metric is introduced to quantify the cell rear- rangements. Cell rearrangements are a key driver of sprout elongation during the early stages of vasculature formation and, as such, are directly related to network growth patterns. Later in this work, we will show how the mixing measure varies for different patterns of vascular net- work formation. Briefly, during simulations, we assign a label to each EC, ι, and record its position in the system at time t, pði; tÞ 2 I. We specify a set of voxels that form a cluster of nearest neighbours in the lattice, I cluster . At the end of the simulation, using recorded cell trajectories, we compute how far away (in a pair-wise fashion) cells that were situated at the voxels in I cluster at time t have moved away from each other during time interval of duration tm. Normalising this quan- tity by the cardinal of I cluster , j I cluster j, and the maximum possible travel distance in the system, dmax, we obtain the mixing measure at time t, MðtÞ, defined as 1 X � � MðtÞ ¼ dði1 ; i2 ; t þ tm Þ dði1 ; i2 ; tÞ : j I cluster j dmax i1 ; i2 such that ð22Þ pði1 ; tÞ ¼ i; pði2 ; tÞ ¼ j i; j 2 I cluster ; i 6¼ j Here d(ι1, ι2, t) is the distance between cells with labels, ι1 and ι2, at time t. It is computed as a distance in a manifold of the simulated vascular network (see Fig 6). This is due to the fact that ECs do not migrate randomly but rather within ECM-free vascular guidance tunnels of the generated vascular network. PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1008055 January 7, 2021 20 / 44 PLOS COMPUTATIONAL BIOLOGY A multiscale model of complex endothelial cell dynamics in early angiogenesis Fig 6. An illustration of the mixing measure. Two cells, labelled by ι1 and ι2, are located at the positions corresponding to the set I cluster at time t. We track their trajectories in the simulated vascular network (dashed black lines) during time tm. The mixing measure is defined as the difference between the distances between these cells at times t and t + tm, d(ι1, ι2, t) and d(ι1, ι2, t + tm), respectively, normalised by the number of cells considered and the maximum distance it is possible to travel in the simulated network. The distance function is defined as a distance within the simulated network (see S1 Text for details). https://doi.org/10.1371/journal.pcbi.1008055.g006 The distance, dmax, in Eq (22), is the maximum distance in the simulated vascular network, defined as follows dmax ¼ max dði1 ; i2 ; Tmax Þ; ð23Þ i1 ; i2 where Tmax is the final simulation time. Detailed descriptions of all the metrics and computational algorithms that we used are given in S1 Text. Results Emergent qualitative features: Branching and VEGF sensitivity Our model exhibits two characteristic features of functional angiogenic structures, namely, branching and chemotactic behaviour. A novel aspect of our multiscale model is that these fea- tures are emergent properties of its dynamics rather than being hardwired into the model. Specifically, branching is a direct consequence of the phenotype-dependent polarity of indi- vidual cells. When a stalk cell within a sprout undergoes a phenotype switch and assumes a tip identity, its exploratoriness, increases (i.e. κ, decreases). This enables the cell to develop a polarity angle that departs from the mean elongation direction of the sprout, μ. As a conse- quence, a new branch forms. In our model, new branches are typically initiated by cells exhib- iting the tip cell phenotype (see S1 Movie). This behaviour is characteristic of ECs observed in biological experiments [7, 9]. Figs 7 and 8 illustrate the branching phenomenon and stabilisa- tion of the network structure (due to accumulation of BM components deposited by ECs) in single realisations of numerical simulations of the model for uniform VEGF distribution at concentrations of 5 and 50 ng/ml (see also S1 and S4 Movies). Note that in all our simulations PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1008055 January 7, 2021 21 / 44 PLOS COMPUTATIONAL BIOLOGY A multiscale model of complex endothelial cell dynamics in early angiogenesis Fig 7. An example of an individual vascular network generated during simulation of our model with uniform VEGF = 5 ng/ml. (A) Temporal evolution of a simulated vascular network captured at 0 mins, 400 mins and 800 mins in an individual realization of our model for uniform VEGF = 5 ng/ml. The colour bar indicates level of Delta, D, (green colour corresponds to tip cells, red—to stalk cells). Arrows indicate the configuration of the orientation landscape, l. (B) The final configuration of the simulated vascular network at 1250 mins (corresponds to the final simulation time, t = Tmax = 2.5). The leftmost panel shows the concentration of Delta, D. Higher concentration (green colour) corresponds to tip cell phenotype, low concentration (red colour)—to stalk. On this plot, arrows correspond to the orientation landscape configuration, l. The central panel shows the final concentration of the ECM, c. The rightmost panel shows the final distribution of the polarity angle, μ, variable. Numerical simulation was performed using Setup 1 from S4 Table. Parameter values are listed S1 and S2 Tables for subcellular and cellular/tissue scales, respectively. For a movie of the numerical simulation, see S4 Movie, left panel. https://doi.org/10.1371/journal.pcbi.1008055.g007 “gaps” within sprouts can arise since we track the positions of cell nuclei and do not account for their true spatial extent. Chemotactic sensitivity in our model is a direct consequence of cell interactions with the ECM. In biological experiments, proteolytic activity of cells was observed to increase as expres- sion levels of Delta rise [83, 84]. In our model, increased levels of extracellular VEGF, V, up- regulate subcellular levels of Delta. As a result, an EC’s ability to degrade the ECM and invade it at a faster rate is enhanced where VEGF levels are high. This can be seen by comparing net- works generated at different uniform VEGF concentrations (see Figs 7 and 8). The network generated at VEGF = 5 ng/ml is small (Fig 7B, leftmost plot), and the vascular guidance tunnels created via proteolysis (middle plot) are not fully formed. By contrast, the simulated network for VEGF = 50 ng/ml (Fig 8B, leftmost plot) has a greater spatial extent, since collagen-free vas- cular tunnels (Fig 8B, middle plot) facilitate cell migration within them, increasing sprouting and cell persistence (see ECM density term in Eq (7)). Our model also exhibits the brush-border phenomenon [11, 85]. We performed a numerical simulation experiment of sprouting initialised from an initial vessel placed in a matrix with PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1008055 January 7, 2021 22 / 44 PLOS COMPUTATIONAL BIOLOGY A multiscale model of complex endothelial cell dynamics in early angiogenesis Fig 8. An example of an individual vascular network generated during simulation of our model with uniform VEGF = 50 ng/ml. (A) Temporal evolution of a simulated vascular network captured at 0 mins, 400 mins and 800 mins in an individual realization of our model for uniform VEGF = 50 ng/ml. The colour bar indicates level of Delta, D, (green colour corresponds to tip cells, red—to stalk cells). Arrows indicate the configuration of the orientation landscape, l. (B) The final configuration of the simulated vascular network at 1250 mins (corresponds to the final simulation time, t = Tmax = 2.5). The leftmost panel shows the concentration of Delta, D. Higher concentration (green colour) corresponds to tip cell phenotype, low concentration (red colour)—to stalk. On this plot, arrows correspond to the orientation landscape configuration, l. The central panel shows the final concentration of the ECM, c. The rightmost panel shows the final distribution of the polarity angle, μ, variable. Numerical simulation was performed using Setup 1 from S4 Table. Parameter values are listed in S1 and S2 Tables for subcellular and cellular/tissue scales, respectively. For a movie of the numerical simulation, see S1 Movie. https://doi.org/10.1371/journal.pcbi.1008055.g008 linearly increasing VEGF gradient. Fig 9 shows the evolution of the network at different times. The brush-border effect is evident at later times and characterised by increased cell numbers and branches in the top regions of the domain where VEGF levels are high. Model calibration Having established that our model exhibits the essential features of branching and chemotactic behaviour observed in experiments, we next compared our simulations with experimental results from [2, 8, 38]. This enabled us to estimate baseline parameter values for processes at the cellular and tissue scales (estimated values of parameters associated with processes acting at the subcellular scale are taken from previous works, see S1 Table). We ran 100 model simulations for uniform VEGF concentrations of 0, 5 and 50 ng/ml (Setup 1, see S4 Table, with final simulation time, Tmax = 2.5, equivalent to 1250 mins, i.e. � 20.8 h) and computed three metrics, namely, displacement, orientation and directionality (see Materials and methods for definitions) as was done in [2, 8]. The results presented in PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1008055 January 7, 2021 23 / 44 PLOS COMPUTATIONAL BIOLOGY A multiscale model of complex endothelial cell dynamics in early angiogenesis Fig 9. Sprouting in static VEGF gradient. Snapshots at different times of a vascular network growing in linear VEGF gradient increasing from VEGF = 0 ng/ml to VEGF = 5 ng/ml. (A) Time = 0 mins, initial setup. (B) Time = 400 mins, new branches appear from the initial sprout mostly in the upper half of the domain (higher VEGF concentration). ECs with lower positions have lower Delta level. (C) Time = 800 mins, the effect of the VEGF gradient can be seen clearly. (D) Time = 1250 mins, the final configuration of the simulated V-shaped (opening towards higher VEGF concentrations) network. Colour bar shows the level of Delta, D. Numerical simulation was performed using Setup 2 from S4 Table with final simulation time, Tmax = 2.5. Parameter values are listed in S1 and S2 Tables for subcellular and cellular/tissue scales, respectively. https://doi.org/10.1371/journal.pcbi.1008055.g009 Fig 10 show that our simulations reproduce general trends of angiogenic sprouting reported in [2, 8]. Specifically, regarding the displacement statistic (Fig 10A), agreement is very good, except for the inconsistency at displacement = 0 μm, i.e. cells that did not move during the considered time interval (15 minutes). This inconsistency arises because we do not include the vascular bed from which cells migrate (we account for it via a boundary condition; see Fig 2A and S1 Appendix). By contrast, in [8], cell displacements from the embedded aortic ring assay were included in the sample (these ECs are mostly quiescent). Similarly, results regarding the orientation statistic (Fig 10B and 10C) are in good agreement with the experimental results. In PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1008055 January 7, 2021 24 / 44 PLOS COMPUTATIONAL BIOLOGY A multiscale model of complex endothelial cell dynamics in early angiogenesis Fig 10. Statistics extracted from simulations of our model. (A) Histograms of cell displacements during a 15 minute time period for (1.) VEGF = 0 ng/ml, (2.) VEGF = 5 ng/ml and (3.) VEGF = 50 ng/ml. Black histograms correspond to the experimental data taken from the Supplementary Material of [8], red lines correspond to displacement curves for each VEGF concentration extracted from our model simulations. (B) A cartoon illustrating the orientation statistic, which is defined as a ratio between the net trajectory and the actual trajectory of a cell during simulation. (C) Box plots of the orientation statistic extracted from model simulations with VEGF = 0, 5, 50 ng/ml. Red crosses indicate box plot outliers. Orientation statistics obtained from experimental data from [2] are shown by blue stars on each box plot. (D) A cartoon illustrating the directionality statistic. (E) The directionality statistics for model simulations with VEGF = 0, 5, 50 ng/ml. (F) The directionality statistics extracted from experimental data in [2]. Numerical simulations were performed using Setup 1 from S4 Table and Tmax = 2.5. Parameter values are listed in S1 and S2 Tables for subcellular and cellular/tissue scales, respectively. All statistics were computed for 100 realisations. https://doi.org/10.1371/journal.pcbi.1008055.g010 PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1008055 January 7, 2021 25 / 44 PLOS COMPUTATIONAL BIOLOGY A multiscale model of complex endothelial cell dynamics in early angiogenesis particular, in our model ECs are more oriented, i.e. more persistent, in higher VEGF concen- trations, which is a feature also observed in [2]. Concerning the directionality statistic (Fig 10D–10F), we note that when VEGF = 0 ng/ml the numbers of anterograde and retrograde cells in the experiments [2] (Fig 10F) and numerical simulations (Fig 10E) are approximately equal. In this scenario, ECM proteolysis is slow and cell migration is mostly constrained to existing sprouts. Thus, any anterograde movement is an overtaking event in which the over- taken cell has to perform a retrograde displacement. As the VEGF concentration increases, ECM proteolysis (see Eqs (18) and (19)) increases and more cells at the leading edge of sprouts can invade the surrounding ECM, elongating the sprouts. This leads to an increase in the ratio of anterograde to retrograde moving cells with the VEGF concentration. Further evidence of agreement with experimental data is found by performing numerical simulations imitating the experimental setup of [38]. We performed simulations of sprouting from a cell bead embedded into the ECM (see Setup 3 from S4 Table) with varying collagen density (which corresponds to the initial ECM concentration, cmax, in our model) and a static linear VEGF gradient. Results from single realisations of different values of cmax are presented in Fig 11 (see also S2 Movie). These results show free cell migration with no preferred direction for low cmax values, typical angiogenic morphology for intermediate values of cmax, and poorly elongating sprouts for higher values of cmax. These findings are consistent with the experimen- tal observations reported in [38]. Furthermore, we note that the “sweet spot” of ECM concen- tration is related to EC ability to form typical angiogenic sprouting structures rather than to their ability to invade the ECM which decreases as the ECM concentration, cmax, increases (see Fig 11). We note that the results presented in this section were generated using a fixed set of param- eter values, except for the concentration of VEGF, V, and the concentration of collagen, cmax. Henceforth, we use these values as baseline parameter values (see S2 Table). Further model validation In this section we validate our model by comparing its predictions with experimental results detailing the behaviour of certain VEGF receptor mutant cells (VEGFR2+/- and VEGFR1+/- mutants with halved gene expressions of VEGFR2 and VEGFR1, respectively), studied by Jakobsson et al. in [7] and described in S1 Appendix. In order to ascertain whether our model can quantitatively reproduce competition between cells of different lineages (wild type (WT) and mutant cells) for the position of the leading cell in a sprout, we designed a series of numerical experiments which mimic the biological experi- ments reported in [7]. We start by simulating of EC competition within linear sprouts that are devoid of collagen matrix (Setup 4 in S4 Table), to ensure proteolysis-free random shuffling of cells within the sprout. We randomly initialise the sprout with cells of two chosen types with probability 50% (50% of WT cells and 50% of a specific mutant cell type) (see Fig 12, left col- umn). Cells are then allowed to shuffle within the sprout, overtaking each other. For each reali- sation we record the total amount of time for which WT and mutant cells occupy the position of the leading cell. As a control, we ran simulations in which two identical cell lines with parameters corre- sponding to WT lineage were mixed in a 1:1 ratio. As expected, the contribution of each WT cell to the leading cell position was approximately 50% (see Table 3). We then performed competition simulations in which WT cells were mixed with the differ- ent mutant cell lines in a 1:1 ratio (i.e. mixing 50% of WT cells with 50% of mutant cells). We repeated these numerical experiments in the presence of DAPT inhibitor which abolishes Notch signalling. The results of individual realisations presented in Fig 12 (see also S3 Movie), PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1008055 January 7, 2021 26 / 44 PLOS COMPUTATIONAL BIOLOGY A multiscale model of complex endothelial cell dynamics in early angiogenesis Fig 11. Cell migration from a cell bead in substrates of different collagen density. Final configurations of simulated vascular networks at time Tmax = 1.0 (corresponding to 500 minutes) of individual simulations used in reproducing the results of the polarisation experiment in [38]. Maximum collagen density (A) cmax = 0.1, (B) cmax = 1.0, (C) cmax = 1.7, (D) cmax = 3.0. The VEGF linear gradient starts with 0 ng/ml at y = 0 and increases up to 5 ng/ml at y = 125 μm. Central bead initial and basement membrane conditions, I BM ¼ I init , are outlined by a black thick line on each plot. Colour bars indicate Delta ligand concentration. Numerical simulations were performed using Setup 3 from S4 Table. Parameter values are listed in S1 and S2 Tables for subcellular and cellular/tissue scales, respectively. For a movie of the numerical simulation, see S2 Movie. https://doi.org/10.1371/journal.pcbi.1008055.g011 illustrate features of the different competition scenarios. For the WT:VEGFR2+/- scenario (see Fig 12A) when the mutant cells compete with WT cells, they almost never acquire the tip phe- notype. Consequently, they are rapidly overtaken by WT cells and accumulate far from the leading edge of the sprout (outlined in cyan on each plot). The leading cell positions are thus occupied predominately by WT cells. By contrast, in the WT:VEGFR1+/- scenario (Fig 12B), PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1008055 January 7, 2021 27 / 44 PLOS COMPUTATIONAL BIOLOGY A multiscale model of complex endothelial cell dynamics in early angiogenesis Fig 12. Initial and final configurations of single realisations of cells shuffling within a linear sprout when two given cell lines are mixed 1:1 (50% to 50%). Left column corresponds to the initial (random, 1:1) distribution of cells, right column—to the final one. The colour bar for Delta level of the WT goes from red colour (stalk cell) to green (tip cell), whereas for the mutant cells the bar goes from purple colour (stalk cell) to yellow (tip cell). (A) 50% of WT cells mixed with 50% of VEGFR2+/- mutant cells, no DAPT treatment. (B) 50% of WT cells mixed with 50% of VEGFR1+/- mutant cells, no DAPT treatment. (C) 50% of WT cells mixed with 50% of VEGFR2+/- mutant cells, both DAPT-treated. (D) 50% of WT cells mixed with 50% of VEGFR1+/- mutant cells, both DAPT-treated. Voxels corresponding to the leading edge of a sprout are outlined by thick cyan lines on each plot. Numerical simulations were performed using Setup 4 from S4 Table. Parameter values are listed in S1 and S2 Tables for subcellular and cellular/tissue scales, respectively, except for the changed parameters for the mutant cells listed in S1 Appendix. Final simulation time, Tmax = 50.0. For a movie of the numerical simulations, see S3 Movie. https://doi.org/10.1371/journal.pcbi.1008055.g012 Table 3. Contribution of WT cells to the leading cell position when mixed 1:1 (50%: 50%) with another type of cell (equivalent WT or specified mutant). Since as an initial setup of simulations we considered a sprout of width 2 (two voxels), there are two equivalent leading positions (Position 1 and Position 2) (outlined on each plot by cyan lines in Fig 12). The results are reported as mean value ± standard deviation for samples obtained from 100 realizations for each experimental scenario. Numerical setup is as specified in Fig 12. Experiment Position 1 Position 2 Ref. value [7] WT:WT 51.1±16.4% 53.3±20.6% 45.8% +/- WT:VEGFR2 93.6±7.1% 90.4±14.8% 87.0% WT:VEGFR1+/- 19.5±9.8% 20.5±17.8% 30.0% WT:VEGFR2+/- +DAPT 51.5±13.7% 52.9±17.8% 47.0% WT:VEGFR1+/- +DAPT 50.3±14.3% 47.6±16.7% 40.6% https://doi.org/10.1371/journal.pcbi.1008055.t003 PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1008055 January 7, 2021 28 / 44 PLOS COMPUTATIONAL BIOLOGY A multiscale model of complex endothelial cell dynamics in early angiogenesis mutant cells acquire the tip phenotype more often than WT cells, and thus contribute more significantly to the leading cell position. Treatment with DAPT forces all ECs (WT, VEGFR2+/- and VEGFR1+/-) to acquire the tip cell phenotype (which is the default when Notch signalling is abolished [12, 13]). Consequently, in these cases both treated mutants have a 50% likelihood of occupancy of the leading position (Fig 12C and 12D). We have also collected statistics from 100 realisations for each scenario and compared the results to the quantitative estimates provided in [7]. The results reported in Table 3 show that the contribution of WT cells to the leading cell position in each scenario is in good agreement with experimental values from [7]. In our simulations, VEGFR2+/- cells are less likely to stay at the leading position than WT cells: they occupy the leading position approximately 7% of the time. By contrast, VEGFR1+/- cells occupy the leading cell position approximately 78% of the time. DAPT restored the balance between the cells of different lineages so that they were on average equally mixed. Since the only parameters that we have modified are those used to mimic mutant cell gene expression and DAPT inhibition of Notch signalling (see S1 Appen- dix) our model provides possible explanation for overtaking dynamics of ECs in angiogenesis. Sensitivity analysis To ascertain how variation in the baseline parameter values of our model affects the behaviour of the system, we have performed an extensive sensitivity analysis. Since the subcellular VEGF-Delta-Notch model has already been calibrated and validated independently [76], we have focused our analysis on the cellular and tissue scale parameters (see S2 Table). Briefly, we have performed our analysis by fixing all the parameters except one at their baseline values, and then vary the focal parameter by ±0.1%, ±5%, ±10%, ±15%, and ±20%. This procedure is repeated for each of the tissue and cellular scale parameters. In order to quantify the impact of the variation of each parameter on both EC behaviour and network structure we have mea- sured the following quantities: • anterograde cell proportion (directionality metric); • orientation; • displacements; • number of branching points per 100 μm2 of vascular network area; • number of vessel segments; • vessel segment lengths. Each of these metrics has been measured over 100 realisations of the multiscale system. For a full account of the details, see S1 Text. The results of our sensitivity analysis are summarised in Fig 13 and S2 Fig. Our analysis shows that system behaviour is robust to variations in most of the model parameters consid- ered. This is indicated by the central cluster in Fig 13 and S2 Fig, highlighted in magenta, which represents those scenarios that exhibit very small deviations from the baseline behav- iour. By contrast, variations of a small number of parameters produce significant deviation from the baseline behaviour (see Table 4 and Fig 13 and S2 Fig). Specifically, we observe that an increase in Dm and a decrease in both Dc and K induce excessive branching, with shorter average vessel length (see Fig 13, hyper-branching region highlighted in brown). By contrast, a decrease in Dm and an increase in both Dc and K induce less branched networks, with longer average vessel length (see Fig 13, hypo-branching region highlighted in grey). These results are in agreement with well-known features of tumour PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1008055 January 7, 2021 29 / 44 PLOS COMPUTATIONAL BIOLOGY A multiscale model of complex endothelial cell dynamics in early angiogenesis Fig 13. Sensitivity analysis: Branching vs. vessel elongation. The results of the sensitivity analysis are represented as scatter plots of the mean number of branching points per 100 μm2 of vascular network area vs. mean vessel segment length, with colouring indicating mean (A) cell number; (B) vascular network area; (C) anterograde cell proportion; (D) orientation and (E) displacements. On these plots, dashed magenta lines indicate the point corresponding to the default parameter values (see S2 Table); magenta highlights the region of the main point clustering. The grey-coloured outlier region corresponds to vascular networks with few branching points and long vessel segments (hypo-branching), whereas the brown outlier region is characterised by short vessel segments and greater number of branching points (hyper- branching). Variations of the parameters that push the system towards one of the outlier regions are indicated on each plot. Panel (F) provides a general summary of these results. Simulation setup as in Setup 1, S3 Table, with Tmax = 2.5. The results are averaged over 100 realisations for each scenario. The subcellular parameters are fixed at their default values in all experiments (see S1 Table). https://doi.org/10.1371/journal.pcbi.1008055.g013 PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1008055 January 7, 2021 30 / 44 PLOS COMPUTATIONAL BIOLOGY A multiscale model of complex endothelial cell dynamics in early angiogenesis Table 4. Sensitivity analysis: Parameters producing significant deviation in system behaviour from the baseline scenario. # stands for decrease of the focal parameter, "—increase of the focal parameter. Par. Description Ref. equation Metrics affected Effect Dc Threshold of level of Delta for initiation of ECM proteolysis Eq (19) All Hyper-branching: Dc#; hypo-branching: Dc " Dm Threshold of level of Delta for initiation of BM assembly Eq (21) All Hypo-branching: Dm#; hyper-branching: Dm" K Cell exploratoriness, i.e. controls the variance of the von Mises Eq (14) All Hyper-branching: K#; hypo-branching: distribution K" EF1 Threshold of the level of cell-cell contact necessary for cell movement Eq (9) Orientation, Hyper-branching: EF1# initiation directionality EF2 Threshold of the inhibitory effect of crowding on cell movement Eq (9) Orientation, Hypo-branching: EF2# directionality https://doi.org/10.1371/journal.pcbi.1008055.t004 vasculature, where the tumour microenvironment inhibits vessel stabilisation; specifically, it hinders the formation of the basal membrane in tumour vasculature, yielding aberrant, exces- sively branched, networks [86]. This phenomenon can be realised by increasing Dm. Further- more, proteolysis is also up-regulated during tumour-induced angiogenesis due to secretion of MMPs by cancer cells. Proteolysis reduces the resistance experienced by the ECs as they migrate towards the tumour [86]. This effect can be accounted for phenomenologically by a reduction in Dc. Increased stimulation with growth factors that can also bind to VEGF recep- tors (as in pathological angiogenesis [86, 87]) can reduce the response of the ECs to chemotac- tic stimuli, due to high occupancy of receptors all over the cell membrane, thus shifting cell behaviour to chemokinesis (non-directional cell migration) [88, 89]. In our model, this transi- tion is controlled by the cell exploratoriness, κ (see Eq (14)): for high values of κ (i.e. higher values of K) cell migration is directed along the sprout elongation vector, whereas for small values of κ (i.e. smaller values of K) cells exhibit exploratory behaviour corresponding to chemokinesis. Our model thus predicts that changes in Dc, Dm and K are likely to occur in tumour- induced angiogenesis. This prediction is supported by current knowledge regarding the effects of the presence of a tumour in the microenvironment [86–88]. Model predictions: Network structure and cell mixing We also simulated the growth of vascular networks formed by a single mutant cell line (VEGFR2+/- or VEGFR1+/-) in the presence/absence of DAPT and compared our findings with the results from Jakobsson et al. who observed that in the absence of DAPT mutant cells mix with WT cells to form normal networks [7]. By contrast, the addition of DAPT leads to unstructured growth [7]. This is consistent with our simulation results (see S3–S6 Figs). Specif- ically, simulations with VEGFR2+/- mutant cells in the absence of DAPT (see S3 Fig and S4 Movie) suggest that the rate of network growth of VEGFR2+/- is slower than for their WT counterparts (see Figs 7 and 8 and S7 Fig). Since Delta levels in VEGFR2+/- cells are lower than in WT cells, they are less able to degrade and invade ECM, and deposit BM components than WT cells. This results in slower sprout elongation and increased branching (see S7 Fig). By contrast, VEGFR1+/- mutant cells possess higher levels of Delta and thus degrade the ECM more efficiently and invade the matrix more quickly than the WT cells (see S5 Fig and S4 Movie). Likewise, since the rate of segregation of BM components, γm, increases with Delta (Eq (21)), VEGFR1+/- cells are more persistent than WT cells. As a result, VEGFR1+/- ECs form less branched networks, with longer sprouts (see S7 Fig). Regarding networks grown with DAPT, since DAPT treatment abolishes all Notch signalling, all cells in the simulations PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1008055 January 7, 2021 31 / 44 PLOS COMPUTATIONAL BIOLOGY A multiscale model of complex endothelial cell dynamics in early angiogenesis Fig 14. Temporal evolution of mixing measure, tip cell proportion and branching structure in a simulated vascular network formed by WT cells. (A) The mixing measure, M(t), as a function of time (the mean value is indicated by a thick line and standard deviation is shown by a colour band). The results are averaged over 100 realisations. (B) Evolution of tip cell proportion as a function of time. The results are averaged over 100 realisations. (C) Snapshots from a single realization of our model simulating a vascular network formed by WT cells at 300, 600 and 1000 minutes. Colour bar indicates the level of Delta. The numerical simulation setup used is Setup 1 from S4 Table with final simulation time Tmax = 2.5. VEGF distribution was fixed uniformly at 5 ng/ml. Parameter values are listed in S1 and S2 Tables for subcellular and cellular/tissue scales, respectively. https://doi.org/10.1371/journal.pcbi.1008055.g014 with DAPT acquire the tip cell phenotype (S5 and S6 Figs), which produces unstructured growth. Since a cells’ ability to compete for the leading position, or, equivalently, cell shuffling, is altered in mutant cells, we sought to understand how cell rearrangements influence the struc- ture of a growing vascular network. To quantify cell rearrangements, we introduce a metric, which we refer to as mixing measure, MðtÞ (see Eq (22)). In Fig 14A, we plot the dynamics of the mixing measure, MðtÞ, obtained by averaging over 100 WT simulations. As the vascular network grows and new sprouts form (see Fig 14C), MðtÞ increases over time. The time evolution of the mixing measure varies for different cell lines. For VEGFR2+/- mutant cells it increases more slowly (Fig 15A), than for the WT cells (Fig 14A). VEGFR2+/- mixing arises more from branching than sprout elongation (see Fig 15C). A similar trend is seen for the VEGFR1+/- cells (compared to WT cells) (Fig 15A). However, contrary to VEGFR2+/-, this is due to high migration persistence of VEGFR1+/- ECs (see Fig 15D). A PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1008055 January 7, 2021 32 / 44 PLOS COMPUTATIONAL BIOLOGY A multiscale model of complex endothelial cell dynamics in early angiogenesis Fig 15. Temporal evolution of mixing measure, tip cell proportion and branching structure in simulated vascular networks formed by VEGFR2+/-and VEGFR1+/- mutant cells. (A) The mixing measure, M(t), as a function of time (the mean value is indicated by a thick line and standard deviation is shown by a band of the corresponding colour). The results are averaged over 100 realisations. (B) Evolution of tip cell proportion as a function of time. The results are averaged over 100 realisations. (C) Snapshots from a single realization of our model simulating a vascular network formed by VEGFR2+/- mutant cells at 300, 600 and 1000 minutes. Colour bar indicates the level of Delta. (D) Snapshots from a single realization of our model simulating a vascular network formed by VEGFR1+/- mutant cells at 300, 600 and 1000 minutes. Colour bar indicates the level of Delta. The numerical simulation setup used is Setup 1 from S4 Table with final simulation time Tmax = 2.5. VEGF distribution was fixed uniformly at 5 ng/ml. Parameter values are listed in S1 and S2 Tables for subcellular and cellular/tissue scales, respectively, except of those changed for mutant cells (see S1 Appendix). https://doi.org/10.1371/journal.pcbi.1008055.g015 PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1008055 January 7, 2021 33 / 44 PLOS COMPUTATIONAL BIOLOGY A multiscale model of complex endothelial cell dynamics in early angiogenesis slower increase in the mixing measure for mutant cells correlates with slower stabilisation of tip cell proportion around its steady state in vascular networks formed by mutant cells (see Figs 14B and 15B for WT and mutant cells, respectively). This supports our hypothesis that cell shuffling is directly related to the phenotypic specifications of ECs. Thus, we note that the cell rearrangement phenomenon is not a cell-autonomous decision but rather a result of con- tact-dependent EC cross-talk, which leads to EC phenotype specification. In all cases, regardless of the cell type and VEGF concentration, the mixing measure increases towards a steady state value (see S8 Fig, left panels). Furthermore, the steady state value of the mixing measure is independent of cell type (MðtÞ � 0:385 as t ! 1, see S9 Fig). This is a result of the fact that, as the vascular network reaches a sufficient size, cells perform proteolysis-free random shuffling within the manifold of the developed sprouts [52]. The rate of this proteolysis-free random shuffling does not depend on the cell line but rather on the tip- to-stalk ratio in the vasculature, which evolves to a steady state regardless of cell type and VEGF concentration (see S8 Fig, right panels). We conclude that the temporal evolution of the mixing measure characterises the resulting network structure to a larger degree than its steady state. Discussion and conclusions Angiogenic sprouting has been extensively investigated from both experimental [1–9] and the- oretical [16–33] perspectives. Integrating the variety of approaches and results from all these fields has allowed researchers to disentangle to a great extent mechanisms of complex EC behaviour such as coordinated EC migration [3, 7, 14] and EC-ECM interactions [33]. None- theless, as new experimental results emerge, biological assumptions of mathematical and computational models must be revisited, the models validated with new experimental data and used to test new hypotheses and generate predictions. As such, this multidisciplinary effort becomes a potential tool in designing biological experiments (e.g. for identifying new drug tar- gets [90, 91], finding new mechanisms for abnormal cell behaviour [92], etc.). Recently, it has become clear that cell rearrangements play a key role in driving vascular network growth in angiogenesis [1, 2, 7]. Defects in cell rearrangements lead to anastomosis failure and smaller vascular networks characterised by superimposed aberrant layers [1]. How- ever, the functional role of this phenomenon remains to be explained [1, 2, 7]. Quantifying cell rearrangements and relating them to vascular patterning for ECs with varying gene expres- sion patterns could be a starting point for understanding the functionality of cell mixing in angiogenesis. In this work, we developed a multiscale model which integrates individual cell gene expres- sion, EC migration and interaction with the local ECM environment. Our model exhibits char- acteristic EC behaviour, such as branching and chemotactic sensitivity, as emergent properties instead of being encoded via ad hoc rules (as has been traditionally done in the literature). The vascular networks generated by our model are capable of reproducing the general traits of sprouting angiogenesis: the networks exhibit branching patterns (see Figs 7 and 8), sprout elongation is enhanced in higher VEGF stimulation (see Fig 8) and the brush-border effect can be observed in networks grown in VEGF gradients (see Fig 9). Our simulation results are in good quantitative agreement with the characteristic trends of angiogenesis observed in experi- ments [2, 7, 8, 38] (see Figs 10 and 11 and Table 3). We then used our model to quantify the phenomenon of cell rearrangement. We defined and introduced a mixing measure, MðtÞ (see Fig 6), for networks formed by WT cells and VEGFR2+/- and VEGFR1+/- mutant cells with impaired gene expression of VEGFR2 and VEGFR1, respectively, used in [7] (see Figs 14 and 15). PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1008055 January 7, 2021 34 / 44 PLOS COMPUTATIONAL BIOLOGY A multiscale model of complex endothelial cell dynamics in early angiogenesis In all cases, in agreement with experimental observations, the mixing measure increases over time, although the specific details of its temporal evolution vary for different cell lines (see S8 Fig). In particular, for mutant cells, we find that mixing is lower due to either poor sprout elongation (VEGFR2+/- lineage, see Fig 15A and 15C and S7 Fig) or elevated cell persistency (VEGFR1+/- lineage, see Fig 15A and 15D and S7 Fig). WT cells form more functional net- works, in the sense of more effective coverage of the domain (and thus future delivery of oxy- gen/nutrients). This is achieved by a balance between branching and sprout elongation which increases the mixing measure for WT cells. We thus showed that the time evolution of the mix- ing measure is directly correlated to the generic features of the vascular pattern. This result supports the claim that shuffling and cell mixing are essential for network formation and structure. We also observe that the mixing measure reaches a steady state (see S9 Fig). We hypothesise that this is directly related to the proportion of tip cells in the network since they are the main driver of cell overtaking. This is substantiated by our results that, while the gene expression of VEGFR1+/- and VEGFR2+/- mutant cells exhibits variations, the steady state of the tip cell pro- portion is the same for all cell lineages (see S8 Fig, right column). Thus, although the branching pattern and effective sprout elongation vary for mutant cells, they generate adequate vascular networks in our simulations (see S3 and S4 Figs). We suggest that pathological network forma- tion is directly related to the imbalance in the tip cell proportion (for example, treatment with DAPT, which abolishes Notch signalling and forces all ECs to acquire tip phenotype, leads to hyper-sprouting [7], see S5 and S6 Figs). Analysing this proportion and what triggers its change might help to understand better what leads to malformations in sprouting angiogene- sis. Furthermore, the results of our sensitivity analysis suggest that variations in the parameters that control ECM remodelling (ECM proteolysis and BM assembly) and cell exploratoriness significantly modify vascular network structure (see Fig 13 and S2 Fig). This is in agreement with experimental evidence of aberrant vessels with excessive branching in tumour-induced angiogenesis [86–88]. We calibrated and validated our model against in vitro experiments carried out over a time scale on which cell proliferation and cell death are negligible (only � 5% of cells were undergo- ing mitosis in the observed time of � 22.4 hours [2]). Whilst cell turnover is neglected in the present study, it will need to be incorporated in any future study that simulates larger vascular networks. To do so, we must first reduce the computational complexity of the model since, in its current implementation, the runtime of a single simulation is up to several hours (computa- tional complexity increases with the number of cells in the system). One way to achieve this is to coarse-grain the subcellular model to a two-state (tip and stalk cell) Markovian system, omitting the dynamics of the intermediate variables [93]. We expect such an approach to reduce the computational complexity of the subcellular model, thus allowing us to run larger scale simulations. We also plan to scale up the model in order to obtain a continuum PDE limit description of it [81, 94–96]. This will further reduce the computational complexity of the simulations. Furthermore, a continuum PDE description of the model will allow us to per- form a more systematic model calibration and sensitivity analysis of the model parameters. To conclude, the model we have developed, with naturally emerging branching and EC che- motactic sensitivity to VEGF, allows us to investigate how changes in intracellular signalling and local cell environment influence cell mixing dynamics and to study the impact cell mixing has on the overall network structure. To our knowledge, it is the first attempt to quantify the cell mixing phenomenon in a theoretical model of angiogenesis. Since only individual cell tra- jectories are required for the computation of this statistic, it can also be extracted from experi- mental data. This, together with our predictions that cell mixing intensity is directly related to the vascular network structure, makes the mixing measure a potential marker for pathological PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1008055 January 7, 2021 35 / 44 PLOS COMPUTATIONAL BIOLOGY A multiscale model of complex endothelial cell dynamics in early angiogenesis angiogenesis. Furthermore, although we used a specific formulation for our subcellular model, the same modelling approach can be applied to a more/less detailed system (e.g. recent works [97, 98]) as long as it reproduces typical phenotype patterning of cells within the vasculature. This flexibility allows us to use our model to test various experimental hypotheses or make pre- dictions, for example, regarding pathological network formation. Supporting information S1 Appendix. Computational simulations. The appendix provides details on setups of our numerical simulations and brief discussion on their computational implementation. (PDF) S1 Fig. Examples of steady state patterns of the VEGFR-Delta-Notch subcellular model for different interaction radii. Final steady state patterns established during single stochastic sim- ulations of the system described by the kinetic reactions outlined in Fig 3D for a uniform hex- agonal lattice of 10 × 12 voxels. (A) Rs = 1.0h, (B) Rs = 2.0h, (C) Rs = 3.0h and the rest of the parameter values as in S1 Table. (TIF) S2 Fig. Sensitivity analysis: Orientation vs. anterograde cell proportion. We performed sim- ulations of our model by varying one of the parameters of the cellular and tissue scales at a time by a fixed per cent and keeping default values for the rest of the parameters (as in S2 Table). Each parameter was varied by ±0.1%, ±1%, ±5%, ±10%, ±15% and ±20%. For each numerical experiment, several quantitative metrics were computed. The results are repre- sented as scatter plots of mean cell trajectory orientation vs. mean anterograde cell proportion with colouring indicating mean (A) cell number; (B) vascular network area; (C) vessel segment length; (D) number of branching points per 100 μm2 of vascular network area and (E) dis- placements. On these plots, dashed magenta lines indicate the point corresponding to the default parameter values; magenta highlights the region of the main point clustering. The grey-coloured outlier region corresponds to vascular networks with less persistent, twisted ves- sels, whereas the brown outlier region is characterised by longer straight vessel segments. Vari- ations of the parameters that push the system towards one of the outlier regions are indicated on each plot. Panel (F) provides a general summary of these results. Simulation setup as in Setup 1, S3 Table, with Tmax = 2.5. The results are averaged over 100 realisations. The subcellu- lar parameters were fixed at their default values in all experiments (see S1 Table). (TIF) S3 Fig. Individual simulations of vascular networks generated by VEGFR2+/- mutant cells. Final configurations of simulated vascular networks of VEGFR2+/- mutant cells growing in uniform concentration of VEGF = 5 ng/ml, plot (A), and VEGF = 50 ng/ml, plot (B). The left- most panels show the amount of Delta, D. Higher values (green colour) correspond to tip cell phenotype, low values (red colour)—to stalk. On these plots arrows correspond to the orienta- tion landscape configuration, l. The central panels indicate the final concentration of the ECM, c. The rightmost panels—final distribution of the mean polarity angle, μ, variable. Numerical simulations were performed using Setup 1 from S4 Table and Tmax = 2.5. Parame- ter values are listed in S1 and S2 Tables for subcellular and cellular/tissue scales, respectively, except of those changed for VEGFR2+/- mutant cells (see S1 Appendix). (TIF) S4 Fig. Individual simulations of vascular networks generated by VEGFR1+/- mutant cells. Final configurations of simulated vascular networks of VEGFR1+/- mutant cells growing in PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1008055 January 7, 2021 36 / 44 PLOS COMPUTATIONAL BIOLOGY A multiscale model of complex endothelial cell dynamics in early angiogenesis uniform concentration of VEGF = 5 ng/ml, plot (A), and VEGF = 50 ng/ml, plot (B). The left- most panels show the amount of Delta, D. Higher values (green colour) correspond to tip cell phenotype, low values (red colour)—to stalk. On these plots arrows correspond to the orienta- tion landscape configuration, l. The central panels indicate the final concentration of the ECM, c. The rightmost panels—final distribution of the mean polarity angle, μ, variable. Numerical simulations were performed using Setup 1 from S4 Table and Tmax = 2.5. Parame- ter values are listed in S1 and S2 Tables for subcellular and cellular/tissue scales, respectively, except of those changed for VEGFR1+/- mutant cells (see S1 Appendix). (TIF) S5 Fig. Individual simulations of vascular networks generated by VEGFR2+/- mutant cells treated with DAPT. Final configurations of simulated vascular networks of VEGFR2+/- mutant cells treated with DAPT growing in uniform concentration of VEGF = 5 ng/ml, plot (A), and VEGF = 50 ng/ml, plot (B). The leftmost panels show the amount of Delta, D. Higher values (green colour) correspond to tip cell phenotype, low values (red colour)—to stalk. On these plots arrows correspond to the orientation landscape configuration, l. The central panels indicate the final concentration of the ECM, c. The rightmost panels—final distribution of the mean polarity angle, μ, variable. Numerical simulations were performed using Setup 1 from S4 Table and Tmax = 2.5. Parameter values are listed in S1 and S2 Tables for subcellular and cel- lular/tissue scales, respectively, except of those changed for VEGFR2+/- mutant cells and DAPT treatment (see S1 Appendix). (TIF) S6 Fig. Individual simulations of vascular networks generated by VEGFR1+/- mutant cells treated with DAPT. Final configurations of simulated vascular networks of VEGFR1+/- mutant cells treated with DAPT growing in uniform concentration of VEGF = 5 ng/ml, plot (A), and VEGF = 50 ng/ml, plot (B). The leftmost panels show the amount of Delta, D. Higher values (green colour) correspond to tip cell phenotype, low values (red colour)—to stalk. On these plots arrows correspond to the orientation landscape configuration, l. The central panels indicate the final concentration of the ECM, c. The rightmost panels—final distribution of the mean polarity angle, μ, variable. Numerical simulations were performed using Setup 1 from S4 Table and Tmax = 2.5. Parameter values are listed in S1 and S2 Tables for subcellular and cel- lular/tissue scales, respectively, except of those changed for VEGFR1+/- mutant cells and DAPT treatment (see S1 Appendix). (TIF) S7 Fig. Quantification of vascular network structure for WT and mutant cells at VEGF = 5 and 50 ng/ml. (A) Number of vessel segments. (B) Vessel segment length (μm). (C) Vascular network area (μm2) at the end of the numerical simulation. (D) Number of branching points per 100 μm2 of vascular network area. Details of definitions of these metrics can be found in S1 Text. In each boxplot, the central line indicates the median, and the horizontal edges of the box represent the 25th and 75th percentiles (for the bottom and top edges, respectively). The outliers are indicated by red cross symbols. Numerical simulation setup used is Setup 1 from S4 Table with final simulation time Tmax = 2.5. Parameter values are listed in S1 and S2 Tables for subcellular and cellular/tissue scales, respectively, except of those changed for mutant cells (see S1 Appendix). Results are averaged over 100 realisations for each experimental scenario. (TIF) S8 Fig. Temporal evolution of mixing measure and tip cell proportion in simulated vascu- lar networks. Left column plots show the mixing measure, M(t), as a function of time (the mean value is indicated by a thick line and standard deviation is shown by a band with PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1008055 January 7, 2021 37 / 44 PLOS COMPUTATIONAL BIOLOGY A multiscale model of complex endothelial cell dynamics in early angiogenesis corresponding colour). Right column plots demonstrate the evolution of tip cell proportion. Simulations were done for networks formed by (A) WT cells; (B) VEGFR2+/- mutant cells; and (C) VEGFR1+/- mutant cells. Numerical simulation setup used is Setup 1 from S4 Table with final simulation time Tmax = 2.5. Parameter values are listed in S1 and S2 Tables for sub- cellular and cellular/tissue scales, respectively, except of those changed for mutant cells (see S1 Appendix). Results are averaged over 100 realisations for each experimental scenario. (TIF) S9 Fig. Mixing measure steady state for VEGF = 0 ng/ml. Plots of mixing measure, M(t), as a function of time for WT, VEGFR2+/- and VEGFR1+/- mutant cells for VEGF = 0 ng/ml (the mean value is indicated by a thick line and standard deviation is shown by a band of the corre- sponding colour). At this concentration of external VEGF, there is no effective sprout elonga- tion, thus cells perform proteolysis-free random shuffling within already existing sprouts [52]. This leads to a steady state of the mixing measure for all cell lineages. Mean values are 0.39, 0.39 and 0.38 for WT, VEGFR2+/- and VEGFR1+/- cells, respectively. Numerical simulations were performed using Setup 1 from S4 Table and Tmax = 2.5. Parameter values are listed in S1 and S2 Tables for subcellular and cellular/tissue scales, respectively, except of those changed for mutant cells (see S1 Appendix). (TIF) S1 Table. Baseline parameter values for the VEGF-Delta-Notch subcellular model. Descrip- tion and reference values used in simulations of the subcellular VEGF-Delta-Notch signalling. (PDF) S2 Table. Parameter values of the cellular and tissue scales used in our simulations. (PDF) S3 Table. Initial conditions for numerical simulations. Here I is the set of all voxels; S is the set of all possible migration directions. DUnif[a, b] is a discrete uniform distribution over all integer numbers lying within the interval [a, b]; Unif[a, b] is the uniform distribution on the interval [a, b]. Baseline gene expression parameters for the VEGF-Delta-Notch signalling are listed in S1 Table. Δinit = 1.0 for all numerical simulations (this value, as, in general, for the value of the OL variable, is non-dimensional). The fluctuation parameter, ξ, is set to 0.1 in all numerical simulations. The exact values for cinit and minit are given for each numerical experi- ment in S4 Table, as well as the set of initial cell positions, I init . For the description of model variables see Table 1 in the main text. (PDF) S4 Table. Setups of simulation experiments. For each setup of numerical simulation we spec- ify the lattice dimensions, NIx and NIy ; the set of indices corresponding to the vascular plexus, I VP ; the initial cell nuclei positions, I init ; the initial polarisation direction, sinit; the initial ECM and BM concentrations, cinit and minit, respectively; the VEGF distribution over the lattice, V; and cell line used in simulations. (PDF) S1 Movie. An example of an individual vascular network generated by wild type ECs dur- ing simulation of our model with uniform VEGF = 50 ng/ml. The leftmost panel shows the concentration of Delta, D. The colour bar indicates level of Delta, D, (green colour corre- sponds to tip cells, red—to stalk cells). Arrows indicate the configuration of the orientation landscape, l. The central panel indicates the concentration of the ECM, c. The rightmost panel —the polarity angle, μ, variable. A circular colour bar indicates the value of μ. The simulation was performed using Setup 1 from S4 Table with final simulation time, Tmax = 2.5. Parameter PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1008055 January 7, 2021 38 / 44 PLOS COMPUTATIONAL BIOLOGY A multiscale model of complex endothelial cell dynamics in early angiogenesis values are listed in S1 and S2 Tables for subcellular and cellular/tissue scales, respectively. (MP4) S2 Movie. Cell migration from a cell bead in substrates of different collagen density. Single realisations of angiogenic sprouting from a cell bead in substrates of different collagen densi- ties (reproducing the results of the polarisation experiment in [38]). Maximum collagen den- sity (A) cmax = 0.1, (B) cmax = 1.0, (C) cmax = 1.7, (D) cmax = 3.0. The VEGF linear gradient starts with 0 ng/ml at y = 0 and increases up to 5 ng/ml at y = 125 μm. Central bead initial and basement membrane conditions, I BM ¼ I init , are outlined by a black thick line on each plot. Colour bars indicate the level of Delta ligand. The simulations were performed using Setup 3 from S4 Table. Parameter values are listed in S1 and S2 Tables for subcellular and cellular/tis- sue scales, respectively. (MP4) S3 Movie. Single realisations of cells shuffling within a linear sprout when two given cell lines are mixed 1:1 (50% to 50%). The cell lines used in each realization are indicated in the titles. In the top row, no treatment with DAPT inhibitor was applied to cells; in the bottom row, all ECs were treated with DAPT. The leading edge corresponds to two rightmost voxels of each sprout. The colour bar for Delta level of the WT goes from red colour (stalk cell) to green (tip cell), whereas for the mutant cells the bar goes from purple colour (stalk cell) to yellow (tip cell). The simulations were performed using Setup 4 from S4 Table. Parameter values are listed in S1 and S2 Tables for subcellular and cellular/tissue scales, respectively, except for the changed parameters for the mutant cells listed in S1 Appendix. Final simulation time, Tmax = 50.0. (MP4) S4 Movie. Examples of an individual vascular networks generated by wild type and mutant (VEGFR2+/- and VEGFR1+/-) ECs during simulation of our model with uniform VEGF = 5 ng/ml. The cell line is indicated in the title of each panel. The colour bar indicates level of Delta, D, (green colour corresponds to tip cells, red—to stalk cells). Arrows indicate the config- uration of the orientation landscape, l. Numerical simulation was performed using Setup 1 from S4 Table with final simulation time, Tmax = 2.5. Parameter values are listed in S1 and S2 Tables for subcellular and cellular/tissue scales, respectively, except for the changed parameters for the mutant cells listed in S1 Appendix. (MP4) S1 Text. Supplementary material. The file contains a more detailed description of the subcel- lular VEGF-Delta-Notch model, metric definitions and algorithms for their implementation, further specifications on simulation of our model and sensitivity analysis of model parameters. (PDF) Author Contributions Conceptualization: Daria Stepanova, Helen M. Byrne, Philip K. Maini, Tomás Alarcón. Investigation: Daria Stepanova, Helen M. Byrne, Philip K. Maini, Tomás Alarcón. Methodology: Daria Stepanova, Helen M. Byrne, Philip K. Maini, Tomás Alarcón. Supervision: Helen M. Byrne, Philip K. Maini, Tomás Alarcón. Writing – original draft: Daria Stepanova, Helen M. Byrne, Philip K. Maini, Tomás Alarcón. PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1008055 January 7, 2021 39 / 44 PLOS COMPUTATIONAL BIOLOGY A multiscale model of complex endothelial cell dynamics in early angiogenesis Writing – review & editing: Daria Stepanova, Helen M. Byrne, Philip K. Maini, Tomás Alarcón. References 1. Angulo-Urarte A, Casado P, Castillo SD, Kobialka P, Kotini MP, Figueiredo AM, et al. Endothelial cell rearrangements during vascular patterning require PI3-kinase-mediated inhibition of actomyosin con- tractility. Nature Communications. 2018; 9(1):4826. https://doi.org/10.1038/s41467-018-07172-3 PMID: 30446640 2. Arima S, Nishiyama K, Ko T, Arima Y, Hakozaki Y, Sugihara K, et al. Angiogenic morphogenesis driven by dynamic and heterogeneous collective endothelial cell movement. Development. 2011; 138 (21):4763–4776. https://doi.org/10.1242/dev.068023 PMID: 21965612 3. Bentley K, Franco CA, Philippides A, Blanco R, Dierkes M, Gebala V, et al. The role of differential VE- cadherin dynamics in cell rearrangement during angiogenesis. Nature Cell Biology. 2014; 16(309):309– 321. https://doi.org/10.1038/ncb2926 PMID: 24658686 4. De Bock K, Georgiadou M, Schoors S, Kuchnio A, Wong BW, Cantelmo AR, et al. Role of PFKFB3- driven glycolysis in vessel sprouting. Cell. 2013; 154(3):651–663. https://doi.org/10.1016/j.cell.2013.06. 037 PMID: 23911327 5. Cruys B, Wong BW, Kuchnio A, Verdegem D, Cantelmo AR, Conradi LC, et al. Glycolytic regulation of cell rearrangement in angiogenesis. Nature Communications. 2016; 7:12240. https://doi.org/10.1038/ ncomms12240 PMID: 27436424 6. Franco CA, Jones ML, Bernabeu MO, Geudens I, Mathivet T, Rosa A, et al. Dynamic endothelial cell rearrangements drive developmental vessel regression. PLoS Biology. 2015; 13(4). https://doi.org/10. 1371/journal.pbio.1002125 7. Jakobsson L, Franco CA, Bentley K, Collins RT, Ponsioen B, Aspalter IM, et al. Endothelial cells dynam- ically compete for the tip cell position during angiogenic sprouting. Nature Cell Biology. 2010; 12 (10):943. https://doi.org/10.1038/ncb2103 PMID: 20871601 8. Sugihara K, Nishiyama K, Fukuhara S, Uemura A, Arima S, Kobayashi R, et al. Autonomy and non- autonomy of angiogenic cell movements revealed by experiment-driven mathematical modeling. Cell Reports. 2015; 13(9):1814–1827. https://doi.org/10.1016/j.celrep.2015.10.051 PMID: 26655898 9. Ubezio B, Blanco RA, Geudens I, Stanchi F, Mathivet T, Jones ML, et al. Synchronization of endothelial Dll4-Notch dynamics switch blood vessels from branching to expansion. Elife. 2016; 5:e12167. https:// doi.org/10.7554/eLife.12167 PMID: 27074663 10. Vega R, Carretero M, Travasso RD, Bonilla LL. Notch signaling and taxis mechanims regulate early stage angiogenesis: A mathematical and computational model. PLoS computational biology. 2020; 16 (1):e1006919. https://doi.org/10.1371/journal.pcbi.1006919 PMID: 31986145 11. Blanco R, Gerhardt H. VEGF and Notch in tip and stalk cell selection. Cold Spring Harbor Perspectives in Medicine. 2013; 3(1):a006569. https://doi.org/10.1101/cshperspect.a006569 PMID: 23085847 12. Hellström M, Phng LK, Hofmann JJ, Wallgard E, Coultas L, Lindblom P, et al. Dll4 signalling through Notch1 regulates formation of tip cells during angiogenesis. Nature. 2007; 445(7129):776–780. https:// doi.org/10.1038/nature05571 PMID: 17259973 13. Phng LK, Gerhardt H. Angiogenesis: a team effort coordinated by notch. Developmental Cell. 2009; 16 (2):196–208. https://doi.org/10.1016/j.devcel.2009.01.015 PMID: 19217422 14. Bentley K, Mariggi G, Gerhardt H, Bates PA. Tipping the Balance: Robustness of Tip Cell Selection, Migration and Fusion in Angiogenesis. PLOS Computational Biology. 2009; 5:1–19. https://doi.org/10. 1371/journal.pcbi.1000549 PMID: 19876379 15. Chen W, Xia P, Wang H, Tu J, Liang X, Zhang X, et al. The endothelial tip-stalk cell selection and shuf- fling during angiogenesis. Journal of Cell Communication and Signaling. 2019; p. 1–11. https://doi.org/ 10.1007/s12079-019-00511-z PMID: 30903604 16. Balding D, McElwain DLS. A mathematical model of tumour-induced capillary growth. Journal of Theo- retical Biology. 1985; 114(1):53–73. https://doi.org/10.1016/S0022-5193(85)80255-1 PMID: 2409405 17. Byrne HM, Chaplain MAJ. Mathematical models for tumour angiogenesis: numerical simulations and nonlinear wave solutions. Bulletin of Mathematical Biology. 1995; 57(3):461–486. PMID: 7537141 18. Anderson ARA, Chaplain MAJ. Continuous and discrete mathematical models of tumor-induced angio- genesis. Bulletin of Mathematical Biology. 1998; 60:857–899. https://doi.org/10.1006/bulm.1998.0042 PMID: 9739618 19. Owen MR, Alarcón T, Byrne HM, Maini PK. Angiogenesis and vascular remodelling in normal and cancerous tissues. J Math Biol. 2009; 58:689–721. https://doi.org/10.1007/s00285-008-0213-z PMID: 18941752 PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1008055 January 7, 2021 40 / 44 PLOS COMPUTATIONAL BIOLOGY A multiscale model of complex endothelial cell dynamics in early angiogenesis 20. Macklin P, McDougall S, Anderson ARA, Chaplain MAJ, Cristini V, Lowengrub J. Multi-scale modelling and non-linear simulation of vascular tumour growth. J Math Biol. 2009; 58:765–798. https://doi.org/10. 1007/s00285-008-0216-9 PMID: 18781303 21. Bartha K, Rieger H. Vascular network remodeling via vessel cooption, regression and growth in tumors. Journal of Theoretical Biology. 2006; 241(4):903–918. https://doi.org/10.1016/j.jtbi.2006.01.022 PMID: 16545398 22. Chaplain MAJ, McDougall SR, Anderson ARA. Mathematical modelling of tumor-induced angiogenesis. Annual Review of Biomedical Engineering. 2006; 8(1):233–257. https://doi.org/10.1146/annurev. bioeng.8.061505.095807 PMID: 16834556 23. Jackson T, Zheng X. A Cell-based Model of Endothelial Cell Migration, Proliferation and Maturation Dur- ing Corneal Angiogenesis. Bulletin of Mathematical Biology. 2010; 72(4):830–868. https://doi.org/10. 1007/s11538-009-9471-1 PMID: 20052558 24. Lowengrub JS, Frieboes HB, Jin F, Chuang YL, Li X, Macklin P, et al. Non-linear modelling of cancer: bridging the gap between cells and tumours. Nonlinearity. 2010; 23:R1–R91. https://doi.org/10.1088/ 0951-7715/23/1/R01 PMID: 20808719 25. Perfahl H, Byrne HM, Chen T, Estrella V, Alarcón T, Lapin A, et al. Multiscale Modelling of Vascular Tumour Growth in 3D: The Roles of Domain Size and Boundary Conditions. PLoS One. 2011; 6: e14790. https://doi.org/10.1371/journal.pone.0014790 PMID: 21533234 26. Travasso RDM, Corvera Poiré E, Castro M, Rodriguez-Manzaneque JC, Hernandez-Machado A. Tumour angiogenesis and vascular patterning: a mathematical model. PLoS One. 2011; 6:e19989. https://doi.org/10.1371/journal.pone.0019989 PMID: 21637756 27. Spill F, Guerrero P, Alarcon T, Maini PK, Byrne HM. Mesoscopic and continuum modelling of angiogen- esis. Journal of Mathematical Biology. 2014; p. 1–48. https://doi.org/10.1007/s00285-014-0771-1 PMID: 24615007 28. Perfahl H, Hughes BD, Alarcón T, Maini PK, Lloyd MC, Reuss M, et al. 3D hybrid modelling of vascular network formation. Journal of Theoretical Biology. 2016; 414:254–268. https://doi.org/10.1016/j.jtbi. 2016.11.013 PMID: 27890575 29. Vilanova G, Colominas I, Gomez H. A mathematical model of tumour angiogenesis: growth, regression and regrowth. Journal of The Royal Society Interface. 2017; 14(126):20160918. https://doi.org/10. 1098/rsif.2016.0918 PMID: 28100829 30. Xu J, Vilanova G, Gomez H. A mathematical model coupling tumour growth and angiogenesis. PLoS One. 2016; 11:e0149422. https://doi.org/10.1371/journal.pone.0149422 PMID: 26891163 31. Pillay S, Byrne HM, Maini PK. Modeling angiogenesis: A discrete to continuum description. Phys Rev E. 2017; 95:012410. https://doi.org/10.1103/PhysRevE.95.012410 PMID: 28208423 32. Moreira-Soares M, Coimbra R, Rebelo L, Carvalho J, Travasso RDM. Angiogenic factors produced by hypoxic cells are a leading driver of anastomoses in sprouting angiogenesis—a computational study. Scientific Reports. 2018; 8:8726. https://doi.org/10.1038/s41598-018-27034-8 PMID: 29880828 33. van Oers RFM, Rens EG, LaValley DJ, Reinhart-King CA, Merks RMH. Mechanical Cell-Matrix Feed- back Explains Pairwise and Collective Endothelial Cell Behavior In Vitro. PLoS Computational Biology. 2014; 10(8):1–14. https://doi.org/10.1371/journal.pcbi.1003774 PMID: 25121971 34. Scianna M, Bell C, Preziosi L. A review of mathematical models for the formation of vascular networks. Journal of theoretical biology. 2013; 333:174–209. https://doi.org/10.1016/j.jtbi.2013.04.037 PMID: 23684907 35. Vilanova G, Colominas I, Gomez H. Computational modeling of tumor-induced angiogenesis. Archives of Computational Methods in Engineering. 2017; 24(4):1071–1102. https://doi.org/10.1007/s11831- 016-9199-7 36. Mantzaris NV, Webb S, Othmer HG. Mathematical modeling of tumor-induced angiogenesis. Journal of mathematical biology. 2004; 49(2):111–187. https://doi.org/10.1007/s00285-003-0262-2 PMID: 15293017 37. Boas SEM, Merks RMH. Tip cell overtaking occurs as a side effect of sprouting in computational models of angiogenesis. BMC Systems Biology. 2015; 9(1):86. https://doi.org/10.1186/s12918-015-0230-7 PMID: 26589386 38. Shamloo A, Heilshorn SC. Matrix density mediates polarization and lumen formation of endothelial sprouts in VEGF gradients. Lab on a Chip. 2010; 10(22):3061–3068. https://doi.org/10.1039/c005069e PMID: 20820484 39. Du Y, Herath SC, Wang Qg, Wang Da, Asada HH, Chen PC. Three-dimensional characterization of mechanical interactions between endothelial cells and extracellular matrix during angiogenic sprouting. Scientific Reports. 2016; 6:21362. https://doi.org/10.1038/srep21362 PMID: 26903154 PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1008055 January 7, 2021 41 / 44 PLOS COMPUTATIONAL BIOLOGY A multiscale model of complex endothelial cell dynamics in early angiogenesis 40. Kirkpatrick ND, Andreou S, Hoying JB, Utzinger U. Live imaging of collagen remodeling during angio- genesis. American Journal of Physiology-Heart and Circulatory Physiology. 2007; 292(6):H3198– H3206. https://doi.org/10.1152/ajpheart.01234.2006 PMID: 17307995 41. Lamalice L, Le Boeuf F, Huot J. Endothelial cell migration during angiogenesis. Circulation Research. 2007; 100(6):782–794. https://doi.org/10.1161/01.RES.0000259593.07661.1e PMID: 17395884 42. McCoy MG, Wei JM, Choi S, Goerger JP, Zipfel W, Fischbach C. Collagen Fiber Orientation Regulates 3D Vascular Network Formation and Alignment. ACS Biomaterials Science & Engineering. 2018; 4 (8):2967–2976. https://doi.org/10.1021/acsbiomaterials.8b00384 43. Sage EH, Vernon RB. Regulation of angiogenesis by extracellular matrix: the growth and the glue. Jour- nal of Hypertension Supplement: Official Journal of the International Society of Hypertension. 1994; 12 (10):S145–52. PMID: 7539492 44. Michaelis UR. Mechanisms of endothelial cell migration. Cellular and Molecular Life Sciences. 2014; 71 (21):4131–4148. https://doi.org/10.1007/s00018-014-1678-0 PMID: 25038776 45. Senger DR, Davis GE. Angiogenesis. Cold Spring Harbor Perspectives in Biology. 2011; 3(8):a005090. https://doi.org/10.1101/cshperspect.a005090 PMID: 21807843 46. Sottile J. Regulation of angiogenesis by extracellular matrix. Biochimica et Biophysica Acta (BBA)- Reviews on Cancer. 2004; 1654(1):13–22. https://doi.org/10.1016/j.bbcan.2003.07.002 PMID: 14984764 47. Myer NM, Myers KA. CLASP1 regulates endothelial cell branching morphology and directed migration. Biology Open. 2017; 6(10):1502–1515. https://doi.org/10.1242/bio.028571 PMID: 28860131 48. Potente M, Carmeliet P. The link between angiogenesis and endothelial metabolism. Annual review of physiology. 2017; 79:43–66. https://doi.org/10.1146/annurev-physiol-021115-105134 PMID: 27992732 49. Rousseau S, Houle F, Kotanides H, Witte L, Waltenberger J, Landry J, et al. Vascular endothelial growth factor (VEGF)-driven actin-based motility is mediated by VEGFR2 and requires concerted acti- vation of stress-activated protein kinase 2 (SAPK2/p38) and geldanamycin-sensitive phosphorylation of focal adhesion kinase. Journal of Biological Chemistry. 2000; 275(14):10661–10672. https://doi.org/10. 1074/jbc.275.14.10661 PMID: 10744763 50. Watanabe T, Noritake J, Kaibuchi K. Regulation of microtubules in cell migration. Trends in cell biology. 2005; 15(2):76–83. https://doi.org/10.1016/j.tcb.2004.12.006 PMID: 15695094 51. Reinhart-King CA, Dembo M, Hammer DA. The dynamics and mechanics of endothelial cell spreading. Biophysical Journal. 2005; 89(1):676–689. https://doi.org/10.1529/biophysj.104.054320 PMID: 15849250 52. Stratman AN, Saunders WB, Sacharidou A, Koh W, Fisher KE, Zawieja DC, et al. Endothelial cell lumen and vascular guidance tunnel formation requires MT1-MMP–dependent proteolysis in 3-dimen- sional collagen matrices. Blood. 2009; 114(2):237–247. https://doi.org/10.1182/blood-2008-12-196451 PMID: 19339693 53. Sopher RS, Tokash H, Natan S, Sharabi M, Shelah O, Tchaicheeyan O, et al. Nonlinear Elasticity of the ECM Fibers Facilitates Efficient Intercellular Communication. Biophysical Journal. 2018; 115(7):1357– 1370. https://doi.org/10.1016/j.bpj.2018.07.036 PMID: 30217380 54. Gerhardt H, Golding M, Fruttiger M, Ruhrberg C, Lundkvist A, Abramsson A, et al. VEGF guides angio- genic sprouting utilizing endothelial tip cell filopodia. The Journal of Cell Biology. 2003; 161(6):1163– 1177. https://doi.org/10.1083/jcb.200302047 PMID: 12810700 55. Stundzia AB, Lumsden CJ. Stochastic simulation of coupled reaction–diffusion processes. Journal of Computational Physics. 1996; 127(1):196–207. https://doi.org/10.1006/jcph.1996.0168 56. Howard M, Rutenberg AD. Pattern formation inside bacteria: fluctuations due to the low copy number of proteins. Physical Review Letters. 2003; 90(12):128102. https://doi.org/10.1103/PhysRevLett.90. 128102 PMID: 12688905 57. Moro E. Hybrid method for simulating front propagation in reaction-diffusion systems. Phys Rev E. 2004; 69:060101(R). https://doi.org/10.1103/PhysRevE.69.060101 PMID: 15244531 58. Bernstein D. Simulating mesoscopic reaction-diffusion systems using the Gillespie algorithm. Physical Review E. 2005; 71(4):041103. https://doi.org/10.1103/PhysRevE.71.041103 PMID: 15903653 59. Engblom S, Ferm L, Hellander A, Lötstedt P. Simulation of stochastic reaction-diffusion processes on unstructured meshes. SIAM Journal on Scientific Computing. 2009; 31(3):1774–1797. https://doi.org/ 10.1137/080721388 60. Flegg MB, Chapman SJ, Erban R. The two-regime method for optimizing stochastic reaction–diffusion simulations. Journal of the Royal Society Interface. 2012; 9(70):859–868. https://doi.org/10.1098/rsif. 2011.0574 PMID: 22012973 PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1008055 January 7, 2021 42 / 44 PLOS COMPUTATIONAL BIOLOGY A multiscale model of complex endothelial cell dynamics in early angiogenesis 61. Hellander A, Hellander S, Lotstedt P. Coupled mesoscopic and microscopic simulation of reaction-diffu- sion processes in mixed dimensions. Multiscale Model Sim. 2012; 10:585–611. https://doi.org/10.1137/ 110832148 62. Franz B, Flegg MB, Chapman SJ, Erban R. Multiscale reaction-diffusion algorithms: PDE-assisted Brownian dynamics. SIAM Journal on Applied Mathematics. 2013; 73(3):1224–1247. https://doi.org/10. 1137/120882469 63. Guerrero P, Byrne HM, Maini PK, Alarcón T. From invasion to latency: Intracellular noise and cell motil- ity as key controls of the competition between resource-limited cellular populations. J Math Biol. 2015; 72:123–156. https://doi.org/10.1007/s00285-015-0883-2 PMID: 25833187 64. Yates CA, Flegg MB. The pseudo-compartment method for coupling partial differential equation and compartment-based models of diffusion. Journal of The Royal Society Interface. 2015; 12(106). https:// doi.org/10.1098/rsif.2015.0141 PMID: 25904527 65. Gillespie DT, Petzold LR, Seitaridou E. Validity conditions for stochastic chemical kinetics in diffusion- limited systems. The Journal of Chemical Physics. 2014; 140(5):02B604_1. https://doi.org/10.1063/1. 4863990 PMID: 24511926 66. Isaacson SA. A convergent reaction-diffusion master equation. The Journal of Chemical Physics. 2013; 139(5):054101. https://doi.org/10.1063/1.4816377 PMID: 23927237 67. Elf J, Ehrenberg M. Spontaneous separation of bi-stable biochemical systems into spatial domains of opposite phases. Systems Biology. 2004; 1(2):230–236. https://doi.org/10.1049/sb:20045021 PMID: 17051695 68. Gillespie DT. A general method for numerically simulating the stochastic time evolution of coupled chemical reactions. Journal of Computational Physics. 1976; 22(4):403–434. https://doi.org/10.1016/ 0021-9991(76)90041-3 69. Gardiner CW. The escape time in nonpotential systems. J Stat Phys. 1983; 30:157–177. https://doi.org/ 10.1007/BF01010872 70. Maier RS, Stein DL. A scaling theory of bifurcations in the symmetric weak-noise escape problem. J Stat Phys. 1996; 83:291–357. https://doi.org/10.1007/BF02183736 71. Guantes R, Poyatos JF. Multistable Decision Switches for Flexible Control of Epigenetic Differentiation. PLoS Computational Biology. 2008; 4(11). https://doi.org/10.1371/journal.pcbi.1000235 PMID: 19043543 72. Perez-Carrasco R, Guerrero P, Briscoe J, Page KM. Intrinsic noise profoundly alters the dynamics and steady state of the morphogen controlled bistable genetic switches. PLoS Comp Bio. 2016; 12: e1005154. https://doi.org/10.1371/journal.pcbi.1005154 PMID: 27768683 73. Sprinzak D, Lakhanpal A, LeBon L, Santat LA, Fontes ME, Anderson GA, et al. Cis-interactions between Notch and Delta generate mutually exclusive signalling states. Nature. 2010; 465(7294):86. https://doi.org/10.1038/nature08959 PMID: 20418862 74. Sprinzak D, Lakhanpal A, LeBon L, Garcia-Ojalvo J, Elowitz MB. Mutual inactivation of Notch receptors and ligands facilitates developmental patterning. PLoS Computational Biology. 2011; 7(6):e1002069. https://doi.org/10.1371/journal.pcbi.1002069 PMID: 21695234 75. Boareto M, Jolly MK, Lu M, Onuchic JN, Clementi C, Ben-Jacob E. Jagged–Delta asymmetry in Notch signaling can give rise to a Sender/Receiver hybrid phenotype. Proceedings of the National Academy of Sciences. 2015; 112(5):E402–E409. https://doi.org/10.1073/pnas.1416287112 PMID: 25605936 76. Boareto M, Jolly MK, Ben-Jacob E, Onuchic JN. Jagged mediates differences in normal and tumor angiogenesis by affecting tip-stalk fate decision. Proceedings of the National Academy of Sciences. 2015; 112(29):E3836–E3844. https://doi.org/10.1073/pnas.1511814112 PMID: 26153421 77. Venkatraman L, Regan ER, Bentley K. Time to Decide? Dynamical Analysis Predicts Partial Tip/Stalk Patterning States Arise during Angiogenesis. PLoS One. 2016; 11:1–23. https://doi.org/10.1371/ journal.pone.0166489 PMID: 27846305 78. Lu M, Jolly MK, Levine H, Onuchic JN, Ben-Jacob E. MicroRNA-based regulation of epithelial–hybrid– mesenchymal fate determination. Proceedings of the National Academy of Sciences. 2013; 110 (45):18144–18149. https://doi.org/10.1073/pnas.1318192110 PMID: 24154725 79. Codling EA, Plank MJ, Benhamou S. Random walk models in biology. Journal of the Royal Society Interface. 2008; 5(25):813–834. https://doi.org/10.1098/rsif.2008.0014 PMID: 18426776 80. Jansen APJ. An introduction to kinetic Monte Carlo simulations of surface reactions. vol. 856. Springer; 2012. 81. Painter KJ, Hillen T. Navigating the flow: individual and continuum models for homing in flowing environ- ments. Journal of the Royal Society Interface. 2015; 12(112):20150647. https://doi.org/10.1098/rsif. 2015.0647 PMID: 26538557 PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1008055 January 7, 2021 43 / 44 PLOS COMPUTATIONAL BIOLOGY A multiscale model of complex endothelial cell dynamics in early angiogenesis 82. Hillen T, Painter KJ, Swan AC, Murtha AD. Moments of von Mises and Fisher distributions and applica- tions. Mathematical Biosciences & Engineering. 2017; 14:673. https://doi.org/10.3934/mbe.2017038 PMID: 28092958 83. del Toro R, Prahst C, Mathivet T, Siegfried G, Kaminker JS, Larrivee B, et al. Identification and func- tional analysis of endothelial tip cell–enriched genes. Blood. 2010; 116(19):4025–4033. https://doi.org/ 10.1182/blood-2010-02-270819 PMID: 20705756 84. Siemerink MJ, Klaassen I, Van Noorden CJ, Schlingemann RO. Endothelial tip cells in ocular angiogen- esis: potential target for anti-angiogenesis therapy. Journal of Histochemistry & Cytochemistry. 2013; 61(2):101–115. https://doi.org/10.1369/0022155412467635 PMID: 23092791 85. Muthukkaruppan V, Auerbach R. Angiogenesis in the mouse cornea. Science. 1979; 205(4413):1416– 1418. https://doi.org/10.1126/science.472760 PMID: 472760 86. Papetti M, Herman IM. Mechanisms of normal and tumor-derived angiogenesis. American Journal of Physiology-Cell Physiology. 2002; 282(5):C947–C970. https://doi.org/10.1152/ajpcell.00389.2001 PMID: 11940508 87. Gupta MK, Qin RY. Mechanism and its regulation of tumor-induced angiogenesis. World journal of gastroenterology: WJG. 2003; 9(6):1144. https://doi.org/10.3748/wjg.v9.i6.1144 PMID: 12800214 88. Plank M, Sleeman B. Tumour-induced angiogenesis: a review. Journal of Theoretical Medicine. 2003; 5 (3-4):137–153. https://doi.org/10.1080/10273360410001700843 89. Barkefors I, Le Jan S, Jakobsson L, Hejll E, Carlson G, Johansson H, et al. Endothelial cell migration in stable gradients of vascular endothelial growth factor a and fibroblast growth factor 2 effects on chemo- taxis and chemokinesis. Journal of Biological Chemistry. 2008; 283(20):13905–13912. https://doi.org/ 10.1074/jbc.M704917200 PMID: 18347025 90. Finley SD, Chu LH, Popel AS. Computational systems biology approaches to anti-angiogenic cancer therapeutics. Drug discovery today. 2015; 20(2):187–197. https://doi.org/10.1016/j.drudis.2014.09.026 PMID: 25286370 91. Tang L, Van De Ven AL, Guo D, Andasari V, Cristini V, Li KC, et al. Computational modeling of 3D tumor growth and angiogenesis for chemotherapy evaluation. PloS one. 2014; 9(1). https://doi.org/10. 1371/journal.pone.0083962 PMID: 24404145 92. Bentley K, Philippides A, Regan ER. Do endothelial cells dream of eclectic shape? Developmental cell. 2014; 29(2):146–158. https://doi.org/10.1016/j.devcel.2014.03.019 PMID: 24780735 93. Dykman MI, Horita T, Ross J. Statistical distribution and stochastic resonance in a periodically driven chemical system. The Journal of chemical physics. 1995; 103(3):966–972. https://doi.org/10.1063/1. 469796 94. Othmer HG, Hillen T. The diffusion limit of transport equations derived from velocity-jump processes. SIAM Journal on Applied Mathematics. 2000; 61(3):751–775. https://doi.org/10.1137/ S0036139999358167 95. Painter KJ, Hillen T. From Random Walks to Fully Anisotropic Diffusion Models for Cell and Animal Movement. In: Cell Movement. Springer; 2018. p. 103–141. 96. Plaza RG. Derivation of a bacterial nutrient-taxis system with doubly degenerate cross-diffusion as the parabolic limit of a velocity-jump process. Journal of mathematical biology. 2019; 78(6):1681–1711. https://doi.org/10.1007/s00285-018-1323-x PMID: 30603994 97. Antfolk D, Sjöqvist M, Cheng F, Isoniemi K, Duran CL, Rivero-Muller A, et al. Selective regulation of Notch ligands during angiogenesis is mediated by vimentin. Proceedings of the National Academy of Sciences. 2017; 114(23):E4574–E4581. https://doi.org/10.1073/pnas.1703057114 PMID: 28533359 98. Kang TY, Bocci F, Jolly MK, Levine H, Onuchic JN, Levchenko A. Pericytes enable effective angiogene- sis in the presence of proinflammatory signals. Proceedings of the National Academy of Sciences. 2019; 116(47):23551–23561. https://doi.org/10.1073/pnas.1913373116 PMID: 31685607 PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1008055 January 7, 2021 44 / 44