The formation of new blood vessels from the existing vasculature, called angiogenesis, is regulated by different cells and processes (Carmeliet, 2005; Eelen et al., 2020; Zhang et al., 2022). It is generally classified as sprouting angiogenesis (sprouting of endothelial cells from a mother vessel based on environmental cues) or intussusceptive angiogenesis (vessel splitting internally into two daughter vessels), with the former being better characterized and studied than the latter, and both essential for regeneration of microvasculature (Arpino et al., 2021). Dysregulation in angiogenesis has been associated with the development and complication of different diseases, such as cancer, diabetic retinopathy, age-related macular degeneration, and cardiovascular diseases, including peripheral arterial disease (PAD) (Zhang et al., 2022). Many studies focus on modulating and controlling angiogenesis as a therapeutic approach for treating ischemic diseases and cancer (Teleanu et al., 2019; Annex and Cooke, 2021; Zhang et al., 2023). PAD, for instance, is an atherosclerotic disease characterized by lower limb ischemia; the greater the extent of the endogenous angiogenic response, the lesser the patients’ symptoms. Often accompanied by severe microvascular disease, PAD affects more than 200 million people worldwide and is associated with 53%–90% of all major amputations of the lower limb (Londero et al., 2019; Aday and Matsushita, 2021). Currently, there is no cure for this disease, and medical therapies have limited efficacy in treating it. Therapeutic angiogenesis remains an approach that aims to alter the course of the disease (Annex and Cooke, 2021; Han et al., 2022). Among the main signaling pathways that regulate angiogenesis are the Notch signaling and the oxygen-sensing hypoxia-inducible factor (HIF) pathways. Under hypoxia, VEGFA (and possibly other angiogenic factors) is released in the hypoxic tissue, guiding endothelial cells (ECs) to the ischemic site to restore blood flow and oxygen availability. Once attracted, ECs compete for the tip position, a dynamic process guided by the Notch signaling pathway and affected by hypoxia (Figure 1) (Naito et al., 2020; Naiche et al., 2022). Additionally, EC survival under hypoxia is stimulated through autocrine VEGFA signaling, although it does not significantly affect angiogenesis (Lee et al., 2007).

The Notch signaling pathway is a major regulator of cell fate, cell differentiation, vascular patterning, and intercellular communication in processes such as angiogenesis and tumor growth. The Notch pathway is highly conserved across vertebrates, and it affects different cells, such as macrophages (Lin et al., 2018), endothelial cells (ECs), smooth muscle cells (SMCs) (Fouillade et al., 2012), and pericytes (Tefft et al., 2022). In EC differentiation, the Notch ligands Delta-like ligand 4 (Dll4) and Jagged 1 render different effects when interacting with Notch receptors. While Dll4 is a negative regulator of tip cell formation, Jagged 1 positively regulates tip cell formation and sprouting (Benedito et al., 2009; Pedrosa et al., 2015). Dll4 and Jagged also have a role in the spatial control of sprouting (Tiemeijer et al., 2022). This competition between Notch ligands to bind to the Notch receptor is required to promote functional revascularization and tissue reperfusion.

Different cells in ischemic tissues release pro-angiogenic factors such as vascular endothelial growth factor (VEGFA) (Couffinhal et al., 1998; Rissanen et al., 2002). VEGF induces endothelial cell migration and sprouting in a chemotaxis-dependent fashion (Lee et al., 2022). This pro-angiogenesis mechanism depends on environmental cues such as low oxygen levels, driven by the HIF signaling pathway. Once low oxygen levels are sensed in the tissue, three types of hypoxia-inducible factors (HIF-1, HIF-2, and HIF-3), are activated, each composed of two subunits (α and β). The α-subunits are regulated by changes in oxygen concentration based on proteolytic degradation and transcriptional regulation. Under ischemia, HIF1-α upregulation relates to tissue inflammation, and it has been an important therapeutic target and molecule of interest in angiogenesis models (Qutub and Popel, 2006; Qutub and Popel, 2007; Qutub and Popel, 2008; Cavadas et al., 2013; Nguyen et al., 2013; Fábián et al., 2016).

