Computational Modeling to Quantify the Contributions of VEGFR1, VEGFR2, and Lateral Inhibition in Sprouting Angiogenesis

Sprouting angiogenesis is a necessary process in regeneration and development as well as in tumorigenesis. VEGF-A is the main pro-angiogenic chemoattractant and it can bind to the decoy receptor VEGFR1 or to VEGFR2 to induce sprouting. Active sprout cells express Dll4, which binds to Notch1 on neighboring cells, in turn inhibiting VEGFR2 expression. It is known that the balance between VEGFR2 and VEGFR1 determines tip selection and network architecture, however the quantitative interrelationship of the receptors and their interrelated balances, also with relation to Dll4-Notch1 signaling, remains yet largely unknown. Here, we present an agent-based computer model of sprouting angiogenesis, integrating VEGFR1 and VEGFR2 in a detailed model of cellular signaling. Our model reproduces experimental data on VEGFR1 knockout. We show that soluble VEGFR1 improves the efficiency of angiogenesis by directing sprouts away from existing cells over a wide range of parameters. Our analysis unravels the relevance of the stability of the active notch intracellular domain as a dominating hub in this regulatory network. Our analysis quantitatively dissects the regulatory interactions in sprouting angiogenesis. Because we use a detailed model of intracellular signaling, the results of our analysis are directly linked to biological entities. We provide our computational model and simulation engine for integration in complementary modeling approaches.

Introduction

Angiogenesis is a pivotal process in various aspects of vertebrate life. In development (Breier, 2000) as well as in regenerative processes like wound healing (Flegg et al., 2015) and bone fracture healing (Checa and Prendergast, 2009; Carlier et al., 2015a), angiogenesis is necessary to support newly forming tissue with oxygen and nutrients. Tumors abuse sprouting angiogenesis to direct vascularization toward them (MacGabhann and Popel, 2006). A quantitative understanding of the molecular mechanisms that shape vascular network structure can hence lead to improved treatment options in regenerative medicine and beyond, such as in oncology.

In sprouting angiogenesis, quiescent endothelial cells (ECs) in an existing vessel adopt a motile tip cell phenotype, release matrix metalloproteases (MMPs) to degrade the extracellular matrix (ECM) around the vessel, and lead a sprout followed by proliferative stalk cells (Logsdon et al., 2014). Eventually, the sprouts extend to neighboring vessels and undergo anastomosis and lumenization, followed by pruning and further maturation of the vessel network (Chappell et al., 2011; Tung et al., 2012). Direct cell-cell signaling and signaling via biochemical gradients establish different cell types in angiogenesis. Endothelial cells in emerging vessels can be roughly divided into three different phenotypes: motile tip cells that lead a vessel sprout, stalk cells that are proliferative but are not motile, and phalanx cells that are neither motile nor proliferative but form strong cell-cell junctions in maturating vessels (Carmeliet and Jain, 2011). The distinction between tip cells, stalk cells, and phalanx cells is usually achieved via imaging of different marker genes (Ubezio et al., 2016). Recent data show that the expression of these marker genes is highly dynamic and results in dynamic cell state changes (Venkatraman et al., 2016). Different microenvironments along a sprout can lead to the emergence of tip-like cells behind the sprout tip that eventually overtake the actual tip cell (Jakobsson et al., 2010). The dynamic expression patterns, on the single cell level, can also lead to oscillatory expression patterns of tip cell associated genes on the tissue level that modify the topology of the emerging network in developmental retinal angiogenesis (Ubezio et al., 2016).

VEGF-A induced Delta-Notch signaling drives tip cell selection (Liu et al., 2003; Suchting et al., 2007; Bentley et al., 2009). VEGF-A activates VEGF-receptor 2 (VEGFR2, also called Flk-1), which leads to the expression of tip cell markers and sprouting. VEGFR2 activity also leads to increased expression of dll4, activating Notch1 in neighboring cells through cleavage of the notch intracellular domain (NICD) which in turn represses VEGFR2 expression (and subsequently tip cell markers and dll4 expression) through regulation of transcription factors HES2 and HEY2 (Suchting et al., 2007). This lateral inhibition of tip cell activation induces patterns where tip cells are surrounded by non-tip cells (Hellström et al., 2007; Lobov et al., 2007; Venkatraman et al., 2016). VEGF/Delta-Notch signaling is not isolated, but a wide range of factors can modify tip cell selection (Harrington et al., 2008; Geudens and Gerhardt, 2011; Benn et al., 2017). Different splice variants of VEGF-A and different compositions of the ECM can lead to different diffusion and binding properties and thus influence sprouting (Keyt and Berleau, 1996; Ruhrberg et al., 2002; Ferrara et al., 2003). Additionally, different members of the VEGF family can interact with different receptors (Cao, 2009).

The decoy receptor VEGF-receptor 1 (VEGFR1, also called Flt-1) binds VEGF-A with a higher affinity than VEGFR2, however leads to negligible activation of downstream targets (Waltenberger et al., 1994). In spite of this, VEGFR1 plays an important role in angiogenesis, as disruption of VEGFR1 leads to embryonic lethality through vascular overgrowth (Fong et al., 1995). The presumed role of VEGFR1 in angiogenesis is to reduce VEGFR2 activity by reducing local VEGF-A availability (Hiratsuka et al., 1998; Roberts et al., 2004). A soluble splicing isoform, sVEGFR1, is secreted and can diffuse through extracellular space to reduce VEGF-A availability in a radius around a secreting cell (Kendall and Thomas, 1993; Roberts et al., 2004). Transcriptional regulation of VEGFR1 is contrary to the regulation of VEGFR2, namely activated by Dll4-Notch1 signaling (Harrington et al., 2008; Funahashi et al., 2010) so that VEGFR1 is mainly expressed in stalk cells and inhibited in tip cells while VEGFR2 is a tip cell marker and suppressed in stalk cells. Accordingly, the general mechanisms by which VEGFR1 regulates angiogenesis have been identified, but their interdependencies with activatory VEGFR2 signaling and lateral inhibition via Delta-Notch are difficult to collectively assess experimentally.

VEGF-A has been implied as a target in several clinically relevant conditions, e.g., ischemia, arthritis, or obesity (Carmeliet, 2003), peripheral artery disease (Clegg et al., 2017), or cancer (MacGabhann and Popel, 2006; Zhang et al., 2010). Anti-VEGF drugs have been very successful in the treatment of age-related macula degeneration and diabetic macula oedema (Virgili et al., 2012; Solomon et al., 2014). To maximize the effectiveness of anti-VEGF drugs in cancer therapy, however, they need to be combined with drugs targeting different pathways (Jain et al., 2006; Carmeliet et al., 2009; Incio et al., 2018; Wagner et al., 2018). Optimization of such combinatorial therapies requires detailed knowledge of the dynamics within and between pathways. Computational studies can generate this knowledge efficiently and can also be used to predict optimal treatment regimens (Barros de Andrade e Sousa et al., 2016). Hence an improved understanding of the dynamic interplay of VEGF pathway components can contribute to improved treatments against a wide variety of conditions.

A considerable number of computational models of sprouting angiogenesis exist (Merks et al., 2004; Bentley et al., 2008, 2009; Jakobsson et al., 2010; Carlier et al., 2012, 2014, 2015b; van Oers et al., 2014; Boas and Merks, 2015; Heck et al., 2015; Walpole et al., 2015; Ubezio et al., 2016; Venkatraman et al., 2016; Bentley and Chakravartula, 2017). Most of these models are based on the “memAgent-Spring” (Bentley et al., 2008), Cellular Potts (Merks et al., 2004), or agent based (Carlier et al., 2014) modeling approaches.

The “memAgent-Spring Model” describes the alignment of membrane patches assigned to specific cells on a pre-defined shape, each patch can have its own dynamics concerning signaling (Bentley et al., 2008). Models using this approach have been used to explain, for example, tip cell overtaking (Jakobsson et al., 2010) and oscillations in lateral inhibition (Ubezio et al., 2016). Although this kind of model can be used to describe tip selection and intercellular signaling in a sprout in high detail, the restriction to a static shape on which the cells are aligned makes it impossible to simulate vascular network formation.

In Cellular Potts models, cells are described by a collection of nodes on a lattice. Cell shapes and movement arise dynamically from re-assigning nodes to minimize an energy function (Merks et al., 2004). This energy function can contain various terms, e.g., to constrain cell area or perimeter or bias movement to a certain direction. Cellular Potts models can describe network formation from individual cells (Köhn-Luque et al., 2013) and have also been used to explain tip cell overtaking (Boas and Merks, 2015). A drawback of Cellular Potts models is that terms of the energy function can be added arbitrarily and are often difficult to link to biological mechanisms. In addition, the computational cost renders simulations of vessel networks in 3D currently unpractical.

Agent based models in which single cells are represented by single nodes on a lattice have been used to describe blood vessel formation in bone healing scenarios (Checa and Prendergast, 2010; Carlier et al., 2015b). Although they omit dynamics of cellular shape, they can be used to simulate collective cellular organization (Checa et al., 2015) and computational complexity is low enough to permit simulation of network formation in 3D. Existing agent based models of angiogenesis ignore or strongly simplify biochemical signaling (Checa and Prendergast, 2010; Carlier et al., 2014, 2015b). In these models, the tip cell phenotype is either enforced (Checa and Prendergast, 2010) or emerges as a result of simple rules related to contact inhibition (Carlier et al., 2014, 2015b). In the latter, model predictions overestimate vessel growth.

