Impact Factor 3.394

The world's 3rd most-cited Physiology journal

Original Research ARTICLE

Front. Physiol., 27 November 2017 | https://doi.org/10.3389/fphys.2017.00960

A Network Model to Explore the Effect of the Micro-environment on Endothelial Cell Behavior during Angiogenesis

  • 1ABACUS-Laboratorio de Matemáticas Aplicadas y Cómputo de Alto Rendimiento, Departamento de Matemáticas, Centro de Investigación y de Estudios Avanzados CINVESTAV-IPN, Mexico City, Mexico
  • 2CompBioLab, Departamento de Biología Molecular y Biotecnología, Instituto de Investigaciones Biomédicas, Universidad Nacional Autónoma de México, Mexico City, Mexico
  • 3Departamento de Física, Instituto Nacional de Investigaciones Nucleares, Mexico City, Mexico

Angiogenesis is an important adaptation mechanism of the blood vessels to the changing requirements of the body during development, aging, and wound healing. Angiogenesis allows existing blood vessels to form new connections or to reabsorb existing ones. Blood vessels are composed of a layer of endothelial cells (ECs) covered by one or more layers of mural cells (smooth muscle cells or pericytes). We constructed a computational Boolean model of the molecular regulatory network involved in the control of angiogenesis. Our model includes the ANG/TIE, HIF, AMPK/mTOR, VEGF, IGF, FGF, PLCγ/Calcium, PI3K/AKT, NO, NOTCH, and WNT signaling pathways, as well as the mechanosensory components of the cytoskeleton. The dynamical behavior of our model recovers the patterns of molecular activation observed in Phalanx, Tip, and Stalk ECs. Furthermore, our model is able to describe the modulation of EC behavior due to extracellular micro-environments, as well as the effect due to loss- and gain-of-function mutations. These properties make our model a suitable platform for the understanding of the molecular mechanisms underlying some pathologies. For example, it is possible to follow the changes in the activation patterns caused by mutations that promote Tip EC behavior and inhibit Phalanx EC behavior, that lead to the conditions associated with retinal vascular disorders and tumor vascularization. Moreover, the model describes how mutations that promote Phalanx EC behavior are associated with the development of arteriovenous and venous malformations. These results suggest that the network model that we propose has the potential to be used in the study of how the modulation of the EC extracellular micro-environment may improve the outcome of vascular disease treatments.

1. Introduction

The circulatory system allows for the existence of large multicellular organisms, ensuring adequate oxygen and nutrient supply. Blood vessels are composed of three main layers. The outermost layer—the tunica adventitia—contains elastic fibers, collagen, and connective tissue. The middle layer—the tunica media—is comprised of smooth muscle cells, collagen, and elastin, and the innermost layer—the tunica intima—, which is exposed to the vessel lumen, is a single-cell layer of endothelium. The circulatory system is not a static structure, it adapts to the changing requirements of the body by means of vasculogenesis, arteriogenesis, and angiogenesis (Betz et al., 2016).

Vasculogenesis is a process that allows for the de novo formation of blood vessels. The formation of the first blood vessels in the embryo involves the differentiation of cells from the mesodermal blood islands into angioblasts, also called endothelial precursor cells (EPCs). During later development, angioblasts may differentiate from hematopoietic stem cells, multipotent bone marrow progenitor cells, myeloid cells (specifically monocytes and macrophages), side population cells, and pluripotent stem cells (Kässmeyer et al., 2009). After the differentiation of EPCs, the cells must migrate and aggregate to form a primitive vascular blood plexus. Then, for the vascular network assembly, three mechanisms have been proposed: (a) Extracellular matrix contact guidance, where the ECs are guided by collagen fibers present in the extracellular matrix and each cell may change the tension and orientation of the collagen fiber network to guide other cells, (b) Autocrine chemotaxis, where the ECs follow a morphogen (such as VEGFA) gradient and then secrete the morphogen altering the gradient to guide other cells, and (c) Cell-to-cell contact, where sprout expansion is guided by contact with multicellular elongated structures or projections of other cells (Czirok, 2013).

Arteriogenesis increases the diameter of existing blood vessels and remodels large blood vessels creating natural bypasses when necessary. Whenever, blood flow is redirected into preexisting arterioles, it creates mechanical forces. Elevated shear stress and circumferential wall stress during a long time period are strong inducers of arteriogenesis (Heil et al., 2006). The endothelium of the arteriolar connections is activated by the mechanical forces, causing monocytes to promote arteriogenesis by secreting growth factors and cytokines that increase the mitosis rate of endothelial and smooth muscle cells (Deindl and Schaper, 2005). Perivascular mast cells mediate shear stress-induced arteriogenesis by coordinating the action of T cells, neutrophils, monocytes, macrophages, and other innate immune cells by means of the secretion of cytokines and MMPs. The activation of perivascular mast cells is achieved by the increase of Nox2-derived reactive oxygen radicals, caused by neutrophil extravasation (Chillo et al., 2016).

Angiogenesis extends, maintains, and remodels existing networks of thin blood vessels, mostly capillaries. There exist two main mechanisms for angiogenesis, namely, sprouting angiogenesis (SA), and splitting angiogenesis, also known as intussusceptive angiogenesis (IA) (Gianni-Barrera et al., 2011). Alterations in blood flow and local changes in the concentration of angiogenic factors such as VEGF may trigger angiogenesis. Laminar shear stress inhibits tubule formation and migration of endothelial cells and favors intussusceptive angiogenesis, while turbulent shear stress causes an increase in cell migration and proliferation, and favors sprouting angiogenesis (Makanya et al., 2009). In skeletal muscle, VEGFA164 over-expression induces vascular growth by intussusception rather than sprouting (Gianni-Barrera et al., 2013).

IA occurs during physiological adaptation i.e., exercised muscles, embryonic development, and pathological situations such as tumor growth. During IA, endothelial cells extend processes into the vascular lumen from opposing walls. Once these processes contact each other, the endothelial cell junctions at the contact site are reorganized. Then, the bilayer is perforated by invading interstitial tissue, pericytes, and myofibroblasts, forming a transluminal pillar. Subsequently, pericytes, fibroblasts, and other supporting cells deposit additional collagen and other stabilizing fibers into the extracellular matrix of the pillar (Makanya et al., 2009), that may increase in girth, until it splits the blood vessel into two independent vascular segments (Patan et al., 1996, 1997). Additionally, several transluminal pillars may fuse to split a vessel or improve local hemodynamic behavior (Kurz et al., 2003). IA has three main advantages over SA: first, IA is achieved with minimal tissue degradation and reduced vascular permeability caused by mural cell detachment, second, a relatively short period of time is sufficient to achieve it, and third, only a relatively low rate of endothelial proliferation is needed (Kurz et al., 2003; Makanya et al., 2009; Gianni-Barrera et al., 2011). IA is necessary for the formation of organ-specific angioarchitecture (intussusceptive microvascular growth), the formation of vascular trees (intussusceptive arborization), angioadaptation and vascular pruning (intussusceptive branching remodeling) (Makanya et al., 2009).

SA is a developmental process that results in a new connection between two existing thin blood vessels (Figure 1) and involves eight related cellular processes: (1) Secretion of angiogenic factors. Shear stress, or an insufficient local supply of oxygen or nutrients, may cause the cells within a tissue to secrete angiogenic factors (Forsythe et al., 1996; Song and Munn, 2011; Kumar et al., 2014). Relevant angiogenic factors include growth factors, chemokines, angiopoietins, endostatin, interferons, and NO among other molecules (Logsdon et al., 2014). (2) Vessel destabilization. Before a new sprout may form, pericytes, myofibroblasts, and other supporting cells must be cleared from the area of the blood vessel where the sprout will form. Also, the ECM surrounding the area must be remodeled. Blood vessel destabilization is mediated by VEGFA, ANG2, NO, and the absence of blood flow (Scharpfenecker et al., 2005; Qin et al., 2013; Korn and Augustin, 2015). (3) Tip and stalk cell differentiation. When certain ECs are exposed to a VEGF gradient some respond to VEGFA and shear stress to become tip cells (TCs), growing filopodia toward the VEGFA gradient. TCs induce neighbor cells to become stalk cells (SCs) by Notch-mediated lateral signaling (Blanco and Gerhardt, 2013). TCs become less sensitive to Notch signaling and SCs become less sensitive to VEGF signaling (Weinstein et al., 2015; Glass et al., 2016). (4) Sprout elongation. The sprout is initially formed by the TC and one or two adjacent SCs. The subsequent proliferation of both the TC and SCs together with SC elongation and rearrangement support stalk elongation toward the VEGFA source resulting in stalk growth (Betz et al., 2016). (5) Lumen formation and expansion. Lumen formation may occur through cord hollowing, cell hollowing, trans-cellular lumen formation, and lumen ensheathment. Hemodynamic forces shape the apical membrane of SCs to form and expand new lumenized vascular tubes (Betz et al., 2016). (6) Anastomosis. Vascular anastomosis is the process that allows angiogenic sprouts and blood vessels to connect. Anastomosis can occur between two sprouts, or between a sprout and a functional blood vessel. The first step in an anastomosis is the formation of a stable contact between two ECs forming a new adherens junction with two layers of apical membrane and a small luminal volume in between. The mechanism that allows the formation of a new multicellular, perfused tubes depends on the presence or absence of blood pressure (Betz et al., 2016). (7) Vessel stabilization. Once a lumenized new blood vessel has formed, ECs release platelet-derived growth factor B (PDGFB). PDGFB attracts pericytes, which incorporate into the vessel wall. S1P, S1PR1, ANG1, TIE2, Ephrin-B2, EPH, and TGFβ regulate blood vessel stabilization and maturation and are regulated by shear stress (Scharpfenecker et al., 2005; Qin et al., 2013; Korn and Augustin, 2015). And (8) Pruning. Vessel pruning is basically the process of sprout formation in reverse. The absence of blood flow, or a higher anti-angiogenic (ANG2) to angiogenic (VEGFA) factor ratio, induces small blood vessel pruning by reabsorption of ECs into the remaining vasculature. Regression of larger blood vessels involves apoptosis (Korn and Augustin, 2015; Betz et al., 2016).

FIGURE 1
www.frontiersin.org