The HIF and the Notch pathways are essential in the response to ischemia, and understanding the mechanistic interactions between them and other pathways that drive EC differentiation, proliferation, and stability can help develop new therapeutic strategies for diseases such as PAD. Even for converging pathways, differences in temporal regulation affect their signaling outcomes (Jin et al., 2014). Computational models have been developed and provided important insights into understanding signaling pathways and cell-cell interactions (Subramanian et al., 2022). For instance, Zhao et al. designed an ordinary differential equation (ODE)-based computational model to understand hypoxia-responsive miRNAs control of the HIF-VEGF pathway in EC (Zhao and Popel, 2015). Also, Venkatraman et al. formulated an ODE-based mathematical model to better understand how local and intracellular conditions influence EC patterning (tip or stalk EC) (Venkatraman et al., 2016). Their model showed a partial and stable tip/stalk stage with duration determined by intracellular signaling. Recently, Kuhn and Checa proposed an elegant model of the contributions of VEGF receptors to lateral inhibition during sprouting angiogenesis, integrating agent-based and ODE-based modeling (Kühn and Checa, 2019). Despite advances in experimental biology methods, observing signaling interactions between two microvascular endothelial cells under specific conditions can be challenging. Thus, in silico models can provide insights to help understand the mechanisms of EC patterning, signaling, and behavior under pathological conditions.

To date, most computational models of endothelial cell patterning during angiogenesis are based on Cellular Pots or agent-based modeling strategies (Kühn and Checa, 2019). Previous agent-based models on Notch signaling report how VEGF concentration relates to the elongation, migration, and proliferation of ECs (Qutub and Popel, 2009). Recently, Koon and others presented a computational model of EC patterning focusing on Notch signaling heterogeneity to explain a greater variety of cell patterning; they report hybrid endothelial cells with tip and stalk characteristics (Koon et al., 2018). A more comprehensive Boolean model, including different signaling pathways in EC during angiogenesis, was presented by Weinstein et al., incorporating the molecular regulatory network, extracellular microenvironments, and loss- and gain-of-function mutations (Weinstein et al., 2017). Mechanistic models of signaling networks in EC have also been developed. Bazzazi et al. designed a rule-based computational model implemented in BioNetGen to assess the effects of thrombospondin-1 (TSP1)-CD47 signaling through VEGF on ERK1/2 and calcium (Bazzazi et al., 2017). Using rule-based modeling, the group also investigated TSP1 inhibition of VEGF signaling to Akt-endothelial nitric oxide synthase (eNOS) (Bazzazi et al., 2018). A hybrid multi-scale model of endothelial cells during angiogenesis has also been presented by Stepanova et al., including features such as branching, cell mixing, and the brush border effect; their work indicated a dependence on the time evolution of cell mixing and the newly formed branching structure (Stepanova et al., 2021). Mechanistic computational models of intracellular and intercellular events should include the major and most relevant pathways known to the field studied (Zhang et al., 2022). Additionally, the methodology used for model building should be consistent and follow guidelines to achieve good modeling practices (Mitra and Hlavacek, 2019).

The complex molecular signaling networks that drive angiogenesis and determine EC fate are largely affected by hypoxia; note that the effect may be context-dependent and differ in different tissues such as ischemic skeletal muscle, tumor, or ischemic retina. Exogenous stimulation through VEGF and signaling between interacting ECs determine cell fate and guide downstream signaling, contributing to functional angiogenesis. Although the basic mechanisms that regulate EC patterning during angiogenesis are known, mechanistic details of the interaction between pathways and their effects and how abnormal conditions affect them need further studies (Zhou et al., 2022). In this work, we assess changes in the dynamics of such networks under hypoxia conditions to understand how VEGF, Notch, and HIF signaling pathways interact with each other and regulate downstream pathways, as well as how they affect EC patterns through time. Here we also present a comprehensive methodological approach for model development, based on recent reviews on good modeling practices, and we apply this methodology to the system of interest. Using a mechanistic model of two interacting ECs, we evaluate their signal exchange under varying oxygen conditions, investigate the effects of different interventions on cell pattern determination, and discuss their potential application in the context of angiogenesis modulation, considering the numerous effects of hypoxia on angiogenesis-involved signaling pathways (Rodriguez et al., 2021).