Computational models have also been used to explore the role of VEGFR1 in angiogenesis. Most of these models focus on the establishment of VEGF-A and VEGF receptor gradients and binding kinetics (MacGabhann and Popel, 2004; Wu et al., 2010; Hashambhoy et al., 2011; Chappell et al., 2016). Walpole et al. (2015) described cellular movement, but not sprouting as a whole, in a Cellular Potts-like framework. To date, there is no computer model that integratively investigates the signaling interactions of VEGFR2, VEGFR1, and Delta-Notch during sprouting, their interdependencies and their effects on the forming vascular network.

Here, we present a computational agent based model to quantitatively dissect the interrelation of the regulatory mechanisms in tip selection and sprouting. We focus on Delta-Notch and VEGF signaling including VEGFR1 and VEGFR2. Our model allows us to assess the experimentally inaccessible interdependencies of intracellular signaling and extracellular conditions. Our analysis shows that VEGFR1 is efficient in guiding sprouts away from existing vessels and it also highlights the importance of Delta-Notch signaling, specifically the degradation of NICD, for angiogenesis. We provide our model in an open and reproducible way, thus facilitating integration into different contexts.

Materials and Methods

Agent Based Model Simulator

We use a custom simulator, AngioABM, for our agent based modeling, implemented in C++11 using the boost libraries (The Boost Community, 2017). The source code is available under the Apache License 2.0 at https://gitlab.com/ModularABM/AngioABM/tree/VEGFR1. Scripts for visualization and analysis are available at https://gitlab.com/ModularABM/ABMTools/tree/VEGFR1.

In AngioABM, an agent based model consists of global variables, local variables and agents placed on a discrete 2D grid. Global variables are single numerical values defined for all positions of the grid. Local variables are 2D matrices of the same dimensions as the grid. Agents contain an id, coordinates, internal variables, and update rules. For clarity, we will refer to local variables using monospace font and to agents internal variables in italics. AngioABM reads model descriptions in XML.

AngioABM uses discrete time steps for time course simulations. At each time step, update and output functions of all local variables and all agents are called iteratively as shown in Figure 1.

FIGURE 1

Figure 1. Flowchart of the agent based model simulations. Colors indicate processes/decisions associated with simulation time (blue), local variables (yellow), agents' intracellular dynamics (turquoise), and agent movement and proliferation (green).

Local variables represent molecules (in particle numbers per lattice site). Basic functions for diffusion, degradation, and adding molecules are implemented in a base class Localvar. These functions can be used to implement derived classes of local variables that exhibit custom dynamics. AngioABM uses forward integration to compute local variables' temporal and spatial evolution. The simulator assumes Dirichlet boundary conditions.

Agents represent cells. Basic functions for agent movement and interaction with the environment are implemented in a base class Agent. Derived classes representing specific cells can use these functions, lists of internal variables and parameters and can additionally use custom rules or formulas to define cell dynamics. Agents are updated asynchronously in random order. Hence, changes in one agents' state can influence another agents' dynamics immediately. This reduces numerical errors arising from the discrete time intervals between agents' information exchange and induces stochasticity into the simulation. To reduce errors in integration, especially in the computation of local variables, numerical stability is checked and internal steps are fine grained accordingly. Agents' parameters can be specified to refer to a fixed value for all agents or to be sampled from a random distribution for each agent. Drawing agents' parameters from random distributions induces cellular heterogeneity.

Agent Based Model of Sprouting Angiogenesis

The model described here consists of the local variables VEGF-A, sVEGFR1, and sVEGFR1b (sVEGFR1 bound to VEGF-A) and agents representing endothelial cells (ECs). VEGF-A, sVEGFR1, and sVEGFR1b diffuse, are degraded and can associate and dissociate.

We assume the 2D grid considered in simulations to be a subarea of the bottom of a well plate containing culture medium as described in Kappas et al. (2008). We simulate vascular endothelial cells on this grid, originating from two initial cell aggregates representing embryoid bodies. We introduce a VEGF-A influx to each grid cell that depends on the difference between a global reference value (representing the average value in the whole volume of culture medium) and the local value. The reference value is maintained constant as the culture medium in experiments was renewed every 48 h (Kearney and Bautch, 2003).

We model the EC agents according to the intracellular signaling model described in the next section. The EC agents extend or retract filopodia based on their active VEGFR2 (VEGFR2a, see Equations S10, S36, and S37). The filopodia determine the radius in which an EC agent senses and binds to VEGF-A. VEGF-A binds to mVEGFR1 and VEGFR2 according to the respective association and dissociation reactions (Equations S28, S32), sVEGFR1 is secreted (Equation S25). An EC cell becomes a tip cell when it satisfies tip cell criteria analogous to criteria described in Blanco et al. (2013):

• VEGFR2_mRNA>VEGFR1_mRNA,

• filopodia>A*,

• $dll{4}_{mRNA}>D^{*}$

with threshold parameters A* and D*. A tip cell determines the direction of the VEGF-A gradient and attempts to move into that direction with a chance of deviating clockwise or counter-clockwise. If no gradient exists, tip cells move into a random direction. The integrity of the sprout is maintained by choosing a random non-tip neighbor when moving. As stalk cells are proliferative, a copy of this neighbor is placed on the tip cells' old position.

Intracellular Signaling Dynamics

The model describes interactions relevant for tip cell selection and lateral inhibition. It consists of 21 variables and 43 parameters that describe rate laws for binding, transcription, translation, and degradation of molecules, see Figure 2 for an overview. In the agent based model, we use a forward Euler algorithm to solve the corresponding differential equations. The complete model of intracellular interactions and a list of model equations are available as Supplementary Material. The model is provided in the Systems Biology Markup Language (SBML) (Hucka et al., 2003), an XML based format for the exchange of computational models of biological processes. The SBML model has been submitted to the BioModels Database (Chelliah et al., 2015) and assigned the identifier MODEL1804030001.

FIGURE 2

Figure 2. Model topology indicating signaling systems in different colors (green: Dll4-Notch1, turquoise: VEGFR1, purple: VEGFR2). dll4₂ indicates Dll4 from neighboring cells. Symbols adhere to the Systems Biology Graphical Notation format (Novère et al., 2009): filled arrows represent reactions, empty arrowheads represent activation, barred arrows represent inhibition, shaded symbols (dll4, VEGF) indicate the same entity in different processes. Circles crossed by a bar represent sources and sinks of molecule production.

Binding of Notch1 to Dll4 from neighboring cells (Dll4_nb) releases the active Notch1 intracellular domain NICD. NICD induces transcription via different downstream factors (Hey2/Hes2). We omit intermediate steps, so NICD directly activates transcription of notch1_mRNA and VEGFR1_mRNA while repressing expression of VEGFR2_mRNA. We model transcription rates v_{transcription} by using Hill-Equations (reviewed, e.g., in Tummler et al., 2014) with additional terms for the activation by an activator mod:

$\begin{array}{cc}\begin{array}{r}{v}_{transcription}=V_{max}·(\frac{1}{a}+(1-\frac{1}{a})\frac{mo{d}^h}{mo{d}^h+M_0^h})\end{array}& (1)\end{array}$

where a indicates the fold increase in transcription upon activation by mod, V_max is the maximal transcription rate upon full activation, M₀ determines the modifier value at which v_{transcription} is half its maximal value and h determines the slope of the function. Generally we assume h = 2.

All other reactions (translation and degradation, protein binding, and dissociation) follow standard Mass-Action kinetics. VEGFR1 is regarded in two isoforms, membrane-bound VEGFR1 (mVEGFR1) and soluble VEGFR1 (sVEGFR1). mVEGFR1 and sVEGFR1 both bind to extracellular VEGF-A, mVEGFR1 directly, sVEGFR1 only after export. Binding to VEGF-A leads to inactive complexes that can diffuse, degrade, or dissociate. VEGFR2 binding to VEGF-A is converted to VEGFR2bound in the model. After dissociation of VEGF-A, VEGFR2 is converted to VEGFR2a, also representing the activity of the downstream signaling cascade. Both VEGFR2bound and VEGFR2a enhance transcription of dll4_mRNA (omitting intermediate steps in transcriptional regulation). We describe transcription and translation as individual processes for each protein, which also generates time delays.

The model also includes cis-inhibition of Dll4-Notch1 signaling since it has been extensively discussed in literature to be necessary for lateral inhibition (Sprinzak et al., 2010; Shaya and Sprinzak, 2011). Additionally, transcriptional auto-activation of notch1 (Boareto et al., 2015) has been implied in supporting switch-like behavior by establishing a positive feed-forward loop on the single cell level (Tyson et al., 2003) and we have implemented according reactions in our model (see Equations S2, S16, S17).

Model Parameterization

