Falsifying computational models of angiogenesis through quantitative comparison with in vitro models

During angiogenesis, endothelial cells migrate from existing vessels, proliferate and collectively organize into new capillaries. In vitro and in vivo experimentation is instrumental for identifying the molecular players and cell behavior that regulate angiogenesis. Alongside experimental work, computational and mathematical models of angiogenesis have helped to show if the current molecular and cellular understanding of cell behavior is sufficient. As input, the model takes (a subset of) the current knowledge or hypotheses of single cell behavior and captures it into a dynamical, mathematical description. As output, it predicts the multicellular behavior following from the actions of many individual cells, e.g., the formation of a sprout or the formation of a vascular network. Paradoxically, computational modeling based on different assumptions, i.e., completely different, sometimes non-intersecting sets of observed single cell behavior, can reproduce the same angiogenesis-like multicellular behavior, making it practically impossible to decide which, if any, of these models is correct. Here we present dynamic analyses of time-lapses of in vitro angiogenesis experiments and compare these with dynamic analyses of mathematical models of angiogenesis. We extract a variety of dynamical characteristics of endothelial cell network formation using a custom time-lapse video analysis pipeline in ImageJ. We compare the dynamical network characteristics of the in vitro experiments to those of the cellular networks produced by computational models. We test the response of the in silico dynamic cell network characteristics to key model parameters and make related changes in the composition of the in vitro environment. We present comparisons with computational model outcomes and argue how models that fail to reproduce these trends can be rejected. All in all, we show how our dynamic approach helps to clarify key endothelial cell interactions required for angiogenesis, and how the approach helps analyze what key changes in network properties can be traced back to changes in individual cell behavior.