We formulate the model network based on biological knowledge and experimental evidence and obtain the corresponding ordinary differential equations representing each biochemical reaction. Our integrative model is based on a modular structure, each module representing a biological process, such as oxygen sensing, VEGF signaling, calcium, and NO cycling, and the Notch pathway as illustrated in Figure 2. The reactions are based on mass action, Hill, and Michaelis-Menten kinetics. The model is implemented in the Matlab SimBiology software (MathWorks Inc, 2023b). The sbml and Matlab files are provided as Supplementary Material, with the initial conditions and parameter values used (Supplementary Tables ST1–ST7). A description of model components is also provided. The model was simulated using MATLAB R2023a and includes the pathways considered more relevant to our understanding of the interaction between ECs and the effect of hypoxia. The model is composed of two cells, interacting through the Notch/Dll4 pathway (to simulate lateral inhibition) and stimulated by VEGF (similar to bolus dosages). The exogenous VEGFA bolus represents exogenous gradients of VEGFA leading to EC patterning (paracrine effect of VEGFA). Under hypoxia, we also model an increase in exogenous VEGFA and in intracellular VEGFA production by the EC. Given evidence from the literature, we assume that VEGFA produced by EC under hypoxia does not affect EC patterning and angiogenesis but is present to promote cell survival signal through Akt and Nitric Oxide pathways (Lee et al., 2007). Therefore, we do not consider the transport of EC-produced VEGFA between the two interacting cells.

In this section, we present the results of our model construction, methodological approaches, and simulations.

Sprouting angiogenesis as a response to hypoxia is a necessary process tightly regulated by an amalgamate of agents and signaling pathways (Hashimoto and Shibasaki, 2015; Naito et al., 2020; Rodriguez et al., 2021; Zhang et al., 2022). Endothelial cell patterning into the tip or stalk phenotypes guide sprout formation, promoting blood flow restoration to hypoxic tissues, such as seen in the context of diseases like PAD (Gustafsson et al., 2005; Annex and Cooke, 2021; Han et al., 2022). Efforts to characterize, understand, predict, and control cell behavior under abnormal conditions are relevant in this field, as they can help develop new therapeutic strategies and interventions in angiogenesis.

In this study, we propose a structured, state-of-the-art methodology to comply with good modeling practices and reduce the computational cost required for fitting and simulating large models. We present the step-by-step procedure to perform each of the main parts of model design, training, and validation, and we apply this methodology to a new integrative model of the main pathways in EC during angiogenesis. Through this kinetics-based model, we represent and simulate the interaction and patterning of 2 ECs stimulated with VEGFA under normoxia and hypoxia conditions. This model was developed as a first step towards investigating the differential signaling and response of each of the two cells, as they are differently stimulated by VEGFA and interact through Notch signaling under varying oxygen conditions. We calibrated the model using data from experiments performed in endothelial cells for both oxygen conditions and defined a ratio for stalk/tip pattern identification based on known ligands and receptors overexpressed in each of these phenotypes (del Toro et al., 2010; Xu and Li, 2022). Considering recent efforts in designing and characterizing “virtual cells” (Zhao et al., 2021; Lim et al., 2022; Zhang et al., 2022), this work contributes to the field by bringing insights into cell-cell interaction, patterning, and signaling with the potential translation of results to pathological conditions (specifically, hypoxia-characterized diseases such as PAD).