We obtained ranges for parameter values from literature (MacGabhann and Popel, 2006; MacGabhann et al., 2006; Bentley et al., 2008; Hashambhoy et al., 2011; Carlier et al., 2012; Imoukhuede et al., 2013; Boareto et al., 2015; Venkatraman et al., 2016). wherever possible (fold increase in transcription of Dll4 upon VEGFR2 stimulation, the binding affinity of VEGFR1 to VEGF-A relative to the binding affinity of VEGFR2 to VEGF-A, abundance of specific proteins in ECs). Where no specific data was available, we used ranges reported in Schwanhäusser et al. (2011) (ranges for transcription, translation, and degradation rates). In total, 17 initial values have been set, 41 parameters pertaining to intracellular processes have been estimated and 9 parameters pertaining to extracellular processes have been estimated based on 4 data points (relative vessel area and branch points per vessel length for wild type and VEGFR1 -/-) from published experimental results (Kappas et al., 2008). Although parameter boundaries were set based on additional literature, we cannot eliminate the possibility that alternative biologically feasible parameterizations could also reproduce the experimental data.

We parameterized the model to reproduce experimental data from Kappas et al. (2008), namely the relative vessel area ($A_{rel}^{EC}$) from mouse embryoid bodies grown in culture wells after 8 days for WT and VEGFR1 -/- mutants. For an overview of published parameter values, values used in this study, and parameter boundaries used, see Table S3 and the notes in the Supplementary SBML model.

To generate the parameterization for the whole model, we first estimated parameters for the intracellular ODE model to generate tip and stalk phenotypes using Copasi (Hoops et al., 2006), then extended the ODE model to describe three cells next to each other and estimated parameters based on the previous iteration to generate lateral inhibition. We used the resulting parameter values as initial guesses for estimating parameters for the whole agent based model. Automatic parameter estimation using a genetic algorithm (Spiesser et al., 2015) was not successful, so we optimized parameter values manually to reproduce experimental data. For this, we iteratively modified individual parameter values and then computed a set of 5 simulations for WT and VEGFR1 -/- followed by computing $A_{rel}^{EC}$ and comparison to experimental data. Model fitting was performed forcing all parameters to remain within the range of values reported in literature, which are indicated in Table S3.

For initial values and parameters related to the general condition of a cell (V_max of transcription and translation) and the NICD degradation rate (k_nicd_degradation) we modified not only mean values, but also the deviations of the log-normal distributions from which they are sampled. Parameters related to entity specific characteristics (e.g., binding affinities or fold induction in gene regulation) are not sampled since we assume that the general state of cells is more variable than the entity specific characteristics.

Sensitivity Analysis

We analyzed the sensitivity of the simulation results with regard to changes in parameters. We used the percentage of the area occupied by ECs ($A_{rel}^{EC}$) as the model output to analyze. For parameter sensitivity analysis, we consider variations in single parameters only. To analyze the sensitivity of the model output to changes in single parameters, we performed 20 simulations per parameter value over a range of parameter values centered around the value originally used. This has been done for 46 parameter values between 75% and 125% of the original value of each parameter (except for the exponents in the Hill equations). Results are available as Figures S2–S47. For specific parameters, we repeated the analysis using 51 different parameter values from 1 to 199% of the original value (depicted in Figures 4, 5). This analysis depends on the parameterization chosen and the effects of changes in two (or more) parameters can not be predicted from combining the effects of two individual parameter changes.

Branching Analysis

To measure branch points, we saved simulation results as image files and analyzed these with ImageJ-MATLAB (Hiner et al., 2017). Images were smoothened to ensure connectivity of agents, then binarized and skeletonized. The skeletonized path was analyzed for branch points using AnalzyeSkeleton (Arganda-Carreras et al., 2010). Branch points per millimeter vessel length were computed using the total length of all branches from the ImageJ analysis.

Results

In-silico Model Reproduces Experimental Data on Wild Type and VEGFR1 -/- Strain

The constructed model with the described parameterization reproduces experimental data on mouse embryoid bodies of wild type and VEGFR1-/- strains as described in Kappas et al. (2008). Over 20 simulation runs for wild type conditions, we observed distinct vessels in a network topology and ECs covered 17.2% ± 5.8 (mean ± standard deviation) of the simulated area, compared to 16.7% reported for experiments. In 20 simulations for VEGFR1-/- conditions, nearly no distinct structures were apparent and ECs covered 55.1% ± 10.5 area (55.34% in experiments). Branch points/mm are also comparable to experimental data (9.9 ± 2.2 branchpoints/mm predicted, 12.6 ± 3 branchpoints/mm). Figure 3 shows representative simulation results for the wild type and the VEGFR1-/- strain, comparable to the experimental results shown in Kappas et al. Kappas et al. (2008). Movies of the simulations are given in Supplementary Material.

FIGURE 3

Figure 3. Representative simulation results: extracellular VEGF-A, and VEGFR2 mRNA for VEGFR1 -/- strain (top) and wild type (bottom) at 60, 90, 120, and 199 h (from left to right).

Initial endothelial cells in the model are seeded at two locations on the grid, corresponding to the seeding of multiple embryoid bodies (Kappas et al., 2008). Simulations with one initial cell group in the center of the grid yielded similar means for the relative area covered by ECs ($A_{rel}^{EC}$) in the WT and VEGFR1-/- strain but substantially wider standard deviations (see Supplementary Material for details and model files).

The transition between tip and stalk cells is smooth. Although the mean values of intracellular variables strongly differ for tip and stalk cells, they overlap: non-tip cells exist for which some intracellular variables have similar values as generally observed in tip cells. Exemplary distributions are given in Figure S1.

Contributions of VEGF-A Receptors

To analyze the contribution of the different VEGF-A receptors to sprouting, we evaluated the sensitivity of $A_{rel}^{EC}$ to changes in transcription and translation rates of the receptors. The parameters analyzed are the rate parameter of translation and the V_max parameter of Equation (1) for transcription.

Concerning V_max of transcription for VEGFR2, we observe that too little VEGFR2 transcription renders cells insensitive to VEGF-A while too much VEGFR2 leads to a hyperactivation of Dll4 and thus too high NICD levels that suppress sprouting. This is true for both VEGFR2 transcription (Figure 4A) as well as for VEGFR2 translation (Figure 4B). As shown in Figure 4A, the sensitivity of $A_{rel}^{EC}$ to VEGFR1 transcription displays a strong inverse relationship. Increased VEGFR1 transcription causes a reduction in tip cell phenotypes (a tip cell needs to satisfy VEGFR2_mRNA>VEGFR1_mRNA) and a reduction in available VEGF-A through binding of VEGF-A to VEGFR1 isoforms. For a complete knockout of VEGFR1, the $A_{rel}^{EC}$ decreases from 69 to 51%. This is due to the directional cues that VEGFR1 provides. Under complete VEGFR1 knockout, vessels grow in random directions and cannot fill space as effectively as if they would grow more orthogonally.

FIGURE 4

Figure 4. Parameter sensitivity: changes in relative vessel area ($A_{rel}^{EC}$) after 192 h for variation in parameter V_max for VEGFR2 and VEGFR1 transcription (A), and the translation rate of mVEGFR1, sVEGFR1, VEGFR2 (B).

Figure 4B shows the effect of changes in the rate of translation of VEGF receptors. The effect for VEGFR2 is comparable to that of changing transcription rate (compare Figure 4A), only that the range of parameter values that allows for sprouting is narrower. An effect of changing the translation of mVEGFR1, the membrane bound form of VEGFR1, is not detectable. There is a weak positive correlation between the translation of sVEGFR1 and the relative vessel area. The sensitivity of branch points/mm exhibits similar trends as the sensitivity of the relative vessel area (Figures S49–S52).

Lateral Inhibition via Dll4-Notch1 Signaling

We have also analyzed the translation parameters of Dll4 and Notch1 (Figure 5A). The translation of Dll4 is anti-correlated to $A_{rel}^{EC}$ since an increase in Dll4 leads to an increase in trans activation of Notch1 which in term suppresses VEGFR2 transcription and thus prevents sprouting. The sensitivity to changes in the transcription rate of dll4 (Figure 5A) does not resemble the sensitivity for the translation rate. This is due to the tip cell criterion $dll{4}_{mRNA}>D^{*}$, dll4 mRNA has to exceed a threshold for sprouting. Accordingly, there is a positive correlation between $A_{rel}^{EC}$ and dll4 transcriptional rate until the trans activation of Notch1 outweighs the tip cell activating influence of dll4, at around half the original parameter value. For higher values of dll4 transcription, the sensitivity resembles that of Dll4 translation. The same trends hold for branch points/mm (Figures S51),

FIGURE 5

Figure 5. Parameter sensitivity: changes in relative vessel area ($A_{rel}^{EC}$) after 192 h for variation in the translation rates of Notch1 and Dll4, and the transcription rate of Dll4 (A), and the degradation rate of notch1 intracellular domain and the parameters for Dll4-Notch1 association in cis (kc) and trans (kt) (B).