Figure 1. (A) Hypoxia induced angiogenesis: When tissue cells are exposed to a microenvironment with an insufficient concentration of Oxygen, they secrete VEGFA in a process mediated by the Hypoxia-inducible factor 1 (HIF-1). Forming a VEGFA gradient (green), (B) Certain epithelial cells (peach) respond to VEGFA and shear stress to become tip cells (TCs): VEGFA, ANG2, shear stress, and NO lead to endothelial cell matrix degradation, loss of pericytes (brown triangular cells). Certain EC become TCs (turquoise) and grow filopodia toward the VEGFA gradient. TCs inhibit neighboring cells from becoming TCs by Notch mediated lateral signaling and Wnt, (C) Stalk growth and anastomosis: The cells neighboring the TCs are induced by Notch to become Stalk cells (SCs). SCs (orange) secrete VEGFR1, reducing the concentration of VEGFA in their microenvironment, undergo Wnt mediated proliferation and elongate toward the VEGFA source resulting in stalk growth. Once a TC reaches another TC or vessel wall, it undergoes VE-cadherin and Macrophage mediated binding, initiating anastomosis, (D) Lumen formation: Lumen formation may occur through cord hollowing (Intracellular vacuoles fuse intracellularly to hollow out stalk cells and generate an interconnected luminal space), cell hollowing, transcellular lumen formation, and lumen ensheatment. Hemodynamic forces shape the apical membrane of SCs to form and expand new lumenized vascular tubes, (E) Vessel stabilization: Once a lumenized new blood vessel has formed, ECs release platelet-derived growth factor B (PDGFB). PDGFB attracts pericytes which incorporate into the vessel wall. S1P, S1PR1, ANG1, TIE2, Ephirin-B2, EPH, and TGF regulate blood vessel stabilization and maturation and are regulated by shear stress.

Due to the enormous biological and medical importance of angiogenesis, many computational and mathematical models have been proposed to explore the molecular mechanism involved in angiogenesis control (Peirce, 2008; Qutub et al., 2009; Scianna et al., 2013; Logsdon et al., 2014; Heck et al., 2015; Qutub and Popel, 2015). Some of the relevant previous models are the following: (a) A computational model exploring the relationship between hemodynamics and angiogenesis in 2D (Gödde and Kurz, 2001). (b) A computational model of oxygen transport in skeletal muscle for sprouting and splitting modes of angiogenesis (Ji et al., 2006). (c) A model that describes and explores the progression of angiogenesis during the healing process (Vermolen and Javierre, 2012). (d) A multicellular model of the early stages of angiogenesis using finite element integration that includes chemotaxis and the interaction between tip cells and stalk cells (Bookholt et al., 2016). (e) Two Boolean models that explore the relationship between the Wnt and VEGF signaling pathways, mechanoreceptors, apoptosis, cell proliferation, and lumen formation during angiogenesis (Bauer et al., 2010; Bazmara et al., 2016). And (f) a multilevel model based on the previously mentioned Boolean models (Bazmara et al., 2015). Notably, the authors of Bauer et al. (2010) and Bazmara et al. (2015) included apoptosis in their model. We did not include apoptosis in our model because thin blood vessel pruning usually occurs by the reabsorption of ECs into the remaining vasculature and seldom involves apoptosis.

Previous models of angiogenesis focused on the role played by a few of the canonical signaling pathways. However, recent discoveries have emphasized the role of TGF signaling and its interaction with the WNT, NOTCH, VEGF/NRP1, HIF, AKT, ERK, mTOR, and TIE signaling pathways, as well as the role of HIFs, Ca2+, NO/eNOS, and cytoskeletal mechanoreceptors during angiogenesis. As a result, none of the previous models explore the interaction among all the aforementioned canonical pathways. It was not possible to know, then, if the biological system was sufficiently well characterized from the point of view of the molecular regulation. Hereby we present a model that integrates the largest set of canonical signaling pathways, thus allowing for a comprehensive characterization of the effect of the extracellular micro-environment on EC behavior during differentiation of ECs angiogenesis. The model presents a qualitative agreement with a large set of experimental published results, showing that the regulatory network is a faithful reconstruction of the central molecular mechanism controlling the cell behavior of endothelial cells during angiogenesis. This characteristic permits the use of the model not only to describe the wild-type development and adaptation process but also to propose targets for intervention in certain diseases. Specifically, our model suggests that favoring a micro-environment that induces Phalanx EC behavior may suffice to improve the treatment of vascular retinal disorders and vascular malformations. Thus, our model can be considered as a platform to study several molecular scenarios affecting angiogenesis.

2. Materials and Methods

2.1. Molecular Basis of the Regulatory Network

To assemble our model of the molecular network involved in angiogenesis control, we first explored how each one of the main stages angiogenesis is regulated and then explored how the molecules involved in the control of each stage interact with those that regulate the other stages. We started by exploring how the ANG/TIE signaling pathway acts as a gatekeeper of EC quiescence. Next, we inquired into the mechanisms that allow lack of sufficient oxygen or nutrients to destabilize blood vessels and trigger the angiogenic process. Then, we probed the mechanism that allows certain EC to be more sensitive to angiogenic signals by regulating VEGFR activity. Later, we analyzed how VEGF signaling activates the signaling pathways ERK1/2, PI3K-AKT, SRC, and p38 MAPK, and additionally phosphorylates STATs. After that, we inquired into the mechanisms that allow mechanoreceptors to respond to shear stress and radial stress to regulate VEGF signaling. Our ensuing action was to scrutinize the mechanism that allows the VEGF, NOTCH, WNT, and TGF signaling pathways to interact and regulate tip and stalk EC behavior. Last, we explored the mechanism that allows NOTCH and WNT to regulate EC proliferation. All those molecular mechanisms, their interactions and some of the most relevant references that describe the experimental evidence are discussed in detail in the first section of the Supplementary Material.

2.2. The Regulatory Network as a Discrete Dynamical System

Boolean networks are discrete dynamical systems, whose simplicity allows the attainment of biologically meaningful results, after a systematic exploration of its dynamical behavior (Dubrova and Teslenko, 2011; Azpeitia et al., 2017). In our model, most variables represent genes or proteins, some represent small molecules, and one represents a mechanical force. Each variable has an activation state, that may be active, represented by a 1, or inactive, represented by a 0. Furthermore, we use a synchronous update approach where the states of all the variables are updated simultaneously. We decided to use a synchronous update scheme in our boolean model because the computational analysis of the asynchronous update is extremely time-consuming, and it is mostly required to explore race conditions and cyclic patterns of molecular activation (Garg et al., 2008; Saadatpour et al., 2010). However, neither race conditions nor cyclic behaviors are explored with our current model.

We use definitions and notation for Boolean networks based on Azpeitia et al. (2017). Let 𝔹 = {0, 1} and Nn+={1,2,,n}, a set of labels. A synchronous Boolean network with n components is a function f : 𝔹n→𝔹n, where the i-th component of f is a function fi : 𝔹n→𝔹n such that fi(x) = f(x)i. A state of the network is an n-tuple x = (x1, x2, …, xn) such that x ∈ 𝔹n. To relate a synchronous Boolean network with a molecular network, we interpret that each component of a state x represents the activation state of a variable denoting a molecule included in the network. The dependency of the state on the discrete time parameter t is denoted as x(t) and obeys the update rule given by f. That is for all t ∈ ℤ

x(t+1)=f(x(t))=(f1(x(t)),f2(x(t)),,fn(x(t))),

where

xi(t+1)=fi(x(t)).

Our Boolean network model is deterministic, and any given initial state of the network reaches an attractor. A fixed point attractor is a state s ∈ 𝔹n such that f(s) = s. We define fol as the l-th iterate of the function f such that fol = f(fo(l−1)). An attractor is a set of states A ⊆ 𝔹n, such that fol(x) = x for any state xA, in other terms, there exist a positive natural number l ∈ ℕ+ = {1, 2, …} such that f(x(t + l)) = f(x(t)) for all x(t) ∈ A. Furthermore, l is the size of the attractor and for any j<l+, f(x(t + j)) ∈ A. Fixed point attractors represent stationary patterns of molecular activation, while larger attractors represent cyclic patterns of molecular activation. Additionally, we assume that each attractor represents an EC behavior.

For simplicity, we refer to the variable xi by its position i in the n-tuple x. For a state x ∈ 𝔹n and one of its components, say the one with label i, we denote by x ~ i the n-tuple resulting from replacing the value of the variable xi by its complement. Given two variables i and j and the update function of variable i, namely fi, j activates i if there exists a pair of network states x, y that differ only in the state of activation of variable j, that is y = x ~ j, xj = 0 and yj = 1, such that fi(y)−fi(x) > 0. Conversely, j inhibits i if there exists a pair of network states x, y that differ only in the state of activation of variable j, that is y = x ~ j, xj = 0 and yj = 1, such that fi(y)−fi(x) < 0. An interaction denoted as the pair (i, j), i, j ∈ ℕn is functional if variable j activates or inhibits variable i. Note that according to this definition, it is possible for variable j to both activate and inhibit variable i depending on the functional context. For instance, let C(t+1) = (A(t)∧¬B(t))∨(¬A(t)∧B(t)). A activates C because if we focus on the cases where B is not active; if A is active, then C is active. A also inhibits C because if we focus on the cases where B is active; C is active only when A is not active.

2.3. Model Assembly

Using the information described in the subsection Molecular basis of the network, we assembled a model of the molecular network involved in angiogenesis control. Then we encoded the model using the standardized text file format required by BoolNet (Müssel et al., 2010), an R package for the analysis of Boolean networks. The models in BoolNet format, and the R scripts we used to simulate and analyze the dynamic behavior of the model are available at https://github.com/NathanWeinstein/Angiogenesis-Model/. During the development of our model, we ensured the existence of stable or cyclic patterns of molecular activation that correspond to Phalanx (AKT+, JAGa−, NRP1−), Stalk(JAGa+, NRP1−), and Tip (NRP+, DLL4a+, AKT−) EC behavior and their reachability under certain micro-environmental conditions (Figure 4A); specifically:

1. (VEGFC_Dp−, VEGFAxxxP−, ANG1+, Oxygen+, ShearStress+, JAGp−, DLL4p−, WNT5a−, WNT7a−, FGF−, IGF−, BMP9−, BMP10−, TGFB1−, VEGFC_D−, and AMPATP−) induces Phalanx EC behavior.

2. (VEGFC_Dp+, VEGFAxxxP−, ANG1+, Oxygen+, ShearStress+, JAGp−, DLL4p−, WNT5a−, WNT7a−, FGF−, IGF−, BMP9−, BMP10−, TGFB1−, VEGFC_D−, and AMPATP−) induces Tip EC behavior.

3. (VEGFC_Dp−, VEGFAxxxP+, ANG1+, Oxygen+, ShearStress+, JAGp−, DLL4p−, WNT5a−, WNT7a−, FGF−, IGF−, BMP9−, BMP10−, TGFB1−, VEGFC_D−, and AMPATP−) induces Tip EC behavior.