As new experimental results become available, revisiting previous assumptions and models is essential, making adjustments that bring us closer to in vivo conditions. Endothelial cell patterning and interaction is a significant field of study in therapeutic angiogenesis since these cells are important determinants of angiogenesis and blood flow restoration (Qiu and Hirschi, 2019; Lim et al., 2022). Considering the many aspects that influence their heterogeneity (Gifre-Renom et al., 2022), a computational approach provides the flexibility required to implement further changes and test hypotheses. Our model reproduces the characteristic EC patterning behavior under normoxia given differential VEGF stimulation (Figure 10) and predicts a distinct behavior under hypoxia (Figure 11). It also supports the requirement of a paracrine differential VEGF stimulation for cell patterning. Additionally, it indicates differential signaling and patterning under hypoxia between the two cells, which is supported by previous studies showing that hypoxia regulates and promotes tip/stalk differentiation through overexpression of Dll4, indirect regulation of VEGFR2 and NRP1 established by Notch signaling, and upregulation of Hes1 (Carmeliet et al., 2009; Rodriguez et al., 2021).

Tip and stalk endothelial cells interact and respond differently to signals and targeting them has a therapeutic potential in terms of tip-stalk cell control and selection (Chen et al., 2019). Through global sensitivity analysis (Figures 8, 9; Supplementary Table ST9), our model allowed us to investigate this distinct behavior in a system of two interacting cells, focusing on which parameters and reactions are more influential in defining each cell’s pattern. Although we expected parameters from the Notch pathway to be the most influential in determining cell patterning, our simulations indicate that, although they play a part in modulating the ECPI, parameters in reactions involved in the VEGF, ERK, and Akt pathways have a higher influence on the ECPI. As these pathways’ activation leads to Notch signaling activation, their influence on the ECPI becomes clear in an integrative model such as the one we present. ECs require differential VEGF stimulation to assume a pattern (tip or stalk), and similar stimulation turns cells toward an undifferentiated state. Previously, others have reported that EC can assume tip, stalk, and quiescent phalanx states, where the latter refers to more mature EC, formed once the vessel is perfused (Pasut et al., 2021). Partial active/inactive states during differentiation have also been reported through computational models (Venkatraman et al., 2016). These two states are seen previous to the point where the cell assumes the tip or stalk position, under stimulation by VEGF. In this work by Venkatraman et al., EC quiescence is noted for similar VEGF stimulation (VEGF at cell 1 = VEGF at cell 2 = 0) and led us to question if the parameters that define cell patterning differ and, if so, how they differ when cells are differentially or similarly stimulated or under different oxygen conditions. Through global sensitivity analysis of the ECPI of both cells subject to similar (Figure 8A) or different (Figure 8B) VEGF stimulation, we found that for several parameters, the differential stimulation that guides patterning leads to a change in parameter effect (moving from negative PRCC to positive PRCC or vice-versa), but this effect was parameter-dependent. This simulation advocates for cell patterning as being a highly regulated event that depends on the intensity of VEGF signaling being sensed by each cell. Our model is also able to reproduce the effects of DAPT inhibition of the Notch-induced cell patterning post VEGF stimulation (Supplementary Figure S14) seen experimentally (Takeshita et al., 2007).