The translation of Notch1 shows a different pattern: For V_max of Notch1 translation below half the original value, no trans activation occurs and hence no inhibition of VEGFR2 occurs and vessels fill the complete area. Between half the original parameter value to almost the original value, lateral inhibition is sufficient to lead to VEGFR1 expression, but not sufficient to induce patterning and sprouting. The increase in relative area beyond the original value indicates functional lateral inhibition and patterning that reach a maximum of about 30% $A_{rel}^{EC}$. Further increase of Notch1 production makes it more unlikely to overcome tip cell inhibition, so that the relative EC area decreases. The sensitivity to transcriptional regulation of notch1 is similar to that of translational regulation (see Figures S40, S45), because notch1 mRNA is not involved in tip cell phenotype conditions, as are VEGFR1, VEGFR2, and dll4 mRNA. Again, the branch points/mm over the parameter ranges follow the trends in relative vessel area (Figure S51).

Concerning the interactions in Notch1-Dll4 signaling, we have analyzed the sensitivity to parameters kc, kt, and the degradation rate of NICD, depicted in Figure 5B. Parameter kc resembles the binding rate of Dll4 to Notch1 on the same cell, leading to cis inhibition of Notch1. Parameter kt describes the binding rate of Dll4 to Notch1 on a neighboring cell, leading to trans activation of Notch1. The trans activation determined through kt is necessary for lateral inhibition and its absence leads to hypersprouting. Conversely, a value of kc above 0.154 is necessary to reduce lateral inhibition and induce sprouting.

The amount of NICD, determined by Notch1 activation and NICD degradation, is crucial for balancing VEGFR1/VEGFR2 levels and thus determine the degree of lateral inhibition (as indicated in Figure 2). As shown in Figure 5, NICD degradation needs to surpass a threshold (0.385 of the estimated value) to allow sufficient production of VEGFR2 and thus allow sprouting at all. If NICD degradation exceeds another threshold (1.192 of the estimated value), VEGFR2 transcription is not reduced sufficiently, lateral inhibition ceases and no tip cells are formed. For values of NICD degradation above 1.692 times the estimated value, the transcriptional activation of VEGFR1 is sufficiently reduced to induce sprouting again, but in the absence of lateral inhibition.

For the Notch1-Dll4 signaling interactions, parameter variations lead similar trends for branch points/mm as for the relative vessel area (compare Figure S52).

For most parameters in the Dll4-Notch1 signaling system, the maximal value of $A_{rel}^{EC}$ is about 30%, because changes in Dll4-Notch1 signaling affect the timing of the onset of sprouting in our simulations, not the overall number of sprouts, which is more strongly driven by the availability of active VEGFR2 (compare also Figures 6B,C).

FIGURE 6

Figure 6. Parameter distributions strongly influence model dynamics. The standard deviation of the parameter distributions were varied from 0.01 to 2.56·μ. Sampled parameters are: V_max of transcription (tc), translation rates (tl), and the degradation rate of NICD (nicd_deg) or all of the above. Initial values for new agents were sampled accordingly under all scenarios. Model output for 20 simulations of each setting have been recorded as $A_{rel}^{EC}$ (A), the maximal number of active tip cells (B), and the time of the first tip cell (C). If no tip cells occur during a simulation run (199 h), we set the time to the first tip cell to 1,000 h. Plotted values are given in Table S2.

Influence of Cell-to-Cell Variability

To assess the influence that the sampling of parameters has on simulation results, we have analyzed simulation runs using distributions with similar means but different standard deviations, as depicted in Figure 6. Generally, low variability leads to low $A_{rel}^{EC}$ and late emergence of tip cell phenotypes, often without any tip cells during the simulated time course. Here, the high homogeneity of the cells and the symmetric initial conditions impose highly similar conditions on each agent so that emergence of tip cells depends solely on the stochasticity induced through the asynchronous agent updates. Increasing the width of the random distributions for transcriptional and translational parameters leads to an increase in sprouting up to a ratio $\frac{SD}{Mean}=0.64$, beyond which both Dll4-Notch1 and VEGF-A signaling do not function properly anymore and tip cell phenotypes emerge randomly. This leads to an almost instantaneous emergence of the first tip cell for $\frac{SD}{Mean}>1$ (Figure 6C) and $A_{rel}^{EC}$ at or below the maximal value reached (Figure 6A).

For the degradation of NICD, an increased variability leads to a quite stable $A_{rel}^{EC}$ as well as a stable maximal number of tips. This might seem surprising because the previous analysis showed that sprouting only occurs in a very tight range of values for the NICD degradation rate, as shown in Figure 5B. The results of the sensitivity analysis however impose the same parameter value on all cells in the model while the variability analyzed here is between cells. The degradation rate of NICD critically determines a cells chance of becoming a tip cell, especially if values of neighboring cells differ and thereby reinforce lateral inhibition.

The variability of the model output does not increase additively when sampling transcriptional rate parameters, translational rate parameters and the NICD degradation rate (“all” in Figure 6). Instead, the model output becomes more stable against increasing parameter variability compared to the “tc” and “tl” groups for intermediate ratios of $\frac{SD}{Mean}$, because pattern formation is reinforced by the stronger lateral inhibition caused by the variability in the NICD degradation rate. For higher rates of $\frac{SD}{Mean}$, the signaling fidelity breaks down and tip cells arise almost randomly as in the “tc” and “tl” groups, slightly stabilized by the influence of the variability in NICD degradation that has a stable influence over all ratios of $\frac{SD}{Mean}$.

Again, the number of branch points by vessel length follows the trends observed in the $A_{rel}^{EC}$.

Discussion

Angiogenesis is a multi-scale process where multiple signals must be integrated to enable cellular decision making leading to the sprouting of new vessels. Here, we present an agent based computational model of sprouting angiogenesis that focuses on intracellular signaling. We show that the release of the decoy receptor sVEGFR1 has two distinct effects. By reducing VEGF-A concentrations in the vicinity of stalk cells, VEGFR1 (1) reduces sprouting and (2) directs sprouts to grow orthogonal to existing structures. The second effect leads to a faster vascularization of an avascular area under homogeneous VEGF-A conditions. We observe this impact of VEGFR1 under a wide range of parameters. Our results are consistent with previous findings on the interaction of VEGF and sVEGFR1 in a static environment (Hashambhoy et al., 2011; Chappell et al., 2016), but our study extends on this by considering a dynamic environment in which production rates as well as cell positions are not fixed but emerge from model dynamics.

Our analysis also reveals that the system dynamics, especially the timing of sprouting, strongly depend on the parameterization of the Dll4-Notch1 signaling system that determines lateral inhibition. Our model predicts that the stability of the notch intracellular domain strongly influences the timing of sprouting without a strong impact on vessel network architecture.

We furthermore show that cellular heterogeneity is an important driver of pattern formation. The case of NICD in our model highlights that varying cell-to-cell variation of a parameter has different effects than varying the magnitude of the same parameter. This also has implications for experimental research, because it is accordingly not sufficient to know the mean value of some parameter in an ensemble of multiple cells but we also need to know its variability specified by the distribution and its moments.

Limitations of the Approach

The model presented here is a 2D lattice based model with fixed cell shapes. This means that its predictions are limited to 2D environments and cannot be directly transferred to clinical conditions, e.g., in tumor growth.

Finding the appropriate balance of model detail or resolution is a fundamental challenge in modeling. Any increase in model complexity increases the amount of experimental data necessary to find a reliable parameterization. Here, we have not included all molecules involved in VEGF signaling in sprouting angiogenesis (see e.g., Ferrara et al., 2003; Cao, 2009 for reviews of involved mechanisms), but have attempted to limit our model to the most important mechanisms with respect to the interplay between VEGF, VEGFR1 and VEGFR2. Nevertheless, the model contains 58 parameters, some of which are not fixed values but sampled from random distributions. For none of the parameters, we could find a reliable value in literature that is certainly applicable to the experimental conditions the model describes. Rather, we use the available experimental data (relative vessel area and branch points per vessel length for WT and VEGFR1 -/-) to generate a parameterization that is within sensible boundaries based on literature data.

This parameterization is not guaranteed to be equal to the physiological parameterization. Because we have restricted parameter boundaries based on experimental data and our data reproduces findings reported in literature (Hashambhoy et al., 2011; Chappell et al., 2016), we expect the model analysis to yield at least qualitatively correct predictions.

VEGFR1 in Sprouting Angiogenesis

Kearney et al. (2002) have shown an increase of relative EC area and a loss of vascular network properties in VEGFR1-/- mutants in in-vitro cultures of mouse embryonic bodies. These defects can be restored by additional expression of sVEGFR1 (Roberts et al., 2004), but not mVEGFR1 Kappas et al. (2008), and are VEGF-dependent (Chappell et al., 2009). The proposed mechanism of action is that sVEGFR1 generates a tight corridor of VEGF-A to guide the tip cell away from stalk cells.

The analysis of our model indicates that VEGFR1 suppresses sprouting by reducing the global VEGF-A concentration with increasing cell density (compare the upper and lower row of simulation results in Figure 3) and that VEGFR1 provides spatial cues by reducing the local VEGF-A concentration around stalk cells, which improves vessel network topology. These spatial cues lead to nearly orthogonal sprout growth instead of random sprout growth that we observe when VEGFR1 transcription is removed completely. This also explains the decrease in relative vessel area when removing VEGFR1 transcription completely in Figure 4A. These spatial cues are also the reason for the weak positive correlation between sVEGFR1 translation and relative vessel area in Figure 4B, indicating that sVEGFR1 has a stronger effect than mVEGFR1. Although our model implies that these effects contribute under a wide range of parameters, the importance of either of these mechanisms depends on the external conditions. Under wound healing conditions with cellular sources of soluble VEGF-A, the importance of VEGFR1 will be higher than in the developing mouse retina, where new vessels follow the tracks of astrocytes supplemented with ECM-bound VEGF-A (Gerhardt et al., 2003).

