Effects of endothelial cell proliferation and migration rates in a computational model of sprouting angiogenesis

Angiogenesis, the recruitment of new blood vessels, is a critical process for the growth, expansion, and metastatic dissemination of developing tumors. Three types of cells make up the new vasculature: tip cells, which migrate in response to gradients of vascular endothelial growth factor (VEGF), stalk cells, which proliferate and extend the vessels, and phalanx cells, which are quiescent and support the sprout. In this study we examine the contribution of tip cell migration rate and stalk cell proliferation rate on the formation of new vasculature. We calculate several vascular metrics, such as the number of vascular bifurcations per unit volume, vascular segment length per unit volume, and vascular tortuosity. These measurements predict that proliferation rate has a greater effect on the spread and extent of vascular growth compared to migration rate. Together, these findings provide strong implications for designing anti-angiogenic therapies that may differentially target endothelial cell proliferation and migration. Computational models can be used to predict optimal anti-angiogenic therapies in combination with other therapeutics to improve outcome.

modeling sprouting angiogenesis and for the prediction of potent anti-angiogenic treatments for reducing tumor size, inhibiting or slowing growth. Modeling can help elucidate the contributions of endothelial cell proliferation and migration to vascular coverage, thus enabling the prediction of which mechanism would be the most effective for drug targeting.
There are several in vitro ways to access proliferation and migration of endothelial cells. A common way of measuring proliferation is with BrdU, which living cells incorporate into their DNA and allows them to be counted 10 . Colorimetric proliferation assays are also common, such as WST-1, cell counting kit-8 assay, which causes the reduction of formazan dye in proportion to the number of living cells and can be measured with a fluorescence plate reader 11,12 . Migration can be measured in real time using an RTCA reader based on electrical impedance 13 . Wound healing type assays are also used, where cells are plated with a stopper in the center of the chamber, which is then removed. After a certain number of hours, the cells that have migrated within that region can be counted 11,12 .
A significant class of studies have been performed in the area of computational modeling of angiogenesis (for reviews see 14,15 ) and as it relates specifically to tumor growth (for reviews see [16][17][18] ). The migration of tip cells, both chemotactic and haptotactic, was examined in response to differences in angiogenic factor distributions 19 . A cellular Potts model was used to understand the growth of vasculature under different extracellular matrix (ECM) conditions 20 . A significant number of studies have focused on the tumor vasculature in regards to drug delivery for cancer treatment 21 . Computational models have examined the effects of interstitial and vessel pressure on drug delivery 22 , transport of different drugs on tumor treatment 23 , delivery of anti-angiogenic drugs on tumor growth 24 , as well as examining the role of vasculature and drug delivery on drug resistance 25 . These models help us understand the overall process of angiogenesis but in this study we are interested in the specific contributions of endothelial cell migration and proliferation.
The importance of endothelial cell migration and proliferation has been a topic of several modeling studies. Burke et al. used a 2D lattice-based model to study endothelial cell proliferation and migration under mechanical stress 26 . They determined that the changes seen under mechanical stress could only be reproduced by increasing migration, increasing proliferation and having biased endothelial cell migration perpendicular to the direction of the strain. Another model used partial differential equations to study angiogenesis as an interplay between cell adhesion, traction force (elongation), and proliferation; the authors found that straight (non-tortuous) vasculature was found when proliferation was triggered by endothelial cell strain 27 . This model focuses on mechanical strain under in vitro conditions, rather than in a tumor microenvironment. In a cellular Potts model, migration and growth were considered to be modulated by blood flow, without which the vasculature would collapse 28 ; the model did not examine the individual contributions of endothelial cell elongation, proliferation and migration.
Popel and his colleagues in a series of studies 29,30 have formulated a class of 3D models of angiogenesis at multiple scales. In parallel with these 3D simulations, whole-body compartmental models have been used to predict VEGF distributions and VEGF receptor occupancy 31,32 as well as pharmacokinetics and pharmacodynamics of anti-angiogenic drugs, such as bevacizumab and aflibercept 17,31,33,34 . Other models focused on the effects of vascular endothelial growth factor (VEGF) and Delta-like 4 Notch ligand on angiogenic tip sprouting 29,35 . In this paper, we develop a modified version of the angiogenic tip sprouting module, including tip cell migration, endothelial cell proliferation, and sprouting. This model is then used to understand the interplay between endothelial cell proliferation and migration rates in tumor angiogenesis. Specifically, the model investigates the therapeutic effects of increasing or decreasing endothelial cell migration and proliferation. Thus, we can use this model to predict angiogenic response to therapeutics affecting migration or proliferation. The problem of the relative contributions of proliferation and migration is very important for the general understanding of the regulation of angiogenesis, and also for translational applications since there are pharmacological agents that can differentially affect these phenomena. For example, there are anti-angiogenic peptide agents that preferentially affect blood or lymphatic endothelial cell proliferation and migration 36,37 . It should also be noted that these results would be important in the field of therapeutic angiogenesis where the growth of blood vessels is stimulated to provide blood flow and oxygen delivery to ischemic tissues, e.g. in the case of coronary and peripheral artery diseases, wound healing, and regenerative medicine 38 . Thus, the results of this study could also provide guidance to pro-angiogenic therapeutic treatments.