Introduction
Angiogenesis, the formation of new blood vessels from existing ones, is crucial for physiological and pathological mechanisms such as embryonic development, wound healing and tumor growth.Proangiogenic factors, including vascular endothelial growth factor A (VEGF-A), stimulate sprouting in endothelial cells (ECs), which make up the inner lining of blood vessels (1)(2)(3).The ECs proliferate and migrate into the extracellular environment and form a vascular network that supplies the surrounding tissue with nutrients and oxygen and removes waste products.A key regulatory factor determining the final structure of the new blood vessel network is the composition of the extracellular matrix (ECM), the scaffold of proteins that surrounds the blood vessels and provides structural support to the newly formed multicellular network (4).Specific integrin binding sites on ECM proteins like collagens, fibronectin and laminins allow cells to exert stresses on the ECM and sense the mechanical properties of their environment (5).Cells also secrete enzymes which digest the ECM, releasing ECMretained growth factors, including VEGF-A, which stimulate EC migration (6).Cell-cell interactions, like Notch signaling and vascular endothelium cadherin (VE-cadherin) adhesions, regulate individual cell behavior within the network, separating the cells into leading tip-and proliferating stalk cells (2).
Similarly, VE-cadherin inhibits VEGFR-2 signaling, making the cells less susceptible to VEGF, resulting in less angiogenic potential of the cells (8).The combination of extracellular, intercellular and intracellular regulation makes angiogenic sprouting a complex and multiscale process.An understanding of how each regulatory mechanism contributes to the final network formation is necessary for biomedical applications where angiogenesis needs to be stimulated, steered or reduced, like tissue engineering (9) or vascular normalization (10).
One of the more commonly used in vitro models to test the effect of environmental changes on angiogenesis is the tube formation assay (TFA) (11).In this assay ECs are placed on top of a basement membrane-like substrate, where the ECs, under the right conditions, spontaneously organize into networks of interconnected tubes.With this assay the influence of various chemical and mechanical stimuli on EC networks has been investigated by comparing the final network between conditions (12).In the analysis of TFA outcomes it is often assumed that a network that contains a higher density of branches indicates a higher angiogenic potential of the cells (i.e., more sprouting).However, the temporal development of the network over time may contain key information on the dynamic processes regulating angiogenesis.Merks et al. (2006), therefore, measured the density of branch points in developing vascular networks in vitro at regular intervals and observed that the network density dropped slowly over time (13).Parsa et al. (2011), analyzed time-lapse videos of in vitro angiogenesis in more detail, and described the development of endothelial networks as a carefully orchestrated sequence of five subevents: [1] rearrangement and aggregation, [2] spreading, [3] elongation and formation of cell-cell contacts, [4] plexus stabilization, and [5] plexus reorganization (14).The first two stages describe the attachment of the cells to the substrate directly after cell seeding.The rate and the extent of cell spreading on the substrate have been linked to substrate stiffness and composition (15) and are assumed to be influenced by the traction forces cells are able to exert on the substrate (16).Later stages in the tube formation assay rely on a combination of cell-cell and cell-ECM interactions.Rüdiger et al. (2020) observed a collapse of EC networks on soft, laminin matrices, confirming the importance of the ECM to provide support for the network to stabilize (4).
Computational modeling is a helpful tool for proposing and testing whether hypothetical developmental mechanisms are sufficient to explain biological observations (17).Different computational models of in vitro angiogenesis all consider angiogenesis as the assembly of ECs into a vascular-like network structure (18)(19)(20).For example, in Manoussaki et al. (1996), the in vitro matrix remodeling behavior of ECs on gelled matrix is described with a continuum model, in which ECs are described as local densities using a system of differential equations (21).Similarly, in Serini et al. (2003), ECs are modeled to migrate upwards a chemoattractant gradient with a continuum model (22).
In Palachanis et al. (2015), ECs are considered as identical 2D ellipses, using a single-particle model (23).Köhn-Luque et al., (2013) used a hybrid continuum and cell-based model to model the effect of VEGF retention in the ECM on EC network formation (6).In cell-based models of angiogenesis, ECs are described as discrete areas of variable shapes, where multicellular behavior arises from the response of single cells to inputs from their microenvironment (24,25).Similarly, Vega et al. (2020) used a hybrid continuum and cell-based model to model the interaction between Notch-and VEGF signaling and cell-ECM adhesion (26).
Interestingly, both continuum models, particle-based models, as well as cell-based models are successful in mimicking EC network formation, and different hypotheses, or model inputs, can in some cases lead to a qualitatively similar output.In this work we focus on three cell-based models as previously developed in our group (13,27,28).The "cell elongation model" assumes that ECs are attracted towards one another through autocrine/paracrine signaling by secreting a chemoattractant (13).In this model it was observed that cells self-organize into vascular-like structures if the cells have an elongated shape.Observations by Parsa et al. (2011) support the necessity for ECs to elongate to form networks, as cells elongate before they start forming a network and continue to elongate as the network develops (14).The "contact-inhibition model" also assumes that ECs are attracted to one another through a secreted chemoattractant, but here it was assumed that VE-cadherin binding mediates contact-inhibition of chemotaxis, following observations by Dejana et al. (2004) that VEcadherin-VEGFR2 complex forming limits cell proliferation (29) and that in the yolk sac of VEcadherin double-knockout mice and in in vitro mouse allantois explants ECs failed to form vascular networks and aggregated into isolated islands (28,30).In disagreement with VEGF guided cell-cell attraction-based models, Rüdiger et al. (2020) observed persistent network formation in the absence of a VEGF-A gradient and they argue that the in vitro network formation is driven by mechanical communication (4).The importance of cell-ECM interactions is further supported by observations by Stephanou et al. (2007), where they observed faster and more lacunae in network formation on fibrinogen gels of intermediate rigidity (31).Finally, the "mechanical cell-cell communication model" assumes that ECs are attracted to one another solely through mechanical interaction with the ECM (27).In this cell-ECM interaction model ECs induce strain in the ECM through contractile forces and migrate up the strain gradient (27).With these many models that successfully reproduce multicellular network formation based on different observed EC behavior, we need to critically exam these models to see which assumptions, if any are correct.
Thus, for a more complete understanding of the integration of different mechanisms driving angiogenesis, we must quantitatively compare these computational models to in vitro experiments by systematically changing model parameters and in vitro experimental setup and iteratively refining the models based on the outcomes of the comparison.As an initial step, we will critically examine three models originating from our own research group: [1] The cell elongation model ( 13); [2] the contact inhibition model (28); and [3] the mechanical cell-cell communication model (27).For this comparison, dynamical analyses of in vitro experiments are required.Here we present new dynamic time-lapse videos of tube formation assays and a new analysis pipeline based on ImageJ and demonstrate its use for model selection for network formation.There are commercial and open-source image analysis tools available for quantitative analyses of angiogenic networks (32)(33)(34)(35), however these tools are designed for single frame end-point analyses and therefore less suitable for highthroughput dynamical analyses of large quantities of data.Our novel image analysis pipeline is suitable to dynamically analyze and compare in vitro and simulated networks.The pipeline allows us to reliably extract network features over time, such as the branch density and length and the average size and number of lacunae.In this study we use this pipeline to compare the output of three cellular Potts models of 2D angiogenesis (13,27,28) to in vitro experiments over time at different cell densities and find that the distance at which cells are able to communicate with each other determines the features of the final, stabilized network, but the speed and manner at which the networks stabilize depends on other characteristics of the cell communication, like the speed at which the cell attracting signal spreads.