Diffusion properties and morphogen gradients are notoriously difficult to quantify in tissues experimentally (Bothma et al., 2010; Köhn-Luque et al., 2011). Our computational analysis of the role of sVEGFR1 provides important insights into the regulatory capabilities of this mechanism and opens up new questions: In contrast to the static environment considered in our simulations, how does the dynamic environment provided by the ECM influence the effects of VEGFR1? This question is of specific importance to clinical settings, e.g., regeneration or tumor growth, where multiple dynamic sources of VEGF-A might be present and the ECM is undergoing active rearrangement.

Also, our analysis indicates that only the complete removal of VEGFR1 transcription abolishes the spatial cues provided by sVEGFR1 (Figure 4A). This is an artifact of our model, because our in-silico cells sense any difference in external VEGF concentrations, unlike in-vitro or in-vivo situations, where cells sense the VEGF-A gradient by differential binding of VEGF-A to receptors on their filopodia.

Lateral Inhibition and Cellular Heterogeneity

In contrast to the gradual changes in simulation output caused by changes in VEGFR1 or VEGFR2 production rates, changes in parameters involved in Delta-Notch1 signaling cause stronger and more abrupt changes over certain parameter value ranges (Figures 4, 5). This indicates that the fine tuning of parameters in lateral inhibition is more important in angiogenesis than the actual VEGF-A receptor kinetics, according to our analysis.

In the model we present, the base state of a cell is a quiescent stalk-cell like state, and the basic activity of lateral inhibition is relatively high. For sprouting, cells have to overcome this basal lateral inhibition. Cellular heterogeneity through sampled parameters is an important driver of the differentiation of individual cells in our model. Although the sampling is random, a specific cell can be primed for sprouting by the sampled parameter values. Neighboring cells might be assigned parameter values that make them less likely to sprout and the difference between these cells is then further amplified through lateral inhibition. The most important parameter to generate heterogeneity for sprouting is the degradation rate of NICD. Even small stochasticity in NICD degradation increases cellular heterogeneity in the lateral inhibition pathway and significantly improves sprouting efficiency. Because the stability of NICD does not directly affect sprouting but the heterogeneity between neighboring cells, larger stochasticity in NICD degradation does not lead to a breakdown of the system, as observed when sampling transcription or translation parameters from wide distributions (Figure 6). The importance of the regulation of NICD stability for signaling has been discussed before (see e.g., Kopan and Ilagan, 2009; Andersson et al., 2011; Herbert and Stainier, 2011; Bray, 2016), mechanisms include different cleavage products and active targeting of NICD for proteosomal degradation, acetylation by SIRT1 or interaction with other signaling pathways via YAP or SMAD3.

Beyond this specific finding, our analysis on the role of parameter distributions in computational models of multicellular ensembles also has a general relevance. Assuming a fixed value for NICD degradation for all cells and varying this value showed only a narrow band of parameter values that lead to efficient sprouting (Figure 5B). Sampling the value of NICD degradation for each cell from a random distribution improves sprouting efficiency (Figures 6A,B). Our simulations show that pattern formation is facilitated by cell-cell heterogeneity. It is not surprising that this effect is most pronounced for a parameter at the core of lateral inhibition, the main driver of pattern formation under the conditions used here. Generally, pattern formation can be reinforced or damped by cellular heterogeneity and this heterogeneity needs to be considered, in computational and experimental studies alike. Beyond transcriptional noise and asymmetric cell division (Costa et al., 2016), we need to consider that cell-cell and cell-matrix interactions also modify and influence cellular signaling (Qazi et al., 2017) and even the overall cellular state (Battich et al., 2015).

Computational Modeling of Sprouting Angiogenesis

The agent based model we present here omits cell shape dynamics and ensures vessel connectivity by dividing stalk cells upon tip movement. This neglects cell shape changes and their effect on network formation (Boas and Merks, 2015). Nevertheless, the model presented here considers biochemical signaling with a higher level of detail than most models of sprouting angiogenesis (Merks et al., 2004; Bentley et al., 2008, 2009; Checa and Prendergast, 2010; Jakobsson et al., 2010; Carlier et al., 2012, 2014; Köhn-Luque et al., 2013; van Oers et al., 2014; Boas and Merks, 2015; Walpole et al., 2015; Ubezio et al., 2016; Venkatraman et al., 2016; Bentley and Chakravartula, 2017). The introduction of separate transcription/translation steps in our model—although it still omits intermediate steps in transcriptional regulation—introduces a dynamic delay between transcriptional activation and protein production that is omitted in most models. Because of the more faithful biochemical model, experimental perturbations can be more directly linked to model entities. The coarse spatial resolution of our model can be interpreted as a different level of magnification from which the biological dynamics are observed.

The complete model code can be downloaded and reused and we hope that further development of the model by us and other groups will help to further understand the dynamics of sprouting angiogenesis. Future development of the model could focus on the interactions between chemical and mechanical signals or in the cross-talk between different signaling pathways.

In order to generate reliable and clinically relevant predictions of angiogenesis to improve regeneration, enable fine-tuned tissue engineering or ablate tumor-dependent angiogenesis, we need to integrate various mechanisms in a dynamic environment. Moreover, these models need to be parameterized based on dedicated experimental data. By providing our model in an open and reproducible manner, we enable its integration into models that consider cell shape dynamics, interactions with the ECM, or network formation in 3D.

Conclusion and Outlook

We present an agent based computer model of sprouting angiogenesis focusing on intracellular signaling and its effects on pattern formation. Specifically, we focus on the role of VEGFR1 and its interrelation with VEGFR2 and Delta-Notch signaling. Our findings do not replace published modeling efforts but are complementary because our focus differs from existing literature. Our simulations go beyond comparable models concerning level of detail of intracellular signaling, cell numbers considered, and simulation time.

Our results support the hypothesis that soluble VEGFR1 provides spatial cues to guide sprouts away from their origin and it predicts that this effect is robust over a wide range of parameters. Our model analysis shows that lateral inhibition, and especially the regulation of NICD stability and its variability across cells, is critical for the regulation of sprouting angiogenesis.

This leads us to argue that complementary modeling approaches need to combined and supplemented with dedicated experimental data to generate reliable predictions, to generate models that do not only reproduce training data but are also capable of predicting validation data correctly and function over a wide range of conditions. The advantage of incorporating dedicated signaling models into cell-based models of angiogenesis is that these models provide direct links between experimental perturbations and model dynamics.

Data Availability

Source Code for the simulator and analysis can be found at https://gitlab.com/ModularABM/AngioABM/tree/VEGFR1 and https://gitlab.com/ModularABM/ABMTools/tree/VEGFR1. The model of intracellular dynamics has been submitted in SBML format to the biomodels repository (Chelliah et al., 2015) and assigned the identifier MODEL1804030001. All other data for this study is included in the manuscript and Supplementary Material.

Author Contributions

CK and SC designed the study, wrote and revised the manuscript, read the manuscript, and approved its content. CK implemented simulator and analysis software, constructed and parameterized models, and analyzed simulation results.

Funding

This work was supported by the German Research Foundation (Deutsche Forschungsgemeinschaft): CH 1123/4-1 and CH 1123/4-2 to SC. This work was further supported by the Einstein Stiftung Berlin (project How do cellular dynamics shape vascular network structure? of the Einstein Center for Regenerative Therapies to CK).

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank Edda Klipp for providing computers at the Theoretical Biophysics Department of the Humboldt-Universität zu Berlin, and Ivo Maintz for helping with them.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphys.2019.00288/full#supplementary-material

References

Andersson, E. R., Sandberg, R., and Lendahl, U. (2011). Notch signaling: simplicity in design, versatility in function. Development 138, 3593–3612. doi: 10.1242/dev.063610

PubMed Abstract | CrossRef Full Text | Google Scholar

Arganda-Carreras, I., Fernández-González, R., Muñoz-Barrutia, A., and Ortiz-De-Solorzano, C. (2010). 3D reconstruction of histological sections: application to mammary gland tissue. Microsc. Res. Tech. 73, 1019–1029. doi: 10.1002/jemt.20829

PubMed Abstract | CrossRef Full Text | Google Scholar

Barros de Andrade e Sousa, L. C. L., Kühn, C., Tyc, K. K. M., and Klipp, E. (2016). Dosage and dose schedule screening of drug combinations in agent-based models reveals hidden synergies. Front. Physiol. 6:398. doi: 10.3389/fphys.2015.00398

PubMed Abstract | CrossRef Full Text | Google Scholar

Battich, N., Stoeger, T., and Pelkmans, L. (2015). Control of transcript variability in single mammalian cells. Cell 163, 1596–1610. doi: 10.1016/j.cell.2015.11.018