4. (VEGFC_Dp−, VEGFAxxxP−, ANG1+, Oxygen+, ShearStress+, JAGp−, DLL4p+, WNT5a+, WNT7a−, FGF−, IGF−, BMP9−, BMP10−, TGFB1+, VEGFC_D−, and AMPATP−) induces Stalk EC behavior.

Importantly, the expected patterns of molecular activation and EC behavior transitions are based on the literature (del Toro et al., 2010; Blancas et al., 2012; Glaser et al., 2016).

2.3.1. Simulation of an EC Behavior Transition

To simulate the transitions in EC behavior, we started with one of the states of an attractor that represents the initial EC behavior. Then, we modified the variables that represent the extracellular micro-environment (VEGFC_Dp, VEGFAxxxP, ANG1, Oxygen, ShearStress, JAGp, DLL4p, WNT5a, WNT7a, FGF, IGF, BMP9, BMP10, TGFB1, VEGFC_D, and AMPATP) without changing the other variables related to the internal state of the cell, to simulate a change of micro-environment that should lead to another EC behavior. We then calculated all state transitions until reaching another attractor that represents a new EC behavior.

2.3.2. Boolean Network Simplification

The size of the state space of a boolean molecular network represented as a directed graph with n nodes —one node for each variable—, grows exponentially as 2n. Simulating and analyzing the dynamic behavior of networks containing more than a hundred nodes requires considerable computational resources. Recently, certain algorithms that use boolean satisfiability (SAT) methods capable of finding the attractors of networks with hundreds of nodes have been developed and implemented (Dubrova and Teslenko, 2011). Nonetheless, analyzing the effects of mutations, summarizing trap spaces, and analyzing the robustness of large networks is still a challenging task. However, it has been proved that it is possible to remove inputs and nodes with both an indegree and an outdegree equal to one without affecting the attractors (Saadatpour et al., 2013). Accordingly, we simplified the model by removing input nodes (nodes with an indegree equal to zero) that are either active, or inactive in all ECs, and are not part of the parameters that specify an extracellular micro-environment. Additionally, we removed output nodes (nodes with outdegree equal to zero). Further, we used edge contraction to merge intermediary nodes (nodes that have either an indegree or outdegree equal to one) that are not transcription factors. The edge contraction operation involves the removal of an edge e (from u to v) and the merger of its two incident vertices, u and v, into a new vertex w. We assigned to w the name of u if v was only regulated by u, in this case we substituted v(t) for u(t) if e was positive or ¬u(t) if e was negative in the components of f that correspond to the variables originally regulated by v. When u only regulated v, we assigned to w the name of v and in fv we substituted u(t) with fu, that is, fv(…, u(t), …) becomes fv(,fu,). These operations allowed us to simplify our model without eliminating feedback circuits. This is relevant because to a large extent, feedback circuits determine the dynamic behavior of a boolean network (Azpeitia et al., 2017). The authors of Veliz-Cuba (2011) and Naldi et al. (2011) studied when the attractors are preserved after similar simplification processes. Additionally, we verified that the EC behaviors and transitions based on the literature were preserved after the simplification process. Further, we also verified that in both the detailed and the simplified model, all single gain and loss of function mutations have a similar effect on the EC behaviors and transitions based on the literature (Supplementary Figures 15, 16). Note that for this verification we only simulated the effect of 4 micro-environments. We only simulated the full effect of the mutations on our simplified model as part of the model validation process.

2.4. Analysis of the Dynamic Behavior of Our Model

First, we obtained all the attractors using the exhaustive SAT-based search available as part of BoolNet that uses an adaption of the algorithm proposed by Dubrova and Teslenko (2011). The exhaustive SAT-based search formulates the attractor search as a boolean satisfiability (SAT) problem that is solved using the PicoSAT solver (Biere, 2008). Then, we classified the attractors based on extracellular micro–environment. After that, for each micro-environment, we inferred the EC behavior represented by each attractor. If all EC behaviors associated to one micro-environment where of the same kind, we added that micro-environment to the set of micro-environments that induce that EC behavior. If not all EC behaviors associated with one micro-environment where of the same kind, we added the micro-environment to the set of micro-environments that induce atypical EC behavior. Finally, we summarized the four sets of micro-environments by grouping them into disjoint subsets using their shared characteristics. To validate our model, we simulated all single gain and loss of function mutations. We then compared the simulated effect of each mutation with its experimentally observed effect as reported in the literature (when available). Biological organisms need to be resilient to mutations and fluctuations in the concentration or level of molecular activation, we refer to this property as robustness. For clarity, it is necessary to indicate which trait is robust, to which perturbation and a method to quantify the resilience to define a robust feature (Félix and Barkoulas, 2015). We measured three robust features: (1) The robustness of Phalanx, Stalk, and Tip EC behaviors to single gain and loss-of-function mutations. This was measured as the percentage of mutations that prevent the existence of any stable or cyclic patterns of molecular activation that correspond to said EC behavior. (2) The robustness of attractor determination to molecular activation noise. First, we generated a set of 1,000,000 aleatory initial states. For each initial state, we created a perturbed copy with a Hamming distance of one by reversing the activatory state of one random variable. We quantified attractor determination robustness to molecular activation noise, as the fraction of the initial states that reached the same attractor as their perturbed copies. (3) The robustness of the trajectories that lead to Phalanx, Stalk, and Tip EC behaviors to molecular activation noise. First, we generated a set of 1,000,000 aleatory initial states. For each initial state, we created a perturbed copy with a Hamming distance of one by reversing the activatory state of one random variable. We quantified the robustness of the EC behaviors to molecular activation noise, as the fraction of the initial states that reached an attractor that represents the same EC behavior as that of their perturbed copies. Additionally, we calculated the sensitivity of each component of the update rule to molecular activation noise. For each update rule component fif, we first generated a set of 500,000 aleatory initial states. For each initial state, we created a perturbed copy with a Hamming distance of one by reversing the activatory state of one random variable. Then we applied the update rule once to each initial state and to its perturbed copy. The fraction of initial states, where after update the activatory state of the variable xi(t + 1) is different for the initial state then it is for its perturbed copy is our estimation of the sensitivity for update rule fi.

3. Results

3.1. The Network Model

The model includes 143 nodes and 256 edges (Figure 2) the update rules of the network are included in Supplementary Section 2. To enable a more thorough analysis of the dynamic behavior of our model, we simplified our model and obtained a network composed of 64 nodes and 163 interactions, a diagram of our simplified model is shown in Figure 3. The update rules that define the dynamic behavior of our model are included as Equations 1–64. The EC behavior transitions integrated into both our detailed and simplified models are summarized in Figure 4A, and Supplementary Figures 1–14. Single gain- and loss-of-function mutations have a similar effect on the behaviors and transitions integrated into both models (Supplementary Figures 15, 16).