Quantification of in vitro endothelial network dynamics
To understand and compare the dynamic behavior during in vitro and in silico EC network formation, we developed an image analysis pipeline to extract network characteristics from 2D angiogenesis timelapses and applied it to in vitro and in silico networks.To capture in vitro EC network formation ECs were seeded on a Matrigel matrix and imaged in phase contrast for 24 hours with an interval of 15 minutes (Fig 1A).The ECs in the timelapse images were segmented from the background using a combination of Gaussian denoising and local variance extraction, followed by Huang thresholding (S1 For the in silico networks we selected three cell-based computational models of 2D angiogenesis (13,27,28) and extracted the network characteristics directly from the output images using the same image analysis pipeline as the segmented in vitro timelapses.In the in vitro culture under standard conditions, the ECs started forming cell-cell connections within the first hour after seeding on a Matrigel matrix and thereafter self-organized into fully connected networks within six hours (Fig 1C).We were able to divide the development of the network characteristics over time in three distinct parts, which we describe following Parsa et al. (14), as: [1] The first hour resembled stages 1 and 2, "rearrangement and aggregation", followed by "spreading"

In silico models can be adjusted to resemble in vitro networks
To compare the dynamics of in silico networks to the in vitro networks we selected model parameters  In the elongated cell model, ECs secreted a chemoattractant, migrated upwards the chemoattractant gradient and elongated towards a target length (13).In the contact inhibition model ECs secreted a chemoattractant and migrated upwards the chemoattractant gradient, but locally inhibited chemotaxis in their neighbors (Fig 2A) (28).To find the parameters that describe the chemoattractant in the two chemotaxis models for which the models most resemble the in vitro networks (i.e., the diffusion coefficient, the decay rate and the secretion rate) we investigated the effect of the diffusion length of the chemoattractant on the network formation.Diffusion length (  ) depends on the diffusion coefficient () and the decay rate () of the chemoattractant according to   = √   .Therefore we altered the diffusion length from   = 10  to   = 100  with steps of 10 µm by varying the decay rate of the chemoattractant, while keeping the diffusion coefficient constant at  = 5 ⋅ 10 −13  2  −1 (Fig 2C).For a full list of the selected parameters, see S1 Table .In vitro, the average number of lacunae per mm² after 24 hours was 4.23 ± 0.81, with an average lacuna area of 0.17 ± 0.03 mm² (Fig 1D and E).When we increased the diffusion length in the two chemotaxis model, we observed that the average number of lacunae decreased (Fig 2C ), in agreement with previous work (13,28).
Based on the outcome of the comparison of the number of lacunae, we selected a decay rate of  = 1.02 ⋅ 10 −4  −1 , which corresponds to a diffusion length of   = 70 .This resulted in an average number of lacunae after 2880 MCS of 3.27 ± 0.47 per mm², with an average area of 0.15 ± 0.02 mm² for the elongated cell model and 4.88 ± 0.35 per mm², with an average lacuna area of 0.10 ± 0.01 mm² for the contact inhibition model (Figs 2D and E).With the selected diffusion length of   = 70 , we observed some similarities between the dynamics of the model outcomes and the in vitro networks: During the initial relaxation time (100 MCS) the ECs in the models reached their target size and shape, which mimics stages 1-3 of the 5 stages described in ( 14 In the mechanical model, the ECs are assumed to stiffen the ECM through the exertion of traction forces on their environment and move upwards the stiffness gradient (Fig 2B) (27).We selected a Young's modulus of 12 kPa, with a Poisson's ratio of 0.45, based on the original work.For a full list of the selected parameters, see S2 Table .For the selected parameters, cells formed a dense network of branches with small lacunae.The composition this network was similar to an early-stage endothelial network.This suggests that in the present form of the mechanical model, the mechanical model successfully reproduces the network formation phases (phases 1-3), but it did not describe the remodeling phases well.
Another noticeable difference between the three models and the in vitro networks was the composition of the branches (Figs 1B and 2B).In vitro the cell-covered area decreased over time as cells move on top of each other, whereas in the in silico models the cell-covered area remained constant (S2 Fig) .This is likely caused by the fact that in vitro, the ECs are able to crawl on top of each other, thus forming thin branches and thick nodes, which is impossible in the two-dimensional cellular Potts model.

Increasing cell density results in an initial increase in number of branches and lacunae in vitro
To investigate whether the models respond similarly to a change in the initial conditions as in vitro networks, we next analyzed network formation at different cell densities.Previous in vitro studies on primary endothelial cells observed an inability for the cells to form tubes under or over a threshold density, similar to a percolation threshold (4,22).Consistent with these previous results, we observed that HMEC-1 cells retained the ability to form cell-cell connections and short branches at 1000 ECs per well (80 cells/mm²), however, they were not able to form a network spanning the full well (Fig 3A).Over 2000 ECs per well (160 cells mm²), we observed the development of a fully connected network spanning the entire well (Fig 3A).With more cells per well a higher number of lacunae and branches formed within the first hours, but for more than 4000 cells per well (320 cells mm²), the networks converged onto the same number of lacunae, branch number and branch length (Figs 3B-D).

The cell elongation model response resembles that of the in vitro networks to a change in cell density
Whereas a higher cell density caused an initial increase in the number of branches for all three in silico models, this difference diminished as the network remodeled in the two chemotaxis-based models (Figs 4A-F), similar to what we observed in vitro (Fig 3B).The mechanical model had a much higher number of branches and the difference between branch number for different cell densities persevered for much longer than in the chemotaxis-based models and the in vitro timelapses (Figs 3B and 4C).
Out of the three models, the elongated cell model most resembled the in vitro network response to a change in cell density, even when subjected to a range of diffusion lengths.

Discussion
In biomedical applications, like tissue engineering (9) or tumor vascular normalization (10), angiogenesis is stimulated, redirected or reduced to properly vascularize the tissue of interest.
Computational models can help understand how different environmental factors and cellular cues affect various aspects of EC migration and network formation and help identify effective regulatory targets.Many computational models of angiogenesis, based on assumptions like cell-cell attraction, cell shape effects, and cell-ECM interaction, have been successful at creating EC networks that resemble those observed in vitro qualitatively (18)(19)(20).To test whether the manner in which these simulated ECs form networks is also indicative of in vitro EC behavior, we selected three cell-based models of 2D angiogenesis based on different assumptions, to examine their dynamical behavior.To We used the outcome of our initial analysis of the in vitro and simulated networks to find a parameter space in which the network features of the model outcomes were similar to those of the in vitro EC networks (Fig 2C).Within this parameter regime, we conducted further comparisons between the models and the in vitro networks to assess if systematic, concurrent parameter changes in vitro and in silico could help assess the correctness of the models.Interestingly, both chemotaxis-based models showed similar dynamic behavior as the in vitro networks.They show an initial cell-cell contact forming phase, followed by a remodeling phase, similar to stages 3 to 5 of in vitro network formation: elongation and formation of cell-cell contacts, plexus stabilization, and plexus reorganization.In the mechanical model, however, the ECs quickly form connections, mimicking early in vitro network formation, but there is little network reorganization in the model outcomes.
We found that the dynamic behavior of the lacunae in the elongated cell model is most similar to those of the in vitro networks (Figs 2D and E).In the elongated cell model the ECs actively elongate before they migrate upwards a chemoattractant gradient (13).Because of their shape, the elongated ECs can form long-range connection to other ECs, and thus form branches and close lacunae more rapidly than rounded ECs.This improves the similarity of the model to early stages of in vitro network formation (Figs 2E and F).Due to the slower movement along the short axis for elongated cells, branches are not very mobile and merge slowly as the network stabilizes (Fig 2F).In the contact inhibition model, the ECs are rounded and less sensitive to the chemoattractant when in contact with other ECs (28).The  4D-F).This might be due to an overall lower sensitivity to the chemoattractant due to the contact inhibition.When conducting our search for the optimal diffusion length for the parameter regime in which we compare the networks, we also observed differences in response to changes in decay and secretion rate between the two chemotaxis models.We observe that the contact inhibited ECs require a stronger chemoattractant signal to form a network than the elongated cells.The cell elongation also allows for network formation at lower cell density than the contact inhibition (Figs 4 and S4).From this we can conclude that elongation can help ECs to form networks under circumstances where contact inhibition cannot.The elongated shape of ECs allows for the ECs to sense a chemoattractant gradient at a larger distance from their center of mass, whereas due to contact inhibition ECs become less sensitive to the chemoattractant overall.
It is important to note that the chemoattractant in the simulations can represent any sort of mid-range cell communication signal.Rüdiger et al. (2020) argue that mechanical forces are the main driver behind in vitro endothelial network formation (4).In this study they continuously perfused the ECs during in vitro network formation to wash away any possible chemoattractant gradient, and still observe network formation.To further argue for mechanical regulation, they show that ECs compress the substrate up to a distance of 30 µm, with measurable ECM displacement up to 100 µm.These distances fall within the range of what we observed in the two chemotaxis models similar to the in vitro networks.In the mechanical model we observed strains up to 40 µm from the cell boundary in the investigated parameter regime (S2 Fig) .However, this was not enough to form networks similar to those observed in vitro.In the mechanical model, as it is implemented here, the strain field resets every MCS and is recalculated based on the current cell shape.In vitro, we and others (4,45) observed that ECs remodel the substrate over time, altering its material properties and response to stress.
Resetting the strain field each timestep omits any changes in the composition of the environment from the model.This could explain why the EC network in the mechanical model stabilizes quickly, whereas in vitro and in the two chemotaxis-based models the EC networks continuously evolve as the simulated chemoattractant gradient gradually stabilizes over time.The effect of matrix remodeling, in the form of strain-stiffening, on EC network formation dynamics could be investigated in vitro with synthetic gels with tunable stiffness.With these gels the travel distance of mechanical perturbations on the system can be altered systematically, to measure its effect on EC network formation.With the use of synthetic gels, it is also possible to wash away any chemoattractant from the system, without the possible interference of ECM retention, which might have influenced the network formation in

Rüdiger et al. (2020).
In the two chemotaxis-based models we assumed each MCS represents 30 seconds in real time and each pixel is 2 µm.This was chosen such that the mean velocity of the cells became 5  ℎ −1 , which, it was previously argued, complies with in vivo observations (13,38).This allows us to translate the physical processes like chemoattractant diffusion and decay rate to physical units ( 2  −1 ;  −1 ).We fixed the diffusion coefficient at  = 5 ⋅ 10 −13  2  −1 .The main driver behind angiogenesis is generally assumed to be VEGF-A (39).The diffusion coefficient of VEGF-A is in the order of 10 −11  2  −1 in collagen-I, 10 −10  2  −1 in Matrigel and 10 −9  2  −1 in water (40)(41)(42), which is 100 to 10000-fold higher than what we implemented in the model.However, ECM and EC membrane receptor retention of VEGF-A could affect VEGF-A diffusion, resulting in a smaller effective diffusion coefficient and decay rate (6).Köhn-Luque et al., 2013 show an aggregation of fluorescently labeled VEGF in the substrate surrounding HUVECs on Matrigel within 5-10 minutes after VEGF addition (6).Additionally, it is yet unclear whether VEGF-A diffusion is the main driver behind angiogenic network formation in vitro.Especially since Rüdiger et al. (2020) observed endothelial network formation under constant perfusion to avoid the emergence of a chemical gradient (4).With a diffusion coefficient of  = 5 ⋅ 10 −13  2  −1 we observe simulated networks most like in vitro at a diffusion length of   = 70  (Fig 2C).With a diffusion coefficient of  = 10 −10  2  −1 , the decay rate would need to be  = 2 ⋅ 10 −2  −1 for a diffusion length of   = 70 , since   = √   , which is much faster than the 90 minute half-life of VEGF-A measured in vitro (43).In vivo, VEGF-A has been shown to have a half-life of 33.7 minutes (44).
There are some steps in this image analysis pipeline that require optimalization.For example, we implemented a watershed step in the pipeline to correct for breaks in EC branches where thin stretched-out ECs are misinterpreted as background during segmentation.This watershed step can result in an overrepresentation of connections between branches, but overall, we observe an accurate and consistent result (S1 Fig) .To circumvent the need for segmentation corrections and other suboptimal parameter settings a neural network could be trained to perform the segmentation steps necessary for the network quantification (37).This network can be trained to recognize thin stretches of ECs where conventional image analysis methods, which rely mostly on contrast, might mistake them for background.However, the training of a neural network requires a lot of correctly annotated training data and undertrained networks will give suboptimal results.More interestingly, a neural network could be trained to recognize network features that are more difficult to extract from segmented images automatically, like cell division events.
Conventional tube formation assays to assess the sprouting ability of ECs are usually only examined for a single characteristic (e.g.number of branches) at a single timepoint (32).However, we observed that network features converged over time independent of their initial observed value (Fig 3B).This indicates that an endpoint observation might not suffice to examine the influence of a compound or knockout on collected cell behavior before a stable network has been formed.We also observed a high biological variation between different samples of the same experimental conditions.This can be explained by variations in fitness of different passage number of ECs and small fluctuations in cell seeding density.Therefore, it is important that the exact number of seeded ECs is also considered when comparing samples within 6 hours after seeding or at low cell densities.Passage number has previously been observed to affect tube formation in mouse 3B-11 endothelial ECs (12).Another effect on the fitness of the cells could be the duration for which ECs spent outside of the incubator prior to the experiment.
In this work we have shown the importance of dynamic analyses of EC network formation by showing how EC network features converge over time.We paved the way for a more systematic comparison of computational models to biological experiments, and we have indicated areas of improvement for existing cellular Potts models of 2D angiogenesis.In the future we will expand our analysis of EC networks and find new parameter regimes to explore both in vitro and in silico.
Cells were cultured in a humidified incubator at 37 °C and 5 % CO2 and split twice per week.Passages used were between 9 and 17.

Tube formation assay
Cells at 70-80% confluency were plated on growth factor reduced (GFR) Matrigel (Corning) on Angiogenesis µ-slides (Ibidi) at a density of 10,000 cells per well unless stated otherwise.For fluorescent plasma membrane labeling cells were incubated with PKH67, PKH26 or CellVue Claret (Sigma-Aldrich) in Diluent C as described in (14).

Imaging
All images were captured at a magnification of 10X using a Nikon Ti inverted phase contrast microscope equipped with a DS-Ri2 camera or a Zeiss Axio Observer inverted phase contrast microscope equipped with an AxioCam 705 mono camera with a motorized stage for 24 hours with 15-minute intervals.Fluorescent images were acquired using the Zeiss Colibri multicolor LED light source.Full well images were created by stitching multiple fields together with the NIS-elements and ZEN Blue software.

Network characterization
Time-lapse images of the full wells were processed and quantified using a custom Fiji/ImageJ pipeline (v1.53f) (Fig. S1).For the segmentation of cells from background the following steps were performed for each slice in the timelapse images: [1] High-frequency features were amplified by first convoluting the phase-contrast image with a Gaussian filter with sigma value of 10 µm, then subtracting the resulting low-frequency image from the original image and finally adding the resulting high-frequency image back to the original image.[2] Hereafter, the local variance in the sharpened images was computed with a radius of 8.75 µm.The local variances were convoluted with a Gaussian filter with a sigma value of 3.5 µm.[3] From the resulting image a threshold was computed with automatic Huang thresholding to create a segmented mask of the network.[4] The mask was then corrected for broken branches by producing a distance map of the lacunae larger than 2000 µm² and watershed segmenting the distance map with a threshold set to 50.Hereafter the smaller lacunae were added back to the watershed network for the final mask.
The branches and nodes within the binary masks were analyzed with the Skeletonize 2D/3D and Analyze Skeleton (2D/3D) plugins applied on a 10x scaled down image to avoid a very noisy skeleton.
The lacunae were analyzed with the Analyze Particles plugin.The resulting output includes the number, size, position and shape of the lacunae, and the number, length, position and connectivity of the branches.
All results were exported to CSV files for further analysis in RStudio.This analysis pipeline was automated to process large quantities of time-lapse images sequentially without user input to maximize reproducibility and minimize personal bias in the data analysis.For the validation of the pipeline, we had six independent participants encircle three lacunae and compared the average lacuna area of their selections to automatically segmented images.The automatic segmentation underestimated the lacunae as selected by the participants by 4.2 ± 11.7 % of the average area (n = 15).The difference between two different participants encircling the same lacunae was 4.5 ± 2 % of average area (n = 3) (S1 Fig).

Computational model
The three computational models of angiogenesis (13,27,28) are all based on hybrid cellular Potts models (CPMs).Hybrid CPMs couple a dynamical description of cell behavior based on the CPM with a detail model of the cellular micro-environment, specifically the concentration of extracellular signaling molecules or strains and stresses in the ECM.The cellular Potts model is a cell-based model that focusses on modeling cell-cell interactions and predicting collective cell behavior (36).Cells are represented as collections of lattice sites  ⃗ with an identical index   on a grid, where each index labels a single cell.Each lattice site represents an area of 2 µm x 2 µm in the chemotaxis models and 5 x 5 µm in the mechanical model.To mimic the circular in vitro dishes, we labeled all lattice sites outside of a circle with a diameter of 3800 µm sites as boundary sites   = −1 in the elongated cell and contact inhibition models and as a fixed node in the mechanical model.
Apart from the shapes of the simulation domains, models were used as described previously (13,27,28).Briefly, during a Monte Carlo step (MCS), a time-step in the simulation representing 30 seconds, cells change their shape and move through the lattice by attempting to copy its index to a neighboring lattice site.There are N copy attempts per MCS, with N the number of lattice sites.If the copy attempt results in a lower effective energy of the system, the copy attempt is accepted., where    ⃗ ⃗⃗ ,  ⃗ ⃗⃗ ′ represents the bond energy between a lattice site  ⃗ and its eight second-order neighbors  ⃗′, and   and   represent energy penalties for deviation from a designated target area (   ) or target length (  ), which we set at 50 and 5 respectively.If the copy attempt does not result in a lower effective energy, the copy is accepted with Boltzmann probability: where T is the "cellular temperature" of the system, a measure of cell motility, which we set at 50, and Δ the change in the effective energy function.
In the elongated cell and contact inhibition models (13,28)  Where  is the strength of the chemotactic response, which we set at 1000.For details see (13,28).
In the mechanical cell-cell communication model (27), cells were able to deform the ECM through cell-shape dependent traction forces.The deformations of the ECM are calculated using the finite element method.An additional term is added to the energy function, such that cells respond to the ECM stiffness: Δ  = −( ⃗,  ⃗ ′ )  (ℎ(( 1 ))( ⃗ 1 •  ⃗  ) 2 + ℎ(( 2 ))( ⃗ 2 •  ⃗  ) 2 , where ( ⃗,  ⃗ ′ ) = 1 for cell extensions and ( ⃗,  ⃗ ′ ) = −1 for cell retractions,   is a parameter which we set at 24, ℎ() mimics the influence of the ECM stiffness  as a function of  1 and  2 , the principal strains, on the tension in the cell-ECM connection,  ⃗  gives the copy direction, and  1 and  2 are the strain orientations.
The preference for higher stiffness  is implemented as a sigmoid function: where  sets the strength of the durotactic response,   determines the stiffness where half this strength is reached and  represents the stiffness sensitivity which determines the steepness of the curve.For details see (27).

Fig) .
Fig).Network features were extracted from the segmented timelapses for each timestep using FIJI