PubMed Abstract | CrossRef Full Text | Google Scholar

Benn, A., Hiepen, C., Osterland, M., Schütte, C., Zwijsen, A., and Knaus, P. (2017). Role of bone morphogenetic proteins in sprouting angiogenesis: differential BMP receptor-dependent signaling pathways balance stalk vs. tip cell competence. FASEB J. 31, 4720–4733. doi: 10.1096/fj.201700193RR

PubMed Abstract | CrossRef Full Text | Google Scholar

Bentley, K., and Chakravartula, S. (2017). The temporal basis of angiogenesis. Philos. Trans. R. Soc. B Biol. Sci. 372:20150522. doi: 10.1098/rstb.2015.0522

PubMed Abstract | CrossRef Full Text | Google Scholar

Bentley, K., Gerhardt, H., and Bates, P. A. (2008). Agent-based simulation of notch-mediated tip cell selection in angiogenic sprout initialisation. J. Theor. Biol. 250, 25–36. doi: 10.1016/j.jtbi.2007.09.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Bentley, K., Mariggi, G., Gerhardt, H., and Bates, P. A. (2009). Tipping the balance: robustness of tip cell selection, migration and fusion in angiogenesis. PLoS Comput. Biol. 5:e1000549. doi: 10.1371/journal.pcbi.1000549

PubMed Abstract | CrossRef Full Text | Google Scholar

Blanco, R., Gerhardt, H., Whitehead, K. J., Smith, M. C. P., Li, Y., Lampugnani, M. G., et al. (2013). VEGF and Notch in tip and stalk cell selection. Cold Spring Harb. Perspect. Med. 3, 1–20. doi: 10.1101/cshperspect.a006569

PubMed Abstract | CrossRef Full Text | Google Scholar

Boareto, M., Jolly, M. K., Ben-Jacob, E., and Onuchic, J. N. (2015). Jagged mediates differences in normal and tumor angiogenesis by affecting tip-stalk fate decision. Proc. Natl. Acad. Sci. U.S.A. 112, E3836–E3844. doi: 10.1073/pnas.1511814112

PubMed Abstract | CrossRef Full Text | Google Scholar

Boas, S. E. M., and Merks, R. M. H. (2015). Tip cell overtaking occurs as a side effect of sprouting in computational models of angiogenesis. BMC Syst. Biol. 9:86. doi: 10.1186/s12918-015-0230-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Bothma, J. P., Levine, M., and Boettiger, A. (2010). Morphogen gradients: limits to signaling or limits to measurement? Curr. Biol. 20, 232–234. doi: 10.1016/j.cub.2010.01.040

CrossRef Full Text | Google Scholar

Bray, S. J. (2016). Notch signalling in context. Nat. Rev. Mol. Cell Biol. 9, 722–735. doi: 10.1038/nrm.2016.94

CrossRef Full Text | Google Scholar

Breier, G. (2000). Angiogenesis in embryonic development - A review. Placenta 21(Suppl.1), 11–15. doi: 10.1053/plac.1999.0525

PubMed Abstract | CrossRef Full Text | Google Scholar

Cao, Y. (2009). Positive and negative modulation of angiogenesis by VEGFR1 ligands. Sci. Signal. 2:re1. doi: 10.1126/scisignal.259re1

PubMed Abstract | CrossRef Full Text | Google Scholar

Carlier, A., Geris, L., Bentley, K., Carmeliet, G., Carmeliet, P., and van Oosterwyck, H. (2012). MOSAIC: A multiscale model of osteogenesis and sprouting angiogenesis with lateral inhibition of endothelial cells. PLoS Comput. Biol. 8:e1002724. doi: 10.1371/journal.pcbi.1002724

PubMed Abstract | CrossRef Full Text | Google Scholar

Carlier, A., Geris, L., van Gastel, N., Carmeliet, G., and Van Oosterwyck, H. (2015a). Oxygen as a critical determinant of bone fracture healing-a multiscale model. J. Theor. Biol. 365, 247–264. doi: 10.1016/j.jtbi.2014.10.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Carlier, A., Lammens, J., Van Oosterwyck, H., and Geris, L. (2015b). Computational modeling of bone fracture non-unions: four clinically relevant case studies. In Silico Cell Tissue Sci. 2:1. doi: 10.1186/s40482-015-0004-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Carlier, A., van Gastel, N., Geris, L., Carmeliet, G., and Van Oosterwyck, H. (2014). Size does matter: an integrative in vivo-in silico approach for the treatment of critical size bone defects. PLoS Comput. Biol. 10:e1003888. doi: 10.1371/journal.pcbi.1003888

PubMed Abstract | CrossRef Full Text | Google Scholar

Carmeliet, P. (2003). Angiogenesis in health and disease. Nat. Med. 9, 653–660. doi: 10.1038/nm0603-653

PubMed Abstract | CrossRef Full Text | Google Scholar

Carmeliet, P., De Smet, F., Loges, S., and Mazzone, M. (2009). Branching morphogenesis and antiangiogenesis candidates: tip cells lead the way. Nat. Rev. Clin. Oncol. 6, 315–326. doi: 10.1038/nrclinonc.2009.64

PubMed Abstract | CrossRef Full Text | Google Scholar

Carmeliet, P., and Jain, R. K. (2011). Molecular mechanisms and clinical applications of angiogenesis. Nature 473, 298–307. doi: 10.1038/nature10144

PubMed Abstract | CrossRef Full Text | Google Scholar

Chappell, J. C., Cluceru, J. G., Nesmith, J. E., Mouillesseaux, K. P., Bradley, V., Hartland, C., et al. (2016). Flt-1 (VEGFR-1) coordinates discrete stages of blood vessel formation. Cardiovasc. Res. 1, 84–93. doi: 10.1093/cvr/cvw091

CrossRef Full Text | Google Scholar

Chappell, J. C., Taylor, S. M., Ferrara, N., and Bautch, V. L. (2009). Local guidance of emerging vessel sprouts requires soluble Flt-1. Dev. Cell 17, 377–386. doi: 10.1016/j.devcel.2009.07.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Chappell, J. C., Wiley, D. M., and Bautch, V. L. (2011). Regulation of blood vessel sprouting. Semin. Cell Dev. Biol. 22, 1005–1011. doi: 10.1016/j.semcdb.2011.10.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Checa, S., and Prendergast, P. J. (2009). A mechanobiological model for tissue differentiation that includes angiogenesis: a lattice-based modeling approach. Ann. Biomed. Eng. 37, 129–145. doi: 10.1007/s10439-008-9594-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Checa, S., and Prendergast, P. J. (2010). Effect of cell seeding and mechanical loading on vascularization and tissue formation inside a scaffold: a mechano-biological model using a lattice approach to simulate cell activity. J. Biomech. 43, 961–968. doi: 10.1016/j.jbiomech.2009.10.044

PubMed Abstract | CrossRef Full Text | Google Scholar

Checa, S., Rausch, M. K., Petersen, A., Kuhl, E., and Duda, G. N. (2015). The emergence of extracellular matrix mechanics and cell traction forces as important regulators of cellular self-organization. Biomech. Model. Mechanobiol. 14, 1–13. doi: 10.1007/s10237-014-0581-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Chelliah, V., Juty, N., Ajmera, I., Ali, R., Dumousseau, M., Glont, M., et al. (2015). BioModels: ten-year anniversary. Nucleic Acids Res. 43, D542–D548. doi: 10.1093/nar/gku1181

PubMed Abstract | CrossRef Full Text | Google Scholar

Clegg, L. E., Ganta, V. C., Annex, B. H., and MacGabhann, F. (2017). Systems pharmacology of VEGF165b in peripheral artery disease. CPT: Pharmacometr. Syst. Pharmacol. 6, 833–844. doi: 10.1002/psp4.12261

PubMed Abstract | CrossRef Full Text | Google Scholar

Costa, G., Harrington, K. I., Lovegrove, H. E., Page, D. J., Chakravartula, S., Bentley, K., et al. (2016). Asymmetric division coordinates collective cell migration in angiogenesis. Nat. Cell Biol. 18, 1292–1301. doi: 10.1038/ncb3443

PubMed Abstract | CrossRef Full Text | Google Scholar

Ferrara, N., Gerber, H. P., and LeCouter, J. (2003). The biology of VEGF and its receptors. Nat. Med. 9, 669–676. doi: 10.1038/nm0603-669

PubMed Abstract | CrossRef Full Text | Google Scholar

Flegg, J. A., Menon, S. N., Maini, P. K., and McElwain, D. L. S. (2015). On the mathematical modeling of wound healing angiogenesis in skin as a reaction-transport process. Front. Physiol. 6:262. doi: 10.3389/fphys.2015.00262

PubMed Abstract | CrossRef Full Text | Google Scholar

Fong, G.-H., Rossant, J., Gertsenstein, M., and Breitman, M. L. (1995). Role of the Flt-1 receptor tyrosine kinase in regulating the assembly of vascular endothelium. Nature 376, 66–70. doi: 10.1038/376066a0

PubMed Abstract | CrossRef Full Text | Google Scholar