Results
Simulations were performed on a 500 × 500 × 500 micron grid with the initial capillaries located along the edges of the cube. The proliferation and migration rates were varied for each simulation and the different simulation runs were averaged over 10 different trials. We will use the following metrics to characterize the vasculature similar to those used in describing vascular networks in tumor xenografts 39 . These are: Vascular Length Density, VLD, in mm/mm 3 ; Bifurcation Density, BD, in 1/mm 3 ; Vascular Segment Length, VSL, (length of vascular segment between adjacent bifurcations) in mm; Vascular Segment Tortuosity, VST, (total distance between bifurcations/ Euclidean distance between bifurcations) in mm/mm; Fractal Dimension (FD). We will report mean values and variation of these variables and also where appropriate their distributions. Note that in characterizing vascular tortuosity we follow the morphological definition in ref. 40 rather than an alternative functional definition in ref. 39.
Sprouting Angiogenesis Progression over Time: We performed simulations under different parameter values and calculated the progression of the emerging vasculature over time. The motivation and focus of this investigation is on how the proliferation (PR) and migration rates (MR) affect vasculature, since these rates can be differentially attenuated by selective agents in proliferative diseases like cancer; alternatively the rates can be differentially increased by agents in the case of ischemic diseases, wound healing or regenerative medicine. In these simulations, VEGF is at a constant value of 20 ng/ml. Figure 1 shows the progression at days 1, 10 and 20 with a fixed migration rate, MR = 10 μ /hr and three different values of the proliferation rate, PR = 0.0147 (A), 0.025 (B), 0.083 (C) 1/hr (corresponding to divisions every 68, 40, and 12 hours, respectively). Doubling times for microvascular endothelial cells (MEC) range from about 12 hours to 4 days 41-43 and human umbilical vein endothelial cells (HUVEC) range from 17 to 72 hours [44][45][46] . PR has a major effect on the progression of the tumor vasculature. At day 20, the values of parameters introduced above are: VLD = 72.4, 130.6, 16.9 mm/mm 3 , VSL = 0.193, 0.233, 0.124 mm, and VST = 1.37, 1.41, 2.60 mm/mm, respectively. When PR is low, the vessels reach the interior of the tumor space but the capillaries do not cover much of the area, Fig. 1A; in this case the VLD, VSL, and VST values are intermediate. Thus, migration allows the vasculature to extend into the tumor space but with fewer cells and thus fewer bifurcations covering the tumor space. When PR is medium, the tumor vasculature spreads throughout the entire region and has high VLD and VSL values and an intermediate VST, Fig. 1B. In this case, there is a balance between proliferation and migration, where proliferation is small enough to allow for the extension of the tip cell agent but fast enough for the vessels to grow as well. This leads to extensions into the tumor space due to migration and better coverage due to branching off of the new cell agents. When PR is high, the resulting vasculature is minimal but very tortuous with a low VLD and VSL and high VST level. The vasculature does not fully extend into the interior of the tumor space, Fig. 1C. This is partially due to the fact that the high PR limits the extension of the vasculature and so the capillaries do not spread into the interior. The growth is further limited by anastomosis; since the sprouts are short they end up growing into another tip cell agent and anastomosing, which in turn limits sprouting due to vessel regression. In general, the vasculature is the fullest at medium values of PR and is most tortuous at high values of PR. We have also performed simulations in the presence of a VEGF gradient to explore the effect of the gradient, see Supplemental Data for more information.
Medium Proliferation Values Produce the Most Vascular Coverage: In Fig. 2 we show simulation images with different PR and MR at day 20. Day 20 is chosen here since it is relevant to the timeframe of tumor xenograft growth that will be used in comparison with experimental data below. With low PR and MR values, the vasculature is very minimal. Having high PR also results in very minimal architectures that are very tortuous. With low PR and high MR, the capillaries are long but there are not many sprouts and the tumor space is not well covered. With medium amounts of PR and MR, there is complete coverage of the tumor space. Here it is clear that the vasculature reaches the middle of the space and the density of the vasculature is not too high. Thus, these results suggest that medium proliferation values would be optimal for cancer vascular recruitment.
Most Vascular Metrics are Determined by Proliferation: To better understand the contribution and synergy of PR and MR in the sprouting angiogenesis model, we conducted a parameter space evaluation, varying MR between 0.24 and 40 μ /hr and PR between 0.015 and 0.083 1/hr for 20 days, Fig. 3. We find that the bifurcation density (BD) and the vascular length density (VLD) are mostly dependent on proliferation. There is an optimal PR value of about 0.025 where the BD and VLD are highest except at very low MR. At very low or very high PR, the BD is low. At low PR, there are fewer cells and therefore there are fewer places to sprout, resulting in a reduction of the number of bifurcations. At high PR, there tend to be too many sprouts formed over a short amount of time resulting in anastomosis, reducing the number of bifurcations and the number of vascular segments. Larger or smaller PR values cause a decrease in BD and VLD. In terms of treatment strategy, this would mean that decreasing proliferation would be an effective strategy for limiting vasculature. Surprisingly though, increasing proliferation would also be an effective strategy. This is due to the fact that highly proliferative vessels tend to anastomose quickly limiting their overall growth. Tortuosity (VST) is dependent on both PR and MR. The vasculature is most tortuous at very low MR or very high PR. This is due to the fact that these values lead to the shorter segments that curve due to their random walk.
Comparison with Experimental Data: We compare our model results to in vivo data derived from a breast cancer xenograft model, where the mouse vasculature was perfused with a polymer solution and excised tumors were imaged with microCT for the purpose of visualizing the entire tumor vascular bed 39,47 . A digital map of the vascular network was obtained with a specially developed 3D reconstruction algorithm. Thus, 3D reconstructions of the entire xenograft tumors were able to be developed and vascular metrics were calculated for these tumors (Imaging Data). These metrics were also compared to data from the literature (Literature Values). Several metrics were used to characterize the vasculature, such as vascular length density, vascular segment length, bifurcation density, tortuosity, and fractal dimension, Table 1. We also present the ranges of these parameters for different tumors from the literature. The simulation results generally fall within the range of those seen in vivo. More detailed and focused simulations are necessary to simulate a specific tumor type.
Regression Distance does not Affect Mean Vessel Length: We varied the distance d r, in which vessels growing near a mature vessel would regress, from 50 to 200 microns. Any new sprout within this distance of a mature vessel would retract, the base example being 100 microns. Using a radial distance of 200 microns, we found overall that the mean BD and VLD were lower than normal, 136 vs 194 bifurcations per mm 3 , and 18 vs 27.9 mm/mm 3 . In contrast, the VSLs were comparable: 0.107 and 0.104 mm. There were opposite effects using a regression distance of 50 microns. BD increased from 194 to 315 and VLD increased from 27.9 to 43.9 mm/mm 3 . Once again the VSL were the same 0.104 mm.
Vascular Network Characteristics as Functions of Time: In the previous simulations we presented the results for t = 20 days; this time point was chosen as a typical duration for tumor xenografts, eg 39 . However, it is of interest to investigate how the network characteristics vary with time as the vasculature continues to evolve. Thus, we ran the simulations for 200 days. We plotted the results for BD averaged over 10 runs as functions of time and found that all simulations reached a steady state, Fig. 4A; we also plotted a sample simulation, Fig. 4B, that exhibits fluctuations associated with sprout formation and retraction. Most of the simulation values resulted in high BD (200 or above) but there were a small sample of cases where there were very few bifurcations < 20. The cases with low BD all had high PR values of 0.083, corresponding to proliferating every 6 hours. In conclusion, the limited number of bifurcations at high PR values does not change with more time.