AKT(t+1)=PIP3(t)    (1)
ALK1(t+1)=BMP9(t)BMP10(t)TGFB1(t)    (2)
ALK5(t+1)=BMP9(t)    (3)
AMPATP(t+1)=AMPATP(t)    (4)
AMPK(t+1)=(AMPATP(t)(¬Oxygen(t)))                          (¬AKT(t))    (5)
ANG1(t+1)=ANG1(t)    (6)
ANG2(t+1)=(¬KLF2(t))(HIF1(t)ETS(t)                          AP1(t)FOXO1(t))    (7)
AP1(t+1)=WNT5a(t)    (8)
βcatenin(t+1)=WNT5a(t)WNT7a(t)    (9)
BMP10(t+1)=BMP10(t)    (10)
BMP9(t+1)=BMP9(t)    (11)
Calcium(t+1)=PLCg(t)ShearStress(t)(¬NO(t))    (12)
DLL4a(t+1)=ETS(t)NICD(t)    (13)
DLL4p(t+1)=DLL4p(t)    (14)
ETS(t+1)=MEK(t)VEGFR33(t)    (15)
FAK(t+1)=SRC(t)Integrin(t)    (16)
FGF(t+1)=FGF(t)    (17)
FOXO1(t+1)=(¬AKT)SIRT1(t)    (18)
HEY1(t+1)=NICD(t)((SMAD1(t)SMAD2(t))                        (¬SIRT1(t)))    (19)
HIF1(t+1)=(AMPK(t)¬TSC(t))¬Oxygen(t)                        SIRT1(t)    (20)
IGF(t+1)=IGF(t)    (21)
Integrin(t+1)=ETS(t)(ShearStress(t)TIE2(t))    (22)
JAGa(t+1)=SMAD1(t)βcatenin(t)    (23)
JAGp(t+1)=JAGp(t)    (24)
KLF2(t+1)=ShearStress(t)    (25)
LEF1(t+1)=βcatenin(t)(LEF1(t)NRARP(t))    (26)
MEK(t+1)=(((PLCg(t)Calcium(t))RAS(t))                       (¬AKT(t)))FGF(t)    (27)
NFAT(t+1)=Calcium(t)    (28)
NICD(t+1)=(¬NRARP(t))NOTCH(t)    (29)
NO(t+1)=Calcium(t)AKT(t)SIRT1(t)    (30)
NOTCH(t+1)=(¬JAGp(t))ETS(t)DLL4p(t)    (31)
NRP1(t+1)=(VEGFAxxx(t)VEGFC_Dp(t))                        (¬NICD(t)ETS(t))    (32)
NRARP(t+1)=NICD(t)    (33)
Oxygen(t+1)=Oxygen(t)    (34)
p38MAPK(t+1)=SRC(t)PLCg(t)ALK1(t)    (35)
PECAM1(t+1)=VEGFR22(t)VEGFR23(t)                              ShearStress(t)VEcadherin(t)    (36)
PIP3(t+1)=(¬NICD(t))(¬WNT5a(t))                        (¬WNT7a(t))(¬NRP1(t))                        (SRC(t)KLF2(t)                        VEcadherin(t)TIE2(t))    (37)
PLCg(t+1)=VEGFR22(t)VEGFR33(t)                         WNT5a(t)WNT7a(t)    (38)
RAS(t+1)=PECAM1(t)KLF2(t)ALK1(t)    (39)
ShearStress(t+1)=ShearStress(t)    (40)
SIRT1(t+1)=AMPK(t)(HIF1(t)FOXO1(t))    (41)
SMAD1(t+1)=(¬SMAD6(t))(¬NRP1(t))ALK1(t)    (42)
SMAD2(t+1)=(¬SMAD6(t))(¬NRP1(t))ALK5(t)    (43)
SMAD6(t+1)=NICD(t)    (44)
SRC(t+1)=FAK(t)ShearStress(t)VEGFR22(t)                       VEGFR23(t)    (45)
STAT3(t+1)=VEGFR22(t)    (46)
TGFB1(t+1)=TGFB1(t)    (47)
TIE2(t+1)=(¬ANG2(t))ANG1(t)(ETS(t)                         KLF2(t))    (48)
TSC(t+1)=AMPK(t)(¬AKT(t))    (49)
VEcadherin(t+1)=ETS(t)(SRC(t)FAK(t)                                    ShearStress(t)HIF1(t))    (50)
VegfA(t+1)=(¬Oxygen(t))HIF1(t)STAT3(t)                            FOXO1(t)NFAT(t)KLF2(t)    (51)
VEGFAxxx(t+1)=VEGFAxxxP(t)VEGFAxxxA(t)    (52)
VEGFAxxxA(t+1)=VegfA(t)IGF(t)((¬NICD(t)                                     ¬HIF1(t)¬ETS(t))                                     NFAT(t))(¬AMPK(t))    (53)
VEGFAxxxd(t+1)=p38MAPK(t)VegfA(t)    (54)
VEGFAxxxP(t+1)=VEGFAxxxP(t)    (55)
VEGFC_D(t+1)=VEGFC_D(t)    (56)
VEGFC_Dp(t+1)=VEGFC_Dp(t)    (57)
Vegfr2(t+1)=(ETS(t)(¬HEY1(t)))¬Oxygen(t)    (58)
VEGFR22(t+1)=Vegfr2(t)(PECAM1(t)                                ((VEGFC_Dp(t)VEGFAxxx(t))                                ¬(VEGFAxxxd(t)HIF1(t)))    (59)
VEGFR23(t+1)=Vegfr2(t)Vegfr3(t)(PECAM1(t)                                  VEGFAxxx(t)VEGFC_Dp(t))    (60)
Vegfr3(t+1)=NICD(t)    (61)
VEGFR33(t+1)=Vegfr3(t)(PECAM1(t)                                 VEGFC_D(t)VEGFC_Dp(t))    (62)
WNT5a(t+1)=WNT5a(t)    (63)
WNT7a(t+1)=WNT7a(t)    (64)
FIGURE 2
www.frontiersin.org

Figure 2. A diagram of our extended model: The ANG/TIE signaling pathway is shown in gray, Shear Stress in white, Oxygen and Energy in blue, NO in turquoise, VEGF in yellow, AKT/SRC in light blue, TGF in pink, NOTCH in orange, WNT in purple, RAS/PLCγ in violet, CyclinD1 in light green, and FGF in green. Ligands are represented as hexagons, other micro-environment variables as octagons, receptors as right arrows, transcription factors as ellipses, and signal transducers as rounded rectangles. Intracellular signaling is represented in black arrows, extracellular signaling is represented with blue arrows. Activatory interactions are shown as regular arrows and inhibitory interactions are shown as blunt arrows.

FIGURE 3
www.frontiersin.org

Figure 3. A diagram of our reduced network model: The ANG/TIE signaling pathway is shown in gray, Shear Stress in white, Oxygen and Energy in blue, VEGF in yellow, AKT/SRC in light blue, TGF in pink, NOTCH in orange, WNT in purple, RAS/PLCγ in violet, and FGF in green. Ligands are represented as hexagons, other micro-environment variables as octagons, receptors as right arrows, transcription factors as ellipses, and signal transducers as rounded rectangles. Intracellular signaling is represented in black arrows, extracellular signaling is represented with blue arrows. Activatory interactions are shown as regular arrows and inhibitory interactions are shown as blunt arrows. The self activatory feedback circuits required to keep the micro-environment constant during the simulation are shown in red.

FIGURE 4
www.frontiersin.org

Figure 4. Endothelial cell behavior: Phalanx EC behavior shown in yellow, Stalk EC behavior shown in orange, Tip EC behavior shown in green, and other EC behavior is shown in gray: (A) Expected EC behavior in an extracellular micro-environment with normal oxygen concentration, ATP to ADP ratio and shear stress, (B) The extracellular micro-environments that cause Phalanx, Stalk, and Tip EC behavior according to the simulation of the dynamic behavior of our simplified model, (C) The extracellular micro-environments that cause other EC behavior, (D) Summary of EC behavior according to our model, the numbers on the edges represent the groups of micro-environments shown as columns in panel (B).

3.2. The Effect of the Extracellular Micro-environment on EC Behavior

One of the main goals of this work is to understand how the concentration of several molecules in the extracellular micro-environment combines with the mechanical forces sensed by the mechano-receptors connected to the cytoskeleton of ECs controls EC behavior. We propose that the presence (1) or absence (0) of sufficient VEGFC_Dp, VEGFAxxxP, ANG1, Oxygen, ShearStress, JAGp, DLL4p, WNT5a, WNT7a, FGF, IGF, BMP9, BMP10, TGFB1, VEGFC_D, and AMPATP in the micro-environment of an EC determines its behavior. Further, we propose that Phalanx, Stalk, and Tip ECs retain the ability respond to the micro-environment in a similar manner and that explains the plasticity in EC behavior that has been experimentally observed (Blancas et al., 2012; Glaser et al., 2016). To investigate the effect of the extracellular micro-environment on EC behavior, we first found all the attractors that can be reached through the simulation of the dynamic behavior of our model. Then, we classified them according to their extracellular micro-environment. After that we interpreted the EC behavior represented by the attractors in each micro-environment. If all the attractors that correspond to a certain micro-environment represent the same kind of EC behavior, then we can state that the micro-environment causes that EC behavior. If most micro-environments cause either Tip, Stalk, or Phalanx EC behavior, then to a large extent the extracellular micro-environment controls EC behavior.

Notably, there are 216 = 65536 possible micro-environments. From these, according to our model, under wild-type conditions 50,572 (77.16675%) micro-environments cause Tip EC behavior, 12,096 (18.45703%) cause Stalk EC behavior, and 96 (0.1464844%) cause Phalanx EC behavior. The characteristics of the groups of micro-environments that lead to Phalanx, Stalk, and Tip EC behavior are summarized in Figure 4B and Table 1. The intracellular molecules that are active or inactive in all the patterns of molecular activation in each group are also summarized in Table 1. The other 2,772 micro-environments (4.229736%) cause atypical dynamical patterns, including attractors that cycle between the Tip, Stalk, and/or Phalanx EC behaviors. This means that 62,764 (95.770374%) of the micro-environments induce a certain EC behavior regardless of the internal pattern of molecular activation (Figure 4D). Therefore, according to our model, in most cases, the extracellular micro-environment controls EC behavior.

TABLE 1
www.frontiersin.org

Table 1. Phalanx, Stalk, and Tip EC behavior: The groups correspond to those in Figure 4B, active molecules shown in blue, inactive molecules shown in red.

Tip ECs are localized at the leading edge of vessel sprouts forming numerous long dynamic filipodia. Additionally, Tip cells migrate toward angiogenic stimuli, do not contribute to lumen formation, and seldom divide. Tip ECs are characterized by expressing high levels of DLL4, CXCR4, ANG2, PDGFB, receptors for axon guidance cues, such as the Netrin receptor UNC5B, APLN, various proteases like uPAR and NRP1, (del Toro et al., 2010; Blancas et al., 2012). We use NRP1 activity as a Tip EC-specific marker, and also require DLL4 expression, because DLL4 directly inhibits neighboring cells from becoming Tip ECs. Additionally, AKT must be inactive in Tip ECs. It is known that an increase above a certain threshold on the concentration of VEGFA or proteolytically processed VEGFC or D in the micro-environment surrounding an EC triggers Tip EC behavior (sections 1.2 and 1.10 in the Supplementary Material). According to the simulated dynamic behavior of our model, the micro-environments that include VEGFAxxxP or VEGFC_Dp and induce Tip EC behavior, also include either ShearStress, WNT5a, WNT7a, FGF, BMP9, BMP10, or TGFB1. Alternatively, the model also allows for the possibility that two groups of micro-environments that lack paracrine VEGF activity may cause Tip EC behavior, achieved by inducing autocrine VEGFA activity.

Stalk ECs trail Tip ECs, proliferate rapidly and contribute to lumen formation. While TIE2 is constitutively expressed in all ECs, its protein is detectable by antibody staining on Stalk ECs but not on Tip ECs. Stalk cells also express the Apelin receptor APJ and JAG1 (del Toro et al., 2010; Blancas et al., 2012). We use autocrine JAG1 as a Stalk EC marker due to the specificity of its expression and its function, which is to suppress Notch signaling in neighboring Tip ECs, further, Stalk ECs, are characterized by their lack of NRP1 activity. A sufficient concentration of WNT, TGF and NOTCH ligands, as well as an absence of VEGF in the extracellular micro-environment of an EC, is known to cause Stalk EC behavior (section 1.10 in the Supplementary Material). According to the simulated dynamic behavior of our model, it is possible to obtain the Stalk EC behavior in a micro-environment that complies with either of the following three lists of requirements: (a) WNT activity, lack of VEGF activity, and low Oxygen or IGF; (b) WNT activity, no VEGF activity, Oxygen, IGF, and sufficient energy; and (c) lack of VEGF, NOTCH, WNT, and IGF ligands that includes one of the TGF ligands.

Phalanx ECs form strong EC–EC bonds to compose the tunica intima in stable blood vessels. The Pericytes and SMCs that cover stable blood vessels secrete ANG1 to maintain the integrity of the layer of Phalanx ECs. Phalanx ECs are characterized by a high level of VEGFR1 (FLT1) and TIE1 expression (Blancas et al., 2012), even though neither is a Phalanx EC specific marker. We use active AKT (Kerr et al., 2016) as well as inactive NRP1 and JAGa as specific Phalanx EC markers. Changes in the extracellular concentration of VEGFs, a decrease in the availability of oxygen or energy within the cell, and shear stress cause ANG2-mediated activation of the ECs that line blood vessels (section 1.8 in the Supplementary Material). According to our model, the lack of VEGF, NOTCH, WNT and TGF pathway activity is necessary to observe a stable Phalanx EC behavior. The simulated Phalanx ECs do not divide and do not recruit mural cells.

3.2.1. Atypical EC Behavior

We performed with our model a systematic study of the dynamical behavior of a regulatory network under all possible combination of the micro-environments. Apart from the clearly identifiable phenotypes mentioned in the previous paragraphs, we observed some atypical responses. If the attractors that correspond to a certain micro-environment represented different EC behaviors, or any of the attractors represented an EC behavior that was different from Tip, Stalk, or Phalanx EC behavior, we considered that the micro-environment causes atypical EC behavior. For completeness, we describe such atypical behaviors in Table 2.

TABLE 2
www.frontiersin.org

Table 2. Atypical EC behavior: The groups correspond to those in Figure 4C, active molecules shown in blue, inactive molecules shown in red.

3.2.2. EC Proliferation

EC proliferation allows the number of ECs to increase during sprout elongation. We describe the effect of the micro-environment on EC proliferation according to the simulated dynamic behavior of our model in Table 3. Note that in accordance with what has been reported in the literature, only Tip and Stalk ECs proliferate.

TABLE 3
www.frontiersin.org

Table 3. EC proliferation: Cyclin D1–mediated activation of the cell cycle requires βcatenin and LEF1 activity. Active molecules shown in blue, inactive molecules shown in red.

3.3. Model Validation

Certain diseases exhibit abnormal angiogenesis, because the affected tissue or organ secretes abnormal amounts of angiogenic signals. Simulating the effect of a pathological extracellular micro-environment on EC behavior can be used to understand how a disease is causing abnormal vascular remodeling, the insights are only valid if the dynamic behavior of the model can reproduce the relevant experimental observations. If an experimental observation includes a sufficiently well-defined extracellular micro-environment and an observed EC behavior. Then the extracellular micro-environment fits only one column in Figure 4B or Figure 4C. If the EC behavior according to our model (shown at the bottom row of the column that corresponds to the micro-environment) is the same as the observed EC behavior, then our model fits that experimental observation.

3.3.1. Tumor Angiogenesis

The micro-environment inside many tumors is hypoxic, containing a high concentration of VEGFA and FGF. This state causes the formation of many leaky blood vessels (Nussenbaum and Herman, 2010). Our model can describe this state, as shown in Figure 4B group 27. The results indicate that the mentioned micro-environment induces Tip EC behavior, and inhibits Phalanx EC behavior, which is consistent with experimental observations.

3.3.2. Pathological Ocular Angiogenesis

Diabetic retinopathy, age-related macular degeneration, retinopathy of prematurity, and other irreversible causes of blindness involve pathological angiogenesis. The capillaries of the retina are unique, the inner layer of the blood-retinal barrier is like that of other capillaries, and is composed of a single layer of ECs. However, the outer layer of the blood-retinal barrier is formed by retinal pigment epithelial cells instead of pericytes and SMCs. Pathological ocular angiogenesis is triggered by hypoxia from neuronal metabolism, inflammatory signals, and oxidative stress. Those micro-environmental conditions cause retinal pigmented epithelium, astrocytes, Müller cells, ECs, ganglion cells to secrete VEGFA (Siemerink et al., 2010). According to our model, the Tip ECs that secrete VEGFA during pathological ocular angiogenesis are likely exposed the extracellular micro-environments in groups 34–35 in Figure 4B, and are affected by oxidative stress, lack of shear stress and have sufficient oxygen. The other Tip ECs involved in pathological ocular angiogenesis and induced by paracrine VEGFA correspond to groups 17–30 in Figure 4B.

Other angiogenic pathologies are caused by mutations that affect how an EC responds to changes in the extracellular micro-environment. We used our simplified model to simulate the effect of all single gain- and loss-of-function mutations on EC behavior. Specifically, we analyzed how each mutation affects the groups of extracellular micro-environments that cause Tip, Stalk, and Phalanx EC behaviors in our simplified model. The effect of some of the mutations has been observed experimentally and it should be possible to simulate the observed behavior using our model. The expected effect of reducing, or enlarging the number of extracellular micro-environments that cause each EC behavior depends on the likelihood of appearance of each micro-environment. Only when almost all or none of the micro-environments lead to a certain EC behavior, and the mutation has been observed in–vitro or in–vivo it is possible to compare the simulated effect of a certain mutation (Supplementary Table 13) with its experimentally observed effect.

Simulated loss of autocrine function of DLL4, ETS, MEK, or NRP1, leads to the loss of functional Tip EC behavior, strongly favoring Stalk EC behavior. Importantly, all four mutations have been observed to cause severe vascular defects in vivo and in vitro (Supplementary Tables 6, 10, and 12). The loss of autocrine DLL4 leads to the formation of a higher number of Tip ECs that do not inhibit their neighbor ECs from becoming Tip ECs (del Toro et al., 2010).

Simulated gain-of-function mutations for proteolytically active VEGFA, VEGFC, and VEGFD as well as NRP1, prevent Stalk EC behavior and cause more than 99% of the extracellular micro-environments to induce Tip EC behavior. In vivo and in vitro, proteolytically active VEGFA, VEGFC, and VEGFD increase blood vessel branching, angiogenesis, and permeability (Supplementary Tables 11, 12).

Simulations indicate that the Phalanx EC behavior is prevented by a loss of AKT, PIP3, or ShearStress function, or alternatively by constitutive ALK1, βcatenin, BMP10, BMP9, IGF, autocrine JAG, NICD, NOTCH, NRP1, SMAD1, TGFβ1, proteolytically active VEGFA, VEGFC, or VEGFD, WNT5a, or WNT7a activity. In vitro and in vivo, loss of AKT, PIP3, or ShearStress leads to mural cell loss, blood vessel destabilization and regression (Supplementary Tables 2, 3). Constitutive βcatenin, IGF, NOTCH, NRP1, SMAD1, proteolytically active VEGFA, VEGFC, or VEGFD, WNT5a, or WNT7a activity induces EC migration, proliferation, survival, or angiogenesis (Supplementary Tables 4, 8–12).

3.4. Robustness Analysis

Molecular regulatory networks must balance the need to ignore noise perturbations with the need to respond adequately to stimuli. A Boolean network can be classified as ordered, critical, or chaotic. Ordered Boolean networks resist most perturbations without any important changes in their dynamic behavior and are not sufficiently sensitive to stimuli. Chaotic Boolean networks tend to magnify perturbations and do not resist enough noise. Critical Boolean networks are selectively sensitive to certain perturbations and are sufficiently resilient to noise to be adequate models of molecular regulatory networks (Lloyd-Price et al., 2012). Additionally, the robustness of each trait has specific implications.

3.4.1. The Robustness of Tip, Stalk, and Phalanx EC Behavior to Single Gain and Loss-of-Function Mutations

The resilience of a functional phenotype to changes in the genotype allows the accumulation of genetic variation in a population, and needs to be achieved without limiting excessively the ability of a species to adapt by evolving different traits (Kirschner and Gerhart, 1998; Jiménez et al., 2015). The simulations showed that 23/128 = 17.96875% of all single gain- and loss-of-function mutations did not affect EC behavior at all. Furthermore, 82/128 = 64.0625% of mutations only cause changes in the response of an EC to certain extracellular micro-environments. The other 23/128 = 17.96875% of the mutations led to the loss of an EC behavior. Then, 4/128 = 3.125% of all mutations cause loss of Tip EC behavior. The same number of mutations cause Stalk EC behavior loss and strongly favor Tip EC behavior. Finally, 18/128 = 14.0625% of the mutations cause loss of Phalanx EC behavior. This set of results imply that our model of the network is robust to the complete loss of any of the main EC behaviors, however many mutations change the number of micro–environments that cause Tip, Stalk, and Phalanx EC behaviors (Supplementary Tables 13, 14).

3.4.2. The Robustness of Attractor Determination and EC Behavior to Molecular Activation Noise

Only 33.0538% of the trajectories followed by the perturbed copies of 1,000,000 random initial states reached the same attractor as the original state. In contrast, when we used 1,000,000 random initial states to test the robustness of EC behavior to molecular activation noise in 98.90088% of the relevant experiments the perturbation did not affect Tip EC, in 95.30536% of the relevant experiments, the perturbation did not affect Stalk EC behavior, and in 86.58824% of the relevant experiments, the perturbation did not affect Phalanx EC behavior. In general, 97.91060% of the random initial states reached the same EC behavior as the one reached by their perturbed copies.

3.4.3. The Sensitivity of Each Component of the Update Rule to Molecular Activation Noise

To understand which variables are more sensitive to stimuli and which ones tend to be more resilient to molecular activation noise. We estimated the sensitivity of each component of the update rule as described in the methods section. The results are shown in Figure 5 and the sensitivity values in section 2.1 in the Supplementary Material. The nodes with the six most sensitive update rules in our network are NRP1, MEK, Integrin, HEY1, SIRT1, AMPK. Even the update rule of NRP, the most sensitive in our model, has a relatively low sensitivity of 0.023716.

FIGURE 5
www.frontiersin.org

Figure 5. Update rule component sensitivity: (A) A darker shade of blue indicates a higher sensitivity in the update rule. Values range from VegfA = 0.002926 to NRP1 = 0.023716. (B) The sensitivities of the components of the update rule arranged from smallest to largest compared to the average sensitivity (0.01515947) which is shown as a red line.

4. Discussion

We presented in this work a reconstruction of the regulatory network involved in the control of angiogenesis, integrating the largest set of canonical signaling pathways to date. The dynamical behavior of the network, simulated as a Boolean network model, recovered the qualitative patterns of molecular activation observed in Phalanx, Tip, and Stalk ECs. Furthermore, the simulated behavior of the model corresponded to what has been reported in the literature regarding the high degree of behavioral plasticity between Phalanx, Stalk, and Tip EC behaviors in response to specific molecular micro-environments. Moreover, the model was also able to describe the effect of gain- and loss-of-function mutations.

4.1. Insights and Predictions Based on the Simulated Dynamic Behavior of Our Model

The qualitative agreement between our model and published data shows that the model is a useful framework to understand the mechanisms that underly normal angiogenesis. Furthermore, it allows generating hypotheses on the mechanisms by which a disruption in the system might lead to deviation in EC behavior, which might eventually lead to a pathogenic phenotype. The qualitative agreement between our model and published results cannot be attributed to some sort of model fitting. This is evidenced by the high robustness observed in the model against the complete loss of any of the main EC behaviors (Supplementary Tables 13, 14), despite the perturbations introduced in the update rules. Nonetheless, when we analyzed the effect of single gain- and loss-of-function mutations, the simulations recovered the observed effects of such mutations under certain micro-environments.

Our micro-environment EC behavior map allows us to put forward the following hypotheses about the requirements for Tip Stalk and Phalanx EC behaviors: (1) In a micro-environment with an active, paracrine VEGF ligand, the presence of either ShearStress, WNT5a, WNT7a, FGF, BMP9, BMP10, or TGFB1 is necessary to induce Tip EC behavior. (2) A micro-environment without VEGF can induce Tip EC behavior if it includes Oxygen, nutrients and IGF (Tip II, and Tip III in Table 1). However, the resulting Tip ECs secrete autocrine VEGFA. (3) DLL4 is not required for a micro-environment to induce Stalk EC behavior. (4) Shear stress and the absence of VEGF, TGF, IGF, WNT, and NOTCH ligands in the micro-environment is needed to observe a stable Phalanx EC behavior.

Based on the simulated effect of constitutive NRP1 activity, we predict that it prevents Stalk EC behavior and induces Tip EC behavior. We predict that constitutive ALK1, BMP9, BMP10, autocrine JAG, NICD, NRP1, SMAD1, and TGFβ1 activity inhibits Phalanx EC behavior based on the simulated effect of the corresponding gain-of-function mutations. Therefore, the model helps predict which mutations cause augmented mural cell loss, EC migration, proliferation, and angiogenesis, concomitant with inhibited Phalanx EC behavior.

Knowing the response of endothelial cells under a specific micro-environment is extremely relevant because inhibiting angiogenesis is an important medical goal during the treatment of vascular retinal disorders and cancer. Most of the drugs that are used to inhibit angiogenesis target the VEGF signaling pathway, inhibiting Tip EC behavior (Yadav et al., 2015). Our model suggests alternative ways to eliminate Tip EC behavior. Specifically, by eliminating the function of DLL4, ETS, MEK, or NRP1. Notably, both NRP1 and DLL4 are located on the cell membrane of ECs and are therefore easily reachable by drugs. Furthermore, in vascular retinal disorders, vascular permeability increases and vascular integrity diminishes, that is associated with intra-ocular hemorrhage and invasive potential of cancer. In principle, an extracellular micro-environment conducive to Phalanx EC behavior would help increase vascular integrity. Finally, stimulating angiogenesis is also an important medical goal during wound healing. It would be possible, thus, to use our model to explore one of the micro-environments that lead to Tip EC behavior and therefore, induce the wound healing process.

Arteriovenous malformations are very frequent in patients who suffer from Hereditary Hemorrhagic Telangiectasia (HHT), a disease associated with reduced ALK1, ENG, or SMAD4 function. In addition, Pulmonary Arterial Hypertension (PAH) is associated with reduced BMPRII or SMAD1 function. Furthermore, venous malformations have been observed in mice with constitutive TIE2 activity, as well as in mice with loss of ERK function. According to our model, the simulated effect of the mutations mentioned above includes an increase in the number of micro-environments that lead to Phalanx EC behavior, suggesting that the mentioned diseases are a consequence of ectopic blood vessel stabilization.

4.2. Assumptions and Limitations of Our Model

In this first version of the model of angiogenesis, we focus on the effect of the extracellular micro-environment on the behavior of a single endothelial cell. By using a Boolean model, we assume that all variables can only be active or inactive. Further, we use a synchronous update approach, therefore, we assume that all variables are activated or inhibited simultaneously. The limitations of our model affect the number of sprouting angiogenesis processes that we can reproduce and the extent to which we can simulate them. Some of the processes that are beyond the scope of our model have been studied using other previously published models (Peirce, 2008; Qutub et al., 2009; Scianna et al., 2013; Logsdon et al., 2014; Heck et al., 2015; Qutub and Popel, 2015) while other processes offer opportunities for further research as specified in the following paragraphs.

4.2.1. Secretion of Angiogenic Factors

According to our model, certain conditions cause ECs to secrete vascular growth factors (Figure 4B columns 31–35), the conditions that cause ECs to secrete active VEGFA (VEGFAxxxA) include sufficient oxygen, IGF, and a low AMP to ATP ratio. Normally, ECs are in contact with blood preventing hypoxia and lack of nutrients. The cells that compose other tissues respond to hypoxia or a high AMP to ATP ratio by secreting angiogenic factors; however, those cells are not included in our model. Additionally, Oxygen and then the secreted VEGF form concentration gradients. A continuous model that includes the geometry of the region or organ of interest as a boundary condition is necessary to simulate the gradient. Moreover, VEGFR1s secretion modulates the VEGFA concentration gradient (Chappell et al., 2016).

4.2.2. Vessel Destabilization

ANG2 activity is associated with mural cell detachment and it is possible to reproduce EC behavior during blood vessel destabilization using our model. However, it is not possible to reproduce pericyte and smooth muscle cell detachment because they are not included in our model. Some previous modeling efforts have included blood vessel destabilization (Zheng et al., 2013). However, in our opinion, mural cell behavior during angiogenesis merits a more detailed exploration.

4.2.3. Tip and Stalk Cell Differentiation

We carefully analyzed tip and stalk EC differentiation using our model emphasizing the interaction between the VEGF, WNT, TGF, NOTCH, Calcium, and NO signaling pathways during Tip and Stalk behavior specification. It is noteworthy that while Tip cells induce Stalk behavior in their neighbors by expressing DLL4 (Blanco and Gerhardt, 2013), according to our model NOTCH signaling inhibits Tip EC behavior only in a small group of micro-environments (Figure 4B, columns 34 and 35). A possible explanation for this apparent discrepancy is that active NOTCH signaling induces the secretion of VEGFR1s, which binds VEGFA, effectively raising the extracellular concentration of VEGFA needed to induce Tip EC behavior in the cells with active NOTCH signaling. In our Boolean model, it is not possible to include the changing VEGFAxxxP threshold, this would require a continuous model. Further, at the multicellular level, the chronological order in which ECs are affected by VEGFA and DLL4-mediated lateral inhibition creates a race condition (Bentley and Chakravartula, 2017). The temporal modulation of Tip and stalk EC behavior, including the effect of filipodia on tip cell sensitivity to VEGF, has been explored by previous modeling efforts (Venkatraman et al., 2016). A continuous, asynchronous, multicellular model that includes Matrix metalloproteinase, Apelin signaling (Palm et al., 2016) and VEGFR1s secretion (Chappell et al., 2016) would offer additional valuable insights.

4.2.4. Sprout Elongation

We simulated the micro-environmental conditions that may cause ECs to divide. However, our model does not include cell shape, which also changes during sprout elongation. Further sprout elongation is a multicellular process and our model includes only one EC. Several previous modeling efforts have studied sprout elongation (Logsdon et al., 2014). The authors of Norton and Popel (2016) analyzed the effect of EC proliferation, elongation, and migration during sprout elongation. Mechanical forces regulate both the location of sprout initiation and the rate of sprout elongation (Ghaffari et al., 2015), included in the model proposed by the authors of Vavourakis et al. (2017). A multi-scale model including cytoskeletal dynamics, molecular activation, and mechanical forces would greatly enhance our understanding of sprout elongation.

4.2.5. Lumen Formation and Expansion

PIP3, FAK, and SRC activity has been associated with vacuole secretion that is one of the main processes involved in lumen formation. According to the simulated dynamic behavior of our model, all Phalanx cells secrete vacuoles, additionally, type III Stalk ECs may also secrete vacuoles. Lumen formation is a multicellular process, that involves vacuole secretion and cytoskeletal remodeling. Simulating lumen formation, EC repulsion and flow-mediated lumen formation is beyond the scope of our current model. The authors of Boas and Merks (2014) focused their modeling efforts on the study of lumen formation.

4.2.6. Anastomosis

Is a multicellular process that involves cytoskeletal remodeling including specific shape changes that are beyond the scope of our model. Anastomosis has been included in several 2D and 3D models (Zheng et al., 2013; Norton and Popel, 2016). ECs with a reduced concentration of membrane-localized VEGFR1 are more likely to form stable connections with incoming sprouts (Nesmith et al., 2017). A multicellular model that integrates VEGFR1 regulation, and how it affects anastomosis, may help explain micro–vascular architecture.

4.2.7. Vessel Stabilization

Phalanx EC behavior is expected in stable blood vessels and is recovered by our model. PDGFB-mediated mural cell recruitment is also recovered by our model. Other multicellular effects of vessel stabilization, such as decreased blood vessel permeability, are beyond the scope of our model. Some previous modeling efforts have included blood vessel stabilization (Zheng et al., 2013). However, in our opinion, mural cell behavior during angiogenesis merits a more detailed exploration.

4.2.8. Pruning

Some of the micro-environments that cause atypical EC behavior without VEGF, FGF, IGF, and without Shear Stress (Figure 4C, group 60) may correspond to EC behavior during pruning. However, pruning involves changes in EC shape, EC fusion events, and EC migration, which have not been included in our model. Pruning is mainly regulated by blood flow. Apoptosis is implicated in the regression of large diameter blood vessels. In the small-diameter blood vessels that are remodeled by angiogenesis, pruning involves EC migration, self-fusion, and contraction before reabsorption into the remaining vasculature (Korn and Augustin, 2015; Betz et al., 2016). The model proposed by the authors of Chen et al. (2012) provided valuable insights into the role of hemodynamics during Zebrafish midbrain vascular pruning.

In conclusion, we developed a Boolean model of the network involved in EC behavior control during angiogenesis. The simulated dynamic behavior of our model corresponds with what has been observed experimentally and published about EC behavior and the effect of single gain- and loss-of-function mutations. The dynamical behavior of the model can qualitatively describe a wide variety of physiopathological states during angiogenesis. We believe that this characteristic makes the model a good platform to study the effect of altering the micro-environments and/or molecular backgrounds on endothelial cells.

Author Contributions

NW, LM, IG, and JK planned the research, wrote the article, analyzed, and discussed the results. NW reviewed the literature, composed the model, wrote the update rules, wrote the required scripts, and made the tables and figures. NW and LM carried out the simulations. IG and JK obtained funding for this project.

Funding

This work was partially supported by ABACUS, CONACyT grant EDOMEX-2011-C01-165873.

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 Elisa Domínguez Hüttinger, Marcos Nahmad Bensusan, Víctor Manuel Dávila Borja, Eugenio Azpeitia, Stalin Muñoz, David Rosenblueth, Elena Álvarez Buylla, and Rosa Angélica Castillo Rodríguez for their invaluable advice.

Supplementary Material

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

Abbreviations

AMP, Adenosine Monophosphate; AMPK, AMP-activated Protein Kinase; ANG, Angiopoietin; APLN, Apelin; ASF, Alternative Splicing Factor; ATP, Adenosine Triphosphate; BMP, Bone Morphogenetic Protein; BTrCP, Beta-Transducin Repeat Containing E3 Ubiquitin Protein Ligase; CXCR4, C-X-C motif chemokine Receptor 4; DAAM, Dishevelled Associated Activator Of Morphogenesis; DLL, Delta-Like canonical notch Ligand; DSH, Dishevelled; EC, Endothelial Cell; ECM, Extracellular Matrix; eNOS, Endothelial Nitric Oxide Synthase; EPC, Endothelial Progenitor Cells; EPH, Ephedrin; ERK, Extracellular signal-Regulated Kinase; FAK, Focal Adhesion Kinase; FOXO1, Forkhead Box O1; FGF, Fibroblast Growth Factor; FZD, Frizzled; GSK3β, Glycogen Synthase Kinase 3 Beta; HSC, Hematopoietic Stem Cell; HEY, Hes related family BHLH transcription factor with YRPW motif; HIF, Hypoxia-Inducible Factor; IA, Intussusceptive Angiogenesis; IGF, Insulin-like Growth Factor; JAG, Jagged; KLF, Kruppel Like Factor; LRP, LDL Receptor Related Protein; LEF, Lymphoid Enhancer binding Factor; MAPK, Mitogen-Activated Protein Kinase; MEK, Mitogen-Activated Protein Kinase Kinase (MAPKK); MMP, Matrix Metalloproteinase; mTOR, mechanistic Target Of Rapamycin; NFAT, Nuclear Factor of Activated T-cells; NICD, NOTCH Intracellular Domain; NO, Nitric Oxide; Nox2, NADPH oxidase 2; NRARP, NOTCH Regulated Ankyrin Repeat Protein; NRP1, Neuropilin 1; PA, Plasminogen Activator; PDGF, Platelet-Derived Growth Factor; PECAM, Platelet and Endothelial Cell Adhesion Molecule; PI3K, Phosphatidylinositol-4,5-bisphosphate 3-Kinase; PIP3, Phosphatidylinositol (3,4,5)-trisphosphate; PKC, Protein Kinase C; PTEN, Phosphatase and Tensin homolog; RHO, Rhodopsin; SA, Splitting Angiogenesis; SC, Stalk Cell; SF2, pre-mRNA-Splicing Factor 2; SIRT, Sirtuin, S1P, sphingosine-1-Phosphate; S1PR, sphingosine-1-Phosphate Receptor; TC, Tip Cell; TGF, Transforming Growth Factor; TIE, Tyrosine kinase with domains similar to Immunoglobulin and Epidermal growth factor; TSC, Tuberous Sclerosis; uPAR, Urokinase Receptor; VEGF, Vascular Endothelial Growth Factor; VEGFR, Vascular Endothelial Growth Factor-Receptor.

References

Azpeitia, E., Muñoz, S., González-Tokman, D., Martínez-Sánchez, M. E., Weinstein, N., Naldi, A., et al. (2017). The combination of the functionalities of feedback circuits is determinant for the attractors number and size in pathway-like boolean networks. Sci. Rep. 7:42023. doi: 10.1038/srep42023

CrossRef Full Text | Google Scholar

Bauer, A. L., Jackson, T. L., Jiang, Y., and Rohlf, T. (2010). Receptor cross-talk in angiogenesis: mapping environmental cues to cell phenotype using a stochastic, Boolean signaling network model. J. Theor. Biol. 264, 838–846. doi: 10.1016/j.jtbi.2010.03.025

PubMed Abstract | CrossRef Full Text | Google Scholar

Bazmara, H., Soltani, M., Sefidgar, M., Bazargan, M., Naeenian, M. M., and Rahmim, A. (2015). The vital role of blood flow-induced proliferation and migration in capillary network formation in a multiscale model of angiogenesis. PLoS ONE 10:e0128878. doi: 10.1371/journal.pone.0128878

PubMed Abstract | CrossRef Full Text | Google Scholar

Bazmara, H., Soltani, M., Sefidgar, M., Bazargan, M., Naeenian, M. M., and Rahmim, A. (2016). Blood flow and endothelial cell phenotype regulation during sprouting angiogenesis. Med. Biol. Eng. Comput. 54, 547–558. doi: 10.1007/s11517-015-1341-4

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Betz, C., Lenard, A., Belting, H.-G., and Affolter, M. (2016). Cell behaviors and dynamics during angiogenesis. Development 143, 2249–2260. doi: 10.1242/dev.135616

PubMed Abstract | CrossRef Full Text | Google Scholar

Biere, A. (2008). Picosat essentials. J. Satisfiabil. Boolean Model. Comput. 4, 75–97.

Google Scholar

Blancas, A. A., Wong, L. E., Glaser, D. E., and McCloskey, K. E. (2012). Specialized tip/stalk-like and phalanx-like endothelial cells from embryonic stem cells. Stem Cells Dev. 22, 1398–1407. doi: 10.1089/scd.2012.0376

PubMed Abstract | CrossRef Full Text | Google Scholar

Blanco, R., and Gerhardt, H. (2013). VEGF and Notch in tip and stalk cell selection. Cold Spring Harb. Perspect. Med. 3:a006569. doi: 10.1101/cshperspect.a006569

PubMed Abstract | CrossRef Full Text | Google Scholar

Boas, S. E., and Merks, R. M. (2014). Synergy of cell–cell repulsion and vacuolation in a computational model of lumen formation. J. R. Soc. Interface 11:20131049. doi: 10.1098/rsif.2013.1049

PubMed Abstract | CrossRef Full Text | Google Scholar

Bookholt, F., Monsuur, H., Gibbs, S., and Vermolen, F. (2016). Mathematical modelling of angiogenesis using continuous cell-based models. Biomech. Model. Mechanobiol. 15, 1577–1600. doi: 10.1007/s10237-016-0784-3

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, Q., Jiang, L., Li, C., Hu, D., Bu, J.-W., Cai, D., et al. (2012). Haemodynamics-driven developmental pruning of brain vasculature in zebrafish. PLoS Biol. 10:e1001374. doi: 10.1371/journal.pbio.1001374

PubMed Abstract | CrossRef Full Text | Google Scholar

Chillo, O., Kleinert, E. C., Lautz, T., Lasch, M., Pagel, J.-I., Heun, Y., et al. (2016). Perivascular mast cells govern shear stress-induced arteriogenesis by orchestrating leukocyte function. Cell Rep. 16, 2197–2207. doi: 10.1016/j.celrep.2016.07.040

PubMed Abstract | CrossRef Full Text | Google Scholar

Czirok, A. (2013). Endothelial cell motility, coordination and pattern formation during vasculogenesis. Wiley Interdiscip. Rev. Syst. Biol. Med. 5, 587–602. doi: 10.1002/wsbm.1233

PubMed Abstract | CrossRef Full Text | Google Scholar

Deindl, E., and Schaper, W. (2005). The art of arteriogenesis. Cell Biochem. Biophys. 43, 1–15. doi: 10.1385/CBB:43:1:001

PubMed Abstract | CrossRef Full Text | Google Scholar

del Toro, R., Prahst, C., Mathivet, T., Siegfried, G., Kaminker, J. S., Larrivee, B., et al. (2010). Identification and functional analysis of endothelial tip cell–enriched genes. Blood 116, 4025–4033. doi: 10.1182/blood-2010-02-270819

PubMed Abstract | CrossRef Full Text | Google Scholar

Dubrova, E., and Teslenko, M. (2011). A SAT-based algorithm for finding attractors in synchronous boolean networks. IEEE/ACM Trans. Comput. Biol. Bioinform. 8, 1393–1399. doi: 10.1109/TCBB.2010.20

PubMed Abstract | CrossRef Full Text | Google Scholar

Félix, M.-A., and Barkoulas, M. (2015). Pervasive robustness in biological systems. Nat. Rev. Genet. 16, 483–496. doi: 10.1038/nrg3949

PubMed Abstract | CrossRef Full Text | Google Scholar

Forsythe, J. A., Jiang, B.-H., Iyer, N. V., Agani, F., Leung, S. W., Koos, R. D., et al. (1996). Activation of vascular endothelial growth factor gene transcription by hypoxia-inducible factor 1. Mol. Cell. Biol. 16, 4604–4613.

PubMed Abstract | Google Scholar

Garg, A., Di Cara, A., Xenarios, I., Mendoza, L., and De Micheli, G. (2008). Synchronous versus asynchronous modeling of gene regulatory networks. Bioinformatics 24, 1917–1925. doi: 10.1093/bioinformatics/btn336

PubMed Abstract | CrossRef Full Text | Google Scholar

Ghaffari, S., Leask, R. L., and Jones, E. A. (2015). Flow dynamics control the location of sprouting and direct elongation during developmental angiogenesis. Development 142, 4151–4157. doi: 10.1242/dev.128058

PubMed Abstract | CrossRef Full Text | Google Scholar

Gianni-Barrera, R., Trani, M., Fontanellaz, C., Heberer, M., Djonov, V., Hlushchuk, R., et al. (2013). VEGF over-expression in skeletal muscle induces angiogenesis by intussusception rather than sprouting. Angiogenesis 16, 123–136. doi: 10.1007/s10456-012-9304-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Gianni-Barrera, R., Trani, M., Reginato, S., and Banfi, A. (2011). To sprout or to split? VEGF, Notch and vascular morphogenesis. Biochem. Soc. Trans. 39, 1644–1648. doi: 10.1042/BST20110650

PubMed Abstract | CrossRef Full Text | Google Scholar

Glaser, D. E., Turner, W. S., Madfis, N., Wong, L., Zamora, J., White, N., et al. (2016). Multifactorial optimizations for directing endothelial fate from stem cells. PLoS ONE 11:e0166663. doi: 10.1371/journal.pone.0166663

PubMed Abstract | CrossRef Full Text | Google Scholar

Glass, D. S., Jin, X., and Riedel-Kruse, I. H. (2016). Signaling delays preclude defects in lateral inhibition patterning. Phys. Rev. Lett. 116:128102. doi: 10.1103/PhysRevLett.116.128102

PubMed Abstract | CrossRef Full Text | Google Scholar

Gödde, R., and Kurz, H. (2001). Structural and biophysical simulation of angiogenesis and vascular remodeling. Dev. Dyn. 220, 387–401. doi: 10.1002/dvdy.1118

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

Heil, M., Eitenmüller, I., Schmitz-Rixen, T., and Schaper, W. (2006). Arteriogenesis versus angiogenesis: similarities and differences. J. Cell. Mol. Med. 10, 45–55. doi: 10.1111/j.1582-4934.2006.tb00290.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Ji, J. W., Tsoukias, N. M., Goldman, D., and Popel, A. S. (2006). A computational model of oxygen transport in skeletal muscle for sprouting and splitting modes of angiogenesis. J. Theor. Biol. 241, 94–108. doi: 10.1016/j.jtbi.2005.11.019

PubMed Abstract | CrossRef Full Text | Google Scholar

Jiménez, A., Cotterell, J., Munteanu, A., and Sharpe, J. (2015). Dynamics of gene circuits shapes evolvability. Proc. Natl. Acad. Sci. U.S.A. 112, 2103–2108. doi: 10.1073/pnas.1411065112

CrossRef Full Text | Google Scholar

Kässmeyer, S., Plendl, J., Custodis, P., and Bahramsoltani, M. (2009). New insights in vascular development: vasculogenesis and endothelial progenitor cells. Anat. Histol. Embryol. 38, 1–11. doi: 10.1111/j.1439-0264.2008.00894.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Kerr, B. A., West, X. Z., Kim, Y.-W., Zhao, Y., Tischenko, M., Cull, R. M., et al. (2016). Stability and function of adult vasculature is sustained by Akt/Jagged1 signalling axis in endothelium. Nat. Commun. 7:10960. doi: 10.1038/ncomms10960

PubMed Abstract | CrossRef Full Text | Google Scholar

Kirschner, M., and Gerhart, J. (1998). Evolvability. Proc. Natl. Acad. Sci. U.S.A. 95, 8420–8427.

PubMed Abstract | Google Scholar

Korn, C., and Augustin, H. G. (2015). Mechanisms of vessel pruning and regression. Dev. Cell 34, 5–17. doi: 10.1016/j.devcel.2015.06.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Kumar, V. S., Binu, S., Soumya, S., Haritha, K., and Sudhakaran, P. (2014). Regulation of vascular endothelial growth factor by metabolic context of the cell. Glycoconj. J. 31, 427–434. doi: 10.1007/s10719-014-9547-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Kurz, H., Burri, P. H., and Djonov, V. G. (2003). Angiogenesis and vascular remodeling by intussusception: from form to function. Physiology 18, 65–70. doi: 10.1152/nips.01417.2002

PubMed Abstract | CrossRef Full Text | Google Scholar

Lloyd-Price, J., Gupta, A., and Ribeiro, A. S. (2012). Robustness and information propagation in attractors of random boolean networks. PLoS ONE 7:e42018. doi: 10.1371/journal.pone.0042018

PubMed Abstract | CrossRef Full Text | Google Scholar

Logsdon, E. A., Finley, S. D., Popel, A. S., and Mac Gabhann, 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

Makanya, A. N., Hlushchuk, R., and Djonov, V. G. (2009). Intussusceptive angiogenesis and its role in vascular morphogenesis, patterning, and remodeling. Angiogenesis 12, 113–123. doi: 10.1007/s10456-009-9129-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Müssel, C., Hopfensitz, M., and Kestler, H. A. (2010). Boolnet–an r package for generation, reconstruction and analysis of boolean networks. Bioinformatics 26, 1378–1380. doi: 10.1093/bioinformatics/btq124

PubMed Abstract | CrossRef Full Text | Google Scholar

Naldi, A., Remy, E., Thieffry, D., and Chaouiya, C. (2011). Dynamically consistent reduction of logical regulatory graphs. Theor. Comput. Sci. 412, 2207–2218. doi: 10.1016/j.tcs.2010.10.021

CrossRef Full Text | Google Scholar

Nesmith, J. E., Chappell, J. C., Cluceru, J. G., and Bautch, V. L. (2017). Blood vessel anastomosis is spatially regulated by flt1 during angiogenesis. Development 144, 889–896. doi: 10.1242/dev.145672

PubMed Abstract | CrossRef Full Text | Google Scholar

Norton, K.-A., and Popel, A. S. (2016). Effects of endothelial cell proliferation and migration rates in a computational model of sprouting angiogenesis. Sci. Rep. 6:36992. doi: 10.1038/srep36992

PubMed Abstract | CrossRef Full Text | Google Scholar

Nussenbaum, F., and Herman, I. M. (2010). Tumor angiogenesis: insights and innovations. J. Oncol. 2010, 1–24. doi: 10.1155/2010/132641

PubMed Abstract | CrossRef Full Text | Google Scholar

Palm, M. M., Dallinga, M. G., van Dijk, E., Klaassen, I., Schlingemann, R. O., and Merks, R. M. (2016). Computational screening of tip and stalk cell behavior proposes a role for apelin signaling in sprout progression. PLoS ONE 11:e0159478. doi: 10.1371/journal.pone.0159478

PubMed Abstract | CrossRef Full Text | Google Scholar

Patan, S., Haenni, B., and Burri, P. H. (1996). Implementation of intussusceptive microvascular growth in the chicken chorioallantoic membrane (CAM):: 1. pillar formation by folding of the capillary wall. Microvasc. Res. 51, 80–98.

PubMed Abstract | Google Scholar

Patan, S., Haenni, B., and Burri, P. H. (1997). Implementation of intussusceptive microvascular growth in the chicken chorioallantoic membrane (CAM): 2. pillar formation by capillary fusion. Microvasc. Res. 53, 33–52.

Google Scholar

Peirce, S. M. (2008). Computational and mathematical modeling of angiogenesis. Microcirculation 15, 739–751. doi: 10.1080/10739680802220331

PubMed Abstract | CrossRef Full Text | Google Scholar

Qin, D., Trenkwalder, T., Lee, S., Chillo, O., Deindl, E., Kupatt, C., et al. (2013). Early vessel destabilization mediated by Angiopoietin-2 and subsequent vessel maturation via Angiopoietin-1 induce functional neovasculature after ischemia. PLoS ONE 8:e61831. doi: 10.1371/journal.pone.0061831

PubMed Abstract | CrossRef Full Text | Google Scholar

Qutub, A. A., MacGabhann, F., Karagiannis, E. D., Vempati, P., and Popel, A. S. (2009). Multiscale models of angiogenesis. IEEE Eng. Med. Biol. Mag. 28, 14–31. doi: 10.1109/MEMB.2009.931791

PubMed Abstract | CrossRef Full Text | Google Scholar

Qutub, A. A., and Popel, A. S. (2015). “Angiogenesis, computational modeling perspective,” in Encyclopedia of Applied and Computational Mathematics, ed B. Engquist (Berlin; Heidelberg: Springer Berlin Heidelberg), 58–67.

Google Scholar

Saadatpour, A., Albert, I., and Albert, R. (2010). Attractor analysis of asynchronous boolean models of signal transduction networks. J. Theor. Biol. 266, 641–656. doi: 10.1016/j.jtbi.2010.07.022

PubMed Abstract | CrossRef Full Text | Google Scholar

Saadatpour, A., Albert, R., and Reluga, T. C. (2013). A reduction method for Boolean network models proven to conserve attractors. SIAM J. Appl. Dyn. Syst. 12, 1997–2011. doi: 10.1137/13090537X

CrossRef Full Text | Google Scholar

Scharpfenecker, M., Fiedler, U., Reiss, Y., and Augustin, H. G. (2005). The Tie-2 ligand angiopoietin-2 destabilizes quiescent endothelium through an internal autocrine loop mechanism. J. Cell Sci. 118, 771–780. doi: 10.1242/jcs.0165

PubMed Abstract | CrossRef Full Text | Google Scholar

Scianna, M., Bell, C., and Preziosi, L. (2013). A review of mathematical models for the formation of vascular networks. J. Theor. Biol. 333, 174–209. doi: 10.1016/j.jtbi.2013.04.037

PubMed Abstract | CrossRef Full Text | Google Scholar

Siemerink, M. J., Augustin, A. J., and Schlingemann, R. O. (2010). Mechanisms of ocular angiogenesis and its molecular mediators. Dev. Ophthalmol. 46, 4–20. doi: 10.1159/000320006

PubMed Abstract | CrossRef Full Text | Google Scholar

Song, J. W., and Munn, L. L. (2011). Fluid forces control endothelial sprouting. Proc. Natl. Acad. Sci. U.S.A. 108, 15342–15347. doi: 10.1073/pnas.1105316108

PubMed Abstract | CrossRef Full Text | Google Scholar

Vavourakis, V., Wijeratne, P. A., Shipley, R., Loizidou, M., Stylianopoulos, T., and Hawkes, D. J. (2017). A validated multiscale in-silico model for mechano-sensitive tumour angiogenesis and growth. PLoS Comput. Biol. 13:e1005259. doi: 10.1371/journal.pcbi.1005259

PubMed Abstract | CrossRef Full Text | Google Scholar

Veliz-Cuba, A. (2011). Reduction of boolean network models. J. Theor. Biol. 289, 167–172. doi: 10.1016/j.jtbi.2011.08.042

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:e0166489. doi: 10.1371/journal.pone.0166489

PubMed Abstract | CrossRef Full Text | Google Scholar

Vermolen, F., and Javierre, E. (2012). A finite-element model for healing of cutaneous wounds combining contraction, angiogenesis and closure. J. Mathe. Biol. 65, 967–996. doi: 10.1007/s00285-011-0487-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Weinstein, N., Ortiz-Gutiérrez, E., Muñoz, S., Rosenblueth, D. A., Álvarez-Buylla, E. R., and Mendoza, L. (2015). A model of the regulatory network involved in the control of the cell cycle and cell differentiation in the Caenorhabditis elegans vulva. BMC Bioinformatics 16:1. doi: 10.1186/s12859-015-0498-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Yadav, L., Puri, N., Rastogi, V., Satpute, P., and Sharma, V. (2015). Tumour angiogenesis and angiogenic inhibitors: A review. J. Clin. Diagn. Res. 9, XE01–XE05. doi: 10.7860/JCDR/2015/12016.6135

PubMed Abstract | CrossRef Full Text | Google Scholar

Zheng, X., Koh, G. Y., and Jackson, T. (2013). A continuous model of angiogenesis: Initiation, extension, and maturation of new blood vessels modulated by vascular endothelial growth factor, angiopoietins, platelet-derived growth factor-b, and pericytes. Discrete Continuous Dyn. Syst. Ser. B 18, 1109–1154. doi: 10.3934/dcdsb.2013.18.1109

CrossRef Full Text | Google Scholar

Keywords: sprouting angiogenesis, network model, mechanical stress, cell differentiation, cell polarization, lateral inhibition

Citation: Weinstein N, Mendoza L, Gitler I and Klapp J (2017) A Network Model to Explore the Effect of the Micro-environment on Endothelial Cell Behavior during Angiogenesis. Front. Physiol. 8:960. doi: 10.3389/fphys.2017.00960

Received: 07 September 2017; Accepted: 10 November 2017;
Published: 27 November 2017.

Edited by:

Matteo Barberis, University of Amsterdam, Netherlands

Reviewed by:

Laurence Calzone, Institut Curie, France
Aleksander S. Popel, Johns Hopkins University, United States

Copyright © 2017 Weinstein, Mendoza, Gitler and Klapp. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Nathan Weinstein, nathan.weinstein4@gmail.com
Jaime Klapp, jaime.klapp@inin.gob.mx