Funahashi, Y., Shawber, C. J., Vorontchikhina, M., Sharma, A., Outtz, H. H., and Kitajewski, J. (2010). Notch regulates the angiogenic response via induction of VEGFR-1. J. Angiogen. Res. 2:3. doi: 10.1186/2040-2384-2-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Gerhardt, H., Golding, M., Fruttiger, M., Ruhrberg, C., Lundkvist, A., Abramsson, A., et al. (2003). VEGF guides angiogenic sprouting utilizing endothelial tip cell filopodia. J. Cell Biol. 161, 1163–1177. doi: 10.1083/jcb.200302047

PubMed Abstract | CrossRef Full Text | Google Scholar

Geudens, I., and Gerhardt, H. (2011). Coordinating cell behaviour during blood vessel formation. Development 138, 4569–4583. doi: 10.1242/dev.062323

PubMed Abstract | CrossRef Full Text | Google Scholar

Harrington, L. S., Sainson, R. C. A., Williams, C. K., Taylor, J. M., Shi, W., Li, J. L., et al. (2008). Regulation of multiple angiogenic pathways by Dll4 and Notch in human umbilical vein endothelial cells. Microvasc. Res. 75, 144–154. doi: 10.1016/j.mvr.2007.06.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Hashambhoy, Y. L., Chappell, J. C., Peirce, S. M., Bautch, V. L., and Mac Gabhann, F. (2011). Computational modeling of interacting VEGF and soluble VEGF receptor concentration gradients. Front. Physiol. 2:62. doi: 10.3389/fphys.2011.00062

PubMed Abstract | CrossRef Full Text | Google Scholar

Heck, T., Vaeyens, M., and Van Oosterwyck, H. (2015). Computational models of sprouting angiogenesis and cell migration: towards multiscale mechanochemical models of angiogenesis. Math. Model. Nat. Phenomena 10, 108–141. doi: 10.1051/mmnp/201510106

CrossRef Full Text | Google Scholar

Hellström, M., Phng, L.-K., Hofmann, J. J., Wallgard, E., Coultas, L., Lindblom, P., et al. (2007). Dll4 signalling through Notch1 regulates formation of tip cells during angiogenesis. Nature 445, 776–80. doi: 10.1038/nature05571

PubMed Abstract | CrossRef Full Text | Google Scholar

Herbert, S. P., and Stainier, D. Y. R. (2011). Molecular control of endothelial cell behaviour during blood vessel morphogenesis. Nat. Rev. Mol. Cell Biol. 12, 551–64. doi: 10.1038/nrm3176

PubMed Abstract | CrossRef Full Text | Google Scholar

Hiner, M. C., Rueden, C. T., and Eliceiri, K. W. (2017). ImageJ-MATLAB: a bidirectional framework for scientific image analysis interoperability. Bioinformatics 33, 629–630. doi: 10.1093/bioinformatics/btw681

PubMed Abstract | CrossRef Full Text | Google Scholar

Hiratsuka, S., Minowa, O., Kuno, J., Noda, T., and Shibuya, M. (1998). Flt-1 lacking the tyrosine kinase domain is sufficient for normal development and angiogenesis in mice. Proc. Natl. Acad. Sci. U.S.A. 95, 9349–9354. doi: 10.1073/pnas.95.16.9349

PubMed Abstract | CrossRef Full Text | Google Scholar

Hoops, S., Gauges, R., Lee, C., Pahle, J., Simus, N., Singhal, M., et al. (2006). COPASI - A COmplex PAthway SImulator. Bioinformatics 22, 3067–3074. doi: 10.1093/bioinformatics/btl485

PubMed Abstract | CrossRef Full Text | Google Scholar

Hucka, M., Finney, A., Sauro, H. M., Bolouri, H., Doyle, J. C., Kitano, H. J., et al. (2003). The systems biology markup language (SBML): a medium for representation and exchange of biochemical network models. Bioinformatics 19, 524–531. doi: 10.1093/bioinformatics/btg015

PubMed Abstract | CrossRef Full Text | Google Scholar

Imoukhuede, P. I., Dokun, A. O., Annex, B. H., and Popel, A. S. (2013). Endothelial cell-by-cell profiling reveals the temporal dynamics of VEGFR1 and VEGFR2 membrane localization after murine hindlimb ischemia. Am. J. Physiol. Heart Circul. Physiol. 304, H1085–H1093. doi: 10.1152/ajpheart.00514.2012

PubMed Abstract | CrossRef Full Text | Google Scholar

Incio, J., Ligibel, J. A., McManus, D. T., Suboj, P., Jung, K., Kawaguchi, K., et al. (2018). Obesity promotes resistance to anti-VEGF therapy in breast cancer by up-regulating IL-6 and potentially FGF-2. Sci. Transl. Med. 10:eaag0945. doi: 10.1126/scitranslmed.aag094

PubMed Abstract | CrossRef Full Text | Google Scholar

Jain, R. K., Duda, D. G., Clark, J. W., and Loeffler, J. S. (2006). Lessons from phase III clinical trials on anti-VEGF therapy for cancer. Nat. Clin. Pract. Oncol. 3, 24–40. doi: 10.1038/ncponc0403

PubMed Abstract | CrossRef Full Text | Google Scholar

Jakobsson, L., Franco, C. A., Bentley, K., Collins, R. T., Ponsioen, B., Aspalter, I. M., et al. (2010). Endothelial cells dynamically compete for the tip cell position during angiogenic sprouting. Nat. Cell Biol. 12, 943–953. doi: 10.1038/ncb2103

PubMed Abstract | CrossRef Full Text | Google Scholar

Kappas, N. C., Zeng, G., Chappell, J. C., Kearney, J. B., Hazarika, S., Kallianos, K. G., et al. (2008). The VEGF receptor Flt-1 spatially modulates Flk-1 signaling and blood vessel branching. J. Cell Biol. 181, 847–858. doi: 10.1083/jcb.200709114

PubMed Abstract | CrossRef Full Text | Google Scholar

Kearney, J. B., Ambler, C. A., Monaco, K.-A., Johnson, N., Rapoport, R. G., and Bautch, V. L. (2002). Vascular endothelial growth factor receptor Flt-1 negatively regulates developmental blood vessel formation by modulating endothelial cell division. Blood 99, 2397–2407. doi: 10.1182/blood.V99.7.2397

PubMed Abstract | CrossRef Full Text | Google Scholar

Kearney, J. B., and Bautch, V. L. (2003). In vitro differentiation of mouse ES cells: hematopoietic and vascular development. Methods Enzymol. 365, 83–98. doi: 10.1016/S0076-6879(03)65006-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Kendall, R. L., and Thomas, K. A. (1993). Inhibition of vascular endothelial cell growth factor activity by an endogenously encoded soluble receptor. Proc. Natl. Acad. Sci. U.S.A. 90, 10705–10709. doi: 10.1073/pnas.90.22.10705

PubMed Abstract | CrossRef Full Text | Google Scholar

Keyt, B. A., and Berleau, L. T. (1996). The carboxyl-terminal domain(111-165) of vascular endothelial growth factor is critical for its mitogenic potency. J. Biol. Chem. 271, 7788–7795. doi: 10.1074/jbc.271.13.7788

PubMed Abstract | CrossRef Full Text | Google Scholar

Köhn-Luque, A., de Back, W., Starruß, J., Mattiotti, A., Deutsch, A., Pérez-Pomares, J. M., et al. (2011). Early embryonic vascular patterning by matrix-mediated paracrine signalling: a mathematical model study. PLoS ONE 6:e0024175. doi: 10.1371/journal.pone.0024175

PubMed Abstract | CrossRef Full Text | Google Scholar

Köhn-Luque, A., de Back, W., Yamaguchi, Y., Yoshimura, K., Herrero, M. A., and Miura, T. (2013). Dynamics of VEGF matrix-retention in vascular network patterning. Phys. Biol. 10:066007. doi: 10.1088/1478-3975/10/6/066007

PubMed Abstract | CrossRef Full Text | Google Scholar

Kopan, R., and Ilagan, M. X. G. (2009). The canonical notch signaling pathway: unfolding the activation mechanism. Cell 137, 216–233. doi: 10.1016/j.cell.2009.03.045

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, Z.-J., Shirakawa, T., Li, Y., Soma, A., Oka, M., Dotto, G. P., et al. (2003). Regulation of Notch1 and Dll4 by vascular endothelial growth factor in arterial endothelial cells: implications for modulating arteriogenesis and angiogenesis. Mol. Cell. Biol. 23, 14–25. doi: 10.1128/MCB.23.1.14-25.2003

PubMed Abstract | CrossRef Full Text | Google Scholar

Lobov, I. B., Renard, R. A., Papadopoulos, N., Gale, N. W., Thurston, G., Yancopoulos, G. D., et al. (2007). Delta-like ligand 4 (Dll4) is induced by VEGF as a negative regulator of angiogenic sprouting. Proc. Natl. Acad. Sci. U.S.A. 104, 3219–3224. doi: 10.1073/pnas.0611206104

PubMed Abstract | CrossRef Full Text | Google Scholar

Logsdon, E. A., Finley, S. D., Popel, A. S., and MacGabhann, F. (2014). A systems biology view of blood vessel growth and remodelling. J. Cell. Mol. Med. 18, 1491–1508. doi: 10.1111/jcmm.12164