Discussion
Computational models can be useful in formulating hypotheses, especially when care is taken to justify each abstraction made and to reconcile it with biological reality. This task is made easier by the discrete modeling of endothelial cells as agents capable of independent decisions as well as inter-agent actions. Discrete models like this can serve as a framework for the type of module-based modeling previously described 14,29,30 which attempt to combine independently validated models of cellular processes as modules in a larger integrated model of angiogenesis, or can be incorporated into other models as processes such as tumor growth or metastatic colonization in which angiogenesis is involved. It is well known that tumor vasculature is abnormal and differs from the normal vasculature of the host organ 48,49 . One of the abnormalities is that the vessels are very tortuous 50 in and around the tumor. Endothelial cell proliferation and migration are controlled by a myriad of factors in the tumor microenvironment including a large number of pro-and anti-angiogenic factors that could be used for pro-and anti-angiogenic therapy 51 . Therefore in an experimental setting it is very difficult to modulate or separate out these two processes. In contrast, it is possible to modulate proliferation and migration rates separately in a computational model. Thus, we examined the effects of both of these properties on the characteristics of a growing vascular network.
The overall growth of vasculature is a result of both proliferation and migration. Cell migration in the model is a combination of elongation and motility, such that it includes elongation because the length of the tip cell agent gets longer. It also includes motility such that the leading node of the tip agent can change direction and the rate in which it migrates. It should be noted that the motility of a tip cell does not directly relate to isolated cell migration measured in vitro, such that isolated cells moving on a plate or in 3D matrix may not follow the same dynamics as a tip cell migrating with stalk and phalanx cells behind it. In the model we show that when either proliferation or migration is acting on its own, the vasculature is shunted and ineffective at covering the tumor space. This can be seen by comparing the VLD, at low PR with high MR yields VLD of around 20-30 mm/mm 3 , high PR with low MR rate yields VLD of around 10 mm/mm 3 . whereas an intermediate value of both yields a VLD of around 80-90 mm/mm 3 . Thus, neither proliferation nor migration alone can create organized vascular growth; a balance is necessary between the two processes.
Comparing the contour plots of the parameter space with the representation of the tumor vasculature, we find interesting results. When examining the contour plot for VLD, it is clear that values with high PR and low to medium MR result in large values of total vasculature. When examining the parameter space, simulations with low PR do not have large BD or VLD. High levels of PR actually decrease the BD and VLD. This suggests that tumors that give uncontrolled proliferatory signals might actually cause the vasculature to be sparse in certain areas of the tumor resulting in areas of hypoxia. Thus, controlling proliferation through vascular normalization 52 , may improve the vessel coverage and improve dysfunctional vasculature and drug delivery.
We can relate our results to work done by Lee et al. who analyzed the effects of anti-angiogenic peptides on the proliferation and migration of endothelial cells 11 . For example, somatotropin-derived peptide SP5028 inhibits migration of microvascular endothelial cells (at 50 μ M concentration at 20 h) by 70% and is also a potent proliferation inhibitor; thus, our model predicts that it would be an effective angiogenic inhibitor. Also, another peptide from the same class, SP5031, is a less potent proliferation inhibitor, but it reduces migration by > 90% which is predicted to be effective at inhibiting angiogenesis. In agreement with these data, the peptide SP5031 significantly reduced VEGF-induced angiogenesis, which was measured by staining with CD34 in an in vivo subcutaneous Matrigel plug assay. Using a collagen-derived biomimetic peptide SP2043, a 95% reduction in migration would inhibit bifurcations and vascular segment density enough to inhibit angiogenesis so treating with 25 μ M of SP2043 is predicted to be sufficient to inhibit angiogenesis in agreement with experimental data 53 .
The major prediction of our model is that PR would have a greater effect than MR on reducing vascular growth. In support of this conclusion, Jackson and colleagues found that proliferation is necessary for the extension of the vasculature otherwise the capillary ceases to extend 54 . Very low MR causes a reduction in BD and SLD but overall proliferation has more of an effect on these metrics. Tortuosity is governed by both PR and MR, such that high tortuosity occurs at very high PR or very low MR. This makes sense because migration causes vessels to be straighter 7 and thus low migration reduces this; high proliferation limits BD because it limits the extension of each tip cell agent. Interestingly, our model predicts that increasing PR might also be an effective strategy at reducing vasculature, as it leads to an increase in sprouts which then anastomose to create non-functional vessels. Others have found that increasing branching can lead to tumor growth inhibition by forming non-productive angiogenesis 55,56 . This is in line with the commonly seen effect of tumor vasculature being tortuous and non-functional, possibly due to over-proliferation.
There are several limitations of the current model. For one, we consider a tumor space that is constant and unchanging, whereas in reality the tumor would grow over time with the recruitment of new vasculature. Thus, we are focused on a section of a tumor that is recruiting new vasculature but not focusing on the expansion of the tumor. The interplay between the tumor and vasculature should be the subject of further studies. Another limitation of this model is that we have explicitly included the extracellular matrix and interplay with matrix metalloproteinases (MMPs). This is an important process but beyond the scope of this study. Another limitation is that we only examine VEGF-A with the assumption that it consists of isoform VEGF 164 . In the present analysis we do not distinguish between different VEGF isoforms, e.g. VEGF 121 , VEGF 165 and VEGF 189 . However, it is well established that the patterns of vascular network formation depends on the differential isoform expression 7 , the current model will need to be expanded to incorporate these differences.
In conclusion, this model examines the complex interplay between tip cell migration and stalk cell proliferation. Specifically, we demonstrate that proliferation greatly influences total vascular coverage and, thus, may be a more effective target for anti-angiogenic therapies. These findings highlight the utility of computational modeling in parsing the critical components of complex biological processes, such as angiogenesis, to better optimize therapeutic development and the identification of potentially synergistic drug combinations.