In this work, we assume that hypoxia’s effects on driving angiogenesis are mostly due to the paracrine effect of VEGF released by other cells on endothelial cells. This is based on previous evidence, where the autocrine effect of VEGF (generated by EC under hypoxia) is shown to be mostly for cell survival, with few contributions to angiogenesis (Lee et al., 2007). However, recent work indicates otherwise (Jin et al., 2019), with VEGF production by EC under hypoxia leading to angiogenesis. We hypothesize that this difference is due to the cell environment tested in each experiment, as exogenous VEGF stimulation by other cell types might shadow the angiogenic effects of autocrine VEGF. Additionally, we assume that the effects of hypoxia on regulating VEGF receptors on cell surface depend on Notch signaling, which we pose as a question to be answered in future experiments. Future experiments should also assess the extent to which Notch signaling participates in defining vascular permeability through NO regulation and whether this interaction interferes with vascular permeability and leakiness profile in pathological conditions such as PAD. Many methods are available to estimate NO in vivo and in vitro (Goshi et al., 2019), and previous works have estimated NO release by HUVECs to be within the ranges of nM to μ M (Østergaard et al., 2007; Ugusman et al., 2014; Janaszak-Jasiecka et al., 2018). Additionally, in vivo measurements of NO report values of about 1  μ M, with the values of EC50 for soluble guanylate cyclase (sGC) in the vascular smooth muscle as low as several nM (Chen et al., 2008). Our simulations predict concentrations in the order of 10 − 1 μ M. This difference can be due to the initial concentration of NO assumed (0  μ M), and more precision of absolute values estimates might require fitting the model to absolute values of NO, instead of normalized to their maximum concentration only. This assessment should also be included in future works. Despite this difference in the order of magnitude predicted by our simulations, compared to experimental in vitro data, our model qualitatively indicates the difference in NO production by the two cells, which was among our stated goals.

The methodological approach we used in this work allowed us to develop the model and perform simulations within reasonable computation time, despite the large size of the model. The modular implementation during global optimization aided in speeding up the process of model fitting. Performing SIA and PIA as a way of achieving good modeling practices and limiting the number of parameters for fitting is an additional aspect that helped in this regard. Our method is, therefore, efficient for working with larger mechanistic computational models. Recently, a new framework and steps to build large mechanistic models was proposed, integrating annotated input text files for specific data, python-based platforms for processing the input files and generating Antimony files to be converted to SBML standards and Python for model simulations (Erdem et al., 2022). Another computational framework for parameterization of large-scale mechanistic models has been proposed, with significant advancement on computation time for very large models (>1000 parameters) (Fröhlich et al., 2018). In this framework, the authors include practical identifiability analysis not as a step prior to optimization/fitting, but as an investigation of prediction uncertainty caused by parameter uncertainties performed post-fitting. Their results pointed to most of the calibrated parameters being poorly identifiable, but still allowing them to obtain low-uncertainty predictions. Practical identifiability requires both a sensitivity analysis (to verify that the parameters are influential to the observables in question) and a collinearity analysis (to exclude collinear parameters), and performing these evaluations before the final optimization of the values helps limit the number of parameters being fitted (therefore, the computational cost and time of model optimization) and improve prediction uncertainty caused by parameter uncertainty. Although the method we proposed in this work relies on a simple collinearity analysis (removing parameters with opposing effects before fitting), this method could be automatized by using platforms such as VisID (Gábor et al., 2017). An improved approach that could be tested in future works is to perform the structural identifiability, followed by our simplified version of the practical identifiability (based on initial guesses for parameters within reasonable ranges), then the global optimization, and as an additional analysis to UQ, a new round of practical identifiability (performed with an automatic platform such as VisID). This method might reduce the number of unidentifiable parameters prior to fitting and provide a more accurate and trustable model.

Despite our efforts to calibrate the model with experimental data from the same line of endothelial cells, we were not always able to find the required data, considering the level of detail and ramifications of our model. Data from different cell lines might produce inaccurate simulations and results (Chi et al., 2003). Another limitation of our model is not including Jagged in the Notch signaling pathway, which should be accounted for in future studies, along with a more detailed model of the effects of hypoxia on the calcium signaling pathway (Berna et al., 2002) and of angiopoietin on the Notch signaling pathway (Machado et al., 2019). Our model is based on in vitro data for calibration and assumptions defined. Translating our results to in vivo conditions requires further experiments and mathematical and physiological assumptions not included in this work. Additionally, obtaining experimental data suggested by the present model would further validate the model. Our model is also limited by the pathways considered, as other signaling pathways might be affected by hypoxia and influence cell behavior, such as Ang-Tie, JAK/STAT, and TLR signaling. Those could be considered for future studies. We also note that in the model presented we consider that the effects of HIF1 and HIF2 are summative on inducing VEGFA signaling. Previous works have shown that these HIF1 and HIF2 undergo a switch regarding their time dynamics, to ensure a continuous activation of pathways response to hypoxia and prolong cell survival (Bartoszewski et al., 2019). However, others have shown that HIF1 and HIF2 differentially regulate the hypoxia response, in a context dependent fashion (Branco-Price et al., 2012). In fact, HIF2 in macrophages has been shown to induce an anti-angiogenesis response by inducing soluble VEGFR1, which sequesters free VEGFA from the environment (Eubank et al., 2011). In ECs, HIFs have also shown opposite effects on regulating interleukin 8 (IL-8), a promoter of angiogenesis (Florczyk et al., 2011). Although in out model such opposite effects are not included, the base model here provided can be further developed to represent and simulate additional differential effects of different HIFs.