Figure 1 .
Figure 1.Image analysis pipeline for in vitro phase contrast timelapses.A) Graphic summary of in

(
Fig 1C).In these stages the ECs were randomly dispersed on top of the Matrigel and started to form attachments to the ECM.The next 6 hours matched Parsa's stage 3, "elongation and formation of cellcell contact".During this stage new connections formed between ECs, which is reflected in an increase in the number of lacunae from 2.78 ± 0.34 to 6.84 ± 1.06 per mm² (Fig 1D).Hereafter, the number of lacunae decreased linearly in time (R²=0.91,p-value < 0.001), resembling Parsa's stage 4 and 5: "plexus stabilization" and "reorganization".The average lacuna area displayed a linear increase due to a higher occurrence of lacuna merging and smaller lacunae closing compared to the formation of new lacunae (Fig 1E) (R²=0.93,p-value < 0.001).The number of branches within the network continuously decreased and their average length increased (Figs 1G and H).During the first three stages branches predominantly increased in length due to merging of branch sprouts, but during stages 4 and 5 we observed an increase in merging of existing branches.
Fig).We observed an average cell area of 1052 ± 596 µm² after three hours (S2 Fig), with a cell length

Figure 2 .
Figure 2. Quantitative comparisons between two chemotaxis-based models and a mechanical ): cell spreading and elongation.Next, we looked at the formation of cell-cell connections as ECs started excreting chemoattractant.The increase in cell-cell connections could be indicated by the decrease in number of branches (Fig 2F).The transition towards stage 4, stabilization of the network, was more gradual in the contact inhibition model than in the in vitro networks and in the elongated cell model (Fig 2E).The number of lacunae decreased faster in the elongated cell model than in the contact inhibition model (Fig 2E) and for the first 1680 MCS we observed fewer and larger lacunae in the contact inhibition model than in the elongated cell model.The number of branches was initially larger in the elongated cell model than in the contact inhibition model.In agreement with the in vitro networks, the branch number decreased quasi-exponentially in the elongated cell model, whereas in the contact inhibition model the number of branches decreased more linearly (Fig 2F).Both chemotaxis-based models have longer branches than the in vitro networks (Fig 2G).This difference in branch number and length could partially be explained by the higher number of sprouts detected in vitro than in the models (Fig 2H).Sprouts are branches that are connected to a single node and the Analyze skeleton tool splits an existing branch into two each time a sprout appears, resulting in two shorter branches.

Figure 3 .
Figure 3.A comparison of the response of in vitro networks and model outcomes to higher

Figure 4 .
Figure 4.A comparison of the response of model outcomes to higher seeding density.Quantitative analysis of the number of branches for different cell densities with different diffusion lengths (  ) in the (A-C) elongated cell model and (D-F) the contact inhibition model, and different stiffness sensitivity () (G-I) in the mechanical model.( = 5.0 ⋅ 10 −13  2  −1 ,  = 1.02 ⋅ 10 −4  −1 ,  = 1 ⋅ 10 −3  −1 ).Shaded area represents standard deviation of six simulations.
quantify dynamical features of in vitro EC network formation (Fig 1B) and of simulated EC network formation outcomes of cellular Potts models of 2D angiogenesis we developed a custom image analysis pipeline (Fig 2A) (13,27,28).

S3Fig.
Overview images of simulated networks for different cell densities after 2880 MCS.A) Elongated cell model.B) Contact inhibition model.C) Mechanical model.