PubMed Abstract | CrossRef Full Text | Google Scholar

MacGabhann, F., Ji, J. W., and Popel, A. S. (2006). VEGF gradients, receptor activation, and sprout guidance in resting and exercising skeletal muscle. J. Appl. Physiol. 102, 722–734. doi: 10.1152/japplphysiol.00800.2006

CrossRef Full Text | Google Scholar

MacGabhann, F., and Popel, A. S. (2004). Model of competitive binding of vascular endothelial growth factor and placental growth factor to VEGF receptors on endothelial cells. Am. J. Physiol. Heart Circul. Physiol. 286, H153–H164. doi: 10.1152/ajpheart.00254.2003

CrossRef Full Text | Google Scholar

MacGabhann, F., and Popel, A. S. (2006). Targeting neuropilin-1 to inhibit VEGF signaling in cancer: comparison of therapeutic approaches. PLoS Comput. Biol. 2:e180. doi: 10.1371/journal.pcbi.0020180

CrossRef Full Text | Google Scholar

Merks, R. M. H., Newman, S. A., and Glazier, J. A. (2004). “Cell-oriented modeling of in vitro capillary development,” in Cellular Automata, eds P. M. A. Sloot, B. Chopard, and A. G. Hoekstra (Berlin; Heidelberg: Springer Berlin Heidelberg), 425–434.

Google Scholar

Novère, N. L., Hucka, M., Mi, H., Moodie, S., Schreiber, F., Sorokin, A., et al. (2009). The systems biology graphical notation : abstract : nature biotechnology. Nat. Biotechnol. 27, 735–742. doi: 10.1038/nbt.1558

CrossRef Full Text | Google Scholar

Qazi, T. H., Mooney, D. J., Duda, G. N., and Geissler, S. (2017). Biomaterials that promote cell-cell interactions enhance the paracrine function of MSCs. Biomaterials 140, 103–114. doi: 10.1016/j.biomaterials.2017.06.019

PubMed Abstract | CrossRef Full Text | Google Scholar

Roberts, D. M., Kearney, J. B., Johnson, J. H., Rosenberg, M. P., Kumar, R., and Bautch, V. L. (2004). The vascular endothelial growth factor (VEGF) receptor Flt-1 (VEGFR-1) modulates Flk-1 (VEGFR-2) signaling during blood vessel formation. Am. J. Pathol. 164, 1531–1535. doi: 10.1016/S0002-9440(10)63711-X

PubMed Abstract | CrossRef Full Text | Google Scholar

Ruhrberg, C., Gerhardt, H., Golding, M., Watson, R., Ioannidou, S., Fujisawa, H., et al. (2002). Spatially restricted patterning cues provided by heparin-binding VEGF-A control blood vessel branching morphogenesis. Genes Dev. 16, 2684–2698. doi: 10.1101/gad.242002

PubMed Abstract | CrossRef Full Text | Google Scholar

Schwanhäusser, B., Busse, D., Li, N., Dittmar, G., Schuchhardt, J., Wolf, J., et al. (2011). Global quantification of mammalian gene expression control. Nature 473, 337–342. doi: 10.1038/nature10098

PubMed Abstract | CrossRef Full Text | Google Scholar

Shaya, O., and Sprinzak, D. (2011). From Notch signaling to fine-grained patterning: Modeling meets experiments. Curr. Opin. Genet. Dev. 21, 732–739. doi: 10.1016/j.gde.2011.07.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Solomon, S. D., Lindsley, K., Vedula, S. S., Krzystolik, M. G., and Hawkins, B. S. (2014). Anti-vascular endothelial growth factor for neovascular age-related macular degeneration. Cochr. Database Syst. Rev. 29:CD005139. doi: 10.1002/14651858.CD005139.pub3

CrossRef Full Text | Google Scholar

Spiesser, T. W., Kühn, C., Krantz, M., and Klipp, E. (2015). Bud-localization of CLB2 mRNA can constitute a growth rate dependent daughter sizer. PLOS Comput. Biol. 11:e1004223. doi: 10.1371/journal.pcbi.1004223

PubMed Abstract | CrossRef Full Text | Google Scholar

Sprinzak, D., Lakhanpal, A., Lebon, L., Santat, L. A., Fontes, M. E., Anderson, G. A., et al. (2010). Cis-interactions between Notch and Delta generate mutually exclusive signalling states. Nature 465, 86–90. doi: 10.1038/nature08959

PubMed Abstract | CrossRef Full Text | Google Scholar

Suchting, S., Freitas, C., le Noble, F., Benedito, R., Bréant, C., Duarte, A., et al. (2007). The Notch ligand Delta-like 4 negatively regulates endothelial tip cell formation and vessel branching. Proc. Natl. Acad. Sci. U.S.A. 104, 3225–3230. doi: 10.1073/pnas.0611177104

PubMed Abstract | CrossRef Full Text | Google Scholar

The Boost Community (2017). Boost C++ Libraries.

Tummler, K., Lubitz, T., Schelker, M., and Klipp, E. (2014). New types of experimental data shape the use of enzyme kinetics for dynamic network modeling. FEBS J. 281, 549–571. doi: 10.1111/febs.12525

PubMed Abstract | CrossRef Full Text | Google Scholar

Tung, J. J., Tattersall, I. W., and Kitajewski, J. (2012). Tips, stalks, tubes: notch-mediated cell fate determination and mechanisms of tubulogenesis during angiogenesis. Cold Spring Harb. Perspect. Med. 2:a006601. doi: 10.1101/cshperspect.a006601

PubMed Abstract | CrossRef Full Text | Google Scholar

Tyson, J. J., Chen, K. C., and Novak, B. (2003). Sniffers, buzzers, toggles and blinkers: dynamics of regulatory and signaling pathways in the cell. Curr. Opin. Cell Biol. 15, 221–231. doi: 10.1016/S0955-0674(03)00017-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Ubezio, B., Blanco, R. A., Geudens, I., Stanchi, F., Mathivet, T., Jones, M. L., et al. (2016). Synchronization of endothelial Dll4-Notch dynamics switch blood vessels from branching to expansion. eLife 5:e12167. doi: 10.7554/eLife.12167

PubMed Abstract | CrossRef Full Text | Google Scholar

van Oers, R. F. M., Rens, E. G., LaValley, D. J., Reinhart-King, C. A., and Merks, R. M. H. (2014). Mechanical cell-matrix feedback explains pairwise and collective endothelial cell behavior in vitro. PLoS Comput. Biol. 10:e1003774. doi: 10.1371/journal.pcbi.1003774

PubMed Abstract | CrossRef Full Text | Google Scholar

Venkatraman, L., Regan, E. R., and Bentley, K. (2016). Time to decide? Dynamical analysis predicts partial tip/stalk patterning states arise during angiogenesis. PLoS ONE 11:e166489. doi: 10.1371/journal.pone.0166489

PubMed Abstract | CrossRef Full Text | Google Scholar

Virgili, G., Parravano, M., Menchini, F., and Brunetti, M. (2012). Antiangiogenic therapy with anti-vascular endothelial growth factor modalities for diabetic macular oedema. Cochr. Database Syst. Rev. 12:CD007419. doi: 10.1002/14651858.CD007419.pub3

PubMed Abstract | CrossRef Full Text | Google Scholar

Wagner, J., Kline, C. L., Zhou, L., Khazak, V., and El-Deiry, W. S. (2018). Anti-tumor effects of ONC201 in combination with VEGF-inhibitors significantly impacts colorectal cancer growth and survival in vivo through complementary non-overlapping mechanisms. J. Exp. Clin. Cancer Res. 37, 1–12. doi: 10.1186/s13046-018-0671-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Walpole, J., Chappell, J. C., Cluceru, J. G., MacGabhann, F., Bautch, V. L., and Peirce, S. M. (2015). Agent-based model of angiogenesis simulates capillary sprout initiation in multicellular networks. Integr. Biol 7, 987–97. doi: 10.1039/C5IB00024F

PubMed Abstract | CrossRef Full Text | Google Scholar

Waltenberger, J., Claesson-Welsh, L., Siegbahn, A., Shibuya, M., and Heldin, C. H. (1994). Different signal transduction properties of KDR and Flt1, two receptors for vascular endothelial growth factor. J. Biol. Chem. 269, 26988–26995.

PubMed Abstract | Google Scholar

Wu, F. T. H., Stefanini, M. O., MacGabhann, F., Kontos, C. D., Annex, B. H., and Popel, A. S. (2010). VEGF and soluble VEGF receptor-1 (sFlt-1) distributions in peripheral arterial disease: an in silico model. Am. J. Physiol. Heart Circul. Physiol. 298, H2174–H2191. doi: 10.1152/ajpheart.00365.2009

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Z., Neiva, K. G., Lingen, M. W., Ellis, L. M., and Nör, J. E. (2010). VEGF-dependent tumor angiogenesis requires inverse and reciprocal regulation of VEGFR1 and VEGFR2. Cell Death Differ. 17, 499–512. doi: 10.1038/cdd.2009.152

PubMed Abstract | CrossRef Full Text | Google Scholar