Another limitation of the model is the estimation of parameters that have not been quantified in the literature or in previous models. Our simulations are constrained to the estimated values. As additional experimental data become available, this constraint could be revisited, and the estimates updated. Additionally, ECs under hypoxia act differently from ECs under hypoxia serum starvation (representing PAD conditions), which is a future aspect to be investigated with our modeling strategy. Previous works discuss the role of calcium oscillations on endothelial cell patterning (Debir et al., 2021). In this work, we consider the signaling through calcium in the endoplasmic reticulum and cytoplasm, with the extracellular calcium effect accounted for through CRAC channels calcium influx and the Ca++–Na + exchanger (under hypoxia). In future work, the focus could be more on the calcium effect on EC patterning, using a more comprehensive model of the signaling through calcium in EC under hypoxia.

Our modeling strategy in this work is non-spatial, differing from several other models in the field (Bentley et al., 2008; Reynolds et al., 2019), and focusing mostly on time-dynamics of species involved. Although observing spatial distribution of tip/stalk cells provides important insights on a tissue-level, for example, regarding cell distribution, diffusion of molecules and cell migration, obtaining a more focused view on intra and intercellular dynamics of a pathway as complex as the Notch pathway is also required to design and optimize pro-angiogenesis therapies focused on modulating Notch signaling. Non-spatial, ODE-based models are also more easily calibrated against experimental data, with clear modeling protocols. Given that non-spatial and spatial analysis are both integral parts for understanding the Notch pathway, its downstream effects and upstream regulators, a next step for this model is to integrate spatial analysis. This integration of ODE-based and spatial models has been recently shown on a viral infection and immune response model (Sego et al., 2021).

The Notch signaling has been a therapeutic focus for treatment of conditions related to angiogenesis. Notch inhibition has been applied to treating tumors (Akil et al., 2021; Jiang et al., 2022) and immune and inflammatory disorders (Rizzo and Ferrari, 2015; Allen and Maillard, 2021). Notch signaling is known to affect different cells, as previously discussed, and its therapeutic application requires a detailed understanding of its effects and its effectors. Our model provides a detailed integrative approach of evaluating Notch signaling and its effects in the context of angiogenesis, and can be extended to integrate Notch inhibitors (e.g., DAPT to simulate therapeutic applications of Notch in different pathological scenarios, including pro- and anti-angiogenic contexts (Niu et al., 2022; You et al., 2023), as represented in Supplementary Figure S14. We present this work as a first step in understanding the interconnection of pathways between endothelial cells during hypoxia-induced angiogenesis, as well as differences between tip and stalk ECs under different stimulatory conditions. We believe it expands the knowledge brought by previous studies in the field, and it is our first approach to working with larger, integrative models in ECs with the goal of building a network that can be used for therapeutic angiogenesis purposes. Additionally, it can be used to test hypotheses related to the pathways included (e.g., vascular permeability control through Notch signaling modulation). In summary, our model provides a highly integrative, data-driven platform that can be built on by simply adding pathway inhibitors or stimulators for testing angiogenesis-related therapies and to study differences in signaling between tip/stalk EC at varying oxygen conditions.