Methods
Model Overview. The capillary sprouting model presented here is a hybrid agent-based model that models endothelial cells (EC) as physical agents on a discrete three-dimensional grid. The coordinates of the agents and events such as proliferation and migration occur in continuous space but are mapped to a discrete grid for representation and storage of environmental values. The model takes advantage of the well-developed agent-based and probabilistic approaches to modeling directional angiogenesis in order to produce realistic capillary morphologies under various settings. The general flow of the code is illustrated in Supplementary Figure 1: tip cell migration and anastomosis, stalk cell proliferation and anastomosis, sprouting and vessel regression.
Endothelial Cell Agent. The environment of the in silico model is a discrete 500 × 500 × 500 micron grid.
Each grid coordinate contains environmental variable [VEGF] and may be occupied by an agent. The default VEGF concentration is assumed to be a constant value of 20 ng/ml. The environmental grid represents the 'tumor space. ' We have also examined several different VEGF gradients, see Supplementary Figure 1. Endothelial cells (EC) are the main components of blood vessels and each agent in the model is assigned an endothelial cell type, either tip, stalk, or phalanx, and comprises a forward node (Node 2) and trailing node (Node 1). An agent ahead of another agent would share its trailing node coordinate with the forward node of the agent behind it, thus forming continuous capillaries. This agent definition is an extension of 35 ; along with node coordinates, agents also store values for radius, activation level, phenotype, and cell cycle length. Each endothelial cell is made up of a group of agents. Sprouting can only occur once for each endothelial cell. Initial Capillaries. The initial configuration consists of two capillaries at diagonal ends of the tumor space. Migration. Migration is the primary means by which the growing sprout is able to change directions in response to chemotactic and haptotactic cues. Migration is limited to tip cell agents in this model. Many previous models have modeled angiogenesis as a reinforced random walk or a diffusion process 57 . We model migration the same way as was done previously 35 : where d base is the base endothelial cell migration rate in μ /hr. Thus the tip cell agent migrates by d mig in the direction chosen by the filopodia search. Migration is limited, such that an agent can only extend 50% of its current length. The maximum total length of an agent is limited to d max . The values of parameters are presented in Supplementary Table 1.
The search for the direction of tip cell migration includes persistence. Persistence is defined as where θ is the angle between the current tip cell agent's direction and the direction being considered, and σ is π /6 57 . The direction score of each direction during the search of the local environment is ρ.
where γ is the weight of persistence. Thus is the probability of moving in the direction (i,j,k). In this model, tip cell agents migrate each iteration, attempting to reorient themselves according to equation (3) 58 . It is likely the case that the tip cell blazes a path forward by releasing matrix metalloproteinases (MMP) that degrade the surrounding ECM and by secreting components of basement membrane such as collagens 59 .
Stalk Cell Proliferation. Proliferation is based on the individual cell cycle that must be completed before proliferating. The amount of time it takes to complete a cell cycle is chosen as above. These values are varied based on the cell's proliferation rate.
Proliferation is implemented first by creating a new agent. The forward node of the new tip cell agent is positioned in the same direction as the old tip in the normalized direction of the old tip cell agent. The tip cell agent remains active and the stalk cell agent deactivates. The new agent is assigned a cell cycle duration. Each new agent is part of a larger endothelial cell unit, thus once the total length of the endothelial cell unit is greater than e min the new stalk agent become part of a new endothelial cell unit. The stalk cell agent only proliferates if the new tip agent does not leave the gridspace. After the stalk cell agent proliferates, it determines whether it will anastomose (see section on anastomosis).
Sprouting. Sprouting is determined probabilistically as long as the new sprout will not be on the boundary. A stalk cell agent must be at the beginning of an endothelial cell unit to bifurcate and form a sprout. The trailing node of the new sprout is the position of the forward node of the sprouting capillary agent and the direction of the forward node of the new sprout is determined randomly. This new agent is an activated tip cell agent and the first cell in a new capillary. The tip cell agent then proliferates to form a stalk cell. Once sprouting has occurred the cells attached to the sprouting capillary cell deactivate and are no longer able to sprout. This is due to Delta-Notch signaling in sprouting angiogenesis which has been demonstrated in vivo and in silico 60 .
Anastomosis. Anastomosis can form in one of two ways: 1) when two tip cell agents migrate into one another they can anastomose and form a connection between the two capillaries; and 2) when a tip cell agent migrates into a non-tip cell (either a stalk agent or a phalanx cell) and they fuse. The two different processes are essentially the same with a few minor modifications. Specifically when a migrating tip cell agent approaches another tip cell, the migrating tip agent's forward node is assigned to the location of the forward node of the position of the other tip cell agent, essentially connecting the two tip cells at their ends. This is logical because the filopodia of two tip cells will sense one another and reorient the migrating tip cell to attach to its end. When a migrating tip cell hits any other cell (either a stalk cell agent or a quiescent phalanx cell), the migrating tip cell agent's forward node will be placed at the point in the grid where they hit and the two capillaries will fuse at this point. The tip cells deactivate when they anastomose. Both capillaries that have anastomosed become mature, such that the cells in each capillary are no longer allowed to proliferate or migrate.
Sprout Regression. Endothelial cells can sense the level of VEGF; if a sprout is close to a functional vessel (with blood flow under in vivo conditions), oxygen is delivered locally by the functional vessel that would prevent Scientific RepoRts | 6:36992 | DOI: 10.1038/srep36992 VEGF elevation in the vicinity of the sprout via HIF1α mechanism. We include sprout regression if a newly formed sprout comes close to a functional vessel. Specifically, we assume that if there is a functional vessel within a distance d r from the tip, other than the capillary from which it originated, the sprout will no longer grow and will regress; we will refer to d r as regression distance. The values of d r are varied within a range shown in Supplementary Table 1 to study the sensitivity of the results to this parameter. Two mature capillaries within a distance d r from each other would have overlapping areas of diffusion. So newly formed sprouts within the regression distance d r from another functional capillary would feel a smaller oxygen gradient, causing regression. On the other hand, sprouts that are farther than the distance d r from a mature vessel would sense a stronger gradient and would continue to grow.