University of Southern California Dissertations and Theses

Mechanistic investigation of pro-angiogenic signaling and endothelial sprouting mediated by FGF and VEGF

(USC Thesis Other)

Mechanistic investigation of pro-angiogenic signaling and endothelial sprouting mediated by FGF and VEGF

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content Mechanistic Investigation of Pro-Angiogenic Signaling and Endothelial Sprouting
Mediated by FGF and VEGF

by
Min Song

A Dissertation Presented to the
FACULTY OF THE USC GRADUATE SCHOOL
UNIVERSITY OF SOUTHERN CALIFORNIA
In Partial Fulfillment of the
Requirements for the Degree
DOCTOR OF PHILOSOPHY
BIOMEDICAL ENGINEERING

August 2021

Copyright 2021 Min Song

ii
Dedication
This dissertation is dedicated to my parents, Xiaoe Gao and Liewu Song, for their endless
support, encouragement, and love.

iii
Acknowledgments
First, I would like to express my sincere appreciation to my advisor, Professor Stacey D.
Finley. I am very grateful that she helped me to grow in many different ways, both professionally
and personally. I started my PhD program with a minimal computational background. She patiently
helped me to get started with understanding and reading code and pointed me towards relevant
academic papers to read. Her constant encouragement and support have helped me through the
ups and downs of these five years. Her assistance at every stage of my research, from incubating
ideas to shaping methods and forming results, guided and encouraged me in my academic research.
This work would not have been possible without her persistent help. Her knowledge and insight
have also been an inspiration to me as well, and her guidance prepared me to be confident and
ready for my future. I am grateful for her support, encouragement, and guidance throughout my
journey in graduate school. I would like to say thank you for pushing me to achieve my very best
both while at USC and beyond.
Next, I would like to extend my sincere gratitude to my committee members, Professors
David Z. D’Argenio, Scott E. Fraser, Paul K. Newton, and Keyue Shen, who provided me with
invaluable advice and supported me to grow as a researcher. Their insightful comments and
suggestions from different perspectives were vital inspirations while forming my dissertation.
I would also like to thank my lab members and alumni, Dr. Mahua Roy, Dr. Jennifer A.
Rohrs, Dr. Qianhui Wu, Ding Li, Sahak Makaryan, Colin G. Cess, Patrick Gelbach, Holly Huber,
Vardges Tserunyan, Niki Tavakoli, Ryland Morlock, and Ariella Simoni for all the time we spent
together in the lab and their great support and assistance in my academic studies. I am thankful to
have worked with such great people and will cherish our time together.

iv
I am also thankful to the Department of Biomedical Engineering, Viterbi School of
Engineering, and all its member staff for their considerate guidance, especially William Yang, who
has been very patient and supportive during my PhD program.
I am very grateful to my friends who have supported me all along: Cheryl Pan, Dr.
Caicheng Chris Lu, Dr. Wenbo Chen, Xiwei She, Kai Wang, Hao Zhou, Nannan Ding, Sisi Tian,
Sida Tao, Li Jiang, Hengzhao Bian, Lindsay Dillon, Xinting Han, Wenliang Wang, Yi Deng, Yijie
Zhou. Thank you for always being there for me and making my studies and life at USC a wonderful
time.
I would especially like to thank my parents, Xiaoe Gao and Liewu Song, for their support
and great love. Without their tremendous understanding and encouragement, it would have been
impossible for me to become the person who I am today. I love and appreciate you more than
words can say. You both mean the world to me.

v
Table of Contents
Dedication ...................................................................................................................................... ii
Acknowledgments ........................................................................................................................ iii
List of Figures ............................................................................................................................. viii
List of Tables ................................................................................................................................ ix
Abstract .......................................................................................................................................... x
Chapter 1: Introduction ............................................................................................................... 1
1.1. Angiogenesis: target for tissue engineering and cancer treatment ............................... 1
1.2. Understanding the interactions of pro-angiogenic factors in mediated signaling
pathways quantitatively is critical for improving current therapeutic strategies. ............. 2
1.3. Computational modeling provides a powerful tool to investigate angiogenesis
processes systematically............................................................................................................ 3
1.4. Dissertation Outline ........................................................................................................... 4
Chapter 2: FGF- and VEGF-induced MAPK signaling ............................................................ 8
2.1. Abstract ............................................................................................................................... 8
2.2. Introduction ........................................................................................................................ 9
2.3. Methods ............................................................................................................................. 12
2.3.1. Model construction ..................................................................................................... 12
2.3.2. Sensitivity Analysis .................................................................................................... 15
2.3.3. Data extraction ............................................................................................................ 16
2.3.4. Model parameters........................................................................................................ 16
2.3.5. ERK phosphorylation response................................................................................... 19
2.4. Results ............................................................................................................................... 21
2.4.1 The validated mathematical model captures the main features of FGF- and VEGF-
stimulated ERK phosphorylation dynamics.......................................................................... 21
2.4.2. FGF produces a greater angiogenic response than VEGF when considering the
maximum ERK phosphorylation produced .......................................................................... 26
2.4.3. The combination of FGF and VEGF has greater effects in inducing maximum pERK
than the summation of the individual effects ........................................................................ 29
2.4.4. The combination of FGF and VEGF shows a fast and sustained pERK response ..... 31
2.4.5. Increasing VEGFR2 density can compensate for the relatively low ERK
phosphorylation induced by VEGF ...................................................................................... 34
2.4.6. ERK phosphorylation induced by VEGF can be promoted by decreasing VEGFR2
internalization and degradation rates .................................................................................... 37
2.5. Discussion.......................................................................................................................... 39
2.6. Conclusions ....................................................................................................................... 45
2.7. Acknowledgements .......................................................................................................... 45
2.8. Additional files ................................................................................................................. 45
Chapter 3: FGF- and VEGF-induced MAPK and PI3K/Akt signaling................................. 46
3.1. Abstract ............................................................................................................................. 46
3.2. Introduction ...................................................................................................................... 47
3.3. Methods ............................................................................................................................. 51
3.3.1 Model construction ...................................................................................................... 51

vi
3.3.2 Sensitivity analysis....................................................................................................... 53
3.3.3. Data extraction ............................................................................................................ 54
3.3.4. Model fitting and validation........................................................................................ 54
3.3.5. Monte Carlo simulations ............................................................................................. 56
3.4. Results ............................................................................................................................... 58
3.4.1 The fitted and validated molecular-detailed mathematical model captures the major
characteristics of FGF- and VEGF-induced ERK and Akt phosphorylation dynamics ....... 58
3.4.2. FGF produces greater maximum pAkt and pERK than VEGF at equimolar
concentrations ....................................................................................................................... 62
3.4.3. Akt activation shows a stronger response than ERK in terms of magnitude .............. 64
3.4.4. ERK activation is more responsive to changing the ligand concentration compared to
Akt......................................................................................................................................... 66
3.4.5. The co-stimulation by FGF and VEGF has a greater impact on phosphorylation of ERK
compared to summation of the ligands’ individual effects ................................................... 70
3.4.6. The model identifies potential targets for influencing ERK and Akt activation and
evaluates their efficacy ......................................................................................................... 73
3.5. Discussion.......................................................................................................................... 78
3.6. Conclusion ........................................................................................................................ 83
3.7. Acknowledgements .......................................................................................................... 83
3.8. Additional files ................................................................................................................. 83
Chapter 4: FGF- and VEGF-induced endothelial sprouting .................................................. 84
4.1. Abstract ............................................................................................................................. 84
4.2. Introduction ...................................................................................................................... 86
4.3. Methods ............................................................................................................................. 90
4.3.1. Model construction ..................................................................................................... 90
4.3.2. Sensitivity analysis...................................................................................................... 98
4.3.3. Data extraction ............................................................................................................ 99
4.3.4. Parameterization ......................................................................................................... 99
4.4. Results ............................................................................................................................. 102
4.4.1. The fitted hybrid agent-based model captures the main features of FGF- and VEGF-
induced endothelial sprouting characteristics ..................................................................... 102
4.4.2. The type and concentration of ligand, length of growth factor stimulation, and initial
number of cells impact endothelial sprouting ..................................................................... 105
4.4.3. The cell proliferation and sprout growth of existing sprouts are predicted to be more
important in the sprouting process compared to the effect of rate of forming a new sprout
............................................................................................................................................. 110
4.4.4. The MAPK and PI3K pathways contribute to cell proliferation, sprout growth, and
probability of sprouting in different ways .......................................................................... 112
4.4.5. ERK pathway regulates vessel network mainly via regulating cell proliferation and NS,
while Akt pathway mainly affects vessel network via regulating sprout growth ............... 115
4.5. Discussion........................................................................................................................ 119
4.6. Conclusion ...................................................................................................................... 126
4.7. Acknowledgements ........................................................................................................ 126
4.8. Supplementary Materials .............................................................................................. 126

vii
Chapter 5: Conclusion .............................................................................................................. 136
5.1. Overview ......................................................................................................................... 136
5.2. Summary ......................................................................................................................... 136
5.3. Future Directions ........................................................................................................... 138
5.4. Conclusion ...................................................................................................................... 140
References .................................................................................................................................. 141

viii
List of Figures
Figure 2-1. Schematic of FGF and VEGF signaling network. ..................................................... 14
Figure 2-2. Model comparison to training data for FGF or VEGF stimulation. ........................... 24
Figure 2-3. Model comparison to validation data. ........................................................................ 25
Figure 2-4. Predicted maximum pERK response. ........................................................................ 29
Figure 2-5. Predicted time response of pERK following stimulation by FGF, VEGF, and their
combination................................................................................................................................... 33
Figure 2-6. Predicted pERK response with varied initial VEGFR2 concentrations. .................... 37
Figure 2-7. Effect of varying VEGFR2 trafficking parameters on pERK response. .................... 39
Figure 3-1. Schematic of FGF and VEGF signaling network. ..................................................... 52
Figure 3-2. Model comparison to training data for FGF or VEGF stimulation. ........................... 60
Figure 3-3. Model comparison to validation data. ........................................................................ 61
Figure 3-4. Predicted pERK and pAkt responses stimulated by single agents. ............................ 63
Figure 3-5. Predicted maximum pERK and pAkt responses with co-stimulation. ....................... 66
Figure 3-6. Comparison of mono- and co-stimulation.................................................................. 72
Figure 3-7. Predicted targets for modulating pERK and pAkt responses. .................................... 78
Figure 4-1. Endothelial spheroid sprouting process. .................................................................... 90
Figure 4-2. Flowchart of the endothelial sprouting agent-based model. ...................................... 93
Figure 4-3. Model comparison to training and validation data for FGF or VEGF stimulation. . 105
Figure 4-4. Predicted sprouting responses stimulated by single agents...................................... 108
Figure 4-5. Predicted sprouting responses in response to FGF and VEGF co-stimulation. ....... 110
Figure 4-6. Predicted rcp, rsg, and p in response to mono-and co-stimulation of FGF and VEGF.
..................................................................................................................................................... 112
Figure 4-7. The contributions of MAPK and PI3K/Akt pathways to rcp and rsg in response to FGF
and VEGF mono-stimulation. ..................................................................................................... 114
Figure 4-8. Predicted representative targets for modulating rcp and rsg. ..................................... 117
Figure 4-9. Predicted effects of ERK in modulating TL (μm), NS, and AL (μm). ..................... 118
Figure 4-10. Predicted effects of Akt in modulating TL (μm), NS, and AL (μm). ..................... 119

ix
List of Tables
Table 3-1. Fold change for pAkt and pMEK in response to varying ligand concentration. ......... 70
Table 3-2. The total sensitivity index Sti values. .......................................................................... 75
Table 4-1. Representative low, intermediate, high levels of FGF and VEGF. ............................. 99
Table 4-2. Influential parameters. ............................................................................................... 116

x
Abstract
Angiogenesis, the formation of new blood capillaries from pre-existing blood vessels, plays
an important role in the survival of tissues, because vessels carry blood throughout the body to
provide oxygen and nutrients required by the resident cells. Angiogenesis has received great
attention, as angiogenic-based strategies are beneficial in many contexts, including in both tissue
engineering and cancer treatment. Specifically, promoting new blood vessel formation to ensure
adequate vasculature is critical for the survival of artificial tissues; on the other hand, inhibiting
angiogenesis is an important strategy for cancer treatment, as vessels not only deliver nutrients to
support cancer growth, they also provide routes for metastasis. However, the outcomes of these
pro- and anti-angiogenic therapies are not all effective. Lack of therapeutic response, low efficacy,
or drug resistance have been observed. One major obstacle for angiogenesis-based therapies is that
the angiogenic response is driven by signals integrated from multiple relevant signaling pathways.
Targeting one pathway could be insufficient, as alternative pathways may compensate,
diminishing the overall effect of the treatment strategy. For example, there are multiple proteins
that promote signaling needed for angiogenesis. Fibroblast growth factor (FGF) and vascular
endothelial growth factor (VEGF), which are two major pro-angiogenic factors, mediate mitogen-
activated protein kinase (MAPK) and phosphatidylinositol 3-kinase/protein kinase B (PI3K/Akt)
pathways. Signaling through these pathways leads to activated extracellular regulated kinase
(ERK) and Akt, which are important in angiogenic responses such as endothelial cell survival,
proliferation, and migration. However, there is a limited understanding of how these promoters
combine together to stimulate angiogenesis. Therefore, understanding FGF- and VEGF-mediated
angiogenic signaling systematically can be beneficial in improving current angiogenic strategies.

xi
We seek to achieve optimal angiogenic outcomes and inform the development of effective
angiogenic therapies by constructing experimentally validated mathematical models. In this work,
I present a series of mechanistic models to address questions regarding the interactions of FGF-
and VEGF-mediated angiogenic signaling in endothelial cells, the effects of the crosstalk in
angiogenic cellular responses mediated by these factors, and potential targets for angiogenic-based
therapies.

1
1. Chapter 1
Chapter 1
Introduction

1.1. Angiogenesis: target for tissue engineering and cancer treatment
Angiogenesis, the formation of new blood capillaries from pre-existing blood vessels, plays
an important role in the survival of tissues, as the vessels carry blood throughout the body to
provide oxygen and nutrients required by the resident cells. Angiogenesis also provides a route for
tumor metastasis. Thus, targeting angiogenesis has been a strategy in many contexts, including in
both tissue engineering and cancer treatment. Specifically, in the context of tissue engineering,
researchers have sought to create artificial tissues to substitute damaged tissues in response to a
great shortage of donors for transplant surgery. One challenge for the success of synthetic tissues
is to ensure sufficient transport of nutrients and gases such as oxygen to the cells. Thus, stimulating
new blood vessel formation is an important strategy for the long-term viability of engineered tissue
constructs. On the other hand, inhibiting angiogenesis is an important strategy for cancer treatment
as tumors grow by obtaining nutrients and oxygen from the blood delivered by vessels, and tumor
metastasis is facilitated by the blood circulation.
The processes of blood vessel formation, particularly for capillaries, are mostly initiated
and mediated by endothelial cells responding to the local physiological conditions (1). Many
different pro-angiogenic growth factors, such as fibroblast growth factor (FGF), vascular
endothelial growth factor (VEGF), and platelet-derived growth factor (PDGF), regulate
angiogenesis (2, 3). These factors promote different cellular processes involving endothelial cells

2
leading to new blood vessel formation, including cell survival, proliferation, migration, and vessel
maturation (4, 5). Strategies to promote or inhibit angiogenesis focus on modulating the effects of
these factors to alter the cellular-level processes they induce, with a focus on endothelial cells.
However, not all approaches to promote or inhibit angiogenesis lead to successful
outcomes. For example, clinical trials have shown no effective improvement in FGF- (6) or VEGF-
induced (7) angiogenesis. Also, bevacizumab, an anti-VEGF agent, has limited effects in certain
cancer types, and it is no longer approved for the treatment of metastatic breast cancer due to its
disappointing results (8). Thus, there is a need to better understand the mechanisms of
angiogenesis, specifically the molecular interactions and signaling required for new blood vessel
formation and how they affect cellular behaviors, in order to establish more effective therapeutic
strategies.

1.2. Understanding the interactions of pro-angiogenic factors in mediated signaling pathways
quantitatively is critical for improving current therapeutic strategies.
Many growth factors regulate angiogenic signaling, especially FGF and VEGF. In addition,
there is crosstalk between intracellular signaling pathways. The overall response in endothelial
angiogenesis is dependent on the integrated signals activated by these growth factor-mediated
pathways to influence cellular decisions, such as proliferation, survival, and migration, and further
promote or inhibit vessel formation. However, there is a limited understanding of the integrated
effects of more than one factor on intracellular signaling reactions and further cellular behaviors
at a detailed level. In pro-angiogenic strategies, some synergistic effects between FGF and VEGF
in endothelial angiogenic activities have been shown (9-11). Also, in the case of inhibiting
angiogenesis, tumors often evade the effects of drugs that target a single factor by making use of

3
alternate compensatory pathways to activate signaling species needed for proliferation and
migration. For instance, FGFR activation may play a role in the resistance mechanism of anti-
angiogenic drugs, especially anti-VEGF treatment (12, 13). Additionally, experiments show high
levels of FGFR1 in tumors that continue to progress, even during anti-VEGF therapy (14). Thus,
it is needed to better understand the mechanism of how multiple angiogenic factors act together
quantitatively on regulating molecular and cellular behaviors.
In this study, we are interested in the angiogenic signals required to initiate vessel growth.
Therefore, we focused on endothelial cell molecular signaling and further cellular responses
promoted by FGF and VEGF, as they are particularly important in early stages of angiogenesis,
while PDGF is more important in maturing the vessels. FGF and VEGF bind to their receptors and
initiate signaling through the mitogen-activated protein kinase (MAPK) and phosphatidylinositol
3-kinase/protein kinase B (PI3K/Akt) pathways to phosphorylate extracellular regulated kinase
(ERK) and Akt, respectively. ERK and Akt are important signaling species in the angiogenesis
process that influence cell proliferation (15, 16), survival (17-23), and migration (21, 24-26).
Thus, we quantitatively investigated the combination effects of two major pro-angiogenic factors,
FGF and VEGF, on activating MAPK and PI3K/Akt signaling on a molecular level and further
early stages of angiogenic cellular responses, specifically endothelial sprouting.

1.3. Computational modeling provides a powerful tool to investigate angiogenesis processes
systematically.
Given the complexity of biochemical reactions comprising angiogenesis signaling
networks, a better understanding of the dynamics of these networks quantitatively is beneficial for
current angiogenesis-based strategies in many contexts. Mathematical modeling serves as a

4
powerful tool to investigate molecular and cellular responses systematically and to guide
experimental design. There are many published models that predict molecular (27-29) and cellular
(30, 31) responses mediated by angiogenic factors. However, such models are mostly designed to
predict responses upon single agent stimulation. Targeting more than one growth factor and
exploring their effects in intracellular signaling and cellular responses in detail deserves more
attention. In addition, many models that focus on specific cellular behaviors significantly reduced
the intracellular signaling network such that the output signal is simply linearly proportional to the
fraction of bound receptors (30-32). Therefore, in this work I constructed a series of computational
models to characterize the intracellular signaling interactions of FGF and VEGF leading to ERK
and Akt activation and utilize downstream molecular signals, maximum pERK and pAkt, as inputs
to describe angiogenic cellular responses, specifically endothelial sprouting.

1.4. Dissertation Outline
In this work, I present a series of computational models that quantitatively predict the
molecular signaling dynamics and cellular angiogenic responses mediated by the activation of
FGF- and VEGF-induced MAPK and PI3K/Akt pathways. This model framework can be used to
mechanistically explain experimental results, guide experimental design, and aid the development
of pro- and anti-angiogenic strategies.
In Chapter 2, I develop a mechanistic model for MAPK signaling pathway in endothelial
cells mediated by the stimulation of FGF and VEGF, which is referred to as the ERK model in this
work. There are models that study FGF- (27, 33) or VEGF- induced signaling (29, 34) alone; this
is the first model that studies their interactions on a molecular level. I first construct a molecular-
detailed biochemical reaction network of FGF- and VEGF-induced MAPK pathway (35). This

5
model is built on previous computational work (33, 36). Specifically, I adapted the FGF signaling
network from the model by Kanodia et al. (33) and simplified the model of VEGF-induced ERK
phosphorylation pathways from Tan and coworkers (36). Importantly, I expanded upon these
previous models to capture the major steps of FGF- and VEGF-induced ERK phosphorylation.
The main model includes 70 reactions, 72 species, and 75 parameters. The reactions, initial
conditions, and parameter values are taken from published literature (33, 36, 37) and listed in
Additional file 2-2: Tables S2 to S4. The detail of model construction is provided in Chapter 2. I
train and validate the mathematical model to experimental data for FGF- and VEGF-induced
species activation in the network. I then apply the model to predict the dynamics of phosphorylated
ERK (pERK) in response to the stimulation by FGF and VEGF individually and in combination.
The model predicts that FGF and VEGF have differential effects on pERK. Additionally, the
model predictions show that VEGFR2 density and trafficking parameters significantly influence
the response to VEGF stimulation. Overall, the model provides mechanistic insight into FGF and
VEGF interactions on ERK phosphorylation in endothelial cells.
In Chapter 3, I expand my ERK model to include PI3K/Akt signaling activated via FRS2
and VEGFR2 by FGF and VEGF, respectively, which is referred to as the ERK-Akt model in this
work. The PI3K/Akt signaling module is adapted from previous modeling work (38). We used the
best fit values of the parameters and initial values for ERK activation to experimental data in our
previous model (35) and took the remaining model parameters and initial values from published
literature (38) as baseline values. We provided the detail of model construction in Chapter 3. I first
train and validate the mathematical model to experimental data for FGF- and VEGF-induced
species activation in the pathway. This model includes 97 reactions, 99 species, and 100
parameters (see Additional file 3-3). The reactions, initial conditions and parameter values are

6
listed in Additional file 3-1: Tables S1-S3. The validated mathematical model is applied to
characterize the dynamics of pERK and pAkt in response to the mono- and co-stimulation by FGF
and VEGF. The model predicts that ERK and Akt activation have distinct responses. Also, the
combination of FGF and VEGF indicates a greater effect in increasing the maximum pERK
compared to the summation of individual effects, which is not seen for maximum pAkt levels. In
addition, our model identifies the influential species and kinetic parameters that specifically
modulate the pERK and pAkt responses, which represent potential targets for angiogenesis-based
therapies. Overall, the model predicts the combination effects of FGF and VEGF stimulation on
ERK and Akt quantitatively and provides a framework to mechanistically explain experimental
results and guide experimental design.
In Chapter 4, I develop a hybrid agent-based model to describe endothelial sprouting driven
by the integrated molecular signals, pERK and pAkt, upon the mono- and co-stimulation of FGF
and VEGF. The molecular interactions, the initial concentrations, and parameters of FGF- and
VEGF-induced MAPK and PI3K/Akt signaling in endothelial cells are adapted from our previous
ERK-Akt model (39). The baseline values of all remaining model parameters are adapted from
literature (31, 40-43) and the details are provided in Supplementary Materials and Table S4-1. To
link the molecular signals with the short-term cellular responses, we assumed that the endothelial
cell proliferation, sprout growth, and probability of sprouting are dependent on the maximum pAkt
and pERK levels upon the stimulation of FGF and VEGF within two hours, following Hill
functions. The detail of model construction is elaborated in Chapter 4. The model is trained and
validated experimental data that describe endothelial proliferation and sprouting. The validated
model is then applied to predict endothelial sprouting in response the mono- and co-stimulation of
FGF and VEGF. The model predicts that FGF plays a dominant role in the combination effects in

7
endothelial sprouting. In addition, the model suggests that ERK and Akt pathways and cellular
responses contribute differently to the overall sprouting process. Last, the model predicts that the
strategies to modulate endothelial sprouting are context dependent. Our model can identify
potential effective pro- and anti-angiogenic targets under different conditions and study their
efficacy. The model provides detailed mechanistic insight into VEGF and FGF interactions in
sprouting angiogenesis.
Altogether, my work focuses on using a systems biology approach to study FGF- and
VEGF-induced intracellular signaling and cellular responses in endothelial cells. The models
predict the combination effects of FGF and VEGF stimulation quantitatively and provides a
framework to synthesize experimental data. More broadly, this modeling can be utilized to identify
targets that influence angiogenic signaling leading to endothelial sprouting and to study the effects
of pro- and anti-angiogenic therapies. These models can be expanded to investigate the effects of
other pro- and anti-angiogenic signaling and other cellular responses.

8
2. Chapter 2
Chapter 2
FGF- and VEGF-induced MPAK signaling

Portions of this chapter are adapted from:
Min Song and Stacey D. Finley. BMC Systems Biology (2018)

2.1. Abstract
Angiogenesis is important in physiological and pathological conditions, as blood vessels
provide nutrients and oxygen needed for tissue growth and survival. Therefore, targeting
angiogenesis is a prominent strategy in both tissue engineering and cancer treatment. However,
not all of the approaches to promote or inhibit angiogenesis lead to successful outcomes.
Angiogenesis-based therapies primarily target pro-angiogenic factors such as vascular endothelial
growth factor-A (VEGF) or fibroblast growth factor (FGF) in isolation. However, pre-clinical and
clinical evidence shows these therapies often have limited effects. To improve therapeutic
strategies, including targeting FGF and VEGF in combination, we need a quantitative
understanding of the how the promoters combine to stimulate angiogenesis.
In this study, we trained and validated a detailed mathematical model to quantitatively
characterize the crosstalk of FGF and VEGF intracellular signaling. This signaling is initiated by
FGF binding to the FGF receptor 1 (FGFR1) and heparan sulfate glycosaminoglycans (HSGAGs)
or VEGF binding to VEGF receptor 2 (VEGFR2) to promote downstream signaling. The model

9
focuses on FGF- and VEGF-induced mitogen-activated protein kinase (MAPK) signaling and
phosphorylation of extracellular regulated kinase (ERK), which promotes cell proliferation. We
apply the model to predict the dynamics of phosphorylated ERK (pERK) in response to the
stimulation by FGF and VEGF individually and in combination. The model predicts that FGF and
VEGF have differential effects on pERK. Additionally, since VEGFR2 upregulation has been
observed in pathological conditions, we apply the model to investigate the effects of VEGFR2
density and trafficking parameters. The model predictions show that these parameters significantly
influence the response to VEGF stimulation.
The model agrees with experimental data and is a framework to synthesize and
quantitatively explain experimental studies. Ultimately, the model provides mechanistic insight
into FGF and VEGF interactions needed to identify potential targets for pro- or anti-angiogenic
therapies.

2.2. Introduction
Angiogenesis is the formation of new blood capillaries from pre-existing blood vessels.
The essential role of blood vessels in delivering nutrients makes angiogenesis important in the
survival of tissues, including tumor growth. Angiogenesis also provides a route for tumor
metastasis. Thus, targeting angiogenesis is a prominent strategy in many contexts, for example, in
both tissue engineering and cancer treatment.
In the context of tissue engineering, there is a large demand for organs needed for transplant
surgery, but a great shortage of donors. The long-term viability of engineered tissue constructs
depends on growth of new vessels from host tissue, and stimulating new blood vessel formation is
an important pro-angiogenic strategy for tissue engineering (44). Alternatively, the formation of

10
new blood vessels is important for cancer growth and metastasis. Thus, inhibiting angiogenesis is
an anti-angiogenic strategy for cancer treatment. Unfortunately, not all approaches to promote or
inhibit angiogenesis lead to successful outcomes. For example, clinical trials have shown no
effective improvement in blood flow or perfusion by fibroblast growth factor (FGF)-induced (6)
or vascular endothelial growth factor-A (VEGF)-induced (7) angiogenesis. Specifically, a double-
blinded randomized controlled trial studied recombinant FGF-induced angiogenesis and showed
no symptomatic improvement (exercise tolerance or myocardial perfusion) following 90 or 180
days of treatment (6). Similarly, in a double-blinded placebo-controlled trial to study the effects of
recombinant human VEGF-induced angiogenesis in animal models, there was no improvement in
angina, in comparison with placebo by day 60. Only a high dose of VEGF (50 ng/kg/min) showed
any effect (7). Also, bevacizumab, an anti-VEGF product for cancer treatment, has limited effects
in certain cancer types, and it is no longer approved for the treatment of metastatic breast cancer
due to its disappointing results (8). Thus, there is a need to better understand the molecular
interactions and signaling required for new blood vessel formation, in order to establish more
effective therapeutic strategies.
The established angiogenesis-based therapies primarily target pro-angiogenic factors such
as FGF and VEGF in isolation. However, both FGF and VEGF bind to their receptors to initiate
mitogen-activated protein kinase (MAPK) signaling and phosphorylate ERK, the final output of
the MAPK pathway (36, 45). This signaling pathway promotes cell proliferation in the early stages
of angiogenesis. Additionally, the combined effects of FGF and VEGF have been shown to be
greater than their individual effects (9, 11). A quantitative understanding of how these promoters
combine together to stimulate angiogenesis could greatly benefit the current pro- and anti-
angiogenic therapies.

11
Mathematical modeling is a useful tool to predict the molecular response mediated by
angiogenic factors. For example, Mac Gabhann and Popel studied interactions between VEGF
isoforms, VEGF receptors (VEGFR1, VEGFR2, NRP1), and the extracellular matrix using a
molecular-detailed model. The model predicted that blocking Neuropilin-VEGFR coupling is
more effective in reducing VEGF-VEGFR2 signaling than blocking Neuropilin-1 expression or
binding of VEGF to Neuropilin-1 (27). Stefanini et al. constructed a pharmacokinetic model that
studied VEGF distribution after intravenous administration of bevacizumab, and they found that
plasma VEGF was increased after treatment (28). Filion and Popel explored myocardial deposition
and retention of FGF after intracoronary administration of FGF using a computational model. The
model predicted that the response time is dependent on the reaction time of the binding of FGF to
FGFR rather than the FGF diffusion time. Receptor secretion and internalization have also been
predicted to be crucial in FGF dynamics (29). Wu and Finley characterized the intracellular
signaling of TSP1-induced apoptosis and predicted response of cell population to TSP1-mediated
apoptosis by mathematical modeling (34). Zheng et al. integrated the effects of VEGF,
angiopoietins (Ang1 and Ang2) and platelet-derived growth factor-B (PDGF-B) on endothelial
proliferation, migration, and maturation using mathematical modeling. Their model illustrated that
competition between Ang1 and Ang 2 acts as an angiogenic switch and that combining anti-
pericyte and anti-VEGF therapy is more effective than anti-VEGF therapy alone in inducing blood
vessel regression (46). Such models are useful to predict molecular responses of VEGF or FGF
stimulation; however, surprisingly, the interactions between these two growth factors have not
been investigated in detail. Targeting multiple growth factors simultaneously, exploiting their
overlapping and redundant signaling pathways, may improve angiogenesis-based therapies. Thus,

12
there is a need for a model that provides quantitative insights into combination effects of FGF and
VEGF.
In the present study, we aim to quantitatively characterize the crosstalk between FGF and
VEGF in MAPK signaling leading to phosphorylated ERK (pERK). We focus on pERK because
pERK promotes cell proliferation (15) and is mostly found in active, rather than quiescent,
endothelial cells (47). We constructed a computational model that incorporates the molecular
interactions between FGF, VEGF, and their receptors, leading to MAPK signaling. We apply the
model to explore how FGF and VEGF promote ERK phosphorylation. This is the first model that
studies FGF and VEGF interactions together on a molecular level. Our model predicts the
combination effects of FGF and VEGF stimulation and shows that FGF plays a dominant role in
promoting ERK phosphorylation. Using this model, we also investigated the effects of the VEGF
receptor VEGFR2, including how VEGFR2 density and trafficking parameters influence the ERK
response. By predicting the effect of VEGFR2 density and trafficking parameters, we can get a
better understanding of the role of VEGFR2 under pathological conditions. Additionally,
understanding with quantitative detail the FGF and VEGF interactions helps identify potential
targets for enhancing pro- or anti-angiogenic therapies.

2.3. Methods
2.3.1. Model construction
We constructed a molecular-detailed biochemical reaction network including FGF, VEGF,
and their receptors FGFR1 and VEGFR2 (Figure 2-1). Signaling is induced by the growth factors
binding to their receptors, culminating with phosphorylation of ERK, through the MAPK cascade.

13
MAPK signaling is initiated through the activation of Raf and FRS2 by VEGF and FGF,
respectively. Activated Raf (aRaf) and phosphorylated FRS2 (pFRS2) phosphorylate MEK at two
sites, and doubly phosphorylated MEK (ppMEK) further phosphorylates ERK. The molecular
interactions involved in the network are illustrated in Figure 2-1. The model is a novel
advancement of published computational models (33, 36). Specifically, we adapted the
competition of FGF and HSGAG to the binding of FGFR1 and the feedback loop from pERK to
FRS2 from the model by Kanodia et al., and we expanded the model by including FGFR trafficking
(internalization, recycling, and degradation) and accounting for both singly- and doubly-
phosphorylated MEK and ERK. It is worth noting that the formation of the tertiary signaling
complex of FGF:HSGAG:FGFR only occurs by FGFR binding to the FGF:HSGAG complex. This
is because the affinity of FGF:HSGAG binding is approximately two times stronger than that of
FGF:FGFR, and FGFR and HSGAG binding is more than two orders of magnitude lower (48). In
addition, we simplified the model of VEGF-induced ERK phosphorylation pathways from Tan
and coworkers; specifically, we only include Ras activation either from Shc-independent or Shc-
dependent pathways (36). Thus, we expanded upon previous models to capture the major steps of
FGF- and VEGF-induced ERK phosphorylation and better understand their interactions.

14

Figure 2-1. Schematic of FGF and VEGF signaling network. Signaling is induced by the
growth factors binding to their receptors, culminating with phosphorylation of ERK, through the
MAPK cascade. MAPK signaling is initiated through the activation of Raf and FRS2 by VEGF
and FGF, respectively. The FGF:HSGAG:FGFR1 complex dimerizes and leads to phosphorylation
of FRS2 (pFRS2). VEGF binds VEGFR2 to activate Ras, forming Ras-GTP, which further
activates Raf (aRaf). Both aRaf and pFRS2 are able to phosphorylate MEK at two sites, and doubly
phosphorylated MEK (ppMEK) further phosphorylates ERK. ppERK provides negative feedback
on the FGF pathway, as it promotes ubiquitination of FRS2 (FRS2u).

The model is simulated using a range of concentrations for FGF and VEGF, based on
published experimental studies. Typically, FGF and VEGF concentrations are usually within the
range of 0 - 50 ng/ml (0 – 2.2 nM) (49-51) and 0 – 100 ng/ml (0 – 2.2 nM) (52-56) respectively in
in vitro studies, although some studies utilized high concentrations such 300 ng/ml VEGF (9), 500
ng/ml FGF (33). And it has been reported that 25 ng/ml (0.56 nM) and 50 ng/ml (1.1 nM) VEGF
significantly increase HUVECs tube formation, and 0.1 ng/ml (0.004 nM) and 1.0 ng/mL (0.04
nM) FGF strongly induced tube formation on Matrigel after 24 hours compared to the control
groups (50). Moreover, Pepper et al. showed that the total sprout length formed by bovine
microvascular endothelial cells started to plateau when treated with 30 ng/ml (1.3 nM) FGF and

15
100 ng/ml (2.2 nM) VEGF (9). To account for these findings, we simulated the model with the
concentration of FGF and VEGF ranging of 0.01 to 2 nM.
The network is implemented as an ordinary differential equation (ODE) model using
MATLAB (Mathworks, Inc.). The main model includes 70 reactions, 72 species, and 75
parameters. The reactions, initial conditions, and parameter values are listed in Additional file 2-
2: Tables S2 to S4. All reactions are assumed to follow the law of mass action. Receptor
internalization, recycling, and degradation are considered in the model. Because the simulated time
is within two hours, we do not consider the degradation of the ligands or signaling species. The
final model is available in Additional file 2-3. We also implement a modified model that includes
heparin to validate the estimated model parameters (described below).

2.3.2. Sensitivity Analysis
To identify the parameters and initial concentrations that significantly influence the model
outputs, we performed the extended Fourier Amplitude Sensitivity Test (eFAST) (57). All targeted
parameters and initial values were varied simultaneously within specified bounds (one order of
magnitude above and below the baseline values), and the effects of multiple model inputs (kinetic
parameters or initial conditions) on the pERK concentration were computed (the total sensitivity
indices, “S
ti
”). We studied the S
ti
values for kinetic parameters and initial concentrations. The S
ti
index can range from 0 to 1, where a higher S
ti
index indicates this input is more influential to the
output.

16
2.3.3. Data extraction
Data from published experimental studies (33, 52-54, 58) were used for parameter fitting
and model validation. The Western blot data was extracted using ImageJ. Experimental data from
plots was extracted using the function grabit.

2.3.4. Model parameters
The trafficking parameters for VEGFR2 and the parameters and initial values that are
involved in the overlap of FGF and VEGF signaling pathways were estimated by fitting the model
to experimental data using Particle Swarm Optimization (PSO) implemented by Iadevaia (59). We
used MATLAB to implement the PSO algorithm. A total of 39 parameters and initial values were
estimated in the fitting (Additional file 2-2: Table S1, and also highlighted in red in Additional file
2-2: Tables S3 and S4). All other parameters were taken from published literature (33, 36, 37).
The parameters characterizing the overlapping MAPK pathway were chosen for fitting because
while FGF and VEGF upstream parameters are well documented individually in literature, a
uniform set of parameters for their interactions is needed for this combined model. The VEGFR2
trafficking parameters were fitted, as they have been shown to significantly affect ERK activation
(36). Additionally, many of the kinetic parameters for the overlapping reactions and the VEGFR2
trafficking rates were shown to significantly influence pERK in the sensitivity analysis (Additional
file 2-1: Figure S1).
PSO starts with a population of initial particles (parameter sets). As the particles move
around (i.e., as the algorithm explores the parameter space), an objective function is evaluated at
each particle location. Particles communicate with one another to determine which has the lowest

17
objective function value. The objective function for each parameter set was used to identify
optimal parameter values. Specifically, we used PSO to minimize the weighted sum of squared
residuals (WSSR):
WSSR(θ)=𝑚𝑖𝑛 ∑(
𝑉 𝑝𝑟𝑒𝑑 ,𝑖 (θ)−𝑉 𝑒𝑥𝑝 ,𝑖 𝑉 𝑒𝑥𝑝 ,𝑖 )
2
𝑛 𝑖=1

where Vexp,i is the ith experimental measurement, Vpred,i is the ith predicted value at the
corresponding time point, and n is the total number of experimental data points. The minimization
is subject to θ, the set of upper and lower bounds on each of the free parameters. The bounds were
set to be one order of magnitude above and below the baseline parameter values, which were taken
from literature and listed in Additional file 2-2: Tables S3 to S4. Although PSO is a global
parameter estimation approach, and the parameter values are varied within each run to minimize
the error, we still ran the algorithm multiple times to attempt to identify the optimal parameter
values within the large search space. We were able to obtain a total of 72 fitted parameter sets,
which were ultimately narrowed down to 16 parameter sets that allowed the model to capture the
training and validation data sets.
The model was fitted against three datasets, specifically: 1) normalized pERK induced by
FGF concentrations varying from 0.16 ng/ml to 500 ng/ml, where pERK level was normalized by
the maximum pERK stimulated by FGF across all six concentrations (0.16, 0.8, 4, 20, 100, and
500 ng/ml) in two hours, experiments conducted using the H1703 cell line (33); 2) normalized
pVEGFR2 (pR2) stimulated by 5 ng/ml VEGF, pR2 was normalized by the maximum pR2 induced
by 5 ng/ml VEGF, experiments conducted using HUVECs (53); 3) normalized pERK induced by
50 ng/ml VEGF, where pERK was normalized by the maximum pERK induced by 50 ng/ml
VEGF, experiments conducted using HUVECs (52). We note that the pERK and pR2 in the model

18
simulation include all free and bound forms of pERK and ppERK, and all free and bound forms
of pR2 except the degraded pR2, respectively.
After model training, we validated the fitted model with three datasets not used in the
fitting. First, we simulated the effects following the addition of heparin. For this case, we added
500 μg/ml heparin, which competes with HSGAG and binds to FGF. There are an additional 26
reactions, 25 species, and 3 parameters for heparin perturbation in the model. Without any fitting,
parameters are all taken from Kanodia et al.. The influence of heparin is illustrated in Additional
file 2-1: Figure S3, and details are provided in Additional file 2-2: Tables S5 and S6. We simulated
the pERK dose response with or without heparin to compare with the experiments described by
Kanodia (33). Having a difference in pERK with and without additional heparin that is greater
than zero indicates that the presence of heparin enhances FGF-induced ERK phosphorylation.
Second, we predicted the phosphorylated pERK following stimulation with 10 ng/ml FGF,
mimicking measurements obtained from BAECs (58). Third, we predicted the VEGF-induced pR2
response upon stimulation with 80 ng/ml VEGF, simulating experiments conducted using
HUVECs (54). For all three datasets, we simulated the experimental conditions without any
additional model fitting and compared to the experimental measurements. A total of 16 parameter
sets with the smallest errors were taken to be the “best” sets based on the model fitting and
validation (Additional file 2-1: Figure S2 and Additional file 2-2: Table S1) and were used for
model simulations.
VEGFR2 density.
To study the impact of VEGFR2 expression on VEGF-induced angiogenesis, we varied
VEGFR2 density within ten-fold of the baseline value (1000 molecules/cell) and predicted the
level and dynamics of ERK phosphorylation.

19

VEGFR2 trafficking parameters.
To investigate the effects of VEGFR2 trafficking in VEGF-induced ERK phosphorylation,
we decreased the trafficking parameters (internalization, recycling, and degradation rates) values
for VEGFR2. We changed the parameters one-by-one or together to be the same level as FGFR
trafficking parameters and predicted the VEGF-induced pERK response. For example, we made
the internalization rate of free VEGFR2 to be the same as the rate at which free FGFR is
internalized.

2.3.5. ERK phosphorylation response
We investigated the ERK phosphorylation response by the stimulation of FGF or VEGF
alone, compared to their combination. In this study, we mainly focus on two aspects of pERK
dynamics: magnitude of the response and timescale of signaling.
2.3.5.1. Magnitude of ERK phosphorylation response
a. Maximum pERK. We calculate the maximum ERK phosphorylation level induced by the
stimulation of FGF, VEGF, or their combination.

b. Ratio, R. To compare the combination effects with FGF and VEGF individual effects, we
introduce the ratio below:
𝑅 =
𝑚𝑎𝑥 pERK(FGF and VEGF)
𝑚𝑎𝑥 pERK(FGF)+𝑚𝑎𝑥 pERK(VEGF)

20
When R is greater than one, it indicates that the combination effect in inducing maximal
pERK is greater than the summation of individual effects; when R is equal to one, it implies that
the combination effect is additive; when R is less to one, it suggests an antagonistic effect between
FGF and VEGF.
2.3.5.2. Timescale of the signaling response.
We use two parameters to characterize the timescale of ERK activation: the time to reach
the maximum pERK (T1) and the time duration that pERK level remains greater than half of its
maximal value (T2). T1 indicates how quickly ERK is phosphorylated: the smaller T1 is, the faster
ERK becomes phosphorylated. T2 indicates how long ERK remains in a phosphorylated state: the
larger T2 is, the more sustained the pERK response.
2.3.5.3. Reaction rates.
We specify the rates of each reaction based on the law of mass action, where the rate of a
chemical reaction is proportional to the amount of each reactant. For example, for the
phosphorylation of VEGFR2:
𝑉𝐸𝐺𝐹 +𝑉𝐸𝐺𝐹𝑅 2
𝑘 𝑝𝑅 2
, 𝑘 𝑑𝑝𝑅 2
↔ 𝑝𝑉𝐸𝐺𝐹𝑅 2
The reaction rate is:
𝑅𝑎𝑡𝑒 = 𝑘 𝑝𝑅 2
∙[𝑉𝐸𝐺𝐹 ]∙[𝑉𝐸𝐺𝐹𝑅 2]− 𝑘 𝑑𝑝𝑅 2
∙[𝑝𝑉𝐸𝐺𝐹𝑅 2]
where 𝑘 𝑝𝑅 2
and 𝑘 𝑑𝑝𝑅 2
are rate constants for the forward and reverse reactions, respectively, and
[VEGF], [VEGFR2], and [pVEGFR2] are the species’ concentrations.
We simplified the VEGFR phosphorylation into one step because it has been reported that
the two VEGFR2 monomers phosphorylate each other upon ligation (60). Also, the specific

21
mechanism of VEGFR2 phosphorylation is not our focus in this study. Therefore, we assume
autophosphorylation upon VEGF binding, as it has been included in other papers as well (36, 38).

2.4. Results
2.4.1 The validated mathematical model captures the main features of FGF- and VEGF-
stimulated ERK phosphorylation dynamics.
We constructed a computational model that characterizes FGF and VEGF interactions
leading to ERK phosphorylation (Figure 2-1). Signaling is mediated by FGF and VEGF binding
to their respective receptors, leading to pERK. The model was trained against published
experimental measurements (33, 52, 53) using Particle Swarm Optimization (PSO) (59) for
parameter estimation. We note that pERK response stimulated by FGF was measured using the
non-small cell lung cancer cell line NCI-H1730 (33), while phosphorylated VEGFR2 (pR2) (53)
and pERK (52) responses induced by VEGF were obtained using human umbilical vein endothelial
cells (HUVECs). In this study, we assume the downstream signaling has the same kinetics across
different cell lines. The FGFR1 and heparan sulfate glycosaminoglycan (HSGAG) levels on
various cell types are fairly consistent: approximately 10
4
to 10
5
molecules/cell for Balb/c3T3 (37),
osteoblasts and bone marrow stromal cells (61), Bovine Aortic Endothelial cells (BAECs) (62),
and NCI-H1730 (33). In addition, other models have made similar assumptions that signaling is
consistent across cell types, for example, using HUVEC data to study macrophage signaling
responses (41). That work, which supposes that the unique signaling responses for specific cell
types are due to their different protein levels rather than the kinetics, is able to match experimental
data. Thus, we assume the FGFR1 signaling pathway for NCI-H1730 is the same as in endothelial
cells.

22
To explore the behavior of the model, we performed a global sensitivity analysis, which
identifies the variables that significantly influence the model outputs. Specifically, we performed
the eFAST (see Methods for more details) and computed the total sensitivity indices (S
ti
) for initial
conditions or kinetic parameters on pERK level by the stimulation of FGF or VEGF. The eFAST
results reveal the specific species and kinetic parameters that affect pERK. These results show the
importance of particular species (Additional file 2-1: Figure S1A): FRS2 and Ptase2 for FGF-
induced signaling and Ras, MEK, and Ptase2 for VEGF-induced signaling. Additionally, the rates
of certain reactions involving these species are also shown to be important (Additional file 2-1:
Figure S1B). This includes “ked2” and “k_dpMEK_p”, which are the dephosphorylation rate of
ppMEK by Ptase2 and the association rate of pMEK and Ptase2, respectively. Moreover, the
eFAST results indicate the importance of the VEGFR2 trafficking parameters. Specifically, the
internalization rates of bound VEGFR (“k_intb”, “k_recb”, and “k_degb”) are shown to affect
ppERK with VEGF stimulation (Additional file 2-1: Figure S1B, bottom two panels). These results
provide the foundation to investigate how the VEGFR2 trafficking rates influence pERK
dynamics, which we explore below. In addition, these results can guide parameter fitting. Since
many of the influential parameters and initial values are involved in the overlap of FGF and VEGF
signaling pathways, along with the VEGFR2 trafficking parameters, we selected those values for
model fitting.
The fitted model shows a good match to the training data (Figure 2-2). It can capture the
biphasic pERK response caused by FGF stimulation, which has been reported by Kanodia et al.
(33) (Figure 2-2 A, C): pERK increases as the FGF concentration increases from low to
intermediate levels, and decreases with increasing FGF at high concentrations. The decrease in
pERK is caused by the competitive binding of FGF to HSGAG and FGFR (33). At lower FGF

23
levels, there are enough FGFR molecules to bind with FGF, forming FGF-FGFR complexes that
can interact with HSGAGs to form the pro-angiogenic tertiary complex (FGF-HSGAG-FGFR).
However, at higher FGF levels, formation of the FGF-FGFR complex is limited by the number of
FGFR molecules, which thus limits the formation of the tertiary complex. This is primarily due to
the different expression levels of FGFR and HSGAGs (2 × 10
4
versus 10
5
molecules/cell,
respectively) (33). Our fitted model recapitulates this biphasic response at the simulated time
points. Also, VEGF-induced upstream (pVEGFR2) and downstream (pERK) dynamics have good
agreement with experimental measurements (Figure 2-2B, D) (52, 53). For the best 16 fits, the
weighted errors range from 17.4 to 18.6 (Additional file 2-2: Table S1). In addition, the fitted
parameters have good consistency. We allowed the values of the initial conditions of kinetic
parameters to be fitted to vary up to two orders of magnitude (10-fold above and below the baseline
starting value). We examined the results of the 16 best parameter sets, focusing on the range
spanning two standard deviations above and below the mean for each fitted initial condition and
parameter. For all the fitted initial conditions and 32 of the 34 fitted parameters, this range was
within one order of magnitude (Additional file 2-1: Figure S2). Thus, the fitting reduced the
bounds.

24

Figure 2-2. Model comparison to training data for FGF or VEGF stimulation. A. Normalized
pERK dynamics in response to FGF concentrations ranging from 0.16 to 500 ng/ml. B. Normalized
VEGFR2 phosphorylation time course following stimulation with 5 ng/ml VEGF. C. Dose
response of pERK for FGF stimulation. D. Normalized ERK phosphorylation time course upon
stimulation with 50 ng/ml VEGF. The circles are experimental data. Curves are the mean values
of the 16 best fits. Shaded regions show standard deviation of the fits.

To validate the model, we compared the model predictions to additional experimental data.
We first applied heparin perturbation to the trained model to reproduce another set of data by
Kanodia et al. Heparin is a soluble source of HSGAGs, and it competes with HSGAGs to bind
with FGF, interfering with FGF-induced signaling (Additional file 2-1: Figure S3) (33). It has been
reported that additional heparin increases FGF-induced ERK phosphorylation at high FGF
concentrations, and decreases FGF-induced ERK phosphorylation at low FGF concentrations in
two hours (33). We validated the model by expanding it to include heparin binding and comparing
to the experimental measurements. We added 500 g/ml of heparin, the same concentration used

25
in experiments (33), and found that the difference of predicted pERK responses upon FGF
stimulation with and without heparin within two hours exhibits the trend observed to occur
experimentally: at high FGF concentrations, the difference between in the presence and absence
of heparin is greater than zero, while the difference is less than zero at lower FGF concentrations
(Figure 2-3A). This shows a qualitative agreement with experimental observations.

Figure 2-3. Model comparison to validation data. A. The differences in pERK induced by
stimulation of FGF with and without 500 μg/ml heparin at four time points are predicted. Each dot
represents one fit. The dots are spread horizontally to avoid overlap of similar responses from
different fits. B. Normalized pERK by the stimulation of 10 ng/ml FGF. Circles are BAEC
experimental data. C. Normalized VEGFR2 phosphorylation time course upon stimulation with
80 ng/ml VEGF. Circles are HUVEC experimental data. Curves are the mean values of the 16 best
fits from the model training. Shaded regions show standard deviation of the fits.

A separate set of experimental measurements of phosphorylated ERK by the stimulation
of 10 ng/ml FGF conducted using BAECs (58) was extracted to further validate FGF-induced

26
endothelial signaling. Our model also has a good agreement with this experimental data (Figure 2-
3B), which further confirms that this model can be used to predict endothelial cell signaling.
Finally, we extracted an independent set of experimental measurements of the
phosphorylated VEGFR2 response by the stimulation with 80 ng/ml VEGF in HUVECs (54) to
validate the VEGF-induced signaling. Our model quantitatively matches the experimental data
(Figure 2-3C). Overall, we find that our trained model can generate reliable predictions for the
endothelial intracellular signaling response stimulated by FGF or VEGF. This further supports our
assumption that different data sets can be used to establish a predictive model of signaling in ECs.
Thus, we used the best fits (based on model training and validation) in subsequent predictions and
analyses.

2.4.2. FGF produces a greater angiogenic response than VEGF when considering the
maximum ERK phosphorylation produced
Using the validated model, we first studied the effects of FGF and VEGF individually on
pERK. We simulate a range of FGF and VEGF concentrations, based on several experimental
studies published in literature (see Methods). Specifically, we chose the ligand concentration range
of 0.01–2 nM to investigate the ERK activation in response to typical levels of FGF and VEGF in
in vitro studies.
The model predicts that FGF is more potent in promoting ERK phosphorylation, compared
to VEGF, at equimolar concentrations. When FGF concentration is varied from 0.01 nM to 2 nM,
the maximum pERK ranges from 4 10
5
molecules/cell to 8 10
5
molecules/cell, while VEGF

27
induces a maximum of 2 10
-2
molecules/cell to 1 10
5
molecules/cell pERK for the same
concentration range (Figure 2-4A). For example, on average (across the 16 best fits), the maximum
pERK induced by 0.5 nM FGF is 8 × 10
5
molecules/cell, while 0.5 nM VEGF induces a maximum
pERK of 9 × 10
2
molecules/cell. Thus, FGF produces a maximum ERK phosphorylation that is
approximately three orders of magnitude higher than that induced by VEGF. Furthermore, the
maximum pERK is more sensitive to increasing the VEGF concentration, as compared to FGF.
The maximum pERK increases steadily with increasing VEGF stimulation, while maximum pERK
remains relatively constant as the level of FGF stimulation increases. The measurements from
Kanodia et al. show that FGF stimulation does not significantly change the maximal pERK level
(Figure 2-2A) for concentrations ranging from 0.8 ng/ml (0.03 nM) to 100 ng/ml (4 nM) (33). This
experimental observation agrees with our model predictions for FGF stimulation shown in Figure
2-4A. We can explain this result by examining the levels of the intermediate signaling species. For
the FGF pathway, the phosphorylated trimeric complex of FGF, FGFR, and HSGAG binds to and
phosphorylates FRS2, and phosphorylated FRS2 (pFRS2) leads to the phosphorylation of MEK.
The resulting doubly phosphorylated MEK (ppMEK) further mediates the phosphorylation of ERK
(Figure 2-1). We found that even with 0.01 nM FGF stimulation, FRS2 is rapidly depleted
(Additional file 2-1: Figure S4A). The shortage of FRS2 limits ppMEK level, which is the substrate
for ERK phosphorylation and further limits pERK level. Therefore, FRS2 level limits the FGF-
induced ERK phosphorylation. Increasing FGF concentration 200-fold (from 0.01 nM to 2 nM)
only doubles the maximum pERK (from 4 × 10
5
molecules/cell to 8 × 10
5
molecules/cell), again
due to the shortage of FRS2. The importance of FRS2 is also shown in the sensitivity analysis
(Additional file 2-1: Figure S1).

28
On the other hand, for the VEGF pathway, phosphorylated VEGFR2 produces Ras-GTP,
which activates Raf. The activated Raf (aRaf) mediates phosphorylation of MEK. As in the FGFR
pathway, ppMEK mediates ERK phosphorylation (Figure 2-1). The model predicts that there is
enough Raf and MEK supply even upon stimulation with a high concentration of VEGF (2 nM),
as shown in Additional file 2-1: Figure S4B. Thus, the maximum pERK increases significantly
with increasing VEGF concentration, compared with FGF-induced ERK phosphorylation.
Furthermore, the maximum pERK induced by FGF or VEGF gets closer as VEGF concentration
increases (Figure 2-4A).
The model predicts that one of the main reasons why FGF induces a greater maximum
pERK response compared to VEGF is related to differences at the receptor level. Despite depletion
of FRS2, high FGFR levels enable robust FGF-mediated signaling. FGFR density is much higher
than VEGFR2 density (20,000 molecules/cell compared to 1000 molecules/cell). Additionally, the
trafficking parameters (internalization, recycling, and degradation rates) for FGFR are lower than
the corresponding VEGFR2 trafficking parameters (Additional file 2-1: Figure S5). Although
internalized VEGFR2 molecules are recycled back to the surface more rapidly than FGFR
molecules, VEGFR2 is internalized and degraded more rapidly than FGFR. Additionally, the
dynamics of FGF receptors and VEGFR2 in their signaling, internalized, and degraded forms upon
stimulation with 0.5 nM FGF or 0.5nM VEGF (Additional file 2-1: Figure S6) indicate that more
FGFR is available to signal instead of being internalized or degraded (non-signaling), compared
to VEGFR2.

29

Figure 2-4. Predicted maximum pERK response. A. Maximum pERK in response to FGF
(yellow) or VEGF (blue) concentrations varying from 0.01 nM to 2 nM. B. Ratio, R, of
combination effects to the summation of individual effects in response to FGF and VEGF. Each
dot represents one fit. The dots are spread horizontally to avoid overlap of similar responses from
different fits. Asterisk indicates statistically significant difference compared with one (p < 0.05).
Bars are median ± 95% confidence interval.

2.4.3. The combination of FGF and VEGF has greater effects in inducing maximum pERK
than the summation of the individual effects
We next studied the combination effects of FGF and VEGF in inducing maximum pERK.
Here, we define a ratio comparing the combination effects to the individual effects. Specifically,
this ratio is the maximum pERK obtained with co-stimulation of FGF and VEGF over the
summation of the maximum pERK for FGF and VEGF stimulation individually (see Methods for

30
more details). In Figure 2-4B (left panel), 0.5 nM FGF in combination with intermediate to high
VEGF concentrations (0.1–2 nM) can produce a significantly greater maximal pERK response
than the summation of their individual effects, as indicated by the ratios being significantly greater
than one (p < 0.05). The ratios for combinations of VEGF concentrations at 0.01 or 0.05 nM with
0.5 nM FGF are slightly greater than one; however those differences are not statistically
significant. Stimulation with 0.5 nM VEGF in combination with a FGF concentration as low as
0.01 nM can exhibit greater combined effects than the summation of their individual effects. For
these cases, the ratios are all significantly greater than one (Figure 2-4B, right panel).
Co-stimulation compensates for the limitations observed when only one pathway is
stimulated. The model predicts that although VEGF-mediated MEK phosphorylation is much
lower than FGF at the equimolar concentrations (Additional file 2-1: Figure S7A, reactions R26
and R28 compared with Additional file 2-1: Figure S7B, reactions R35 and R378), there are
sufficient levels of Raf available compared to FRS2 (Additional file 2-1: Figure S4) to promote
MEK phosphorylation. Thus, VEGF co-stimulation provides a way to overcome the limitation of
the FRS2 level. On the other hand, in comparison with VEGF stimulation alone, the presence of
FGF in the co-stimulation provides a high level of pMEK because MEK gets phosphorylated by
pFRS2 much faster than by aRaf (Additional file 2-1: Figure S7C, reactions R26 and R28
compared with reactions R35 and R37). Together, these results explain why the combination of
FGF and VEGF has a greater effect on ERK phosphorylation than the summation of their
individual effect.
The effects of FGF and VEGF co-stimulation are more sensitive to VEGF, as compared to
FGF. That is, increasing the VEGF concentration increases the ratio, while the ratio does not
change with varying FGF concentrations. Additionally, the combination of FGF with VEGF shows

31
a more additive response at low VEGF concentrations (< 0.1 nM). At higher concentrations (> 0.1
nM), increasing VEGF concentration increases the ratio (Figure 2-4B).
Overall, the model predictions show that combinations of FGF and VEGF produce more
ERK phosphorylation, compared to their individual effects. Additionally, the model indicates that
VEGF-induced maximum pERK is more sensitive to varying the ligand concentration than FGF-
induced maximum pERK, both for stimulation with VEGF or FGF alone and for co-stimulation.

2.4.4. The combination of FGF and VEGF shows a fast and sustained pERK response
2.4.4.1 The combination of FGF and VEGF exhibits a fast pERK response
In addition to studying the magnitude of the predicted pERK level upon FGF and VEGF
mono- and co-stimulation, we investigated the timescale of the pERK response. First, we analyzed
the time for the pERK level to reach its maximum value, termed “T1” (see Methods for more
details) in response to the stimulation of FGF or VEGF individually. We found that FGF generally
produces a faster response than VEGF stimulation at the same concentrations. Here, we
characterize the timescale of the response in terms of the time it takes to reach maximum pERK.
At low FGF concentrations (< 0.5 nM), the T1 values for FGF and VEGF are not significantly
different (Figure 2-5A). However, at high FGF concentrations (≥ 0.5 nM), FGF shows a
significantly faster T1 response than VEGF. Specifically, for FGF concentrations ranging from 0.5
to 2 nM, the induced pERK response peaks within six minutes, while for the same range of VEGF
concentrations, pERK reaches its peak value within 8 to 22min (Figure 2-5A). Experimental data
from Kanodia et al. show that the values of T1 are all within eight minutes (33); and for 50 ng/ml
(1.1 nM) VEGF stimulation, T1 is 15 minutes (52). Thus, although we did not explicitly fit the

32
model to the T1 values shown in experiments, our model predictions agree with those data.
Together with Figure 2-4A, the model predicts that FGF can induce a greater amount of ERK
phosphorylation within a shorter period of time, compared to VEGF.
This difference in how T1 is affected by the two pro-angiogenic factors is caused by the
availability of upstream species needed to activate MEK, and subsequently, ERK. As described
above, on average, there is more Raf and aRaf available to promote downstream signaling upon
VEGF stimulation, as compared to FRS2 and pFRS2 available with FGF stimulation. At low FGF
concentrations, FRS2 is not depleted as quickly, and more ERK can be phosphorylated. As the
FGF concentration used for stimulation increases, FRS2 is more rapidly depleted, and the maximal
concentration of pERK happens more quickly. On the other hand, at high VEGF concentrations,
there are still sufficient levels of Raf available to become phosphorylated and lead to ERK
phosphorylation. Thus, the time to reach the maximal pERK concentration continues to increase
as VEGF stimulation level increases.
As for the combination effects, we found that when the VEGF concentration is varied from
0.01 nM to 0.5 nM, co-stimulation with 0.5 nM FGF significantly speeds up ERK phosphorylation,
compared to VEGF stimulation alone (Figure 2-5A). For 0.5 nM VEGF stimulation, increasing
the FGF concentration decreases T1, compared to VEGF stimulation alone. This decrease in T1 is
significantly different than VEGF stimulation alone for FGF concentrations greater than 0.5 nM
(Figure 2-5A). Overall, these results indicate that pERK responds faster with FGF stimulation, as
compared to VEGF stimulation.

33

Figure 2-5. Predicted time response of pERK following stimulation by FGF, VEGF, and their
combination. A. T1, time to reach the maximum pERK in response to growth factor stimulation.
Asterisk indicates statistically significant difference compared to corresponding VEGF
concentration (p < 0.05). B. T2, time that pERK is maintained above half of its maximum value in
response to treatments. Each dot represents one fit. The dots are spread horizontally to avoid
overlap of similar responses from different fits. Asterisk indicates statistically significant
difference compared to corresponding FGF concentration (p < 0.05). Yellow: FGF; Blue: VEGF;
Red: combination. Bars are median ± 95% confidence interval.

2.4.4.2 The combination of FGF and VEGF induces sustained pERK response
We explored how long ERK can remain phosphorylated above its half-maximal value,
termed “T2” (see Methods for more details), as another means of characterizing the timescale of
the ERK response. The values of T2 for FGF and VEGF stimulation alone with concentrations

34
ranging from 0.01 nM to 1 nM are not significantly different (T2 is approximately 9min). However,
a higher VEGF concentration produces a more sustained pERK response. Specifically, 2 nM
VEGF shows significantly higher T2 (18 min on average) than 2 nM FGF stimulation (9 min on
average) (Figure 2-5B).
Regarding the combination effects, we found that with 0.5nM FGF, VEGF concentrations
greater than 0.5 nM are able to maintain ERK phosphorylation above its half-maximal value
significantly longer, compared to FGF stimulation alone at the same concentrations. That is, T2 is
significantly greater for combinations of 0.5 nM VEGF with FGF concentrations ranging from
0.01 nM to 2nM, compared to FGF or VEGF stimulation alone (Figure 2-5B).
To identify the reasons why pERK shows a more transient dynamic in response to FGF
stimulation compared to VEGF stimulation at certain concentrations (> 0.5 nM), we compared the
levels of the intermediate signaling species following stimulation of FGF or VEGF alone. For the
co-stimulation of FGF and VEGF, the depletion of FRS2 still limits production of ppMEK with
FGF stimulation (Additional file 2-1: Figure S8), similar to the case of mono-stimulation of FGF
(Additional file 2-1: Figure S4). However, signaling through VEGFR compensates for this
limitation (Additional file 2-1: Figure S8).

2.4.5. Increasing VEGFR2 density can compensate for the relatively low ERK
phosphorylation induced by VEGF
At its baseline level, the density of VEGFR2 is 20 times lower than FGFR1. This large
difference contributes to the predicted results presented above. However, VEGFR2 upregulation

35
has been observed in tumor growth. Experimental measurements of receptor expression show that
some subpopulations of tumor endothelial cells (ECs) have high receptor levels: 13% of tumor-
derived ECs have 7500 VEGFR2 molecules/cell after three weeks of tumor growth, and 5% of the
tumor-derived ECs have 16,200 VEGFR2 molecules/cell after six weeks of tumor growth (63).
Therefore, we sought to understand the effects of varying the VEGF receptor density on the
predicted pERK response to gain some insights into VEGF-mediated signaling in pathological
conditions.
We found that the maximum pERK induced by VEGF increases when VEGFR2 density
increases (Figure 2-6A). At equimolar concentrations of FGF and VEGF stimulation (0.5 nM),
increasing VEGFR2 density by five-fold can increase the maximum pERK level to the same order
of magnitude as FGF stimulation alone (Figure 2-6B), which makes the effects of VEGF-mediated
signaling sizeable in conditions of increased VEGFR2 density. Specifically, the model predicts
that maximum pERK induced by the combination of FGF and VEGF increases more than 90%
when VEGFR2 is increased by five-fold. In contrast, decreasing VEGFR2 by ten-fold leads to an
11.5% decrease in maximum pERK induced by the combination of FGF and VEGF. Thus,
VEGFR2 density significantly impacts ERK phosphorylation with FGF and VEGF co-stimulation.
In addition, the maximum ppERK level is higher upon stimulation by VEGF, compared to FGF,
when VEGFR2 density is increased (Additional file 2-1: Figure S9A). Moreover, although the
reaction rates for ERK phosphorylation by stimulation of FGF are slightly higher than VEGF
during the first 10min, VEGF induces higher rates between 10 to 60min (Additional file 2-1: Figure
S9B, reactions R42 and R43). In addition, VEGF exhibits faster phosphorylation for pERK than
FGF (Additional file 2-1: Figure S9B, reactions R44 and R45). This indicates that the effect of
VEGF is dominant in the combination effect when VEGFR2 density is increased by five-fold.

36
In Figure 2-6C, as VEGFR2 density decreases, the ratio characterizing the ERK signaling
response with a combination of FGF and VEGF compared to the summation of their individual
effects becomes closer to one. This nearly additive combination effect occurs because reducing
VEGFR2 makes the effect of VEGF stimulation negligible. Thus, the ratio is approximately one.
The ratio increases when VEGFR2 density is increased by two-fold, which indicates a stronger
combination effect. However, the summation of individual effects surpasses the combination
effects (the ratio is less than one) when VEGFR2 density is more than five-fold higher than the
baseline level. The reason for this is due to the competition between FGF and VEGF for
downstream resources. Specifically, when VEGFR2 is increased more than two-fold, ERK is
depleted (Additional file 2-1: Figure S10). This makes the combination effects less than the
individual effects, causing the ratio to be less than one. Finally, by increasing VEGFR2, VEGF
more strongly impacts the dynamics of pERK, as both T1 (Figure 2-6D) and T2 (Figure 2-6E)
increase with increasing VEGFR2 density.

37
Figure 2-6. Predicted pERK response with varied initial VEGFR2 concentrations. A.
Maximum pERK induced by 0.5 nM FGF or 0.5 nM VEGF alone. B. Maximum pERK induced
by the combination of 0.5 nM FGF and 0.5 nM VEGF. C. Ratio, R, of 0.5 nM FGF in combination
with 0.5 nM VEGF. D. T1, time to reach the maximum pERK in response to treatments. E. T2,
time that pERK is maintained above half of its maximum value in response to treatments. Yellow:
FGF; Blue: VEGF; Red: combination. Each dot represents one fit. The dots are spread horizontally
to avoid overlap of similar responses from different fits. Bars are median ± 95% confidence
interval.

2.4.6. ERK phosphorylation induced by VEGF can be promoted by decreasing VEGFR2
internalization and degradation rates
Because FGFR trafficking parameter values are lower than the corresponding VEGFR2
trafficking parameters (Additional file 2-1: Figure S5) and pERK is sensitive to these trafficking
rates (Additional file 2-1: Figure S1), we explored the role of trafficking parameters in pERK
response. That is, we investigated the effects of decreasing the VEGFR2 trafficking parameters
individually or together to the same level as FGFR trafficking parameters (Figure 2-7).
Specifically, we simulated the signaling dynamics with 0.5 nM VEGF when all of the VEGFR2
trafficking parameters are decreased to be the same as the FGFR trafficking rates shown in
Additional file 2-1: Figure S5. We then decreased each VEGFR2 trafficking parameter one-by-
one to be the same as the corresponding FGFR trafficking rate shown in Additional file 2-1: Figure
S5. By performing these simulations, we can determine how the trafficking of each pool of
VEGFR2 molecules influences the response to VEGF stimulation.
We found that decreasing the internalization rates of the free and bound forms of VEGFR2,
individually or in combination, with 0.5 nM VEGF leads to a significant increase in maximum
ERK phosphorylation (Figure 2-7A, “fitted” compared to “k_intf”, “k_intb”, or “k_int”). In fact,

38
by decreasing VEGFR2 internalization (“k_int”) to be the same as the FGFR internalization rate,
the maximal VEGF-mediated ERK phosphorylation reaches the same magnitude as the response
induced by 0.5 nM FGF. This qualitative trend is expected, since decreasing receptor
internalization rates makes more VEGFR2 available on the cell surface to bind to the ligand,
inducing downstream signaling. However, the model provides detail about the quantitative effects
of these changes.
Decreasing the recycling rates of the free and bound forms of VEGFR2 individually or in
combination significantly decreases ERK phosphorylation for 0.5 nM VEGF stimulation (Figure
2-7A, “fitted” compared to “k_recf”, “k_recb”, or “k_rec”). A lower receptor recycling rate makes
more non-signaling internalized VEGFR2 remain inside the cell longer and recycles available
VEGFR2 to the cell surface more slowly. Together, these effects limit ERK phosphorylation.
Interestingly, changing the recycling rate leads to a wider range of responses compared to changing
internalization or degradation rates, especially for recycling of free VEGFR2 (Figure 2-7A).
Lowering the rate at which the free form of VEGFR2 is degraded significantly increases
the maximum ERK phosphorylation induced by 0.5 nM VEGF to be the same magnitude of the
maximal pERK level induced by 0.5 nM FGF alone (Figure 2-7A, “fitted” compared to “k_degf”,
“k_degb”, or “k_deg”). A lower free receptor degradation rate makes more VEGFR2 available to
promote signaling.
We also examined the timescales of the ERK response. The model predicts that decreasing
the trafficking rates slows down the dynamics of VEGF-induced ERK phosphorylation (Figure 2-
7B, C), both in terms of T1 and T2. Furthermore, we studied the effects of changing VEGFR2
trafficking parameters on the combination effects. We found that the combination of 0.5 nM FGF
and 0.5 nM VEGF with lower VEGFR2 trafficking parameters has similar results as 0.5 nM

39
VEGF-induced pERK. That is, decreasing VEGFR2 internalization and degradation rates leads to
greater ERK phosphorylation (Additional file 2-1: Figure S11). Overall, lower VEGFR2
trafficking parameters leads to an increased impact of VEGF in the combination effects.

Figure 2-7. Effect of varying VEGFR2 trafficking parameters on pERK response. A.
Maximum pERK, B. T1, and C. T2. The panels show the effect of 0.5 nM FGF (Yellow) or 0.5
nM VEGF (Blue) predicted using the fitted parameter values (“fitted” x-axis label). We ran the
model with 0.5 nM VEGF when all of the VEGFR2 trafficking parameters are decreased (“all” x-
axis label) to be the same as the FGFR trafficking rates shown in Additional file 2-1: Figure S5.
Finally, we decreased each VEGFR2 trafficking parameter individually to be the same as the
corresponding FGFR trafficking rate shown in Additional file 2-1: Figure S5. We omitted some
points for T1 and T2 when the pERK does not reach the maximum value in two hours. Each dot
represents one fit. Bars are median ± 95% confidence interval.

2.5. Discussion
We developed an intracellular signaling model of the crosstalk between two pro-angiogenic
factors, FGF and VEGF. The molecular-detailed model represents the reaction network of
interactions on a molecular level, based on reactions documented in literature. The kinetic
parameters are taken from experimental measurements, where available. Unknown parameters

40
were estimated by fitting the model to experimental data. Additionally, we validated the model
using three separate sets of data.
This is a novel model of FGF and VEGF interactions, taking into account previous
modeling work (33, 36), with a focus on the MAPK cascade and the pERK response as an indicator
for pro-angiogenic signaling. The fitted model predicts the pERK response upon stimulation by
FGF and VEGF, alone and in combination. We particularly focus on the pERK response since
pERK promotes cell proliferation (15), one aspect of early-stage of angiogenesis. Additionally, it
has been shown that pERK is mostly found in active rather than quiescent endothelial cells (47),
and pERK has been used as way to characterize the pro-angiogenic response in other studies (64-
66).
Overall, FGF is predicted to potently and rapidly promote ERK phosphorylation compared
to VEGF stimulation. VEGF also plays an important role in pERK dynamics. Altogether, the
model shows that the pERK level in response to FGF, VEGF, and their combination is dose-
dependent and that some combinations induce greater maximum ERK phosphorylation than the
summation of their individual effects.
Our results reveal that the strength of VEGF-mediated ERK signaling is a combination of
the absolute receptor expression level, the receptor availability, and some intrinsic characteristic
of the receptors or the structure of the signaling pathway. Firstly, there is imbalanced receptor
expression level (high FGFR1 compared to VEGFR2), which is one of the reasons for VEGF’s
lower ERK activation in HUVECs. Increasing the expression of VEGFR2 by five-fold (from 1000
up to 5000 receptors/cell), without changing VEGFR2 trafficking parameters allows the maximum
VEGF-induced pERK level to be the same as what is achieved through FGF stimulation at

41
equimolar concentrations. Additionally, changing the VEGFR2 trafficking rates to be the same as
those for FGFR, without changing VEGFR2 density, also allows the same maximal pERK to be
achieved for equimolar concentrations of FGF and VEGF. Thus, both the absolute receptor level
and availability of receptors directly affect the signal strength. However, these are not the only
causes of the difference in the ability to drive ERK signaling, because even with a five-fold
increase in VEGFR2 (to 5000 molecules/cell), its expression is still less than that of FGFR, which
is present at 20,000 molecules/cell. Thus, to reach the same maximum pERK level at equimolar
concentrations of FGF and VEGF stimulation, the required VEGFR2 level is much lower than
FGFR level, independent of the high VEGFR2 trafficking rates. This indicates that there is also
some intrinsic ability of VEGFR2 that enables ERK activation. This is particularly relevant, as
endothelial VEGF receptor expression is upregulated in tumors (60), and can reach nearly 2 10
4

molecules/cell in tumor-derived endothelial cells after six weeks of tumor growth (63). The effect
of VEGF-mediated pERK signaling may also be due to the structure of the signaling network and
expression levels of the pathway intermediates (such as Raf, which is not present in the FGF
signaling pathway). Overall, the ability of VEGF to promote ERK signaling is due to a
combination of factors. Excitingly, our model is able to predict the contribution of each of these
factors.
The model predictions are consistent with several experimental studies. Multiple
experiments show that FGF induces the same level of angiogenic response at lower concentration
in comparison to VEGF (50, 51, 67), and their combination induces greater angiogenic responses
(9). Additionally, the model predicts that decreasing VEGFR2 internalization and degradation
rates can increase the impact of VEGF in combination effects. This result complements
experiments showing that receptor trafficking plays a critical role in angiogenic signaling (68).

42
Overall, our molecular-detailed model helps synthesize these experimental data and observations
related to VEGF- and FGF-stimulated signaling.
One application of our work is that the model can also be linked with computational models
that predict events on the cellular scale. Our model culminates with ERK activation,
complementing published models that substantially simplify the intracellular signaling and focus
on specific cellular behavior, such as proliferation (30), the probability of sprout formation and the
speed of vessel growth (31), or tumor growth (32). However, these models reduced the intracellular
signaling network such that the output signal is simply linearly proportional to the fraction of
bound receptors. In comparison, our mechanistic model considers intracellular signaling and
quantitatively analyzes pERK response, which could be a better indicator for these cellular
behaviors. For example, Hendrata and Sudiono constructed a computational model that includes
molecular, cellular, and extracellular scales to study tumor apoptosis (32). Our model can be
utilized in combination with such models to more accurately predict cellular behavior.
Our model can also be used for exploring mechanisms that regulate the magnitude and
dynamics of pERK upon FGF and/or VEGF stimulation, as has been done in other modeling work
(69-71). For example, Edelstein et al. showed that ligand depletion diminishes cooperative
interactions between ligands and binding sites, and that receptor concentration plays an important
role in biological signal transduction (69). Such depletion of the ligand that initiates the signaling
could also be explored using our mechanistic model. In other work, Saucerman and Bers combined
a cardiac myocyte excitation-contraction computational model with biochemical reaction models
to investigate how calmodulin (CaM), calcineurin, and CaM-dependent kinase are spatially and
temporally activated by local calcium signals (70). Our model can be expanded to explore spatial
effects as well. Recently, Romano and coworkers studied the competition of seven proteins for

43
CaM binding and concluded that this competition contributes to synaptic plasticity (71). This
model of binding competition is relevant to our system, for example in the case of competitive
binding the activators (FRS2 or aRaf) and phosphatases to species such as MEK. Such competition
can be examined in detail in future work.
Our model can also aid in studying the efficiency of pro- or anti-angiogenic therapies.
Some pro- or anti-angiogenic treatments have not been very effective, particularly those targeting
only a single signaling family (6, 7). However, targeting both FGF and VEGF may be a promising
strategy, given the potential synergistic effects predicted by our model and demonstrated in
experimental studies. In fact, multiple groups have reported interesting interactions between FGF
and VEGF (9, 11, 72, 73). This crosstalk may be exploited to aid in angiogenesis-based therapies,
and our model can be helpful in understanding their interactions and combination effects. Model
predictions for species’ dynamics and reaction rates provide mechanistic insight into FGF and
VEGF interactions. Our predictions show that the low success in targeting VEGF alone could be
due to low receptor numbers and fast internalization, recycling and degradation. Although FGF
has greater effects in inducing ERK phosphorylation, its effects can be enhanced by the addition
of VEGF. Thus, our model can be used to investigate the efficiency of targeting both FGF and
VEGF as an alternative strategy.
In addition to the amount of FGF or VEGF the cells are stimulated with, other factors can
influence the magnitude of timescale of the pERK response. This includes the growth factor
concentration gradient, heterogeneity in a population of cells, and genetic mutations. Our
computational modeling provides a platform for many interesting and relevant studies that can be
helpful to characterize signaling dynamics that mediate endothelial cell sprouting during the early
stages of angiogenesis, in response to extracellular signals.

44
Our model is the first to combine the signaling networks of FGF and VEGF, providing
novel quantitative insight into the effect of combined FGF and VEGF treatment. However, we
recognize some limitations in our model. Firstly, this model does not include FGF activated Ras-
Raf signaling because the protein-protein interactions in this pathway are still not clear. As more
information becomes available as to the detailed mechanisms of those reactions, we can expand
the model to include FGF-mediated Ras signaling. Second, we assumed that the internalized and
degraded phosphorylated receptors would not signal. Receptor trafficking processes such as
internalization have conventionally been thought to downregulate extracellular signals. However,
data suggest that VEGFR2 may signal even when internalized (74-76). Since the role of
internalized receptors is still somewhat debated and to focus on the signaling dynamics mediated
by cell surface receptors, we chose to exclude the effect of internalized receptors in our model.
This assumption can be relaxed in future studies. Additionally, all bound forms of FGFR1 are
assumed to have the same internalization, recycling, and degradation rates as a simplification and
because there are conflicting values reported in literature (37, 77). We tried various FGFR1
trafficking parameters; however, this did not significantly change the model predictions or our
overall conclusions. Third, this model only includes VEGFR2, although VEGF binds to VEGFR1
and neuropilin-1 (NRP1). These receptors also contribute to angiogenesis and may be incorporated
into the model in future studies. Finally, we studied pERK response in two hours. We omit ligand
secretion and protein degradation during this time and do not predict long-term responses. In the
future, we can expand our model to predict the cellular response over a longer period of time.

45
2.6. Conclusions
In summary, our molecular-detailed model quantifies ERK phosphorylation upon
stimulation by two major pro-angiogenic factors, FGF and VEGF, and provides insights into the
molecular interactions between these proteins. Specifically, the model predicts the combination
effects of FGF and VEGF on ERK phosphorylation and quantitatively shows the magnitude and
time scale of the pERK response. Because of the complexity of this biological system, it may be
challenging to get a comprehensive understanding of the system using experiments that only focus
on a few molecular species. Our computational modeling provides a quantitative framework to
explore the system as a whole, generating novel mechanistic insight and complementing
experimental studies.

2.7. Acknowledgements
The authors thank members of the Finley research group for critical comments and
suggestions. The authors acknowledge the support of the US National Science Foundation
(CAREER Award 1552065).

2.8. Additional files
Additional file 2-1: Supplementary figures. (PDF 2453 kb)
Additional file 2-2: Supplementary tables. (XLSX 32 kb)
Additional file 2-3: MATLAB file containing computational model. (DOCX 15 kb)

46
3. Chapter 3
Chapter 3
FGF- and VEGF-induced MPAK and PI3K/Akt signaling

Portions of this chapter are adapted from:
Min Song and Stacey D. Finley. Cell Commun. Signal. (2020)

3.1. Abstract
Angiogenesis plays an important role in the survival of tissues, as blood vessels provide
oxygen and nutrients required by the resident cells. Thus, targeting angiogenesis is a prominent
strategy in many different settings, including both tissue engineering and cancer treatment.
However, not all of the approaches that modulate angiogenesis lead to successful outcomes.
Angiogenesis-based therapies primarily target pro-angiogenic factors such as vascular endothelial
growth factor-A (VEGF) or fibroblast growth factor (FGF) in isolation, and there is a limited
understanding of how these promoters combine together to stimulate angiogenesis. Targeting one
pathway could be insufficient, as alternative pathways may compensate, diminishing the overall
effect of the treatment strategy.
To gain mechanistic insight and identify novel therapeutic strategies, we have developed a
detailed mathematical model to quantitatively characterize the crosstalk of FGF and VEGF
intracellular signaling. The model focuses on FGF- and VEGF-induced mitogen-activated protein
kinase (MAPK) signaling to promote cell proliferation and the phosphatidylinositol 3-

47
kinase/protein kinase B (PI3K/Akt) pathway, which promotes cell survival and migration. We fit
the model to published experimental datasets that measure phosphorylated extracellular regulated
kinase (pERK) and Akt (pAkt) upon FGF or VEGF stimulation. We validate the model with
separate sets of data.
We apply the trained and validated mathematical model to characterize the dynamics of
pERK and pAkt in response to the mono- and co-stimulation by FGF and VEGF. The model
predicts that for certain ranges of ligand concentrations, the maximum pERK level is more
responsive to changes in ligand concentration compared to the maximum pAkt level. Also, the
combination of FGF and VEGF indicates a greater effect in increasing the maximum pERK
compared to the summation of individual effects, which is not seen for maximum pAkt levels. In
addition, our model identifies the influential species and kinetic parameters that specifically
modulate the pERK and pAkt responses, which represent potential targets for angiogenesis-based
therapies.
Overall, the model predicts the combination effects of FGF and VEGF stimulation on ERK
and Akt quantitatively and provides a framework to mechanistically explain experimental results
and guide experimental design. Thus, this model can be utilized to study the effects of pro- and
anti-angiogenic therapies that particularly target ERK and/or Akt activation upon stimulation with
FGF and VEGF.

3.2. Introduction
Angiogenesis is the formation of new blood capillaries from pre-existing blood vessels.
The essential role of blood vessels in delivering nutrients to tissues makes angiogenesis important
in many different settings, including both physiological and pathological conditions.

48
Physiologically, angiogenesis is involved in the growth of normal blood vessels during
development such as placental vascularization during pregnancy (78, 79) and the wound healing
process (80, 81). Pathological angiogenesis is crucial in many diseases, including cancer (82).
Thus, targeting angiogenesis is a prominent strategy in many contexts, for example, in both tissue
engineering and cancer treatment. In the context of tissue engineering, researchers have sought to
create artificial tissues to substitute damaged tissues in response to a great shortage of donors for
transplant surgery. Implementing strategies that promote the formation of adequate vasculature is
critical for the long-term viability of engineered tissue constructs. Therefore, stimulating new
blood vessel formation is an important strategy for tissue engineering (44). On the other hand,
inhibiting angiogenesis is a strategy for cancer treatment, as the formation of new blood vessels is
important for cancer growth and metastasis. Therefore, understanding the angiogenesis process is
very beneficial to current strategies that target vessel formation.
Many different pro-angiogenic growth factors, such as fibroblast growth factor (FGF),
vascular endothelial growth factor (VEGF), and platelet-derived growth factor (PDGF), mediate
angiogenesis (2, 3). These factors promote different cellular processes involving endothelial cells
leading to new blood vessel formation, including proliferation, migration, survival, and vessel
maturation (4, 5). Strategies to promote or inhibit angiogenesis focus on modulating the effects of
the factors that promote these cellular-level processes.
Unfortunately, not all approaches to promote or inhibit angiogenesis lead to successful
outcomes. For example, clinical trials have shown no effective improvement in angiogenesis upon
stimulation by FGF (6) or VEGF (7). Also, bevacizumab, an anti-angiogenic agent designed to
sequester VEGF extracellularly, inhibiting VEGF-mediated signaling by preventing VEGF from
binding to its receptor (83, 84), has limited effects in certain cancer types, and it is no longer

49
approved for the treatment of metastatic breast cancer due to disappointing results in patients (8).
Thus, there is a need to better understand the molecular interactions and signaling required for new
blood vessel formation, in order to establish more effective therapeutic strategies.
Given the complex set of biochemical reactions comprising angiogenesis signaling
networks, it is essential to apply computational modeling to better understand the dynamics of
these networks. Computational modeling serves as a powerful tool to investigate molecular
responses mechanistically and to guide experimental design. Indeed, many models have been
developed to explore the angiogenic response mediated by growth factors. Models focused on the
extracellular-level interactions (27, 28, 85, 86) enhance our understanding of the distribution of
angiogenic factors, which affects downstream angiogenic signaling. These models can be used to
study strategies that regulate the distribution of angiogenic factors in tumor tissue. For example,
Li and Finley constructed a compartmental whole-body model to study the effect of anti-
angiogenic therapies targeting VEGF and TSP1 signaling in a simulated cancer patient cohort
(85). Also, models that study intracellular signaling (34, 38) can help identify potential targets and
explore their efficacy for pro- or anti-angiogenic therapies.
Modulating angiogenesis signaling networks can involve targeting multiple angiogenic
factors. There are few models that simulate the effects of more than one factor on intracellular
signaling reactions at a detailed level. However, this insight is needed to better understand the
effect of the angiogenic factors, mechanistically study experimental data, and guide new
experiments. Moreover, in the case of inhibiting angiogenesis, tumors often evade the effects of
drugs that target a single factor by making use of alternate compensatory pathways to activate
signaling species needed for proliferation and migration. For instance, FGFR activation may play
a role in the resistance mechanism of anti-angiogenic drugs, especially anti-VEGF treatment (12,

50
13). Additionally, experiments show high levels of FGFR1 in tumors that continue to progress,
even during anti-VEGF therapy (14). FGF and VEGF have been shown to be particularly important
in the early stages of angiogenesis, and we are interested in signaling crosstalk between these
factors required to initiate vessel growth. Thus, we aim to quantitatively investigate the
combination effects of FGF and VEGF on activating signaling in endothelial cells and identify
potential intracellular targets by building a molecular-detailed computational model that
incorporates the crosstalk between these pathways.
Specifically, in our model, FGF and VEGF bind to their receptors and initiate signaling
through the mitogen-activated protein kinase (MAPK) and phosphatidylinositol 3-kinase/protein
kinase B (PI3K/Akt) pathways to phosphorylate ERK and Akt, respectively. Previously, we
studied the response of phosphorylated ERK (pERK) upon stimulation by FGF and VEGF (35),
using mathematical modeling to gain insight into proliferation signaling, one aspect of the early
stage of angiogenesis. However, angiogenesis involves not only proliferation, but also survival
and migration of the endothelial cells. To get a more comprehensive understanding of this process,
we expand the model to now incorporate the PI3K/Akt pathway, which has been shown to play an
important role in cell survival (17-21) and migration (21, 24, 25). Thus, in this study we examine
the responses of pERK and phosphorylated Akt (pAkt) following mono- and co-stimulation by
FGF and VEGF using mathematical modeling. The model predicts that the maximum pERK level
is more responsive to changing the ligand concentration compared to the maximum pAkt level for
certain concentration ranges. Also, co-stimulation with FGF and VEGF indicates a greater effect
in increasing the maximum pERK compared to the summation of individual effects, which is not
seen for maximum pAkt levels. Using this model, we also identified the influential species and
kinetic parameters that specifically regulate the pERK and/or pAkt responses, indicating potential

51
targets for pro- or anti-angiogenic therapies. The model predictions provide mechanistic insight
into FGF and VEGF interactions in angiogenesis signaling. More broadly, this model provides a
framework to study the efficacy of angiogenesis-based therapies.

3.3. Methods
3.3.1 Model construction
We constructed a molecular-detailed model that describes the intracellular network of
FGF- and VEGF-induced ERK and Akt phosphorylation in endothelial cells. The model
significantly expands on previous modeling work studying ERK (35) and Akt (38). In our model,
FGF binding to FGFR1 and HSGAG activates FRS2 and then initiates PI3K/Akt pathway, and
VEGF binding to its receptor, VEGFR2, phosphorylates VEGFR2 and activates PI3K directly. In
addition, activated FRS2 and Raf trigger MAPK pathway upon stimulation by FGF and VEGF,
respectively.
The molecular interactions involved in the network are illustrated in Figure 3-1. This
network is implemented as an ordinary differential equation (ODE) model using MATLAB. The
main model includes 97 reactions, 99 species, and 100 parameters (see Additional file 3-3). The
reactions, initial conditions and parameter values are listed in Additional file 3-1: Tables S1-S3.
All reactions are assumed to follow the law of mass action. Receptor internalization, recycling,
and degradation are considered in the model, as these processes occur on a relatively fast timescale.
However, because the simulated time is within two hours, we do not consider the degradation of
the ligands or signaling species.
We note that the concentrations of extracellular ligands and intracellular species are
considered with different relative volumes. Specifically, the concentrations of extracellular ligands

52
are expressed relative to the volume of the cell culture media, while the concentrations of
intracellular species are usually considered in a volume of a cell. In this study we focused on
endothelial cells, which has been reported to have a mean cell volume of 1009 ± 180 𝜇 m
3
(1.01 ±
0.18 pL) (87). Therefore, we used 1 pL cell volume to convert concentrations of intracellular
species from molecules/cell to nM. This same conversion factor has been used in other
computational work to study endothelial cell signaling responses (36, 38). Details regarding
interconversion between these units are provided in Additional file 3-1.

Figure 3-1. Schematic of FGF and VEGF signaling network. Signaling is induced by growth
factors binding to their receptors, culminating with phosphorylation of ERK and Akt, through the
MAPK and PI3K/Akt cascades, respectively.

Raf aRaf
pFRS2 FRS2 FRS2u
Ras-GTP Ras-GDP
pMEK MEK ppMEK
ERK pERK ppERK
Ptase 2
Ptase 3
HSGAG FGF-2 FGFR-1 VEGF-A VEGFR-2
PI3K
pFRS2:pPI3K
pR2:pPI3K
PIP3
Akt
pAkt ppAkt
PIP2
PTEN
PP2A PP2Aoff
Ptase 1
+PDK1
pF:H:R:R:H:F
pR2
PTP1B

53
3.3.2 Sensitivity analysis
Before fitting the model to experimental data, we first performed a sensitivity analysis to
identify the parameters and initial concentrations that significantly influence the model outputs for
model training, using the extended Fourier Amplitude Sensitivity Test (eFAST) (57) method.
Since the parameters and initial values for ERK activation were fit to experimental data in our
previous model (35), we used the best fit values and held them constant during the sensitivity
analysis. All remaining model parameters and initial values were varied simultaneously within two
orders of magnitude above and below the baseline values, where the baseline values were taken
from published literature (38). In this way, the effects of multiple model inputs (kinetic parameters
or initial conditions) on the pERK and pAkt concentrations were computed (the total sensitivity
indices, “Sti”). The Sti values can range from 0 to 1, where a higher Sti index indicates the input
is more influential to the output. Based on the experimental data that are used for model training,
we calculated the Sti values using eFAST for all the same concentrations and time points as what
was used in the experiments. The highest Sti value (Stimax) across all of the concentrations and time
points was selected to represent the sensitivity index for each variable.
We also performed eFAST for the trained and validated model to identify potential targets
for pro- and anti-angiogenic strategies. All parameters and initial concentrations were varied
simultaneously within two orders of magnitude above and below the baseline values. In this case,
the baseline values for the fitted variables were the median values estimated from model fitting.
We calculated the Sti values to quantify how all the variables affected pERK and pAkt. Based on
the behaviors of maximum pERK and pAkt that reach a plateau as the FGF and/or VEGF
concentration increases, we selected five representative concentrations to capture the low (0.01
nM and 0.05 nM), intermediate (0.1 nM), and high (0.5 nM and 1 nM) levels of responses. We

54
calculated the Sti values using eFAST for these five concentrations of FGF and VEGF stimulation
at 15 time points ranging from zero to 120 minutes. Again, the Stimax across all the concentrations
and time points were compared for all the variables.

3.3.3. Data extraction
Data from published experimental studies (33, 52, 54, 88-90) were used for parameter
fitting and model validation. The Western blot images were analyzed using ImageJ. Experimental
data from plots was extracted using the grabit function in MATLAB.

3.3.4. Model fitting and validation
Nine influential variables with Stimax values greater than 0.8 were identified by performing
eFAST (Additional file 3-1: Table S4). The value of 0.8 was chosen as the cutoff to balance the
fitting results and the computational expense. However, this included two correlated parameters,
the kinetic rates k_pFRS2PI3K and kd_pFRS2PI3K, which are the forward and reverse rates,
respectively, of the reaction pFRS2+PI3K ↔ pFRS2:pPI3K. Thus, we held k_pFRS2PI3K constant
and fitted the rest of the influential variables (Additional file 3-1: Table S4).
Therefore, a total of eight variables (four initial conditions and four kinetic parameters)
were estimated by fitting the model to experimental data using Particle Swarm Optimization (PSO)
(59). PSO starts with a population of initial particles (parameter sets). As the particles move around
(i.e., as the algorithm explores the parameter space), an objective function is evaluated at each
particle location. Particles communicate with one another to determine which has the lowest
objective function value. The objective function for each parameter set was used to identify

55
optimal parameter values. Specifically, we used PSO to minimize the weighted sum of squared
residuals (WSSR):
WSSR(θ)=𝑚𝑖𝑛 ∑(
𝑉 𝑝𝑟𝑒𝑑 ,𝑖 (θ)−𝑉 𝑒𝑥𝑝 ,𝑖 𝑉 𝑒𝑥𝑝 ,𝑖 )
2
𝑛 𝑖=1

where Vexp,i is the ith experimental measurement, Vpred,i is the ith predicted value at the
corresponding time point, and n is the total number of experimental data points. The minimization
is subject to 𝜃 , the set of upper and lower bounds on each of the fitted parameters. The bounds
were set to be two orders of magnitude above and below the baseline parameter values, which
were taken from literature.
The model was fitted using four datasets as shown in Figure 3-2, represented by circles.
We note that the datasets shown in Figure 3-2C and D were also used for fitting pERK levels in
our previous work (35). Here, we wanted to ensure that the model can still match this data, even
upon expanding the model to include the PI3K/Akt pathway.
Model simulations were compared to experimental measurements. Specifically, the
relative change of the responses was calculated as following:
relative change(t)=
𝑟𝑒𝑠𝑝𝑜𝑛𝑠𝑒 (t)−𝑟𝑒𝑠𝑝𝑜𝑛𝑠𝑒 (𝑡 𝑟𝑒𝑓 )
𝑟𝑒𝑠𝑝𝑜𝑛𝑠𝑒 (𝑡 𝑟𝑒𝑓 )

where 𝑟𝑒𝑠𝑝𝑜𝑛𝑠𝑒 (t) is the level of pERK, pAkt, or phosphorylated VEGFR2 (pR2) at time t, and
𝑟𝑒 𝑠𝑝𝑜𝑛𝑠𝑒 (𝑡 𝑟𝑒𝑓 ) is the response (pERK, pAkt, or pR2) at a reference time point tref. Here, the pERK
and pAkt in the model simulation include all free and bound forms of singly- and doubly-
phosphorylated ERK and Akt, respectively.
We note that the FGF-induced pERK response reported by Kanodia et al. was measured
using the non-small cell lung cancer cell line NCI-H1730, while the rest of experimental
measurements mentioned above were obtained using human umbilical vein endothelial cells

56
(HUVECs). As in previous work (35), we assumed the FGFR2 signaling kinetics for NCI-H1730
are the same as in HUVECs, as FGFR1 and HSGAG levels are fairly consistent for various cell
types (33, 37, 61, 62).
We first fitted the model 40 times to experimental data. However, we noticed that for the
parameter sets with the lowest error, the fitted values for the k_pFRS2fPIP3 parameter were all at
the upper bound. To exclude the possibility of arbitrary bounds limiting the parameter searching
space, we took this upper bound (20 s
-1
) as the baseline value and expanded the bounds for this
parameter to be 0.2 – 2,000 s
-1
. The targeted variables were estimated another 40 times with the
new bounds. With this second round of fitting, none of the parameters were estimated to be at one
of the bounds (Additional file 3-1: Table S5).
After model training, we validated the model with three datasets not used in the fitting. We
first predicted the VEGF-induced pR2 relative change upon stimulation with 80 ng/ml (1.78 nM)
VEGF (54). We also simulated the change of pAkt upon stimulation with 10 ng/ml (0.43 nM)
FGF- (89) or 20 ng/ml (0.44 nM) VEGF-induced (90), respectively.
For all three datasets, we simulated the experimental conditions without any additional
model fitting and compared to the experimental measurements. A total of 15 parameter sets with
the smallest errors were taken to be the “best” sets based on the model fitting and validation
(Additional file 3-2: Figure S2 and Additional file 3-1: Table S5) and were used for all model
simulations.

3.3.5. Monte Carlo simulations
To study the robustness of the system, we ran the fitted model 1,000 times by generating
1,000 values for all parameters and non-zero initial concentrations, sampling from a normal

57
distribution. For initial concentrations and parameters that were estimated by fitting to the
experimental data, the mean values ( ) were the median of the fitted values, and we used the
standard deviation ( ) calculated from the fitted parameter sets. For all other model variable
values, we set  to be the baseline values and calculated  to capture 99.7% of the possible values
given the range of   50%  (i.e.,   3). It is worth noting that with this sampling, it is possible
to get negative values, though this is unlikely to occur. However, if any negative values were
selected, we resampled until all the sampled variables are positive.

3.3.5.1. Signaling responses
We investigated the ERK and Akt phosphorylation responses upon stimulation by FGF or
VEGF alone or in combination.
a. Maximum pERK and pAkt. In our model simulations, for simplification, representative values
were used as indicators for the magnitude of pERK or pAkt responses, specifically the
maximum values. We calculated the maximum ERK and Akt phosphorylation levels induced
by the stimulation by FGF, VEGF, or their combination.

b. Ratio, R. To compare the combination effects to the effects of FGF and VEGF individually,
we introduce the ratio below:

𝑅 (𝑟𝑒𝑠𝑝𝑜𝑛𝑠𝑒 )=
𝑚𝑎𝑥 response(FGF and VEGF)
𝑚𝑎𝑥 response(FGF)+𝑚𝑎𝑥 response(VEGF)

When R is greater than one, it indicates that the combination effect in inducing the
maximum response is greater than the summation of individual effects; when R is equal to one,

58
it implies that the combination effect is additive; when R is less than one, it means that the
combination of FGF and VEGF effect does not surpass the summation of individual effects
and implies a competitive effect.

c. Fold change, F. To explore the efficiency of varying the identified influential variables, we
define F as the predicted maximum pERK and pAkt levels when the parameters in Table 3-2
were varied by 0.1- or 10-fold individually compared to the baseline values. When F is greater
than one, varying the parameter enhances the response; when F is equal to one, varying the
parameter has no effect on the response; when F is less than one, it indicates that varying the
parameter inhibits the response.

3.4. Results
3.4.1 The fitted and validated molecular-detailed mathematical model captures the major
characteristics of FGF- and VEGF-induced ERK and Akt phosphorylation dynamics
We developed an intracellular signaling model of the crosstalk between two pro-angiogenic
factors, FGF and VEGF. The signaling is initiated by FGF binding to FGF receptor 1 (FGFR1)
and heparan sulfate glycosaminoglycans (HSGAGs) or VEGF binding to VEGF receptor 2
(VEGFR2), both promoting downstream signaling (Figure 3-1). The model focuses on FGF- and
VEGF-induced signaling through the MAPK and PI3K/Akt pathways, leading to activation of
ERK and Akt, respectively. We consider that activated ERK and Akt include both the singly and
doubly phosphorylated forms of each species (i.e., pERK and ppERK, as well as pAkt and ppAkt).
For simplicity, we collectively refer to these species as phosphorylated ERK and Akt, pERK and

59
pAkt, respectively. The model reactions, initial conditions, and parameter values are given in
Additional file 3-1: Tables S1 to S3.
The parameters and initial concentrations involved in FGF- and VEGF-initiated MAPK
signaling are taken from the “best” fit from our previous model (35). As a starting point, the newly
introduced parameters and initial concentrations involved in the PI3K/Akt pathway are acquired
from literature (38). The influential model parameters and initial conditions were estimated by
fitting the model to experimental data, as described below.
For model training, we aimed to first identify the model variables (kinetic parameters and
initial concentrations) that significantly influence the model outputs, phosphorylated ERK and
Akt. To do so, we performed the eFAST (57) (see Methods for more details) and analyzed the total
sensitivity indices (Sti) for the species’ concentrations and kinetic rates that are involved in the
PI3K/Akt pathway. The highest Sti values (Stimax) across all of the concentrations and time points
for the 32 newly introduced variables were compared, and nine of them (Additional file 3-1: Table
S4) were identified as influential to pERK and/or pAkt induced upon stimulation by FGF and
VEGF. Of these, eight were not correlated (denoted by red text in Additional file 3-1: Table S4),
and we estimated their values by fitting the model to experimental measurements (33, 52, 88) using
PSO (59) (see Methods for more details).
The fitted model shows a good agreement with experimental results (Figure 3-2). It
quantitatively captures the FGF- and VEGF-induced pAkt dynamics from experimental
observations (52, 88) (Figure 3-2A-B). In addition, this expanded model retains the ability to
reproduce the measured pERK levels promoted by VEGF or FGF stimulation (Figure 3-2C-D),
including the biphasic pERK dose response following stimulation with FGF (Additional file 3-2:
Figure S1) reported by Kanodia et al. (33), which our previous model also reproduced. The

60
weighted errors for 15 best fits are all approximately 20.7 (Additional file 3-1: Table S5). Also,
the estimated values of the fitted variables show good consistency (Additional file 3-2: Figure S2).

Figure 3-2. Model comparison to training data for FGF or VEGF stimulation. A. Relative
change of pAkt for 100 ng/ml (4.35 nM) FGF stimulation compared with the reference time point
(10 min). B. Relative change of Akt phosphorylation upon stimulation with 50 ng/ml (1.11 nM)
VEGF compared with reference time point (60 min). C. Relative change of ERK phosphorylation
following stimulation with 50 ng/ml (1.11 nM) VEGF compared to the pERK level at a reference
time point (30 min). D. Normalized pERK dynamics in response to FGF concentrations ranging
from 0.16 to 500 ng/ml (0.007–21.74 nM), where pERK level was normalized by the maximum
pERK stimulated by FGF across all six concentrations in 2 hours. Circles in Panels A-C are
experimental data from HUVECs, and circles in Panel D are experimental data from the NCI-
H1730 cell line. Curves are the mean model predictions of the 15 best fits. Shaded regions show
standard deviation of the fits

In addition to matching data used for fitting, the model is consistent with independent
experimental observations. To validate the model, the predictions were compared to three
Figure 2
A B C
D

61
independent sets of experimental data (Figure 3-3). The model-predicted pAkt dynamics with 10
ng/ml (0.43 nM) FGF- or 20 ng/ml (0.44 nM) VEGF-induced pAkt agree with additional
experimental data from Pisanti et al., 2011 (89) and Schneeweis et al., 2010 (90), respectively
(Figure 3-3A-B). In addition, model predictions match the levels of VEGFR2 phosphorylation
following stimulation with 80 ng/ml (1.78 nM) VEGF extracted from a separate set of data from
Chabot et al., 2009 (54) (Figure 3-3C).

Figure 3-3. Model comparison to validation data. A. Relative change of pAkt upon stimulation
with 10 ng/ml (0.43 nM) FGF compared with reference time point (30 min). B. Relative change
of Akt phosphorylation upon stimulation with 20 ng/ml (0.44 nM) VEGF compared with reference
time point (45 min). C. Relative change of VEGFR2 phosphorylation upon stimulation with 80
ng/ml (1.78 nM) VEGF using the reference time point of 7 min. Circles are experimental data from
HUVECs. Curves are the mean model predictions of the 15 best fits from the model training.
Shaded regions show standard deviation of the fits.

We performed Monte Carlo simulations (see Methods for more details) to investigate the
predicted pERK and pAkt levels given variability in the initial conditions and parameters. The
model predictions with parameters values randomly varied within the range of the estimated values
can still capture pERK, pAkt, and pR2 dynamics stimulated by FGF and VEGF (Additional file 3-
2: Figures S3 and S4). These Monte Carlo simulations suggest that the overall dynamics of the
model outputs, pERK and pAkt, are relatively robust to variability or uncertainty in initial species’
concentrations and parameters in the signaling network.
Figure 3
A B C

62

3.4.2. FGF produces greater maximum pAkt and pERK than VEGF at equimolar
concentrations
We first explored the individual effects of FGF and VEGF on pERK and pAkt using the
trained and validated model. The dynamics of pERK and pAkt stimulated by 0.5 nM FGF and 0.5
nM VEGF are shown in Figure 3-4A. The model predicts transient activation of ERK and Akt
following stimulation by FGF or VEGF. The species’ concentrations are predicted to peak within
30 minutes and return to basal level after 60 minutes, as seen in experimental data used for model
fitting and validation. These predicted time courses show that 0.5 nM FGF stimulation leads to
higher maximum levels of pERK and pAkt, compared to 0.5 nM VEGF stimulation.
We also simulated the pERK and pAkt dynamics for a range of concentrations of FGF or
VEGF. Here, we use the maximum pERK and pAkt levels within the two hours simulated by the
model as indicators for pERK and pAkt responses. Maximum pAkt (green) and pERK (purple)
levels increase with the increase of FGF or VEGF concentrations (Figure 3-4B). The model
predicts that the average levels of the maximum pAkt and pERK across the 15 best fits, given 0.5
nM FGF stimulation, are 1.8 10
3
nM and 1.0 10
3
nM, respectively. In comparison, 0.5 nM VEGF
induces an average maximum pAkt and pERK of 6.3 10
2
nM and 5.2 10
-1
nM, respectively. Thus,
0.5 nM FGF produces averaged maximum Akt and ERK phosphorylation levels that are 3-fold
and nearly 2000-fold higher than that induced by 0.5 nM VEGF, respectively.
We can use the detailed model to explain these results. In our previous work (35), we
described that the main reasons that FGF induces a greater maximum pERK response compared
to VEGF are the relatively high level of FGFR density compared to VEGFR2 (33.2 versus 1.7 nM;
20,000 versus 1,000 molecules/cell, respectively) and lower internalization and degradation rates

63
for FGFR compared to the corresponding VEGFR2 parameters. These differences also make FGF-
induced pAkt higher than VEGF-induced pAkt. Indeed, the model predicts that increasing
VEGFR2 level by 10-fold can increase the 0.5 nM VEGF-induced maximum pAkt to
approximately the same maximum pAkt level induced by 0.5 nM FGF (Additional file 3-2: Figure
S5A). In addition, decreasing VEGFR2 internalization and degradation rates to be the same level
as the corresponding FGFR internalization and degradation rates leads to an increase in VEGF-
induced maximum pAkt level (Additional file 3-2: Figure S5B). This is because lower
internalization and degradation rates lead to more signaling complexes available for signal
transduction. Together, these receptor-related properties (density, internalization, and degradation)
lead to stronger signaling induced by FGF.

Figure 3-4. Predicted pERK and pAkt responses stimulated by single agents. A. Predicted
time courses of pERK and pAkt stimulated by 0.5 nM FGF and 0.5 nM VEGF. Curves are the
mean predictions for the 15 best fits from the model training. Shaded regions show standard

64
deviation of the fits. B. Maximum pERK (Purple) and pAkt (Green) in response to FGF (left) and
VEGF (right) for concentrations varying from 0.01 nM to 1 nM. Bars are mean ± standard
deviation of model predictions. Note that the y-axes are not on the same scale.

3.4.3. Akt activation shows a stronger response than ERK in terms of magnitude
We also compared Akt and ERK activation responses by the mono- and co-stimulation of
FGF and VEGF. The model predictions show that the maximum pAkt level is higher than the
maximum pERK level in response to same ligand concentration, whether considering FGF or
VEGF stimulation (Figure 3-4B and Figure 3-5). As shown in Figure 3-4B, when FGF
concentration is varied from 0.01 to 1 nM, the averaged maximum pAkt and pERK range from
1.610
3
nM to 1.8 10
3
nM and 5.0 10
2
nM to 1.0 10
3
nM, respectively. Similarly, for the same
concentration range of 0.01 – 1 nM, VEGF induces an averaged maximum pAkt of 3.5 10
1
nM to
7.610
2
nM pAkt and maximum pERK of 8.8 10
-5
nM to 4.3 10
0
nM (Figure 3-4B). Thus, for
FGF mono-stimulation, the maximum pAkt level induced by low FGF concentration (0.01 nM) is
even higher than the maximum pERK level activated by a high concentration of FGF (1 nM). The
same holds true for VEGF mono-stimulation – a low concentration of VEGF (0.01 nM) induces a
higher maximum pAkt than the maximum pERK stimulated by high VEGF concentration (1 nM).
We then studied the effects of co-stimulation of FGF and VEGF in inducing maximum
pERK and pAkt. The dynamics of pERK and pAkt stimulated by 0.5 nM FGF and 0.5 nM VEGF
in combination are shown in Figure 3-5A and B, respectively. Similar to mono-stimulation, the
model also predicts transient activation of ERK and Akt by the co-stimulation of FGF and VEGF.
The dynamics of pERK and pAkt are predicted to reach their maximum level within 30 minutes
and return to basal level after 60 minutes. We also predicted the averaged maximum pERK and
pAkt induced by co-stimulation of FGF and VEGF in a range of 0.01 – 1 nM across the 15 best

65
fits (Figure 3-5C-D). Maximum pAkt shows a greater response at all combinations of FGF and
VEGF stimulation compared to the maximum pERK induced by the same combinations (Figure
3-5C-D).
Studying the model mechanistically, we found that the main reason that Akt showed a
higher level of activation compared to ERK is because the initial concentration of PP2A, the
phosphatase that acts on pAkt (2.5 nM), is much lower than the initial concentrations of the
phosphatases that act on pMEK and pERK, Ptase 2 and Ptase 3 (3.7 10
2
and 1.7 10
3
nM,
respectively). In the model, PP2A can be produced by PP2Aoff (see reactions R87 and R88 in
Additional file 3-1: Table S1); however, the initial amount of PP2Aoff (1.1 10
2
nM) is not high
enough to make the PP2A level comparable to that of Ptase2 and Ptase3.
To confirm the effect of phosphatases, we decreased Ptase2 and Ptase3 levels to be 2.5 10
0

nM, which is the same level as PP2A. The model predicts that when Ptase2 and Ptase3 levels are
decreased, stimulation by FGF, VEGF, or their combination in the range of 0.01 – 1 nM induced
5.610
2
– 8.310
2
, 3.4 10
1
– 7.6 10
2
, and 5.6 10
2
– 9.6 10
2
nM for maximum pAkt, respectively.
In comparison, the maximum pERK is predicted to be 3.5 10
3
, 9.5 10
1
– 3.5 10
3
, and 3.5 10
3

nM, respectively, for these three cases, as pERK saturated to reach its maximum level such that
all of the ERK initially present (3.5 10
3
nM) was phosphorylated. Thus, the predicted maximum
pERK level surpassed maximum pAkt when Ptase2 and Ptase3 levels decreased, confirming that
the relatively high levels of Ptase 2 and Ptase3 limit the ERK phosphorylation compared with Akt.

66

Figure 3-5. Predicted maximum pERK and pAkt responses with co-stimulation. Predicted
time courses of pERK (A) and pAkt (B) stimulated by the combination of 0.5 nM FGF and 0.5 nM
VEGF. Curves are the mean predictions for the 15 best fits from model training. Shaded regions
show standard deviation of the fits. Maximum pERK (C) and pAkt (D) in response to co-
stimulation by FGF and VEGF for concentrations varying from 0.01 nM to 1 nM

3.4.4. ERK activation is more responsive to changing the ligand concentration compared to
Akt
We compared pERK and pAkt behaviors in terms of their responsiveness (i.e., sensitivity)
to changes in FGF and VEGF concentrations. The fold change of pAkt for VEGF stimulation of
0.01 nM compared to 1 nM for the baseline model is predicted to be 21.7, while the fold change

67
of pERK is 4.9 10
4
. The fold changes for pAkt and pERK for FGF stimulation at 0.01 nM
compared to 0.4 nM are 1.2 and 2.0. In addition, there appears to be an optimal ligand
concentration required to attain the maximum response for FGF-induced pERK and pAkt, as their
dose-response curves plateau at approximately 0.4 nM and 0.1 nM, respectively, for stimulation
with a ligand concentration in the range of 0.01 nM – 1 nM FGF (Figure 3-4B). The maximum
level of pAkt also plateaus at 0.4 nM VEGF. Before the response saturates, maximum pERK shows
a steeper increase than maximum pAkt for increasing levels of either FGF or VEGF (Figure 3-4).
In addition, maximum pERK following stimulation by VEGF continues to increase for the
concentration range simulated here (Figure 3-4B). Altogether, these results indicate that pERK is
more responsive to changing the ligand concentrations, as compared to pAkt in the range of 0.01
nM to the saturation concentration.
We next studied the effects of FGF and VEGF in combination on pERK and pAkt. The
maximum pERK attained upon co-stimulation shows an increase with increasing FGF and VEGF
(Figure 3-5C). In contrast, maximum pAkt induced by co-stimulation with FGF and VEGF is
approximately the same for all of the combinations examined (Figure 3-5D). Thus, the model
suggests that with FGF and VEGF co-stimulation, phosphorylation of ERK is more dose-
dependent compared to Akt activation, similar to modeling predictions for mono-stimulation.
Given its mechanistic detail, we can use the model to explain the reason for the greater
sensitivity of maximum pERK compared to pAkt. We found that the reason Akt activation is less
responsive to varying ligand concentrations is due to the interaction between pAkt and the
phosphatase PP2A, both in terms of the binding rates and negative feedback. Regarding the
binding rates between activated (phosphorylated) species and their phosphatases, we compare
doubly phosphorylated Akt and doubly phosphorylated MEK, as these species are each in the first

68
layer of phosphorylation reactions that lead to formation of doubly phosphorylated MEK and Akt
following FGF or VEGF binding to their corresponding receptors. This layered structure of the
phosphorylation reactions enables nonlinear signaling amplification. The association rate of ppAkt
and PP2A (k_aPP2A = 1.0 10
-3
1/nM/s) is two orders of magnitude lower than the association rate
of ppMEK and the phosphatase Ptase2 (k_dpMEK_pp = 1.4 10
-1
1/nM/s). This relatively low
association rate of ppAkt and PP2A leads to accumulation of both pAkt and ppAkt, since the
phosphorylated species bind slowly to the phosphatase. The accumulation of the phosphorylated
Akt due to this lower association rate more strongly influences the levels of phosphorylated Akt
compared to the effect of increasing the concentrations of the ligands. Therefore, the total pAkt is
relatively stable in response to stimulation. In comparison, since the association rate of
phosphorylated MEK and Ptase2 is faster, less phosphorylated MEK accumulates, and the system
remains responsive to increases in the ligand concentrations.
Another contributor to the differences in activation of MAPK versus the PI3K/Akt pathway
is negative feedback induced by phosphorylated Akt. Production of PP2A, the phosphatase that
acts on pAkt and ppAkt, is promoted by ppAkt itself (38, 91) (see reactions R87 and R88 in
Additional file 3-1: Table S1). This feedback highly regulates pAkt and ppAkt levels such that
when more ppAkt is produced, the level of PP2A also increases. The presence of this negative
feedback loop makes phosphorylated Akt levels less responsive to increased ligand concentration.
In contrast, the amounts of the phosphatases Ptase2 (which dephosphorylates MEK) and Ptase3
(which dephosphorylates ERK) do not depend on upstream species. Thus, varying the amount of
FGF or VEGF does not affect the phosphatases’ concentrations.
We confirmed the influence of ppAkt-PP2A binding and PP2A negative feedback by
systematically altering the network for Akt activation to be more like the network for ERK

69
activation (Table 3-1). We first increased the association rate of doubly phosphorylated Akt and
PP2A to be the same as the rate for doubly phosphorylated MEK and Ptase2 binding. For this case,
the fold change of ppAkt comparing two levels of VEGF stimulation (1 nM versus 0.01 nM) is
6.010
3
. When the network is modified such that PP2A is not activated by doubly phosphorylated
Akt, the fold change is predicted to be 3.4 10
3
. If the ppAkt-PP2A association rate is increased
and negative feedback is removed, the fold change in ppAkt is 1.0 10
5
. Thus, in all three of these
cases where the network is modified, the fold change for ppAkt is much higher than the baseline
case. The fold change of ppAkt in response to two levels of FGF stimulation (0.4 nM versus 0.01
nM) is also predicted to be higher than the baseline case when the association rate of doubly
phosphorylated Akt and PP2A is increased (Table 3-1). However, when the negative feedback is
removed, Akt is used up, even by the stimulation of 0.01 nM FGF. Thus, the fold change of ppAkt
in response to FGF stimulation is even lower than the baseline model because of the shortage of
Akt. These simulations confirm that the unresponsiveness of pAkt to changes in ligand
concentration is due to the particular properties of the Akt activation pathway. Overall, the model
predictions and analysis provide quantitative insight that helps to better understand how the
signaling response can be modulated to achieve a desired effect.

70

Table 3-1. Fold change for pAkt and pMEK in response to varying ligand concentration.

3.4.5. The co-stimulation by FGF and VEGF has a greater impact on phosphorylation of
ERK compared to summation of the ligands’ individual effects
To explore how FGF and VEGF influence pERK and pAkt responses together, we
compared the combination effects to the summation of individual effects in inducing maximum
pERK and pAkt. The dynamics of pERK and pAkt stimulated by 0.5 nM FGF and 0.5 nM VEGF
in combination (solid lines) and in summation (dashed lines) are shown in Additional file 3-2:
Figure S6. We observed greater pERK levels induced by the co-stimulation of 0.5 nM FGF and
0.5 nM VEGF compared to the summation of individual effects within before the activation gets
attenuated. On the other hand, the summation of 0.5 nM FGF and 0.5 nM VEGF induced pAkt is
greater than the co-stimulation at all simulated time. To more concisely represent the signaling
response induced by the growth factors, we compared the maximum pERK and pAkt induced by
FGF and VEGF co-stimulation to the summation of individual effects. We defined this ratio as
R(response), where response stands for pERK or pAkt (see Methods for more details).

71
Figure 3-6A shows that R(pERK) is greater than one for combinations of FGF and VEGF
ranging from 0.01 nM to 1 nM, specifically R(pERK) ranges from 1.01 to 1.44, and 53% of
combinations that we simulated have R(pERK) greater than 1.25. This indicates that the
combination effect of inducing maximum pERK is greater than the summation of individual
effects, which is consistent with our previous work (35). Our simulations show that 0.01 nM FGF
induces a maximum of 5.0´10
2
nM pERK, while 1 nM VEGF produces 4.3 nM pERK. Given this
100-fold difference, we believe a 25% increase in maximum pERK by the co-stimulation
compared to FGF stimulation alone is significant. As explained in 3.2, the reason VEGF is not as
potent as FGF in inducing maximum pERK could be due to the low VEGF receptor level and high
trafficking parameters, compared to FGF receptors (35).
In contrast, R(pAkt) for combinations of FGF and VEGF concentrations ranging from 0.01
nM to 1 nM is less than one (Figure 3-6B). This suggests that the combination of FGF and VEGF
is not as effective in inducing maximum pAkt, as compared to the summation of the responses
induced by each ligand individually.
We again applied the model to explain these predicted behaviors. We find that the values
of R(pERK) are greater than one because the co-stimulation compensates for limitations observed
when only one ligand is applied (35). Because there is abundant Raf available in comparison to
limited FRS2 level, VEGF co-stimulation helps to overcome the stoichiometric limitations of
FRS2 for FGF mono-stimulation. Also, because pFRS2 phosphorylates MEK faster than aRaf,
FGF co-stimulation provides a high level of pMEK. Therefore, FGF and VEGF co-stimulation
exhibits a greater effect in phosphorylating ERK than the summation of individual effects.
The reason why the R(pAkt) values are less than one is due to features of the network.
Specifically, phosphorylated Akt is more stable in response to the change of ligand concentration

72
as a result of the low association rate for ppAkt and the phosphatase PP2A and the negative
feedback of ppAkt promoting the production of its own phosphatase, PP2A. As explained in the
previous section, with mono-stimulation, the pAkt level remains relatively constant, even as the
ligands’ concentrations change. Additionally, the pAkt levels induced by combinations of FGF
and VEGF are approximately the same as the pAkt levels induced by FGF alone (Figure 3-4B and
Figure 3-5B). Thus, summing the individual effects to include the VEGF-induced maximum pAkt
level means that the denominator in the ratio R is greater than the numerator, forcing R to be less
than one.
We note that mono- and co-stimulation of FGF and VEGF affects not only the magnitude
of pERK and pAkt levels, but also the time required to reach the maximum responses and the
duration of the responses. The timescale of the pERK response was a major focus of our previous
work (35). Briefly, we found that the combination of FGF and VEGF exhibits a fast and sustained
pERK response compared to mono-stimulation (35).

Figure 3-6. Comparison of mono- and co-stimulation. Ratios, R, comparing the combination
effects to the summation of individual effects in response to FGF and VEGF for maximum pERK
(A) and pAkt (B)

73

3.4.6. The model identifies potential targets for influencing ERK and Akt activation and
evaluates their efficacy
We applied the model to determine the parameters and initial concentrations that
significantly influence pERK and pAkt levels, providing insight for researchers seeking to
effectively modulate the MAPK and PI3K signaling pathways required for angiogenesis. We
identified the influential model variables by evaluating the Stimax calculated in the eFAST global
sensitivity analysis (Additional file 3-2: Figure S7 and Additional file 3-1: Table 2) (see Methods
for more details). As a rule of thumb, we consider variables to not be influential if their Stimax
values are lower than 0.5, and the variables that have Stimax values greater than 0.7 are taken as
influential. The influential variables (Table 3-2) could be potential targets for pro- or anti-
angiogenic strategies. Interestingly, this analysis shows that species’ concentrations and kinetic
parameters of the upstream signaling network are strong regulators for both pERK and pAkt levels
and are shown to have high Stimax values for both pERK and pAkt levels. This includes: initial
concentrations of HSGAG, VEGFR2, and FGFR; kinetic parameters kf5a, k_pR2, and kf0, which
are involved in ligand receptor binding reactions; as well as the parameter k_1PI3K, which is the
association rate of pR2 and PI3K and plays an important role in the competition between the two
pathways. Not surprisingly, species and kinetic parameters involved in intermediate or
downstream signaling leading to ERK are influential and specific modulators of pERK. Similarly,
we identify model variables that specifically influence Akt. For instance, Stimax values of the MEK-
Raf dissociation rate (kd_aMEKRaf), pMEK phosphorylation rate mediated by pFRS2 (kf37) and
the ERK level are shown be greater than 0.7 for pERK but less than 0.1 for pAkt; while Akt and
PI3K levels that are involved in PI3K/Akt pathway are only influential to pAkt.

74
The eFAST analysis tells which model variables the variances in the predicted pERK and
pAkt can be attributed to. However, it is also important to determine how finitely changing those
influential variables affects the output. That is, we aim to understand i) whether the influential
variables promote or inhibit ERK and Akt activation and ii) how much the pERK and pAkt levels
change when the influential variables are changed. Therefore, the ratios of the maximum pERK
and pAkt levels compared to the baseline values were predicted when the parameters in Table 3-2
were varied by 0.1- and 10-fold (Figure 3-7). This ratio is defined as the fold change, F, in the
response (see Methods for more details). We consider the parameters that cause log2(F) to be
greater than |1| (i.e., a two-fold change) as effective targets.

75
Table 3-2. The total sensitivity index Sti values.

We calculated log2(F) for the variables in Table 3-2, identifying effective targets that
modulate pERK and pAkt upon stimulation with 0.5 nM FGF or 0.5 nM VEGF, or with co-
stimulation with 0.5 nM FGF and 0.5 nM VEGF. These predictions complement in vitro studies
that focus on the responses induced by angiogenic agents. The model predicts that increasing ERK
and MEK levels can strongly promote FGF- and VEGF-induced pERK (Figure 3-7A and C), and
increasing Akt promotes FGF- and VEGF-induced pAkt (Figure 3-7B and D). Similarly,
decreasing ERK, MEK and Akt can effectively inhibit pERK and pAkt. These predicted targets
are intuitive, as they are directly related to the signaling species of interest.
Excitingly, the model predicts several other targets. For example, increasing the
phosphatase Ptase2, which dephosphorylates pMEK and ppMEK, and decreasing FRS2, which is
an upstream FGF-mediated signaling species, significantly inhibits FGF-induced pERK (Figure 3-
7A). Also, decreasing PIP2, the substrate for producing PIP3, which further phosphorylates Akt
and pAkt, is another effective means of inhibiting FGF-induced Akt phosphorylation (Figure 3-
7B). In addition, our model predicts that the initial concentrations of VEGFR2 and Ras-GDP
positively regulate VEGF-induced ERK phosphorylation (Figure 3-7C), while increasing the
VEGFR2 level and the PIP3 activation rate (k_fPIP3) are effective strategies to promote VEGF-
induce pAkt (Figure 3-7D). Also, increasing the concentrations of phosphatases Ptase1 and Ptase2,
which deactivate Raf and MEK, respectively, can inhibit VEGF-induced pERK, as they are
negative regulators of ERK phosphorylation (Figure 3-7C). Increasing PTEN (the phosphatase for
PIP3) and the Ras-GDP level can inhibit VEGF-induced pAkt as well (Figure 3-7D). It is
noteworthy that the model predicts that the VEGFR2 level and k_pR2 positively regulate both

76
VEGF-induced pERK and pAkt. Thus, VEGFR2 and k_pR2 are targets for promoting or inhibiting
the activation of both pathways mediated by VEGF.
Interestingly, the model predicts strategies to explore the effect on ERK and Akt activation
individually. There are three potential targets that have opposing effects for the two pathways and
can be utilized to enhance the signal transduction for one pathway and dampen the response of the
other pathway. Also, these targets could play a role in the mechanism of resistance in which
inhibiting one pathway enables greater activation of the other. Specifically, decreasing the PI3K
level and the pVEGFR2 and PI3K association rate (k_1PI3K) can enhance VEGF-induced pERK,
but reduces VEGF-induced pAkt. This opposing effect is because decreasing PI3K level and
k_1PI3K can reduce the signal transduction for PI3K/Akt pathway, however decreasing PI3K level
means there is relatively less PI3K (upstream species for PI3K/Akt pathway) competing against
Ras-GDP (upstream species in the MAPK pathway) for pVEGFR2 induced by VEGF (Figure 3-
1). Also, decreasing pVEGFR2 and PI3K association rate (k_1PI3K) reduces the competition of
PI3K/Akt pathway for pVEGFR2. Therefore, there is relatively more pVEGFR2 utilized for
activating the MAPK pathway when the competition of PI3K/Akt pathway is reduced, and this
further leads to an elevated pERK level induced by VEGF. In addition, increasing Ras-GDP level
promotes VEGF-induced pERK but inhibits VEGF-induced pAkt. This is also caused by the
competition between two pathways. Increased Ras-GDP level consumes more pVEGFR2, which
limits PI3K activation by pVEGFR2 and further reduces pAkt level induced by VEGF.
Our model predicts other potential effective pro- and anti-angiogenic strategies (Figure 3-
7); however, it is of interest to investigate the effects of the crosstalk. Since FGF and VEGF are
typically both present in physiological or pathological conditions, it is relevant to identify variables
that affect activation of ERK and Akt upon co-stimulation with FGF and VEGF (Figure 3-7E-F).

77
The model shows that an increase in ERK, MEK, and VEGFR2 levels, as well as increasing the
kd_RasGDP rate promotes ERK phosphorylation. Decreasing the Ptase1 level also enhances
pERK. In addition, increasing the initial level of Akt is the most effective pro-angiogenic strategy
to enhance Akt phosphorylation. On the other hand, increasing Ptase2, and decreasing the MEK,
ERK, and FRS2 levels inhibit pERK. Lastly, decreasing the Akt, PI3K and PIP2 levels is effective
anti-angiogenic strategies to inhibit pAkt.
Finally, we compared the effect of the identified effective potential targets under different
treatments (FGF-, VEGF-, and FGF/VEGF-stimulation). Interestingly, we found that some
potential targets predicted to have an effect in response to mono-stimulation had only limited
effects in response to co-stimulation. For instance, increasing the initial level of Ptase1 was
predicted to effectively inhibit VEGF-induced pERK (Figure 3-7C); however, increasing Ptase1
leads to a log2(F) value of only -0.3 (an 0.8-fold change) when we simulate FGF and VEGF co-
stimulation. This implies that FGF-mediated signaling can diminish the inhibitory effect of
increasing Ptase1 that occurs with VEGF-induced pERK, and further illustrates the effect of
compensatory pathways in the overall results. Our simulations show that it is critical to study the
network systematically to identify potential effective targets for specific conditions.

78

Figure 3-7. Predicted targets for modulating pERK and pAkt responses. Log2(F) for 0.5 nM
FGF-induced pERK (A) and pAkt (B); 0.5 nM VEGF-induced pERK (C) and pAkt (D); and
combination of 0.5 nM FGF- and 0.5 nM VEGF-induced pERK (E) and pAkt (F). x-axes are
log2(F), y-axes are variables from Table 3-1. Bars are mean ± standard deviation of model
predictions

3.5. Discussion
We developed an intracellular signaling model of the crosstalk between two pro-angiogenic
factors, FGF and VEGF. The model focuses on pERK and pAkt responses as indicators for
signaling promoted by the two pro-angiogenic factors. In this study, we built on our previous
modeling work and incorporated PI3K/Akt pathway to get a more comprehensive understanding
of the angiogenesis process, as the PI3K/Akt pathway is important in regulating cell survival (17-
21) and migration (21, 24, 25).
Figure 7
A
B
C
D
E
F

79
Excitingly, the primary model predictions are supported by experimental students,
providing confidence that the model can be used to examine novel aspects of FGF- and VEGF-
mediated signaling. Our model predicts that the maximal levels of FGF- and VEGF-induced pERK
and pAkt plateau as the ligand concentration increases. This prediction is supported by
experimental studies showing an optimal concentration for FGF- and VEGF-induced human
umbilical vein cells (HUVECs) tube formation on Matrigel for 24 hours (0.1 ng/ml (0.004 nM)
and 25 ng/ml (0.56 nM), respectively), also 0.1 ng/ml (0.004 nM) FGF exhibits approximately
same level of increase in HUVECs proliferation and migration as 25 ng/ml (0.56 nM) VEGF
stimulation (50). In addition, the model predicts that the combination of FGF and VEGF
stimulation induces ERK phosphorylation to a greater extent than the sum of the individual effects
of FGF and VEGF. In contrast, the combination of FGF and VEGF does not promote enhanced
Akt phosphorylation compared to the summation of the response stimulated by FGF and VEGF
individually. These predictions are consistent with experimental observations.
Researchers have shown that endothelial sprouting is FGF and VEGF dose dependent (9,
50), and that the combination of FGF and VEGF induces greater total sprout length than
summation of individual effects (9). Goto et al. also demonstrated a synergistic effect on
endothelial cell proliferation upon co-stimulation by FGF and VEGF (10). In addition, it has been
reported that FGF and VEGF have significantly greater effects in combination, compared to their
individual effects in angiogenesis in vivo (11). Specifically, the systolic pressure ratio of ischemic
limb to healthy limb, the stem artery diameter, as well as the capillary density of New Zealand
White rabbits treated with FGF and VEGF in combination were significantly greater than FGF or
VEGF treated alone (11). These results are consistent with the model prediction that R(pERK) is
greater than one, where pERK is expected to directly influence proliferation. Moreover, Ratajska

80
et al. showed that co-stimulation with FGF and VEGF did not have synergistic effect on migration
distance in E12 embryonic hearts (92), which is consistent with our model prediction that R(pAkt)
is less than one for all combinations of FGF and VEGF simulated, assuming pAkt directly
influences migration. The model predictions for the pERK and pAkt responses following
stimulation by FGF and VEGF mirrors these experimental observations, providing confidence in
the model and its utility.
The molecular-detailed model presented here can be applied in various ways. We can use
the model to increase understanding of the FGF- and VEGF-mediated angiogenic mechanisms and
provide quantitative insight regarding the downstream signaling that mediates a cell’s response.
As such, our work complements models that predict cellular behavior. Norton and Popel
constructed a computational model to study vessel growth in tumor and showed that the
proliferation rate has a greater effect on the spread and extent of vascular growth compared to the
migration rate (42). The simulations from our model are in line with their results, as pERK is more
responsive to changing the ligand concentration (from 0.01 nM to the saturation concentration)
compared to pAkt. And our model provides a detailed mechanistic explanation regarding their
model predictions. Thus, the model can be utilized in combination with other modeling
frameworks that predict cellular behaviors but do not yet take intracellular signaling into account
(30, 31).
This model can also be used to study the efficiency of pro- or anti-angiogenic therapies.
Currently, there are inhibitors targeting the ERK and Akt signaling networks, such as LY294002
and wortmannin (PI3K inhibitors), and PD98059 (MEK inhibitor) (93). These inhibitors reduce
pAkt and pERK levels (93) and further inhibit endothelial migration (94, 95) and proliferation
(96), respectively. These inhibitors also reduce overall tube formation (93). Interestingly, Hoeflich

81
and coworkers showed that the MEK inhibitor PD0325901 upregulates the PI3K pathway
signaling (97). Our model is consistent with this observation and shows that inhibiting MEK
significantly reduces pERK induced by FGF or the combination of FGF and VEGF, but actually
increases the pAkt response (Figure 3-7). The model predicts other instances where targeting
certain parameters leads to opposing effects on pERK and pAkt. For example, decreasing the PI3K
level and the rate of VEGFR2-induced PI3K phosphorylation or increasing Ras-GDP level can
inhibit VEGF-induced pAkt, but promote VEGF-induced pERK (Figure 3-7 C, D). Overall, our
model can predict the important variables that influence pERK and pAkt and how the
concentrations of these signaling species are affected. These predictions can supplement
experimental studies and provide insight into investigating the efficiency of targeting particular
variables as pro- or anti-angiogenic strategies.
We acknowledge some limitations in our model. Firstly, some assumptions were made
during model construction. We simplified certain reactions that occur upstream of activating
MEK/ERK and PI3K/Akt pathways in order to focus on the effects of FGF and VEGF and their
interactions. It has been reported that the PLC𝛾 activation via VEGFR2 and FGFR further leads
to PKC activation (4, 98). However, the molecular detail relating PKC to ERK signaling is not
clear. For example, some studies show that PKC may activate Ras and trigger Raf-MEK-ERK
signaling (99-101), while PKC has been showed to activate ERK via Raf (102-104). Another study
reported that PKCα may activate MEK, independently of Raf and Ras, to further activate ERK
signaling (105). Therefore, we simplified certain reactions, and we can incorporate this detail when
the protein-protein interactions in this pathway becomes available. Also, we excluded VEGFR1
and neuropilin-1 (NRP1) since VEGFR2 is thought to be the main receptor on endothelial cells
(106). While it has been shown that VEGFR1 promotes signal transduction (107), it is largely

82
considered to be a decoy receptor (108). In addition, NPR1 primarily acts as a coreceptor for
VEGFR1 and VEGFR2 (106). We can incorporate the contributions of VEGFR1 and NRP1 into
the model in future studies. If these VEGF binding molecules were included, the effective
concentration of VEGF would be lower and thus the magnitude of VEGF-induced certain
responses may change. Moreover, we modeled FGF-mediated activation of Akt via FRS2 in the
same way that VEGF promotes Akt activation through VEGFR2, and we used the VEGF kinetic
parameters as a starting guess for the parameter fitting. This is because the FRS2-mediated protein-
protein interactions that promote Akt signaling are not fully known, and there is a scarcity of
quantitative data for the kinetics rates of FGF-induced PI3K activation. It has been reported that
phosphorylated Akt deactivates Raf (109, 110); however, experimental and computational studies
have shown that MAPK and PI3K/Akt pathways act independently for a number of different cell
types (111, 112). Since there is lack of quantitative data for the kinetic parameters for Akt-
mediated deactivation of Raf, and those parameters were not shown to influence the main model
outputs (their Sti values are less than 0.22, see detail in Sensitivity analysis section in Methods),
we did not include this feedback in our model. In the future, we can implement this interaction as
more detailed mechanistic information becomes available. Finally, we only studied the pERK and
pAkt responses over two hours in order to understand the initial effects of FGF and VEGF
stimulation. Future work can expand our model to predict the downstream effects of this initial
signaling, which occur on a longer timescale. However, despite these limitations, our model
provides novel insights into angiogenic signaling, complements experimental studies, and is a
platform for a range of future investigations.

83
3.6. Conclusion
In conclusion, we developed a mathematical model to characterize the dynamics of pERK
and pAkt following stimulation with two main pro-angiogenic factors, FGF and VEGF. The model
quantitatively studies particular aspects of FGF and VEGF interactions network in ERK and Akt
phosphorylation, provides mechanistic insight into their signaling network, and identifies specific
potential angiogenic targets that can be altered to modulate ERK and Akt activation. The model
provides a molecularly detailed understanding of the regulation of endothelial cell angiogenesis
signaling in terms of ERK and Akt activation upon stimulation with FGF and VEGF. Thus, our
work can aid in the development of pro- and anti-angiogenic strategies that particularly target ERK
and/or Akt responses induced by FGF and VEGF.

3.7. Acknowledgements
The authors are grateful to members of the Finley research group for critical discussions.
This work is supported by the US National Science Foundation (CAREER Award 1552065 to
S.D.F.).

3.8. Additional files

Additional file 3-1: Table S1. List of model reactions. Table S2. List of species with non-zero
initial concentrations. Table S3. List of model pa- rameters. Table S4. The total sensitivity index
Sti values. Table S5. Fitted initial concentrations and parameters with adjusted bounds.
Additional file 3-2: Supplementary Figures and Legends.
Additional file 3-3: Mathematical model. MATLAB.m file containing model code.

84
4. Chapter 4
Chapter 4
FGF- and VEGF-induced endothelial sprouting

4.1. Abstract
Background: The essential role of blood vessels in delivering nutrients makes
angiogenesis important in the survival of tissues, including both physiological and pathological
conditions such as wound healing and tumor growth, respectively. Angiogenesis also provides a
route for tumor metastasis. Thus, targeting angiogenesis is a prominent strategy in both tissue
engineering and cancer treatment. However, not all approaches to promote or inhibit angiogenesis
lead to successful outcomes. Pro- and anti-angiogenic therapies primarily target pro-angiogenic
factors such as vascular endothelial growth factor (VEGF) and fibroblast growth factor (FGF) in
isolation. There is a limited understanding of how these promoters combine together to stimulate
angiogenesis. We aim to quantitatively characterize the crosstalk between VEGF- and FGF-
mediated angiogenic signaling in endothelial cells and the effects of the interactions on a cellular
level, specifically endothelial sprouting, in order to generate new mechanistic insights and identify
novel therapeutic strategies.
Methods: We constructed a hybrid agent-based mathematical model that characterizes
endothelial sprouting driven by FGF and VEGF-mediated signaling. The molecular interactions of
FGF- and VEGF-induced mitogen-activated protein kinase (MAPK) signaling and the
phosphatidylinositol 3-kinase/protein kinase B (PI3K/Akt) pathway in endothelial cells were
adapted from our previous work. To link the molecular signals with the short-term cellular

85
responses, we assumed that the endothelial cell proliferation, sprout growth, and probability of
sprouting are dependent on the maximum pAkt and pERK levels upon the stimulation of FGF and
VEGF within two hours, following Hill functions. We predicted the total sprout length (TL),
number of sprouts (NS), and average length (AL) by the mono- and co-stimulation of FGF and
VEGF. We fitted these model predictions to experimental data that describe endothelial cell
proliferation and total sprout length upon FGF and VEGF stimulation. The model was validated
against experimental data not used in model fitting.
Results: We apply the trained and validated mathematical model to characterize
endothelial cell proliferation and sprouting in response to the mono- and co-stimulation of FGF
and VEGF. The model predicts that the type and concentration of ligand, length of growth factor
stimulation, and initial number of cells are important in endothelial sprouting. Also, FGF plays a
dominant role in the combination effects in endothelial sprouting. In addition, the model suggests
that cell proliferation and sprout growth of existing sprouts are more important in the sprouting
process compared to the effect of the chance of forming a new sprout. Moreover, the ERK pathway
regulates the vessel network mainly via regulating cell proliferation and NS, while the Akt pathway
mainly affects the vessel network via regulating sprout growth. Last, the model predicts that the
strategies to modulate endothelial sprouting are context dependent, and our model can identify
potential effective pro- and anti-angiogenic targets under different conditions and study their
efficacy.
Conclusions: The model provides detailed mechanistic insight into VEGF and FGF
interactions in sprouting angiogenesis. The model predicts the combination effects of FGF and
VEGF stimulation quantitatively and provides a framework to synthesize experimental data. More

86
broadly, this model can be utilized to identify targets that influence angiogenic signaling leading
to endothelial sprouting and to study the effects of pro- and anti-angiogenic therapies.

4.2. Introduction
Angiogenesis, the formation of new blood capillaries from pre-existing blood vessels, plays
an important role in the survival of tissues, as the vessels carry blood throughout the body to
provide oxygen and nutrients required by the resident cells. Thus, researchers have been putting
effort toward targeting angiogenesis as a strategy in many contexts, including in both tissue
engineering and cancer treatment. Specifically, one challenge for the success of synthetic tissues
is to ensure sufficient transport of nutrients and gases such as oxygen to the cells. Thus, stimulating
new blood vessel formation is an important strategy for the long-term viability of engineered tissue
constructs. On the other hand, inhibiting angiogenesis is an important strategy for cancer treatment
as tumors grow by obtaining nutrients and oxygen from the blood delivered by vessels, and tumor
metastasis is facilitated by the blood circulation.
The processes of blood vessel formation, particularly for capillaries, are mostly initiated
and mediated by endothelial cells responding to the local physiological conditions (1). Many
different pro-angiogenic growth factors, such as fibroblast growth factor (FGF), vascular
endothelial growth factor (VEGF), and platelet-derived growth factor (PDGF), regulate
angiogenesis (2, 3). These factors promote different cellular processes involving endothelial cell
survival, proliferation, migration, and vessel maturation leading to new blood vessel formation (4,
5). Strategies to promote or inhibit angiogenesis focus on modulating the effects of these factors
to alter the cellular-level processes they induce, with a focus on endothelial cells.

87
However, not all approaches to promote or inhibit angiogenesis lead to successful
outcomes. For example, clinical trials have shown no effective improvement in FGF- (6) or VEGF-
induced (7) angiogenesis. Also, bevacizumab, an anti-VEGF agent, has limited effects in certain
cancer types, and it is no longer approved for the treatment of metastatic breast cancer due to its
disappointing results (8). Thus, there is a need to better understand the mechanism of angiogenesis,
specifically the molecular interactions and signaling required for new blood vessel formation and
how they affect cellular behaviors, in order to establish more effective therapeutic strategies.
In addition, there is crosstalk between intracellular signaling pathways. The overall
response in endothelial angiogenesis is dependent on the integrated signals activated by these
growth factor-mediated pathways to influence cellular decisions, such as proliferation, survival,
and migration, and further promote or inhibit vessel formation. However, the integrated effects of
more than one factor on intracellular signaling reactions at a detailed level have not been studied
enough. In pro-angiogenic strategies, some synergistic effects between FGF and VEGF in
endothelial angiogenic activities have been shown (9-11). However, the mechanism of how they
act together quantitatively on regulating molecular and cellular behaviors is still not clear. Also,
in the case of inhibiting angiogenesis, tumors often evade the effects of drugs that target a single
factor by making use of alternate compensatory pathways to activate signaling species needed for
proliferation and migration. For instance, FGFR activation may play a role in the resistance
mechanism of anti-angiogenic drugs, especially anti-VEGF treatment (12, 13). Additionally,
experiments show high levels of FGFR1 in tumors that continue to progress, even during anti-
VEGF therapy (14). Thus, it is necessary to better mechanistically understand the effects of
multiple angiogenic factors and their crosstalk in activating angiogenic signaling and further
cellular responses.

88
In this study, we are interested in the angiogenic signals required to initiate vessel growth.
Therefore, we focused on the molecular signaling and further cellular responses in the process of
endothelial sprouting promoted by FGF and VEGF, as they are particularly important in early
stages of angiogenesis, while PDGF is more important in maturing the vessels. FGF and VEGF
bind to their receptors and initiate signaling through the mitogen-activated protein kinase (MAPK)
and phosphatidylinositol 3-kinase/protein kinase B (PI3K/Akt) pathways to phosphorylate
extracellular regulated kinase (ERK) and Akt, respectively. ERK and Akt are important signaling
species in the angiogenesis process that influence cell proliferation, survival, and migration. Thus,
we aim to quantitatively investigate the combination effects of two major pro-angiogenic factors,
FGF and VEGF, on activating MAPK and PI3K/Akt signaling on a molecular level and further
promoting endothelial sprouting, on a cellular level in endothelial cells.
Given the complexity of biochemical reactions comprising angiogenesis signaling
networks, a better understanding of the dynamics of these networks quantitatively is beneficial for
current angiogenesis-based strategies in many contexts, for example, in both tissue engineering
and cancer treatment. Computational modeling serves as a powerful tool to investigate molecular
and cellular responses systematically. There are many published models that predict molecular
(27-29) and cellular (30, 31) responses mediated by angiogenic factors. However, such models are
mostly designed to predict responses upon single agent stimulation. Targeting more than one
growth factor and exploring their effects in intracellular signaling and cellular responses in detail
deserves more attention. In addition, many models that focus on specific cellular behaviors
significantly reduced the intracellular signaling network such that the output signal is simply
linearly proportional to the fraction of bound receptors (30-32). Therefore, we constructed a hybrid
agent-based model to characterize the intracellular signaling interactions of FGF and VEGF and

89
utilized downstream signals, maximum pERK and pAkt, as inputs to describe angiogenic cellular
responses in the process of endothelial spheroid sprouting (Figure 4-1). Overall, we predict the
contributions of: (1) different cellular processes, (2) pro-angiogenic factors influence sprouting,
and (3) underlying molecular signals to cell sprouting. Thus, we generate novel cellular- and
molecular-level insights related to cell sprouting. In particular, in this study, we described total
sprout length (TL), number of sprouts (NS), and average length (AL), which are the well-studied
metrics in in vitro studies, by the mono- and co-stimulation of FGF and VEGF. The model predicts
that the type and concentration of ligand, length of growth factor stimulation, and initial number
of cells are important in endothelial sprouting and FGF is a dominant factor driving the
combination effects in endothelial sprouting. Also, the model suggests that proliferation of
endothelial cells and the growth of existing sprouts are more important in sprouting process
compared to the effect of the chance of forming a new sprout. In addition, MAPK and PI3K/Akt
pathways regulate the vessel network in different ways. Moreover, the strategies to modulate
endothelial sprouting are context dependent, and our model can identify potential effective pro-
and anti-angiogenic targets under different conditions and study their efficacy.

90

Figure 4-1. Endothelial spheroid sprouting process. (A) Activated endothelial cells become tip
cells and start to migrate, and they secrete proteolytic enzymes to degrade and remodel local
extracellular matrix. The stalk cells behind tip cells are proliferative. Finally, the endothelial cells
sprout into linear cord-like structures. (B) The sprouting process involves tip cell migration, stalk
cell proliferation, and their elongation. (C) Growth factors, FGF and VEGF, binding to their
receptors initiate intracellular signaling and regulate cellular responses.

4.3. Methods
4.3.1. Model construction
We constructed a hybrid agent-based model that describes cellular responses, including
cell proliferation, sprout growth, and the formation of new sprouts. These cellular responses are
driven by molecular signals, pERK and pAkt, upon the mono- and co-stimulation of FGF and
VEGF. The molecularly-detailed biochemical reaction network that characterizes the MAPK and
PI3K/Akt pathways induced by FGF and VEGF is adapted from our previous work (39), which is
referred to as the ERK-Akt model in this study.
Proliferation
Elongation
Sprouting and migration
A
B
C
Endothelial cell
Tip cell

91
To link the molecular signals (pERK and pAkt) produced by the angiogenic factors (FGF
and VEGF) with the short-term cellular responses, we made four assumptions: 1) the endothelial
cell responses (cell proliferation, sprout growth, and probability of sprouting) are dependent on the
maximum pAkt and pERK levels upon the stimulation of FGF and VEGF within two hours; 2) the
intrinsic properties of endothelial cells to grow and form sprouts are stable within three days of
simulated cell culture; 3) daughter cells inherit all the properties from mother cells; and 4) the
maximum pERK and pAkt drive endothelial cell proliferation and sprouting following Hill
functions:
𝑓 (𝑆 )=
𝑉 𝑚𝑎𝑥 1+(
𝐾𝑚
[S]
)
𝑛
where S refers to either pERK or pAkt, [S] is their maximal concentration, Km is the substrate
concentration where the proliferation or sprouting rate is half of its maximum value, 𝑉 𝑚𝑎𝑥 , and n
is the Hill coefficient.

4.3.1.1. Module of endothelial cell proliferation
The endothelial cell proliferation module is a hybrid agent-based model that simulates each
endothelial cell as one agent that has its own cell proliferation rate and divides based on its own
cell doubling time. The pERK and pAkt levels were used to inform the rate of cell proliferation,
as has been done in other computational work (113). The total number of endothelial cells was
quantified to account for cell proliferation, which is a net result of survival and proliferation of
endothelial cells in response to FGF and/or VEGF stimulation. The ERK-Akt model predicts the
dynamics of the molecular signals pAkt and pERK within two hours upon ligand stimulation,
which are the inputs to calculate the rate of cell proliferation.

92
Based on the assumptions above, the average cell proliferation rate (rcp) is given by
𝑟 𝑐𝑝
=𝑟 𝑐 𝑝 𝑏𝑎𝑠𝑎𝑙 +𝑟 𝑐 𝑝 𝑝𝐸𝑅𝐾 +𝑟 𝑐 𝑝 𝑝𝐴𝑘𝑡
=𝑟 𝑐𝑝 _𝑏𝑎𝑠𝑎𝑙 +𝑘 𝑐𝑝 _𝑝𝐸𝑅𝐾 ×𝑓 𝑐𝑝
(max(𝑝𝐸𝑅𝐾 ))+𝑘 𝑐𝑝 _𝑝𝐴𝑘𝑡 ×𝑓 𝑐𝑝
(max(𝑝𝐴𝑘𝑡 ))
=𝑟 𝑐𝑝 _𝑏𝑎𝑠𝑎𝑙 + 𝑘 𝑐𝑝 _𝑝𝐸𝑅𝐾 ×
𝑉 max_𝑝𝐸𝑅𝐾 _𝑐𝑝
1+(
𝐾𝑚
𝑝𝐸𝑅𝐾 _𝑐𝑝
[max (𝑝𝐸𝑅𝐾 )]
)
𝑛 𝑐𝑝
+𝑘 𝑐𝑝 _𝑝𝐴𝑘𝑡 ×
𝑉 max_𝑝𝐴𝑘𝑡 _𝑐𝑝
1+(
𝐾𝑚
𝑝𝐴𝑘 𝑡 _𝑐𝑝
[max (𝑝𝐴𝑘𝑡 )]
)
𝑛 𝑐𝑝

=𝑟 𝑐𝑝 _𝑏𝑎𝑠𝑎𝑙 +
𝑉 𝑝𝐸𝑅𝐾 _𝑐𝑝
1+(
𝐾𝑚
𝑝𝐸𝑅𝐾 _𝑐𝑝
[max (𝑝𝐸𝑅𝐾 )]
)
𝑛 𝑐𝑝
+
𝑉 𝑝𝐴𝑘𝑡 _𝑐𝑝
1+(
𝐾𝑚
𝑝𝐴𝑘𝑡 _𝑐𝑝
[max (𝑝𝐴𝑘𝑡 )]
)
𝑛 𝑐𝑝

where 𝑟 𝑐𝑝 _𝑏𝑎𝑠𝑎𝑙 is the basal cell proliferation rate when no stimuli are applied, 𝑘 𝑐𝑝 _𝑝𝐸𝑅𝐾 and
𝑘 𝑐𝑝 _𝑝𝐴𝑘𝑡 are the cell proliferation rate constants, 𝑉 max_𝑝𝐸𝑅𝐾 _𝑐𝑝
and 𝑉 max_𝑝𝐴𝑘𝑡 _𝑐𝑝
are the maximum
rates of cell proliferation driven by maximum pERK and pAkt, respectively. For simplification,
we use 𝑉 𝑝𝐸𝑅𝐾 _𝑐𝑝
and 𝑉 𝑝𝐴𝑘𝑡 _𝑐𝑝
to represent 𝑘 𝑐𝑝
×𝑉 max_𝑝𝐸𝑅𝐾 _𝑐𝑝
and 𝑘 𝑐𝑝
×𝑉 max_𝑝𝐴𝑘𝑡 _𝑐𝑝
respectively.
And [max (𝑝𝐸𝑅𝐾 )] and [max (𝑝𝐴𝑘𝑡 )] are maximum pERK and pAkt levels within two hours
upon the ligand simulation, respectively. 𝐾𝑚
𝑝𝐸𝑅𝐾 _𝑐𝑝
and 𝐾𝑚
𝑝𝐴𝑘𝑡 _𝑐𝑝
are the maximum pERK and
pAkt levels that produce the half maximal cell proliferation rates, respectively. 𝑛 𝑐𝑝
is the Hill
coefficient for the cell proliferation rate.
The average cell proliferation rate (rcp) indicates that the average doubling time for
endothelial cells is 1/rcp. To account for the cell heterogeneity within a cell population, we assigned
a cell proliferation rate chosen from a normal distribution with a mean (μ) of the calculated average
cell proliferation rate and a standard deviation (σ) to capture 99.7% of the possible values given
the range of μ ± 25%μ (i.e., μ ± 3σ) for each cell.
The total cell number at a particular time T (𝑁 𝑡𝑜𝑡 (𝑇 ) ) is given by
𝑁 𝑡𝑜𝑡 (𝑇 )=∑2
𝑓𝑙𝑜𝑜𝑟 (𝑇 ×𝑟 𝑐𝑝 _𝑖 )
𝑁 𝑖𝑛𝑡 𝑖=1

93
where cell i from an initial cell population that consists of 𝑁 𝑖𝑛𝑡 of HUVECs has a cell proliferation
rate of 𝑟 𝑐𝑝 _𝑖 , T is cell culture time, 𝑓𝑙𝑜𝑜𝑟 (𝑇 ×𝑟 𝑐𝑝 _𝑖) rounds (𝑇 ×𝑟 𝑐𝑝 _𝑖) to the nearest integer less
than or equal to (𝑇 ×𝑟 𝑐𝑝 _𝑖) .

Note that the mean cell proliferation from experimental data is usually calculated for
several replicates. To compare with experimental data, we calculated the average total cell number
for ten simulations per condition.

4.3.1.2. Model of endothelial cell sprouting
Model overview. The endothelial cell sprouting model is also a hybrid agent-based model
that simulates each endothelial cell as one agent that has its own properties and makes its own
cellular decisions. The model utilizes a probabilistic approach to model sprouting. The general
flow of the model is shown in Figure 4-2.

Figure 4-2. Flowchart of the endothelial sprouting agent-based model. The model simulates
each endothelial cell as one agent that has its own properties and makes its own cellular decisions
in every time step.

94

First the model checks each cell to see whether it is a tip cell (i.e., if it is a leading cell in a
trail of cells) and if so, the cell migrates, leading to sprout elongation with an assigned sprout
growth rate. If the cell is not a tip cell, it has a chance to become a tip cell based on a certain
probability: the model generates a random number, and if the random number is greater than the
given probability threshold, then this cell becomes a tip cell and starts migrating instead of
proliferating. If the random number is not greater than the threshold for the probability of becoming
a tip cell, the model then checks whether at this time, the cell is ready to proliferate based on
whether the simulated time has reached the assigned cell doubling time. The cell’s doubling time
is defined by the inverse of its cell proliferation rate (refer to above section on the cell proliferation
module). If the check is yes, the cell divides to generate a daughter cell that is assumed to inherit
all the properties from the mother cell. If the elapsed simulation time is not enough for the cell to
proliferate, it remains quiescent. The model repeats this process at the next time point until the end
of the simulation. We update results every hour, which corresponds to the time scale over which
cellular responses are usually studied in vitro (50, 114, 115).

Initial state. Initially, the model consists of a specific number of endothelial cells (Nint) in
a spheroid. Because we only focus on the number of sprouts and sprout lengths, which are well-
studied metrics in in vitro sprouting assays, we do not consider spatial effects or sprouting
directions in this study.
Each of the initial cells are assigned a cell proliferation rate (refer to above section on the
cell proliferation module), sprout growth rate, and a probability of forming a new sprout. The
average values of the sprout growth rate and a probability of forming a new sprout for a population

95
of the endothelial cells were calculated, as described below in detail. The sprout growth rate and a
probability of forming a new sprout for each cell were chosen from normal distributions, where
the mean and standard deviation for each is based on from the calculated values for the population.
The mean (μ) is taken as the calculated average probability of forming a new sprout and average
sprout growth rate and a standard deviation (σ) is set to capture 99.7% of the possible values given
the range of μ ± 100%μ and μ ± 25%μ (i.e., μ ± 3σ) for each cell, respectively.
During the first iteration, each of the initial cells is assigned a cell proliferation rate, a
sprout growth rate, and a probability to become a tip cell, and then all decisions follow the
flowchart shown in Figure 4-2. Each cell can only become a tip cell once, and it stops proliferating
if it becomes a tip cell.

Sprouting. Endothelial sprouting is dependent on cell proliferation, migration, and
elongation. To quantify endothelial sprouting, we consider the formation of a new sprout and the
growth of existing sprouts, which are determined by the probability of forming a new sprout (p),
and sprout growth rate (rsg), respectively. We note that the sprout growth rate here is the net rate
of sprout elongation caused by cell proliferation, migration and elongation.
Based on the assumptions made above, the average probability of forming a new sprout 𝑝
from a cell during a time period ∆t (∆t = 1 hour) is
p=𝑝 𝑏𝑎𝑠𝑎𝑙 +𝑝 𝑝𝐸𝑅𝐾 +𝑝 𝑝𝐴𝑘𝑡
=𝑝 𝑏𝑎𝑠𝑎𝑙 + 𝑘 𝑝 _𝑝𝐸𝑅𝐾 ×𝑓 𝑝 (max(𝑝𝐸𝑅𝐾 ))+𝑘 𝑝 _𝑝𝐴𝑘𝑡 ×𝑓 𝑝 (max(𝑝𝐴𝑘𝑡 ))
=𝑝 𝑏𝑎𝑠𝑎𝑙 +𝑘 𝑝 _𝑝𝐸𝑅𝐾 ×
𝑉 max_𝑝𝐸𝑅𝐾 _𝑝 1+(
𝐾𝑚
𝑝𝐸𝑅𝐾 _𝑝 [max (𝑝𝐸𝑅𝐾 )]
)
𝑛 𝑝 +𝑘 𝑝 _𝑝𝐴𝑘𝑡 ×
𝑉 max_𝑝𝐴𝑘𝑡 _𝑝 1+(
𝐾𝑚
𝑝𝐴𝑘𝑡 _𝑝 [max (𝑝𝐴𝑘𝑡 )]
)
𝑛 𝑝

96
=𝑝 𝑏𝑎𝑠𝑎𝑙 +
𝑉 𝑝𝐸𝑅𝐾 _𝑝 1+(
𝐾𝑚
𝑝𝐸𝑅𝐾 _𝑝 [max (𝑝𝐸𝑅𝐾 )]
)
𝑛 𝑝 +
𝑉 𝑝𝐴𝑘𝑡 _𝑝 1+(
𝐾𝑚
𝑝𝐴𝑘𝑡 _𝑝 [max (𝑝𝐴𝑘𝑡 )]
)
𝑛 𝑝
where 𝑝 𝑏𝑎𝑠𝑎𝑙 is the basal probability rate of forming a new sprout when no stimuli are applied,
𝑘 𝑝 _𝑝𝐸𝑅𝐾 and 𝑘 𝑝 _𝑝𝐴𝑘𝑡 are the probability constants, 𝑉 max_𝑝𝐸𝑅𝐾 _𝑝 and 𝑉 max_𝑝𝐴𝑘𝑡 _𝑝 are the maximum
probability rates of sprout formation driven by maximum pERK and pAkt in an hour, respectively.
For simplification, we use 𝑉 𝑝𝐸𝑅𝐾 _𝑝 and 𝑉 𝑝𝐴𝑘𝑡 _𝑝 to represent 𝑘 𝑝 _𝑝𝐸𝑅𝐾 ×𝑉 max_𝑝𝐸𝑅𝐾 _𝑝 and
𝑘 𝑝 _𝑝𝐴𝑘𝑡 ×𝑉 max_𝑝𝐴𝑘𝑡 _𝑝 respectively. 𝐾𝑚
𝑝𝐸𝑅𝐾 _𝑝 and 𝐾𝑚
𝑝𝐴𝑘𝑡 _𝑝 are the maximum pERK and pAkt
levels that produce the half maximal probability rates of sprout formation, respectively. 𝑛 𝑝 is the
Hill coefficient for the probability rate of forming a new sprout.
The average sprout growth rate (rsg) is
𝑟 𝑠𝑔
=𝑟 𝑠𝑔 _𝑏𝑎𝑠𝑎𝑙 +𝑟 𝑠𝑔 _𝑝𝐸𝑅𝐾 +𝑟 𝑠𝑔 _𝑝𝐴𝑘𝑡
=𝑟 𝑠𝑔 _𝑏𝑎𝑠𝑎𝑙 +𝑘 𝑠𝑔 _𝑝𝐸𝑅𝐾 ×𝑓 𝑠𝑔
(max(𝑝𝐸𝑅𝐾 ))+𝑘 𝑠𝑔 _𝑝𝐴𝑘𝑡 ×𝑓 𝑠𝑔
(max(𝑝𝐴𝑘𝑡 ))
=𝑟 𝑠𝑔 _𝑏𝑎𝑠𝑎𝑙 +𝑘 𝑠𝑔 _𝑝𝐸𝑅𝐾 ×
𝑉 max_𝑝𝐸𝑅𝐾 _𝑠𝑔
1+(
𝐾𝑚
𝑝𝐸𝑅𝐾 _𝑠𝑔
[max (𝑝 𝐸 𝑅𝐾 )]
)
𝑛 𝑠𝑔
+𝑘 𝑠𝑔 _𝑝𝐴𝑘𝑡 ×
𝑉 max_𝑝𝐴𝑘𝑡 _𝑠𝑔
1+(
𝐾𝑚
𝑝𝐴𝑘𝑡 _𝑠𝑔
[max (𝑝𝐴𝑘𝑡 )]
)
𝑛 𝑠𝑔

=𝑟 𝑠𝑔 _𝑏𝑎𝑠𝑎𝑙 +
𝑉 𝑝𝐸𝑅𝐾 _𝑠𝑔
1+(
𝐾𝑚
𝑝𝐸𝑅𝐾 _𝑠𝑔
[max (𝑝𝐸𝑅𝐾 )]
)
𝑛 𝑠𝑔
+
𝑉 𝑝𝐴𝑘𝑡 _𝑠𝑔
1+(
𝐾𝑚
𝑝𝐴𝑘𝑡 _𝑠𝑔
[max (𝑝𝐴𝑘𝑡 )]
)
𝑛 𝑠𝑔

where 𝑟 𝑠𝑔 _𝑏𝑎𝑠𝑎𝑙 is the basal sprout growth rate when no stimuli are applied, 𝑘 𝑠𝑔 _𝑝𝐸𝑅𝐾 and 𝑘 𝑠𝑔 _𝑝𝐴𝑘𝑡
are the sprout growth rate constants, 𝑉 max_𝑝𝐸𝑅𝐾 _𝑠𝑔
and 𝑉 max_𝑝𝐴𝑘𝑡 _𝑠𝑔
are the maximum rates of
sprout length increase driven by maximum pERK and pAkt, respectively. For simplification, we
use 𝑉 𝑝𝐸𝑅𝐾 _𝑠𝑔
and 𝑉 𝑝𝐴𝑘𝑡 _𝑠𝑔
to represent 𝑘 𝑠𝑔 _𝑝𝐸𝑅𝐾 ×𝑉 max_𝑝𝐸𝑅𝐾 _𝑠𝑔
and 𝑘 𝑠𝑔 _𝑝𝐴𝑘𝑡 ×𝑉 max_𝑝𝐴𝑘𝑡 _𝑠𝑔

respectively. 𝐾𝑚
𝑝𝐸𝑅𝐾 _𝑠𝑔
and 𝐾𝑚
𝑝𝐴𝑘𝑡 _𝑠𝑔
are the maximum pERK and pAkt levels that produce the

97
half maximal sprout growth rates, respectively. 𝑛 𝑠𝑔
is the Hill coefficient for the sprout growth
rate.
Next, we used these three parameters (rcp, p, and rsg) to characterize endothelial cell
sprouting: number of sprouts (NS), total sprout length (TL), and average sprout length (AL) in a
certain period of time. The number of sprouts is determined by counting the number of tip cells
predicted by the model. The total sprout length is the summation of all sprout lengths. The average
sprout length is calculated as the total sprout length divided by the number of sprouts.

Constraint on probability of forming a sprout
It has been reported that tip cells activate Notch signaling and prevent neighboring cells
from becoming a tip cell (116-118). Thus, we applied a constraint in our model to account for the
effects of lateral inhibition. We adapted a rate constant Pmax (5×10
−4
μm
−1
h
−1
), which determines
the maximum probability of sprout formation per unit time and vessel length in a rat corneal assay
(31). To match HUVEC data (119, 120), we adjusted the Pmax to 10
−3
μm
−1
h
−1
. Thus, Pmax is used
to limit the maximum number of sprouts at every time step.

Model outputs
The TL, NS, and AL from experimental data are usually obtained by quantifying the
sprouts on a focal plane. In order to make comparisons to experimental data, we calculated the TL,
NS, and AL on a focal plane by assuming the sprouts are uniformly distributed on a spheroid.
Thus, the TL and NS are scaled to the ratio of the number of cells on the focal plane and the total
number of cells.

98
Also, the mean values of TL, NS, and AL from experimental data are usually calculated
for several randomly selected spheroids per experimental group. To compare with experimental
data, we calculated the average TL, NS, and AL on a focal plane for ten simulations per condition.

4.3.2. Sensitivity analysis
Before fitting the model to experimental data, we first performed a sensitivity analysis to
identify the parameters that significantly influence the model outputs, using the extended Fourier
Amplitude Sensitivity Test (eFAST) (57) method. Since the initial concentrations and parameters
involved in the ERK-Akt model are adapted from previous work (39), we used the median fitted
values and held them constant during the sensitivity analysis. All remaining model parameters
were varied simultaneously within two orders of magnitude above and below the baseline values,
where the baseline values are provided in Table S4-1. In this way, the effects of multiple model
inputs on 𝑟 𝑐𝑝
, 𝑟 𝑐𝑝 _𝑝𝐸𝑅𝐾 , and 𝑟 𝑐𝑝 _𝑝𝐴𝑘𝑡 in the cell proliferation module and 𝑟 𝑠𝑔
, 𝑟 𝑠𝑔 _𝑝𝐸𝑅𝐾 , 𝑟 𝑠𝑔 _𝑝𝐴𝑘𝑡 ,
p, 𝑝 𝑝𝐸𝑅𝐾 , and 𝑝 𝑝𝐴𝑘𝑡 in the sprouting model were computed. Specifically, the eFAST method gives
the total sensitivity indices, “Sti”, which can range from 0 to 1, where a higher Sti index indicates
the input is more influential to the output. We calculated the Sti values using eFAST for all the
same ligand concentrations as the experimental data that were used for model training. The highest
Sti value (Stimax) across all of the concentrations was selected to represent the sensitivity index for
each parameter.
We also performed eFAST for the trained and validated model to identify potential targets
for pro- and anti-angiogenic strategies. All parameters and initial concentrations in the ERK-Akt
model were varied simultaneously within two orders of magnitude above and below the baseline
values. The fitted variables were held constant at the median values estimated from model fitting.

99
We calculated the Sti values to quantify how all the variables affected rates of cell proliferation,
sprout growth, and the probability of forming a new sprout. Based on the effects in influencing
TL, low, intermediate, and high levels of FGF and VEGF (Table 4-1) were selected as
representative ligand concentrations. We calculated the Sti values using eFAST for the nine
possible combinations of low, intermediate, and high levels of FGF and VEGF stimulation. Again,
the Stimax

across all the combinations were compared for all the variables.

Low (ng/ml) Intermediate (ng/ml) High (ng/ml)
FGF 0.03 0.1 10
VEGF 0.1 4 25
Table 4-1. Representative low, intermediate, high levels of FGF and VEGF.

4.3.3. Data extraction
Data from published experimental studies (119-122) were used for parameter fitting and
model validation. Experimental data from plots was extracted using the grabit function in
MATLAB.

4.3.4. Parameterization
4.3.4.1. Cell proliferation module
Fitting. The initial concentrations and parameters involved in the ERK-Akt model are
adapted from previous work (39). Five influential variables with Sti values greater than 0.5 were
identified by performing eFAST in the cell proliferation module (Figure S4-1A). They were
estimated against experimental measurements using Particle Swarm Optimization (PSO) (59) to
minimize the objective function (the difference between model predictions and experimental data).

100
PSO starts with a population of initial particles (parameter sets). As the particles move around (i.e.,
as the algorithm explores the parameter space), an objective function is evaluated at each particle
location. Particles communicate with one another to determine which has the lowest objective
function value. The objective function for each parameter set was used to identify optimal
parameter values. Specifically, we used PSO to minimize the weighted sum of squared residuals
(WSSR):
WSSR(θ)=𝑚𝑖𝑛 ∑(
𝑉 𝑝𝑟𝑒𝑑 ,𝑖 (θ)−𝑉 𝑒𝑥𝑝 ,𝑖 𝑉 𝑒𝑥𝑝 ,𝑖 )
2
𝑛 𝑖=1

where Vexp,i is the ith experimental measurement, Vpred,i is the ith predicted value at the
corresponding time point, and n is the total number of experimental data points. The minimization
is subject to θ, the set of upper and lower bounds on each of the free parameters. The bounds were
set to be two orders of magnitude above and below the baseline parameter values, which are listed
in Table S4-1.
The cell proliferation module was fitted 200 times using two datasets, specifically the
relative proliferation of HUVECs stimulated by 0.03 - 1 ng/ml FGF (121) and 0.1 - 1 ng/ml VEGF
(122) for 48 hours compared with the reference FGF and VEGF concentration points of 1 ng/ml,
respectively (Figure 4-3A-B). Note that the simulated initial number of cells are the same as
experimental data, specifically, West et al. (121) and Jih et al. (122) cultured 10
4
cells and 5000
cells initially respectively. Specifically, the relative change of the HUVEC proliferation was
calculated as following:
𝑟𝑒𝑙𝑎𝑡𝑖𝑣𝑒 𝑝𝑟𝑜𝑙𝑖𝑓𝑒𝑟𝑎𝑡𝑖𝑜𝑛 (𝑐 )=
𝑝𝑟𝑜𝑙𝑖𝑓𝑒𝑟𝑎𝑡𝑖𝑜𝑛 (𝑐 )−𝑝𝑟𝑜𝑙𝑖𝑓𝑒𝑟𝑎𝑡𝑖𝑜𝑛 (𝑐 𝑟𝑒𝑓 )
𝑝𝑟𝑜𝑙𝑖𝑓𝑒𝑟𝑎𝑡𝑖𝑜𝑛 (𝑐 𝑟𝑒𝑓 )

where 𝑝𝑟𝑜𝑙𝑖𝑓𝑒𝑟𝑎𝑡𝑖𝑜𝑛 (𝑐 ) is the HUVEC proliferation in response to ligand concentration c, and
𝑝𝑟𝑜𝑙𝑖𝑓𝑒𝑟𝑎𝑡𝑖𝑜𝑛 (𝑐 𝑟𝑒𝑓 ) is the HUVEC proliferation at a reference concentration point cref.

101

Validation. We then validated the model with two datasets not used in the fitting. We
predicted results for 1 – 30 ng/ml FGF- and 1 – 10 ng/ml VEGF-induced HUVEC relative
proliferation on day 2 (using the reference ligand concentration points of 1 ng/ml) and compared
with experimental data, specifically the relative proliferation of HUVECs stimulated by 1 – 30
ng/ml FGF (121), 0.5 – 10 ng/ml FGF (122), and 1 – 10 ng/ml VEGF (122) for 48 hours.
For all three datasets (121, 122), we simulated the experimental conditions without any
additional model fitting and compared to the experimental measurements. A total of 21 parameter
sets with the smallest errors were taken to be the “best” sets based on the model fitting and
validation (Figures 4-3A-B and S4-2A and Table S4-2) and the median values were used for the
sprouting model.

4.3.4.2. Sprouting model
Seven influential variables with Sti values greater than 0.75 were identified by performing
eFAST in the sprouting model (Figure S4-1B). Due to the lack of experimental data, we first
estimated all the unknown parameters 500 times by fitting the model to experimental observations,
specifically 0 – 64 ng/ml VEGF-induced total sprout length for 24 hours cultured with 500-cell
spheroid initially (119). After model training, we validated the model with another dataset not used
in the fitting. Specifically, we used the fold change of the average length and number of sprouts
induced by 25 ng/ml FGF and 25 ng/ml VEGF for 24 hours cultured with 400-cell spheroid initially
compared to the control (120) for validation. We simulated the experimental conditions without
any additional model fitting and compared to the experimental measurements. A total of 15
parameter sets with the smallest errors were taken to be the “best” sets based on the model fitting

102
and validation (Figures S4-2B and S4-3 and Table S4-3). The non-influential parameters were held
constant at the median of the fitted values, and the seven influential variables were estimated using
PSO 300 times using data from Heiss et al. (0 – 64 ng/ml VEGF-induced total sprout length) (119).
We again compared model predictions to Liebler et al. data (120) without any additional model
fitting. A total of 18 parameter sets with the smallest errors were taken to be the “best” sets based
on the model fitting and validation (Figures 4-3C-D and S4-2C and Table S4-4) and were used for
all model simulations presented below.

4.4. Results
4.4.1. The fitted hybrid agent-based model captures the main features of FGF- and VEGF-
induced endothelial sprouting characteristics
We developed a hybrid agent-based mathematical model that describes angiogenic cellular
responses in the process of endothelial sprouting driven by integrating molecular signals, pERK
and pAkt, upon the mono- and co-stimulation of FGF and VEGF (Figure 4-2). The model focuses
on the endothelial proliferation, sprout growth, and new sprout formation, which are assumed to
be dependent on the maximum pERK and pAkt levels following Hill functions. The model
parameter values are given in Supplementary Tables S4-1-S4-4. The molecular-detailed
biochemical reaction network that characterizes the MAPK and PI3K/Akt pathways induced by
FGF and VEGF is adapted from our previous work (39), and the parameters and initial
concentrations are taken from the median of the fitted values (39). The newly introduced
parameters were estimated by fitting the model to experimental data, as described below.

103
4.4.1.1. Cell proliferation module
To identify the influential parameters to the model outputs, rcp, rcp_pERK, and rcp_pAkt, we
performed the eFAST (57) (see Methods for more details) and analyzed the Stimax for the newly
introduced parameters (Figure S4-1A). All five parameters (Table S4-2) were identified as
influential and were estimated by fitting the model to experimental measurements (121, 122) using
PSO (59) (see Methods for more details).
The fitted model shows a good agreement with experimental results (Figure 4-3A-B). It
quantitatively captures the main features of FGF- and VEGF-induced endothelial cell proliferation
from experimental observations (121, 122). In addition to matching data used for fitting, model
predictions were compared to experimental data not used in the model fitting (121, 122) to validate
the model (Figure 4-3A-B). The model is consistent with experimental observations and can
capture the plateau behavior at high FGF and VEGF concentrations (Figure 4-3A-B). The
weighted errors for 21 best fits are all approximately 4.1 (Table S4-2). Also, the estimated values
of the fitted variables show good consistency (Figure S4-2A) and the median values were used for
the sprouting model.

4.4.1.2. Sprouting model
We again identified seven variables that are influential to rsg, rsg_pERK, and rsg_pAkt, p, p_pERK,
and p_pAkt with Sti values greater than 0.75 in the sprouting model using eFAST (Figure S4-1B).
Due to the lack of experimental data, we first estimated all the unknown parameters by fitting the
model to experimental observations showing VEGF-induced total sprout length (119). We selected
15 “best” fits that showed good match to experimental observations (120) (Figure S4-3). The non-
influential parameters were then set at the median of fitted values (Figures S4-1B and S4-2B) and

104
the influential parameters were estimated 300 times, from which 18 “best” fits were selected
(Figures 4-3C and S4-2C). The fitting results can capture the main features of VEGF-induced
endothelial total length for 24 hours (Figure 4-3C). We again compared model predictions to
independent experimental data (120) without any additional model fitting. It showed good
agreement with the fold change of the average length and number of sprouts induced by 25 ng/ml
FGF and 25 ng/ml VEGF for 24 hours compared to the control (Figure 4-3D). The weighted errors
for 18 best fits are 0.001 – 0.013 (Table S4-4). Also, the estimated values of the fitted variables
show good consistency (Figure S4-2C). The ability to predict experimental data not used for
estimating the model parameters suggests the model is reliable to make new predictions.

105
Figure 4-3. Model comparison to training and validation data for FGF or VEGF stimulation.
(A) Relative change of endothelial cell proliferation for 10
4
cells cultured for 48 hours in response
to 0.01 – 30 ng/ml FGF stimulation compared with the reference FGF concentration points of 1
ng/ml. (B) Relative change of endothelial cell proliferation for 5000 cells cultured for 48 hours by
the stimulation of 0.1 – 10 ng/ml VEGF compared with the reference VEGF concentration points
of 1 ng/ml. (C) Total sprout length induced by 0 – 64 ng/ml VEGF for 24 hours cultured with 500-
cell spheroid initially. (D) The fold change of the average length and number of sprouts induced
by 25 ng/ml FGF and 25 ng/ml VEGF for 24 hours cultured with 400-cell spheroid initially
compared to the control. Circles, squares, and diamonds in Panels A-C are experimental data (119-
122). Circles in Panel A, squares in Panels A and B, and diamonds in Panel C are experimental
data from West et al. (121), Jih et al. (122), and Heiss et al. (119), respectively. The light yellow
circles and light blue squares in Panel A-B are experimental data used for model fitting, while the
oranges circles and squares and dark blue squares are experimental data used for model validation.
Curves in Panels A-B and C are the mean model predictions of the 21 and 18 best fits, respectively.
Shaded regions show standard deviation of the fits. Solid and dashed bars in Panel D are mean ±
standard deviation of Liebler et al. data (120) and model predictions, respectively.

4.4.2. The type and concentration of ligand, length of growth factor stimulation, and initial
number of cells impact endothelial sprouting
4.4.2.1. Mono-stimulation
We first compared the effects of four inputs: type of ligand (FGF and VEGF), concentration
of the ligand (low, intermediate, and high) (Table 4-1), length of growth factor stimulation (1-3
days), and initial number of cells (250, 500, and 750 cells) on sprouting characteristics (TL, NS,
and AL) (Figure 4-4). Note that AL equals TL/NS, providing the relative change in TL compared
to NS. A higher AL suggests that the TL is due to the growth of existing sprouts, whereas a lower
AL indicates that the formation of new sprouts plays a more important role in TL. Also, we take
low, intermediate, and high levels of FGF and VEGF based on their effects in inducing TL as
representative concentrations to study the effects of ligand concentration. Note that for the high
levels of FGF and VEGF, the total sprout length is at a plateau level.

Type of ligand (FGF and VEGF)

106
Generally, FGF induces greater sprouting responses than VEGF in terms of TL (Figure 4-
4A-B) and NS (Figure 4-4C-D) at the same concentrations. On day 1, FGF shows slightly higher
sprouting responses compared to VEGF, and the differences between FGF- and VEGF-induced
TL (Figure 4A-B) and NS (Figure 4-4C-D) increase with time. This indicates that FGF-induced
TL and NS increase faster than VEGF-induced TL and NS, which results in greater sprouting
responses induced by FGF in the long-term, compared to VEGF. In addition, FGF-induced AL is
higher than VEGF for all concentrations on day 1 (Figure 4-4E-F). On day 2 and 3, FGF showed
higher effects in AL at low to intermediate concentrations (Figure 4-4E), while VEGF showed
higher effects in AL at a high concentration (Figure 4-4F). While increasing FGF and VEGF
concentration both lead to greater NS and TL, a higher FGF concentration causes a greater increase
in NS compared to TL, which causes a lower AL (Figure 4-4A, C, and E). The opposite is true for
VEGF, where a higher concentration causes a greater increase in TL compared to NS, producing
a higher AL (Figure 4-4B, D, and F). Thus, increasing concentrations of FGF and VEGF affect the
average length of sprouts in different ways.

Concentration of the ligands (low, intermediate, and high)
FGF- and VEGF-induced TL, NS, and AL are dose dependent (Figure 4-4). Specifically,
FGF- and VEGF-induced TL (Figure 4-4A-B) and NS (Figure 4-4C-D) increase with increasing
FGF or VEGF concentration. Also, VEGF-induced AL on day 1 increases with the increase in
VEGF concentration (Figure 4-4F). Moreover, FGF- and VEGF-induced AL has a biphasic dose
response on days 1-3 and days 2-3 respectively (Figure 4-4E-F). The dose response of the AL for
FGF is U-shaped (Figure 4-4E), while the dose-response for VEGF is an inverted U (Figure 4-4F).
This is caused by the difference in the relative change in TL compared to NS induced by different

107
FGF and VEGF concentrations. As explained above, at high FGF concentrations, the increase in
TL is contributed most by the increase in the formation of new sprouts rather than the growth of
existing sprouts, while VEGF showed opposite effects.

Length of growth factor stimulation (1-3 days)
TL, NS, and AL increase with the increase in length of cell stimulation (Figure 4-4). Note
that since we are studying responses within three days, we assume that anastomosis has not
happened yet or is at a minimum level. Otherwise, we would expect that the AL reaches a plateau
at some point and might even decrease in a long term.

Initial number of cells (250, 500, and 750 cells)
TL and NS increase with the increase in the initial number of cells, but AL is independent of the
initial number of cells (Figure 4-4). Since the characteristics of endothelial sprouting are
qualitatively similar among the groups of different initial number of cells, we take 500 cells as a
representative initial number of cells to investigate the effects of FGF and/or VEGF stimulation to
reduce computational cost in the remainder of this study.

108

Figure 4-4. Predicted sprouting responses stimulated by single agents. Response to FGF
stimulation, left panels: Predicted TL (μm) (A), NS (C), and AL (E) stimulated by low,
intermediate, and high levels of FGF. Response to VEGF stimulation, right panels: Predicted TL
(μm) (B), NS (D), and AL (F) stimulated by low, intermediate, and high levels of VEGF. 250-,
500-, and 750-cell spheroid sprouting responses when simulated for 1, 2, and 3 days. Bars are
mean model predictions of 18 best fits.

4.4.2.2. Co-stimulation
We next studied the effects of FGF and VEGF co-stimulation in endothelial sprouting.
Generally, we found that similar to mono-stimulation, TL, NS, and AL are FGF and VEGF dose
dependent. Specifically, TL and NS increase with the increase in FGF or VEGF concentrations
while AL has a biphasic dose response (Figure 4-5). Also, TL, NS, and AL increase with the
increase in length of cell stimulation (Figure 4-5).
Furthermore, we found that FGF plays a dominant role in the combination effects in
endothelial sprouting as co-stimulation exhibits the features of FGF-induced sprouting (Figures 4-
0
5000
10000
15000
250 500 750 250 500 750 250 500 750
0
5000
10000
15000
250 500 750 250 500 750 250 500 750
0
20
40
60
250 500 750 250 500 750 250 500 750
0
20
40
60
250 500 750 250 500 750 250 500 750
0
100
200
300
400
250 500 750 250 500 750 250 500 750
0
100
200
300
400
250 500 750 250 500 750 250 500 750
Day1 Day2 Day3 Day1 Day2 Day3
Figure 4
A
B
C
E
D
F

109
4, 4-5, and S4-4). Specifically, the effects of co-stimulation on TL, NS, and AL are more sensitive
to FGF concentration change compared to VEGF (Figure 4-4). Also, TL, NS, and AL induced by
co-stimulation are approximately the same level as FGF stimulation alone, as the ratios of
combination effects, relative to FGF mono-stimulation, are approximately equal to one (Figure
S4-4A-C left). In comparison, TL, NS, and AL induced by co-stimulation are greater than VEGF
stimulation alone, as the ratios of combination effects, relative to VEGF mono-stimulation, are
greater than one (Figure S4-4A-C right).
In summary, endothelial sprouting is ligand- and dose-dependent and has different short-
term and long-term responses. In addition, the initial number of cells is important in the sprouting
process. Also, FGF plays a dominant role in the effects of FGF and VEGF co-stimulation on
endothelial sprouting. Moreover, the predicted effects of the co-stimulation by FGF and VEGF on
endothelial sprouting were not significantly greater than FGF mono-stimulation alone.
To understand the contributions of cellular behaviors in the process of endothelial
sprouting, we next investigated the effects of FGF and VEGF on rcp, rsg, and p. NS is a result of
the number of cells and the probability of a sprout formation, which are determined by rcp and p
for a certain number of cells initially present. TL is a result of NS and the growth of the sprouts
determined by rsg. And AL equals TL/NS. Thus, understanding how the rate of cell proliferation,
rate of sprout growth, and the probability of sprouting depend on growth factor concentration gives
insight into the observable sprouting features (number of sprouts, total sprout length, and average
sprout length).

110

Figure 4-5. Predicted sprouting responses in response to FGF and VEGF co-stimulation. Co-
stimulation of FGF- and VEGF- induced TL (μm) on day 1 (A), day 2 (D), and day 3 (G); NS on
day 1 (B), day 2 (E), and day 3 (H); and AL (μm) on day 1 (C), day 2 (F), and day 3 (I).

4.4.3. The cell proliferation and sprout growth of existing sprouts are predicted to be more
important in the sprouting process compared to the effect of rate of forming a new sprout
4.4.3.1. Mono-stimulation
We studied the cell proliferation, sprout growth, and the probability of sprouting in
response to FGF and VEGF mono-stimulation. We found that FGF and VEGF mono-stimulation
both show sigmoidal dose response curves for cell proliferation (Figure 4-6A) and sprout growth
(Figure 4-6B). Additionally, FGF induces faster cell proliferation and sprout growth than VEGF

111
at the same concentrations. This is because FGF-induced rcp and rsg are higher than VEGF-induced
rcp and rsg, respectively, at the same ligand concentrations (Figure 4-6A-B). In addition, the large
variation in p suggests that the exact value for the probability of forming a new sprout does not
significantly affect the model predictions (Figure 4-6C). This indicates that NS is more dependent
on the number of cells compared to the probability of sprout initiation. It also suggests that the
number of cells/cell proliferation and the growth of existing sprouts are more important in
contributing to the TL over time than the chance of forming a new sprout. Thus, rcp and rsg are the
main focus in the remainder of this study.

4.4.3.2. Co-stimulation
Similar to mono-stimulation, we next investigated the effects of FGF and VEGF co-
stimulation in the rate of cell proliferation, rate of sprout growth, and the probability of sprouting.
We found that rcp, rsg, and p increase as the ligand concentration increases (Figure 4-6D-F).
Moreover, FGF is dominant in the combination effects, as the co-stimulation exhibits the dose-
dependent features and approximately the same magnitude in inducing rcp, rsg, and p as in response
to FGF mono-stimulation. Specifically, rcp and rsg are more sensitive to FGF concentration and
relatively independent of the VEGF concentration (Figure 4-6D-E). Also, the values of rcp, rsg, and
p are only slightly higher with FGF and VEGF co-stimulation, compared to FGF mono-
stimulation; while with co-stimulation, rcp is significantly higher, and the values of rsg and p are
also much higher compared to VEGF mono-stimulation (Figure S4-5). It suggests that FGF is
dominant in promoting cell proliferation, sprout growth, and the probability of sprouting in the
combination effects. It is also consistent with the dominant role of FGF in the combination effects
in the sprouting characteristics specifically TL, NS, and AL observed in the previous section.

112
In summary, rcp and rsg are more important in endothelial sprouting compared to p.
Moreover, FGF plays a dominant role in the combination effects.

Figure 4-6. Predicted rcp, rsg, and p in response to mono-and co-stimulation of FGF and
VEGF. Effects of mono-stimulation of FGF (yellow) or VEGF (blue) on rcp (A), rsg (B), and p (C).
Effects of co-stimulation of FGF and VEGF on rcp (D), rsg (E), and p (F). Curves in Panels A-C
are the mean model predictions of 18 best fits. Shaded regions show standard deviation of the fits.

4.4.4. The MAPK and PI3K pathways contribute to cell proliferation, sprout growth, and
probability of sprouting in different ways
We next studied the contributions of MAPK and PI3K pathways in rcp, rsg, and p to gain
insight into the mechanisms of the endothelial sprouting process in response to mono- and co-
stimulation of FGF and VEGF.

4.4.4.1. Mono-stimulation
Cell proliferation

113
FGF and VEGF stimulate cell proliferation in different ways: FGF-induced cell
proliferation is dominated by pERK, while VEGF-induced cell proliferation is promoted by pAkt
(Figure 4-7A-B). The dose dependent feature of cell proliferation induced by FGF is mostly
contributed by ERK phosphorylation, as the rate of cell proliferation influenced by pERK exhibits
the dose-dependent feature of the overall cell proliferation rate, while the rate of cell proliferation
influenced by pAkt is independent of the FGF concentration (Figure 4-7A). In addition, at low
FGF concentrations (< 0.12 ng/ml), Akt phosphorylation plays a more important role in the cell
proliferation rate, compared to its impact at high FGF concentrations (Figure 4-7A). As the FGF
concentration increases, the impact of rcp_ERK increases and eventually surpasses the influence of
rcp_pAkt (or Akt activation) at higher FGF concentrations (Figure 4-7A). In contrast, the cell
proliferation behavior and dose-dependent feature induced by VEGF are mostly contributed by
Akt phosphorylation, while ERK phosphorylation shows a negligible contribution (Figure 4-7B).
Moreover, the effect of basal cell proliferation minimally contributes to the response induced by
FGF or VEGF stimulation (Figure 4-7A-B).

Sprout growth
The MAPK and PI3K pathways contribute differently to FGF- and VEGF-induced sprout
growth (Figure 4-7C-D). First, we found that the basal rsg (6.87 μm/h) plays a major role in the
overall rsg, accounting for 42% and 63% of the FGF- and VEGF-induced rsg at their plateau levels,
respectively (Figure 7C-D). In addition to the basal rsg, at low FGF concentrations (< 0.03 ng/ml),
Akt phosphorylation also plays a substantial role in rsg (Figure 4-7C). As the FGF concentration
increases, the impact of rsg_ERK increases. Eventually, the contributions of pAkt and pERK to the
sprout growth rate plateau at 4.07 μm/h and 5.52 μm/h (Figure 4-7C), which are the mean values

114
of fitted VpAkt_sg and VpERK_sg respectively. When VEGF is lower than 3 ng/ml, basal rsg is dominant
in VEGF-induced rsg (Figure 4-7D). As VEGF concentration increases, rsg_pAkt increases and
plateaus at 4.07 μm/h (VpAkt_sg) and rsg_ERK is negligible in VEGF-induced rsg (Figure 4-7D).

Figure 4-7. The contributions of MAPK and PI3K/Akt pathways to rcp and rsg in response to
FGF and VEGF mono-stimulation. Contributions of ERK (purple), Akt (green), and basal (gray)
for FGF-induced rates of cell proliferation, rcp (A) and sprout growth, rsg (C). Contributions of
ERK (purple), Akt (green), and basal (gray) for VEGF-induced cell proliferation, rcp (B) and sprout
growth, rsg (D).

4.4.4.2. Co-stimulation
We then investigated the contributions of the MAPK and PI3K pathways in rcp and rsg by
the co-stimulation of FGF and VEGF. We found that ERK and Akt activation contribute differently
to cell proliferation and sprout growth upon stimulation with FGF and VEGF in combination, and
FGF plays a dominant role of in combination effects in cell proliferation and sprout growth.

115
Cell proliferation
First, rcp, rcp_pERK, and rcp_pAkt stimulated by a combination of FGF and VEGF mirror the
corresponding responses stimulated by FGF alone (Figures 4-6A, 4-6D, 4-7A and S4-6A-B).
Specifically, the dose-dependent manner of the cell proliferation rate induced by FGF and VEGF
co-stimulation is mostly impacted by rcp_pERK, while rcp_pAkt is independent of the FGF or VEGF
concentrations (Figure S4-6A-B). Also similar to FGF mono-stimulation, Akt phosphorylation
plays a more substantial role in the cell proliferation rate at FGF concentrations lower than 0.12
ng/ml, while rcp_ERK increases as the FGF concentration increases and surpasses the influence of
rcp_pAkt at higher FGF concentrations (>0.12 ng/ml) (Figure S4-6A).

Sprout growth
The dose-dependent manner of sprout growth rate induced by FGF and VEGF co-
stimulation is mostly influenced by rsg_pERK (Figure S4-6C), while rsg_pAkt is independent of FGF
or VEGF concentrations (Figure S4-6D). In addition, Akt phosphorylation is more important in rsg
at FGF concentrations lower than 0.03 ng/ml, and the impact of rsg_ERK increases as the FGF
concentration increases (Figure S4-6C-D). These predictions show the same feature as FGF-
induced sprout growth (Figure 4-7C), indicating the dominant role of FGF on sprout growth.

4.4.5. ERK pathway regulates vessel network mainly via regulating cell proliferation and NS,
while Akt pathway mainly affects vessel network via regulating sprout growth
We applied the model to identify the parameters and initial concentrations in the ERK-Akt
model that are influential to the cell proliferation rate and sprout growth rate. This allows us to
gain mechanistic insight into how to modulate the effects of intracellular signals in endothelial

116
sprouting. We identified the influential variables by performing eFAST and evaluating the
calculated Stimax values. The variables (Table 4-2) that have Stimax values greater than 0.3 are
identified as influential and considered as potential targets for pro- and anti-angiogenic strategies.
Specifically, ERK and MEK are predicted to be influential to rcp and rsg, while Akt and Ptase2 are
influential to rsg. We note that no parameter was identified as influential to the probability of
sprouting, further justifying our focus solely on the rates of cell proliferation and sprout growth.
The model predictions suggest that the MAPK pathway is more influential to rcp, and MAPK and
PI3K/Akt pathways are both influential to rsg.

rcp rsg
ERK
MEK
Akt
ERK
Ptase2
MEK
Table 4-2. Influential parameters.

However, eFAST only tells us that those variables are influential, but the information of
how those variables influence model outputs are limited. We then varied two representative
influential variables, ERK and Akt, by 0.1- and 10-fold and predicted the rcp and rsg (Figure 4-8),
as well as TL, NS, and AL (Figures 4-9 and 4-10) and compared with the baseline model
predictions. We found that the strategies to modulate endothelial sprouting are context dependent,
and ERK and Akt pathways regulate vessel network differently.

117

Figure 4-8. Predicted representative targets for modulating rcp and rsg. Predicted rcp (A) and
rsg (D) from baseline model. Predicted rcp (B) and rsg (E) when ERK is varied by 0.1- (left) and 10-
fold (right). Predicted rcp (C) and rsg (F) when Akt is varied by 0.1- (left) and 10-fold (right).

Varying ERK
In Figure 4-8A and B, upregulating ERK is predicted to be effective in promoting cell
proliferation at low to intermediate but not at high FGF concentrations. In contrast, downregulating
ERK is predicted to be more effective in inhibiting cell proliferation at intermediate to high FGF
concentrations (Figure 4-8A and B). In addition, upregulating ERK is more effective in promoting
sprout growth at low FGF concentrations (Figure 4-8D and E). Downregulating ERK is most
effective in inhibiting sprout growth at intermediate and high FGF, particularly in combination
with low VEGF concentration (Figure 4-8D and E).
We next investigated how finitely varying ERK affects endothelial sprouting, specifically,
TL, NS, and AL on days 1-3 (Figure 4-9). The model showed no obvious effects of increasing or
decreasing ERK by 10-fold on TL, NS, and AL on day 1 (Figure 4-9Ai, Bi, and Ci). In addition,
increasing the ERK level is effective in promoting TL and NS at low to intermediate FGF levels
but has no obvious effects at high FGF concentration on day 2 (Figure 4-9Aii, Bii). and day 3
(Figure 4-9Aiii, Biii). Also, decreasing ERK is effective in inhibiting TL and NS at intermediate

118
and high FGF levels, but this strategy is not very promising at low FGF level on day 2 (Figure 4-
9Aii, Bii) and day 3 (Figure 4-9Aiii, Biii). Furthermore, the increase in TL is less than the increase
in NS, which causes a decrease in AL (Figure 4-9Cii, Ciii) since AL equals TL/NS.

Figure 4-9. Predicted effects of ERK in modulating TL (μm), NS, and AL (μm). Predicted TL
(A), NS (B), and AL (C) when ERK is varied by 0.1- (left) and 10-fold (right), compared with
baseline model predictions (middle) on days 1-3 (i-iii).

Varying Akt
In Figure 4-8A and C, we did not observe obvious effects in rcp when increasing or
decreasing Akt by 10-fold. However, downregulating Akt is predicted to be effective in inhibiting
sprout growth at all FGF and VEGF combinations, while upregulating Akt does not have obvious
effects in rsg (Figure 4-8D and F).
We then investigated how finitely varying Akt affects TL, NS, and AL on days 1-3 (Figure
4-10). We found that increasing Akt level has no obvious effects on TL, NS, or AL at any of the
simulated combinations of FGF and VEGF on days 1-3 (Figure 4-10). However, decreasing Akt
level leads to a decrease in AL and TL at all combinations of FGF and VEGF on days 1-3, but no
obvious effects in NS were observed when the Akt level was decreased (Figure 4-10).

119

Figure 4-10. Predicted effects of Akt in modulating TL (μm), NS, and AL (μm). Predicted TL
(A), NS (B), and AL (C) when Akt is varied by 0.1- (left) and 10-fold (right), compared with
baseline model predictions (middle) on days 1-3 (i-iii).

In summary, targeting the ERK pathway is predicted to control the vessel network mainly
via regulating cell proliferation and NS, while targeting Akt pathway mainly influences the vessel
network via regulating sprout growth. In addition, the effects of the molecular signaling pathways
on sprouting is not obvious after one day of stimulation but is more effective in a long-term (days
2-3). Overall, our model can identify potential effective pro- and/or anti-angiogenic targets and
predicts the effects of perturbing those targets under different conditions.

4.5. Discussion
We developed a hybrid agent-based model characterizing the endothelial sprouting process
driven by integrating molecular signals, pERK and pAkt, upon the mono- and co-stimulation of
two pro-angiogenic factors FGF and VEGF. The intracellular signaling model of ERK and Akt
activation in response to FGF and VEGF stimulation in endothelial cells was adapted from our

120
previous work (39). The endothelial sprouting process was modeled by assuming the cellular
responses, including cell proliferation, sprout growth, and the probability of forming a new sprout
are driven by pERK and pAkt, following Hill functions. Unknown parameters were estimated by
fitting the model to experimental data. Additionally, we validated the model using a separate set
of data.
The fitted model predicts the TL, NS, and AL upon stimulation by FGF and VEGF, alone
and in combination, on days 1-3. We particularly focus on TL, NS, and AL because they are
metrics examined most often in in vitro studies (119, 120, 123-125). The model predicts that the
type and concentration of ligand, length of cell stimulation by the ligand, and initial number of
cells are important in endothelial sprouting. The predicted dose dependent sprouting responses
including TL induced by FGF and VEGF in Figure 4-4 are consistent with the experimental data,
which showed dose-dependent total tubular length induced by FGF and VEGF (50). Also, in the
same study, 0.1 ng/ml FGF exhibits approximately the same level of increase in HUVECs
proliferation and migration as 25 ng/ml VEGF after 24 hours (50), which agrees with our model
prediction that FGF induces greater sprouting responses than VEGF at same concentrations
(Figure 4-4). Furthermore, our model suggests that FGF promotes greater sprouting responses in
the long-term compared to VEGF (Figure 4-4). In addition, the model predicts that FGF plays a
dominant role in the combination effects in endothelial sprouting, and co-stimulation of FGF and
VEGF only slightly increases TL, NS, and AL, compared to FGF simulation alone within three
days (Figures 4-4, 4-5, and S4-4). Also, the effects of low FGF concentration in combination with
high VEGF concentration in TL, NS, and AL shows no obvious difference compared to VEGF
mono-stimulation within three days (Figures 4-4, 4-5, and S4-4). It is consistent with experimental
observation which showed no significant increase in HUVEC proliferation in 72 hours, migration

121
in 8 hours, and total tubular length in 24 hours stimulated by the combination of 0.1 ng/ml FGF
and 25 ng/ml VEGF compare with their mono-stimulation (50). Moreover, our model can
supplement experimental endothelial sprouting assays to differentiate and analyze the
contributions of cellular and molecular responses during the overall sprouting process including
cell proliferation, sprout growth and formation, and the activation of the ERK and Akt pathways.
Specifically, the model predicts that the ERK pathway regulates vessel network mainly via
regulating cell proliferation and NS, while Akt pathway mainly affects vessel network via
regulating sprout growth (Figures 4-8, 4-9, and 4-10), which is in line with the literature results
reporting that ERK is believed to mainly promote cell proliferation (15) and Akt is more important
in cell survival (17-21) and migration (21, 24, 25).
Compared to other models that study cellular behaviors, our mechanistic model considers
the intracellular signaling, and quantitatively analyzes cellular responses driven by integrating
molecular signals, pERK and pAkt. Our model explicitly examines how cellular behaviors are
driven by pERK and pAkt, which are downstream signals that regulate angiogenic cellular
responses. Thus, we can apply our model to mechanistically study the roles of intracellular
signaling species in affecting endothelial sprouting. For example, Tong and Yuan constructed a
computational model to study vessel growth in rat cornea based on assumptions that cellular
responses are only dependent on FGF bound FGFR, where the probability of sprout formation and
the speed of vessel growth are linearly proportional to the fraction of FGFR occupied by FGF (31).
Our model can complement such models to understand intracellular mechanisms that regulate
cellular responses. In another study, Norton and Popel constructed a computational model to study
vessel growth in tumors and showed that the proliferation rate has a greater effect on the spread
and extent of vascular growth compared to migration rate (42). Our model predictions agree with

122
these previous modeling works, as we predict that cell proliferation and the number of cells are
critical factors that contribute to endothelial sprouting (Figure 4-4) and that varying ERK level
seems to be more influential in TL than varying Akt level (Figures 4-9 and 4-10). Our work goes
further, in that the model considers the intracellular signaling and provides mechanistic insight
into the signaling factors driving endothelial sprouting.
Our model can be utilized to study the efficiency of pro- or anti-angiogenic therapies. The
model predicts that potential strategies to modulate endothelial sprouting are context dependent
(Figures 4-8, 4-9, and 4-10) and it can identify potential effective pro- and/or anti-angiogenic
targets under different conditions and study their efficacy. Moreover, the intracellular angiogenic
signals such as pERK and pAkt are believed to play important roles in cellular behaviors,
especially cell survival, proliferation, and migration (15, 17-21, 24, 25). However, there is limited
quantitative understanding of how integrating the intracellular angiogenic signals affect cellular
behaviors, especially in cases where approaches to inhibit angiogenesis have counterintuitive
effects on pERK and pAkt. For example, it has been shown that the MEK inhibitor PD0325901
upregulates the PI3K pathway signaling (97). Our model can investigate such experimental
observations and predict the overall cellular responses driven by the integrated molecular signals
from pERK and pAkt.
This model can be utilized in combination with other modeling framework that predict
intracellular signaling to provide more mechanistic insight into certain cellular responses. For
example, there are models that study pro-angiogenic signaling, including sphingosine kinase 1 and
calcium responses induced by VEGF, which are the downstream signals of ERK1/2 (126), and
also angiopoietin-Tie signaling in endothelial cells (127), which has been shown to be important
in vessel development, permeability, vascular homeostasis and remodeling (127-130). Other

123
models study anti-angiogenic signaling such as TSP1-CD36 signaling that influences endothelial
cell apoptosis (34) and antagonizes VEGF-induced eNOS signaling (131, 132), as nitric oxide is a
major vasodilator (131) and important in angiogenesis and vascular permeability (133, 134). Still
other models characterize sprouting angiogenesis behaviors by considering how Notch signaling
in endothelial cells determines the tip cell and vessel branching as reported in other relevant studies
(117, 135-137). Our model can be combined with these existing models of angiogenic signaling
and applied in various way to provide quantitative insight on both molecular and cellular levels.
We acknowledge some limitations in our model. Although many studies reported that
pERK and pAkt play important roles in cell survival, proliferation, and migration, the upstream
species in the network (for example, pVEGFR2, pFGFR1, and PI3K), can activate other pathways
that are not included in this work, but may also contribute to relevant angiogenic cellular responses
(4, 41). Also, to ease model construction, we excluded VEGFR1 and neuropilin-1 (NRP1) although
their binding with VEGF plays a role in angiogenesis. We can incorporate their contributions into
the model in future studies.
We also made some assumptions to simplify our model. First, the molecular interactions
usually happen on the order of seconds to minutes. However, it usually takes hours or even days
to respond on a cellular level. To bridge the difference in time scale, we assumed that the cellular
behaviors in short term are driven by the molecular signals within two hours, specifically the
maximum pERK and pAkt levels. Other computational models have made similar assumptions.
For example, Adlung et al. correlated the molecular signals, cyclin D2, cyclin G2, p27, and pS6,
which were characterized within 60 minutes, to mCFU-E, BaF3-EpoR, and 32D-EpoR cell
proliferation, which were analyzed after 14-20 hours, 62 hours, and 38 hours, respectively (113).
In addition, Tong and Yuan assumed cellular responses are only dependent on the fraction of FGF

124
bound FGFR, and they are linearly related to investigate FGF-induced angiogenic dose responses
in rat corneal pocket assay (31). They analyzed migration distance, total and average vessel length,
and total number of vessels at time from 0 to 120 hours, while the time for FGF bound FGFR to
reach 90% of steady-state value was reported to be within 25 minutes (31). Also, Padera et al.
constructed a mathematical model to study FGF-mediated cellular response by assuming the signal
that drives proliferation of cultured cells from the F32 cell line is only dependent on ligand bound
receptor signaling complexes (30). The model was validated with F32 cell proliferation data, which
were evaluated by counting the cell number after 72 hours incubation (30).
Moreover, our main focus is to study the intrinsic properties of endothelial cells responding
to FGF and/or VEGF stimulation. Unfortunately, the experimental setup in in vitro studies may
involve environmental factors that influence cell behaviors other than our focus, FGF and VEGF.
For example, to maintain sufficient nutrient supply, cell culture media is typically changed every
three days, which would inevitably induce a change in temperature and air composition, affecting
cellular behaviors. In addition, cell proliferation (121, 122, 138), migration (114, 115, 139), and
sprouting assays (119, 123, 138) are usually conducted in a short period of time following cell
culture, typically within three days. Thus, to reduce the effects of other possible factors, we only
studied short-term cellular behaviors, specifically three days, when nutrients and space are still
sufficient. We assumed the cell intrinsic properties, specifically the cell proliferation rate, sprout
growth rate, and probability rate of forming a new sprout, remain constant within three days
considering a relatively stable experimental environment. Other computational studies have also
made similar assumptions. For example, Norton and Popel investigated the effects of endothelial
proliferation and migration rates on vascular growth by simulating vasculature using various time-
independent proliferation and migration rates for times of up to 200 days (42). Tong and Yuan

125
studied FGF-induced angiogenic dose responses in rat corneal pocket assay for times of up to five
days (31).
In addition, we assumed the daughter cells inherit the same cell proliferation rate, sprout
growth rate, and the probability rate of forming a new sprout from the mother cells, again
considering a relatively stable experimental environment within three days. Similar assumptions
have also been made by other computational work. For example, Roy and Finley built a multiscale
computational model to study pancreatic tumor growth with an assumption that the daughter cell
inherits all properties and the last intracellular metabolic state of the parent cell (140). Letort et al.
built a multiscale agent-based model assuming daughter cell inherits its signaling network state
from the mother cell and applied the model to study the effects of cell heterogeneity in tumor
growth in response to tumor necrosis factor treatment (141).
Last, we assumed that the maximum pERK and pAkt drive endothelial cell proliferation
and sprouting follow Hill functions. There are indeed some computational studies that apply a
threshold function to decide if a cell is eligible to migrate or proliferate (142, 143). However,
HUVECs have been reported to respond to low levels of FGF or VEGF stimulation. For example,
Bai et al. showed that FGF concentration as low as 0.1 ng/ml significantly increased HUVECs
total tubular length compared to control on Matrigel for 24 hours (50). VEGF has been shown to
induce a half-maximal effect on tubule formation on Laponite substrates at 0.01 g/ml (144). Also,
Wolfe et al. showed that the EC50 value for VEGF-induced tube length response is 0.67 ng/ml for
HUVECs co-cultured with normal human dermal fibroblasts (145). In addition, angiogenic cellular
responses, for instance endothelial sprouting (9, 50) and vessel density (146), have been shown to
be FGF and VEGF dose dependent and will reach a plateau if FGF or VEGF concentration is
higher than a certain saturation level. Thus, we decided not to use a threshold function for

126
angiogenic cellular responses. We instead applied a Hill function to account for the observed
angiogenic cellular responses, even for low growth factor concentrations.
In the future, we can incorporate more details as more mechanistic information becomes
available. Despite these limitations, our model provides quantitative insight into angiogenic
signaling and cellular responses and can be utilized as a framework for future mechanistic studies.

4.6. Conclusion
In conclusion, we developed a mathematical model to characterize endothelial sprouting
driven by pERK and pAkt in response to the stimulation of two main pro-angiogenic factors, FGF
and VEGF. The model quantitatively studied FGF- and VEGF-mediated cell proliferation, sprout
growth, and formation of new sprouts, and provided mechanistic insight into endothelial sprouting.
The understanding of the regulation of angiogenesis signals on a molecular scale, and further on a
cellular level, can better aid the development of pro- and anti-angiogenic strategies.

4.7. Acknowledgements
The authors thank members of the Finley research group for critical comments and
suggestions. The authors acknowledge the support of the US National Science Foundation
(CAREER Award 1552065).

4.8. Supplementary Materials
Parameterization
Cell proliferation module

127
It has been reported that HUVECs up to 45 generations showed a minimal doubling time
of 14 hours when cultured with an optimized medium (40). Since ERK activation is believed to
mainly promote cell proliferation, the upper bound of 𝑉 𝑝𝐸𝑅𝐾 _𝑐𝑝
is set at 1/14 1/h. It has been
reported that inhibiting PLCγ, an upstream species of the ERK pathway, leads to 50% proliferation
inhibition, while blocking PI3K results in 28% proliferation inhibition (41). Therefore, we set
𝑉 𝑝𝐴𝑘𝑡 _𝑐𝑝
to be (28/50) * 𝑉 𝑝𝐸𝑅𝐾 _𝑐𝑝
. We set the baseline values of 𝐾𝑚
𝑝𝐸𝑅𝐾 and 𝐾𝑚
𝑝𝐴𝑘𝑡 at 1×10
4

molecules/cell to capture the reported range of maximal pERK and pAkt levels: 2.6×10
−3
to
6×10
5
molecules/cell and 4.7×10
3
to 1.1×10
6
molecules/cell, respectively, when FGF and
VEGF vary from low to high levels. Since cell doubling time for HUVECs range from 17 to 72
hours (42), we set the baseline value of 𝑟 𝑐𝑝 _𝑏𝑎𝑠𝑎𝑙 at 1/72 1/h.

Sprouting model
Norton and Popel studied the effects of proliferation and migration rates on sprouting
angiogenesis, and they varied te migration rate from 0.24 to 40 μm/h (42). Also, Tong and Yuan
investigated angiogenesis in rat corneal pocket assay and chose the maximum rate of sprout length
increase to be 20 μm/h (31). Based on this prior work, we took 16 μm/h as the baseline value of
sprout growth rate. Also, it has been shown that inhibiting PLCγ leads to 65% HUVEC migration
inhibition, while blocking PI3K results in 48% HUVEC migration inhibition (41). Therefore, we
set the baseline values of 𝑉 𝑝𝐸𝑅𝐾 _𝑠𝑔
at 9.2 μm/h (16 μm/h * 65/(65+48)), 𝑉 𝑝𝐴𝑘𝑡 _𝑠𝑔
at (48/65) *
𝑉 𝑝𝐸𝑅𝐾 _𝑠𝑔
, and 𝑟 𝑠𝑔 _𝑏𝑎𝑠𝑎𝑙 at 0.4 μm/h.
In addition, it has been shown that the average rate of sprout formation is approximately
10
−4
1/μm/h (43), we set the baseline of 𝑝 𝑏𝑎𝑠𝑎𝑙 at 10
−4
1/μm/h and the upper bound of 𝑉 𝑝𝐸𝑅𝐾 _𝑝 at
10
−3
1/μm/h. Due to lack of experimental evidence to correlate 𝑉 𝑝𝐸𝑅𝐾 _𝑝 and 𝑉 𝑝𝐴𝑘𝑡 _𝑝 , we set the

128
ratio k at 1 as a starting point to estimate the relative contribution of ERK and Akt pathways to the
overall probability of forming a new sprout. We used the average endothelial cell diameter to
convert p from units of 1/μm/h to 1/cell/h. In this study we focused on endothelial cells, which
have been reported to have a mean cell volume of 1009 ± 180 𝜇 m
3
(87). Therefore, we calculated
the average diameter (d) of an endothelial cell using 1009 𝜇 m
3
cell volume and assuming that
endothelial cells are spherical: (4/3)*pi*(d/2)
3
= 1009 𝜇 m
3
. By this calculation, d is calculated to
be 12.4 μm.
We set all the baseline values of Hill coefficients at 2 as a starting point. We provided Table
S4-1 to summarize all the baseline parameter values.

Figure S4-1. Results of eFAST sensitivity analysis. The highest total sensitivity index (Stimax)
values for (A) parameters in the cell proliferation module and (B) parameters in the sprouting
model. The color bars indicate the Stimax values. Influential parameters are highlighted as red.

129

Figure S4-2. Model values estimated in fitting. (A) Distribution of the parameters in the cell
proliferation module. (B) Distribution of all the parameters in the sprouting model. Influential
parameters are highlighted as red. (C) Distribution of the influential parameters in the sprouting
model. Each dot represents one fit. Bars are median ± 95% confidence interval.

Figure S4-3. Model comparison to training and validation data for FGF or VEGF
stimulation. (A) Total sprout length induced by 0 – 64 ng/ml VEGF for 24 hours cultured with
500-cell spheroid initially. (B) The fold change of the average length and number of sprouts
induced by 25 ng/ml FGF and 25 ng/ml VEGF for 24 hours cultured with 400-cell spheroid initially
compared to the control. Diamonds in Panel A are experimental data from Heiss et al. (119). Curve
in Panel A is the mean model predictions of the 15 best fits. Shaded regions show standard
deviation of the fits. Solid and dashed bars in Panel B are mean ± standard deviation of Liebler et
al data (120) and model predictions.

V
pERK_sg

V
pERK_p

k
n
sg

n
p

Km
pERK_sg

Km
pERK_p

Km
pAkt_sg
Km
pAkt_p

r
sg_basal

p
basal

1 ×10
-5
1 ×10
-3
1 ×10
-1
1 ×10
1
1 ×10
3
1 ×10
5
V
pERK_sg

V
pERK_p

k
n
sg

n
p

Km
pERK_sg

Km
pERK_p

Km
pAkt_sg
Km
pAkt_p

r
sg_basal

p
basal

1 ×10
-5
1 ×10
-3
1 ×10
-1
1 ×10
1
1 ×10
3
1 ×10
5
V
pERK_sg

V
pERK_p

k
n
sg

n
p

Km
pERK_sg

Km
pERK_p

Km
pAkt_sg
Km
pAkt_p

r
sg_basal

p
basal

1 ×10
-5
1 ×10
-3
1 ×10
-1
1 ×10
1
1 ×10
3
1 ×10
5
V
pERK_sg
V
pERK_p
k
n
sg
n
p
r
sg_basal
p
basal
1 ×10
-5
1 ×10
-3
1 ×10
-1
1 ×10
1
1 ×10
3
Km
pERK_cp
Km
pAkt_cp
V
pERK_cp
n
cp
r
cp_basal
1 ×10
-4
1 ×10
-2
1 ×10
0
1 ×10
2
1 ×10
4
1 ×10
6
Figure S2
A B C
Figure S3
(Heiss et al., 2015)
A
B
0
0.5
1
1.5
2
2.5
25 ng/ml FGF 25 ng/ml VEGF 25 ng/ml FGF 25 ng/ml VEGF
Average Length Number of Sprouts
Liebler et al., 2012 Model Predictions
Fold Change

130

Figure S4-4. Comparison of mono- and co-stimulation in TL, NS, and AL. Ratios comparing
the combination effects to the individual effects in response to FGF and VEGF for TL (i), NS (ii),
and AL (iii) on day 1 (A), day 2 (B), and day 3 (C).

Figure S4-5. Comparison of mono- and co-stimulation in rcp, rsg, and p. Ratios comparing the
combination effects to the individual effects in response to FGF for rcp (A), rsg (C), and p (E).
Ratios comparing the combination effects to the individual effects in response to VEGF for rcp (B),
rsg (D), and p (F).

131

Figure S4-6. The contributions of MAPK and PI3K/Akt pathways in rcp, rsg, and p in
response to FGF and VEGF co-stimulation. rcp contributed by ERK (A) and Akt activation (B),
rsg contributed by ERK (C) and Akt activation (D), and p contributed by ERK (E) and Akt
activation (F).

132

Table S4-1. List of parameters.

Table S1. List of parameters.
Baseline Values Units
V_pERK_cp 0.0007143 1/h
V_pAkt_cp (28/50)*V_pERK_cp 1/h
Km_pERK_cp 1.00E+04 molecules/cell
Km_pAkt_cp 1.00E+04 molecules/cell
n_cp 2
rcp_basal 0.01389 1/h
V_pERK_sg 9.2035 μm/h
V_pAkt_sg (48/65)*V_pERK_vg μm/h
Km_pERK_sg 1.00E+04 molecules/cell
Km_pAkt_sg 1.00E+04 molecules/cell
n_sg 2
rsg_basal 0.4 μm/h
V_pERK_p 1.00E-05 1/ μm/h
V_pAkt_p k*V_pERK_p 1/ μm/h
k 1
Km_pERK_p 1.00E+04 molecules/cell
Km_pAkt_p 1.00E+04 molecules/cell
n_p 2
p_basal 1.00E-04 1/ μm/h

133

Table S4-2. Fitted parameters in the cell proliferation module.

134

Table S4-3. All fitted parameters in the sprouting model.

135

Table S4-4. Fitted parameters in the sprouting model.

Table S4. Fitted parameters in the sprouting model.
Set
Error
V_pERK_sg
V_pERK_p k n_vg n_p
1 0.002 5.90E+00 9.87E-04 3.21E-01 3.54E+01 3.55E+01
2 0.009 4.85E+00 7.41E-04 9.67E-02 4.63E+01 1.17E-01
3 0.005 5.36E+00 1.46E-04 3.25E+00 1.01E+02 1.04E+01
4 0.013 5.78E+00 1.60E-04 2.34E+00 1.86E+02 3.23E-01
5 0.005 5.74E+00 2.08E-05 4.74E+01 1.65E+02 2.03E-02
6 0.008 5.48E+00 9.47E-05 5.58E+01 1.06E+02 3.62E+01
7 0.007 5.70E+00 1.25E-04 1.86E+01 1.05E+02 2.54E-01
8 0.006 5.49E+00 7.66E-04 1.31E+00 9.46E+01 7.35E+01
9 0.007 5.68E+00 3.74E-05 7.00E+01 1.02E+02 5.93E-02
10 0.008 5.09E+00 1.95E-05 4.71E+01 1.66E+02 4.48E-01
11 0.004 5.70E+00 8.76E-04 6.65E-01 1.05E+02 1.79E+00
12 0.004 5.57E+00 3.15E-04 1.36E+00 1.76E+01 3.84E-01
13 0.008 5.66E+00 2.05E-04 3.46E+01 9.82E+01 2.68E-01
14 0.007 5.59E+00 1.32E-04 4.15E+01 1.10E+02 4.12E+01
15 0.008 5.53E+00 7.77E-04 4.33E+01 9.80E+01 3.24E-02
16 0.003 5.74E+00 2.64E-05 2.78E+01 1.38E+02 2.18E-01
17 0.007 5.33E+00 9.20E-04 3.16E+00 1.14E+02 8.07E-01
18 0.003 5.13E+00 9.41E-04 2.01E-01 1.80E+02 3.88E-02
Median 5.58E+00 1.82E-04 1.09E+01 1.05E+02 3.54E-01

136
5. Chapter 5
Chapter 5
Conclusion

5.1. Overview
In this work, I have developed experimentally validated models to characterize how two
main pro-angiogenic factors, FGF and VEGF, mediate angiogenic signaling and cellular responses
in endothelial cells. I have constructed molecularly-detailed models of FGF- and VEGF-induced
MAPK and PI3K/Akt pathways, which lead to ERK and Akt activation, respectively. The models
were used to study the pERK and pAkt responses upon the mono- and co-stimulation of FGF and
VEGF. The model simulations provide a mechanistic understanding of the distinct responses of
pERK and pAkt. I then linked the molecular signals to cellular responses to investigate the effects
of FGF and VEGF interactions in endothelial cell proliferation and sprouting. Moreover, I
identified potential targets that can enhance sprouting by performing eFAST and studied their roles
in influencing endothelial sprouting.

5.2. Summary
The modeling work for how pro-angiogenic growth factors mediate intracellular signaling
and endothelial sprouting process described in Chapters 2-4 provides mechanistic understanding
of sprouting angiogenesis regulated by FGF and VEGF. In Chapter 2, I first constructed an
experimentally validated model of FGF- and VEGF-induced MAPK pathway and focused on ERK
activation in endothelial cells, which mainly promotes cell proliferation (15). I then applied the

137
model to predict the dynamics of ERK phosphorylation in response to mono- and co-stimulation
of FGF and VEGF. The model predicts that FGF is more potent in inducing maximum pERK than
VEGF at equimolar concentrations. Additionally, VEGFR2 density and trafficking parameters are
predicted to significantly influence the response to VEGF stimulation. Moreover, the model
suggests that the combination of FGF and VEGF has greater effects in inducing maximum pERK
than the summation of the individual effects.
The angiogenesis process involves not only cell proliferation but also other cellular
responses, including cell survival and migration. Thus, in Chapter 3, I expanded the ERK model
to include PI3K/Akt signaling, which mostly promotes cell survival (17-21) and migration (21, 24,
25). I applied the trained and validated mathematical model to characterize the dynamics of pERK
and pAkt in response to the stimulation of FGF and VEGF individually and in combination. The
model predicts distinct responses of ERK and Akt activation in response to FGF and VEGF. Akt
activation shows a stronger response than ERK in terms of magnitude as the maximum pAkt level
is predicted to be higher than maximum ERK activation in response to FGF and VEGF at the same
concentrations. The model also predicts that for certain ranges of ligand concentrations, the
maximum pERK level is more responsive to changes in ligand concentration compared to the
maximum pAkt level. Also, the combination of FGF and VEGF has a greater impact on increasing
the maximum pERK compared to the summation of individual effects, which is not seen for
maximum pAkt levels. In addition, our model identifies the influential species and kinetic
parameters that specifically modulate the pERK and pAkt responses, which represent potential
targets for angiogenesis-based therapies. This work provides quantitative insight into FGF and
VEGF interactions in promoting ERK and Akt phosphorylation.

138
Since there is limited quantitative understanding of how integrate intracellular angiogenic
signals affect cellular behaviors, investigating the cellular responses quantitatively with
consideration of molecular details is beneficial in understanding the angiogenesis process
mechanistically. Therefore, I developed a predictive, hybrid agent-based model to characterize
endothelial cell proliferation and sprouting by accounting for both pAkt and pERK responses in
Chapter 4. The experimentally validated model predicts that the type and concentration of ligand,
length of growth factor stimulation, and initial number of cells are important in endothelial
sprouting. Also, FGF is predicted to play a dominant role in the combination effects in endothelial
sprouting. In addition, the model investigates the contributions of cellular responses, cell
proliferation, the growth of existing sprouts, and the chance of forming a new sprout in the
sprouting process. Moreover, the model addresses questions regarding the relative roles of pERK
and pAkt in overall endothelial sprouting. Last, the model predicts that the strategies to modulate
endothelial sprouting are context dependent, and our model can identify potential effective pro-
and anti-angiogenic targets under different conditions and study their efficacy.

5.3. Future Directions
The ERK model and ERK-Akt model have provided novel insights into the mechanisms
of the interactions of FGF and VEGF in MAPK and PI3K/Akt signaling in endothelial cells.
However, there are other important intracellular angiogenic signaling in endothelial cells that play
an important role in the angiogenesis process. The balance between pro-angiogenic growth factors
(such as FGF and VEGF) and anti-angiogenic factors (such TSP1) regulate angiogenesis. Our
model can be expanded to include other angiogenic signaling pathways to improve angiogenic-
based strategies, by combining with published models. For example, angiopoietin (Ang)-Tie

139
signaling in endothelial cells has been shown to be important in vessel development, permeability,
vascular homeostasis and remodeling (127-130). Zhang et al. constructed a computational model
to characterize the Ang-Tie signaling pathway in endothelial cells. In addition, TSP1 signaling is
important in endothelial cell apoptosis (34). Wu and Finley built a molecular-detailed model to
describe TSP1-CD36 signaling (34). Also, Bazzazi et al. constructed a mechanistic model to
investigate TSP1-mediated inhibition of VEGF-induced Akt-eNOS signaling in endothelial cells
(132). Moreover, Bazzazi et al. analyzed VEGF-induced TSP1-CD47 signaling leading to ERK1/2
activation and calcium elevation using a detailed rule-based model (147). Models of these
signaling pathways can be combined with our model to better understand angiogenic signaling and
crosstalk between families.
In addition, our sprouting model does not consider the molecular signals that regulate tip
cell formation. There are many models characterize sprouting angiogenesis behaviors by
considering Notch signaling in endothelial cells to determine the tip cell and vessel branching (117,
135-137). For example, Vega et al. constructed a cellular Potts model where notch signaling
governs tip cell selection and vessel branching (135). In another study, Weinstein et al. constructed
a Boolean model including several intracellular pathways such as FGF, VEGF, Ang-Tie,
PLCγ/Calcium, PI3K/Akt, NO, Notch signaling to study the effects of microenvironment in
endothelial sprouting (137). Thus, our model can be expanded to incorporate Notch signaling in
endothelial cells, which is believed to regulate endothelial cell phenotype and important in
endothelial sprouting (116-118, 148), to determine the tip cells and vessel branching and further
characterize sprouting angiogenesis behaviors.
In our model, we excluded VEGFR1 and neuropilin-1 (NRP1) since VEGFR2 is thought
to be the main receptor on endothelial cells (106). While it has been shown that VEGFR1 promotes

140
signal transduction (107), it is largely considered to be a decoy receptor (108, 149). In addition,
NPR1 primarily acts as a coreceptor for VEGFR1 and VEGFR2 (106). We can incorporate the
contributions of VEGFR1 and NRP1 into the model in future studies.
On the cellular level, our model simulates endothelial sprouting, but only considers
endothelial cells. Angiogenesis-based strategies are important in tissue engineering and cancer
treatment. Our model can be expanded to include other cell types for example to model the tumor
microenvironment, fibroblasts, macrophages, and cancer cells to simulate various aspects of the
tissue in different contexts to identify potential targets and investigate the optimal angiogenesis-
based strategies under different conditions.

5.4. Conclusion
Ultimately, this work provides a framework to mechanistically explain experimental
results and guide experimental design to study pro-angiogenic signaling and endothelial sprouting.
This work generates novel insights into the angiogenesis processes and is beneficial to gain
quantitative understanding of molecular and cellular responses. Moreover, this model can be
utilized to explore optimal pro- and anti-angiogenic strategies systematically and guide the
angiogenic-based therapies. Furthermore, this framework can be expanded to include other pro-
and anti-angiogenic signaling and other cellular responses to gain a more comprehensive
understanding of the angiogenesis process.

141
References
1. Shiu YT, Weiss JA, Hoying JB, Iwamoto MN, Joung IS, Quam CT. The role of mechanical
stresses in angiogenesis. Crit Rev Biomed Eng. 2005;33(5):431-510.
2. Ucuzian AA, Gassman AA, East AT, Greisler HP. Molecular mediators of angiogenesis. J
Burn Care Res. 2010;31(1):158-75.
3. Risau W. Angiogenic growth factors. Prog Growth Factor Res. 1990;2(1):71-9.
4. Cross MJ, Claesson-Welsh L. FGF and VEGF function in angiogenesis: signalling
pathways, biological responses and therapeutic inhibition. Trends Pharmacol Sci. 2001;22(4):201-
7.
5. Bai Y, Bai L, Zhou J, Chen H, Zhang L. Sequential delivery of VEGF, FGF-2 and PDGF
from the polymeric system enhance HUVECs angiogenesis in vitro and CAM angiogenesis. Cell
Immunol. 2018;323:19-32.
6. Simons M, Annex BH, Laham RJ, Kleiman N, Henry T, Dauerman H, et al.
Pharmacological treatment of coronary artery disease with recombinant fibroblast growth factor-
2: double-blind, randomized, controlled clinical trial. Circulation. 2002;105(7):788-93.
7. Henry TD, Annex BH, McKendall GR, Azrin MA, Lopez JJ, Giordano FJ, et al. The VIVA
trial: Vascular endothelial growth factor in Ischemia for Vascular Angiogenesis. Circulation.
2003;107(10):1359-65.
8. Sasich LD, Sukkari SR. The US FDAs withdrawal of the breast cancer indication for
Avastin (bevacizumab). Saudi Pharm J. 2012;20(4):381-5.
9. Pepper MS, Ferrara N, Orci L, Montesano R. Potent synergism between vascular
endothelial growth factor and basic fibroblast growth factor in the induction of angiogenesis in
vitro. Biochem Biophys Res Commun. 1992;189(2):824-31.
10. Goto F, Goto K, Weindel K, Folkman J. Synergistic effects of vascular endothelial growth
factor and basic fibroblast growth factor on the proliferation and cord formation of bovine capillary
endothelial cells within collagen gels. Lab Invest. 1993;69(5):508-17.
11. Asahara T, Bauters C, Zheng LP, Takeshita S, Bunting S, Ferrara N, et al. Synergistic effect
of vascular endothelial growth factor and basic fibroblast growth factor on angiogenesis in vivo.
Circulation. 1995;92(9 Suppl):II365-71.

142
12. Dieci MV, Arnedos M, Andre F, Soria JC. Fibroblast growth factor receptor inhibitors as
a cancer treatment: from a biologic rationale to medical perspectives. Cancer Discov.
2013;3(3):264-79.
13. Sobhani N, Ianza A, D'Angelo A, Roviello G, Giudici F, Bortul M, et al. Current Status of
Fibroblast Growth Factor Receptor-Targeted Therapies in Breast Cancer. Cells. 2018;7(7).
14. Semrad TJ, Kim EJ, Tanaka MS, Sands J, Roberts C, Burich RA, et al. Phase II Study of
Dovitinib in Patients Progressing on Anti-Vascular Endothelial Growth Factor Therapy. Cancer
Treat Res Commun. 2017;10:21-6.
15. Chambard JC, Lefloch R, Pouyssegur J, Lenormand P. ERK implication in cell cycle
regulation. Biochim Biophys Acta. 2007;1773(8):1299-310.
16. Yu JS, Cui W. Proliferation, survival and metabolism: the role of PI3K/AKT/mTOR
signalling in pluripotency and cell fate determination. Development. 2016;143(17):3050-60.
17. Franke TF, Yang SI, Chan TO, Datta K, Kazlauskas A, Morrison DK, et al. The protein
kinase encoded by the Akt proto-oncogene is a target of the PDGF-activated phosphatidylinositol
3-kinase. Cell. 1995;81(5):727-36.
18. Yao R, Cooper GM. Requirement for phosphatidylinositol-3 kinase in the prevention of
apoptosis by nerve growth factor. Science. 1995;267(5206):2003-6.
19. Gerber HP, McMurtrey A, Kowalski J, Yan M, Keyt BA, Dixit V, et al. Vascular
endothelial growth factor regulates endothelial cell survival through the phosphatidylinositol 3'-
kinase/Akt signal transduction pathway. Requirement for Flk-1/KDR activation. J Biol Chem.
1998;273(46):30336-43.
20. Alon T, Hemo I, Itin A, Pe'er J, Stone J, Keshet E. Vascular endothelial growth factor acts
as a survival factor for newly formed retinal vessels and has implications for retinopathy of
prematurity. Nat Med. 1995;1(10):1024-8.
21. Chen J, Somanath PR, Razorenova O, Chen WS, Hay N, Bornstein P, et al. Akt1 regulates
pathological angiogenesis, vascular maturation and permeability in vivo. Nat Med.
2005;11(11):1188-96.
22. Mebratu Y, Tesfaigzi Y. How ERK1/2 activation controls cell proliferation and cell death:
Is subcellular localization the answer? Cell Cycle. 2009;8(8):1168-75.
23. Lu Z, Xu S. ERK1/2 MAP kinases in cell survival and apoptosis. IUBMB Life.
2006;58(11):621-31.

143
24. Shiojima I, Walsh K. Role of Akt signaling in vascular homeostasis and angiogenesis. Circ
Res. 2002;90(12):1243-50.
25. Morales-Ruiz M, Fulton D, Sowa G, Languino LR, Fujio Y, Walsh K, et al. Vascular
endothelial growth factor-stimulated actin reorganization and migration of endothelial cells is
regulated via the serine/threonine kinase Akt. Circ Res. 2000;86(8):892-6.
26. Huang C, Jacobson K, Schaller MD. MAP kinases and cell migration. J Cell Sci.
2004;117(Pt 20):4619-28.
27. Mac Gabhann F, Popel AS. Targeting neuropilin-1 to inhibit VEGF signaling in cancer:
Comparison of therapeutic approaches. PLoS Comput Biol. 2006;2(12):e180.
28. Stefanini MO, Wu FT, Mac Gabhann F, Popel AS. Increase of plasma VEGF after
intravenous administration of bevacizumab is predicted by a pharmacokinetic model. Cancer Res.
2010;70(23):9886-94.
29. Filion RJ, Popel AS. Intracoronary administration of FGF-2: a computational model of
myocardial deposition and retention. Am J Physiol Heart Circ Physiol. 2005;288(1):H263-79.
30. Padera R, Venkataraman G, Berry D, Godavarti R, Sasisekharan R. FGF-2/fibroblast
growth factor receptor/heparin-like glycosaminoglycan interactions: a compensation model for
FGF-2 signaling. FASEB J. 1999;13(13):1677-87.
31. Tong S, Yuan F. Dose response of angiogenesis to basic fibroblast growth factor in rat
corneal pocket assay: II. Numerical simulations. Microvasc Res. 2008;75(1):16-24.
32. Hendrata M, Sudiono J. A Computational Model for Investigating Tumor Apoptosis
Induced by Mesenchymal Stem Cell-Derived Secretome. Comput Math Methods Med.
2016;2016:4910603.
33. Kanodia J, Chai D, Vollmer J, Kim J, Raue A, Finn G, et al. Deciphering the mechanism
behind Fibroblast Growth Factor (FGF) induced biphasic signal-response profiles. Cell Commun
Signal. 2014;12:34.
34. Wu Q, Finley SD. Predictive model identifies strategies to enhance TSP1-mediated
apoptosis signaling. Cell Commun Signal. 2017;15(1):53.
35. Song M, Finley SD. Mechanistic insight into activation of MAPK signaling by pro-
angiogenic factors. BMC Syst Biol. 2018;12(1):145.
36. Tan WH, Popel AS, Mac Gabhann F. Computational model of VEGFR2 pathway to ERK
activation and modulation through receptor trafficking. Cell Signal. 2013;25(12):2496-510.

144
37. Filion RJ, Popel AS. A reaction-diffusion model of basic fibroblast growth factor
interactions with cell surface receptors. Ann Biomed Eng. 2004;32(5):645-63.
38. Tan WH, Popel AS, Mac Gabhann F. Computational Model of Gab1/2-Dependent
VEGFR2 Pathway to Akt Activation. PLoS One. 2013;8(6):e67438.
39. Song M, Finley SD. ERK and Akt exhibit distinct signaling responses following
stimulation by pro-angiogenic factors. Cell Commun Signal. 2020;18(1):114.
40. Friedl P, Tatje D, Czpla R. An optimized culture medium for human vascular endothelial
cells from umbilical cord veins. Cytotechnology. 1989;2(3):171-9.
41. Weddell JC, Chen S, Imoukhuede PI. VEGFR1 promotes cell migration and proliferation
through PLCγ and PI3K pathways. NPJ Syst Biol Appl. 2018;4:1.
42. Norton KA, Popel AS. Effects of endothelial cell proliferation and migration rates in a
computational model of sprouting angiogenesis. Sci Rep. 2016;6:36992.
43. Tong S, Yuan F. Numerical simulations of angiogenesis in the cornea. Microvasc Res.
2001;61(1):14-27.
44. Lovett M, Lee K, Edwards A, Kaplan DL. Vascularization strategies for tissue engineering.
Tissue Eng Part B Rev. 2009;15(3):353-70.
45. Korc M, Friesel RE. The role of fibroblast growth factors in tumor growth. Curr Cancer
Drug Targets. 2009;9(5):639-51.
46. Zheng XM, Koh GY, Jackson T. 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 Cont Dyn-B.
2013;18(4):1109-54.
47. Lassoued W, Murphy D, Tsai J, Oueslati R, Thurston G, Lee WM. Effect of VEGF and
VEGF Trap on vascular endothelial cell signaling in tumors. Cancer Biol Ther. 2010;10(12):1326-
33.
48. Ibrahimi OA, Zhang F, Hrstka SC, Mohammadi M, Linhardt RJ. Kinetic model for FGF,
FGFR, and proteoglycan signal transduction complex assembly. Biochemistry. 2004;43(16):4724-
30.
49. Chang LK, Garcia-Cardena G, Farnebo F, Fannon M, Chen EJ, Butterfield C, et al. Dose-
dependent response of FGF-2 for lymphangiogenesis. Proc Natl Acad Sci U S A.
2004;101(32):11658-63.

145
50. Bai Y, Leng Y, Yin G, Pu X, Huang Z, Liao X, et al. Effects of combinations of BMP-2
with FGF-2 and/or VEGF on HUVECs angiogenesis in vitro and CAM angiogenesis in vivo. Cell
Tissue Res. 2014;356(1):109-21.
51. Zheng J, Wen Y, Song Y, Wang K, Chen DB, Magness RR. Activation of multiple
signaling pathways is critical for fibroblast growth factor 2- and vascular endothelial growth
factor-stimulated ovine fetoplacental endothelial cell proliferation. Biol Reprod. 2008;78(1):143-
50.
52. Bruns AF, Herbert SP, Odell AF, Jopling HM, Hooper NM, Zachary IC, et al. Ligand-
stimulated VEGFR2 signaling is regulated by co-ordinated trafficking and proteolysis. Traffic.
2010;11(1):161-74.
53. Lamalice L, Houle F, Huot J. Phosphorylation of Tyr1214 within VEGFR-2 triggers the
recruitment of Nck and activation of Fyn leading to SAPK2/p38 activation and endothelial cell
migration in response to VEGF. J Biol Chem. 2006;281(45):34009-20.
54. Chabot C, Spring K, Gratton JP, Elchebly M, Royal I. New role for the protein tyrosine
phosphatase DEP-1 in Akt activation and endothelial cell survival. Mol Cell Biol. 2009;29(1):241-
53.
55. Zhang Z, Neiva KG, Lingen MW, Ellis LM, Nor JE. VEGF-dependent tumor angiogenesis
requires inverse and reciprocal regulation of VEGFR1 and VEGFR2. Cell Death Differ.
2010;17(3):499-512.
56. Chang YS, Munn LL, Hillsley MV, Dull RO, Yuan J, Lakshminarayanan S, et al. Effect
of vascular endothelial growth factor on cultured endothelial cell monolayer transport properties.
Microvasc Res. 2000;59(2):265-77.
57. Marino S, Hogue IB, Ray CJ, Kirschner DE. A methodology for performing global
uncertainty and sensitivity analysis in systems biology. J Theor Biol. 2008;254(1):178-96.
58. Rikitake Y, Kawashima S, Yamashita T, Ueyama T, Ishido S, Hotta H, et al.
Lysophosphatidylcholine inhibits endothelial cell migration and proliferation via inhibition of the
extracellular signal-regulated kinase pathway. Arterioscler Thromb Vasc Biol. 2000;20(4):1006-
12.
59. Iadevaia S, Lu Y, Morales FC, Mills GB, Ram PT. Identification of optimal drug
combinations targeting cellular networks: integrating phospho-proteomics and computational
network analysis. Cancer Res. 2010;70(17):6704-14.
60. Olsson AK, Dimberg A, Kreuger J, Claesson-Welsh L. VEGF receptor signalling - in
control of vascular function. Nat Rev Mol Cell Biol. 2006;7(5):359-71.

146
61. Dupree MA, Pollack SR, Levine EM, Laurencin CT. Fibroblast growth factor 2 induced
proliferation in osteoblasts and bone marrow stromal cells: a whole cell model. Biophys J.
2006;91(8):3097-112.
62. Zhao B, Zhang C, Forsten-Williams K, Zhang J, Fannon M. Endothelial cell capture of
heparin-binding growth factors under flow. PLoS Comput Biol. 2010;6(10):e1000971.
63. Imoukhuede PI, Popel AS. Quantitative fluorescent profiling of VEGFRs reveals tumor
cell and endothelial cell heterogeneity in breast cancer xenografts. Cancer Med. 2014;3(2):225-
44.
64. Shin M, Beane TJ, Quillien A, Male I, Zhu LJ, Lawson ND. Vegfa signals through ERK
to promote angiogenesis, but not artery differentiation. Development. 2016;143(20):3796-805.
65. Mavria G, Vercoulen Y, Yeo M, Paterson H, Karasarides M, Marais R, et al. ERK-MAPK
signaling opposes Rho-kinase to promote endothelial cell survival and sprouting during
angiogenesis. Cancer Cell. 2006;9(1):33-44.
66. Nagasawa-Masuda A, Terai K. ERK activation in endothelial cells is a novel marker during
neovasculogenesis. Genes Cells. 2016;21(11):1164-75.
67. Zheng J, Bird IM, Melsaether AN, Magness RR. Activation of the mitogen-activated
protein kinase cascade is necessary but not sufficient for basic fibroblast growth factor- and
epidermal growth factor-stimulated expression of endothelial nitric oxide synthase in ovine
fetoplacental artery endothelial cells. Endocrinology. 1999;140(3):1399-407.
68. Gourlaouen M, Welti JC, Vasudev NS, Reynolds AR. Essential role for endocytosis in the
growth factor-stimulated activation of ERK1/2 in endothelial cells. J Biol Chem.
2013;288(11):7467-80.
69. Edelstein SJ, Stefan MI, Le Novere N. Ligand depletion in vivo modulates the dynamic
range and cooperativity of signal transduction. PLoS One. 2010;5(1):e8449.
70. Saucerman JJ, Bers DM. Calmodulin mediates differential sensitivity of CaMKII and
calcineurin to local Ca2+ in cardiac myocytes. Biophys J. 2008;95(10):4597-612.
71. Romano DR, Pharris MC, Patel NM, Kinzer-Ursem TL. Competitive tuning: Competition's
role in setting the frequency-dependence of Ca2+-dependent proteins. PLoS Comput Biol.
2017;13(11):e1005820.
72. Kano MR, Morishita Y, Iwata C, Iwasaka S, Watabe T, Ouchi Y, et al. VEGF-A and FGF-
2 synergistically promote neoangiogenesis through enhancement of endogenous PDGF-B-
PDGFRbeta signaling. J Cell Sci. 2005;118(Pt 16):3759-68.

147
73. Murakami M, Nguyen LT, Hatanaka K, Schachterle W, Chen PY, Zhuang ZW, et al. FGF-
dependent regulation of VEGF receptor 2 expression in mice. J Clin Invest. 2011;121(7):2668-78.
74. Lampugnani MG, Orsenigo F, Gagliani MC, Tacchetti C, Dejana E. Vascular endothelial
cadherin controls VEGFR-2 internalization and signaling from intracellular compartments. J Cell
Biol. 2006;174(4):593-604.
75. Basagiannis D, Zografou S, Murphy C, Fotsis T, Morbidelli L, Ziche M, et al. VEGF
induces signalling and angiogenesis by directing VEGFR2 internalisation through
macropinocytosis. J Cell Sci. 2016;129(21):4091-104.
76. Sun J, Huang W, Yang SF, Zhang XP, Yu Q, Zhang ZQ, et al. Galphai1 and
Galphai3mediate VEGF-induced VEGFR2 endocytosis, signaling and angiogenesis. Theranostics.
2018;8(17):4695-709.
77. Sperinde GV, Nugent MA. Mechanisms of fibroblast growth factor 2 intracellular
processing: a kinetic analysis of the role of heparan sulfate proteoglycans. Biochemistry.
2000;39(13):3788-96.
78. Zygmunt M, Herr F, Munstedt K, Lang U, Liang OD. Angiogenesis and vasculogenesis in
pregnancy. Eur J Obstet Gynecol Reprod Biol. 2003;110 Suppl 1:S10-8.
79. Chen DB, Zheng J. Regulation of placental angiogenesis. Microcirculation. 2014;21(1):15-
25.
80. Tonnesen MG, Feng X, Clark RA. Angiogenesis in wound healing. J Investig Dermatol
Symp Proc. 2000;5(1):40-6.
81. Demidova-Rice TN, Durham JT, Herman IM. Wound Healing Angiogenesis: Innovations
and Challenges in Acute and Chronic Wound Healing. Adv Wound Care (New Rochelle).
2012;1(1):17-22.
82. Kerbel RS. Tumor angiogenesis. N Engl J Med. 2008;358(19):2039-49.
83. Keating GM. Bevacizumab: a review of its use in advanced cancer. Drugs.
2014;74(16):1891-925.
84. Vasudev NS, Reynolds AR. Anti-angiogenic therapy for cancer: current progress,
unresolved questions and future directions. Angiogenesis. 2014;17(3):471-94.
85. Li D, Finley SD. The impact of tumor receptor heterogeneity on the response to anti-
angiogenic cancer treatment. Integr Biol (Camb). 2018;10(4):253-69.

148
86. Li D, Finley SD. Exploring the Extracellular Regulation of the Tumor Angiogenic
Interaction Network Using a Systems Biology Model. Frontiers in Physiology. 2019;10.
87. Rubin DB, Drab EA, Bauer KD. Endothelial cell subpopulations in vitro: cell volume, cell
cycle, and radiosensitivity. J Appl Physiol (1985). 1989;67(4):1585-90.
88. Maffucci T, Raimondi C, Abu-Hayyeh S, Dominguez V, Sala G, Zachary I, et al. A
phosphoinositide 3-kinase/phospholipase Cgamma1 pathway regulates fibroblast growth factor-
induced capillary tube formation. PLoS One. 2009;4(12):e8285.
89. Pisanti S, Picardi P, Prota L, Proto MC, Laezza C, McGuire PG, et al. Genetic and
pharmacologic inactivation of cannabinoid CB1 receptor inhibits angiogenesis. Blood.
2011;117(20):5541-50.
90. Schneeweis C, Grafe M, Bungenstock A, Spencer-Hansch C, Fleck E, Goetze S. Chronic
CRP-exposure inhibits VEGF-induced endothelial cell migration. J Atheroscler Thromb.
2010;17(2):203-12.
91. Schoeberl B, Pace EA, Fitzgerald JB, Harms BD, Xu L, Nie L, et al. Therapeutically
targeting ErbB3: a key node in ligand-induced activation of the ErbB receptor-PI3K axis. Sci
Signal. 2009;2(77):ra31.
92. Ratajska A, Torry RJ, Kitten GT, Kolker SJ, Tomanek RJ. Modulation of cell migration
and vessel formation by vascular endothelial growth factor and basic fibroblast growth factor in
cultured embryonic heart. Dev Dyn. 1995;203(4):399-407.
93. Matsui J, Wakabayashi T, Asada M, Yoshimatsu K, Okada M. Stem cell factor/c-kit
signaling promotes the survival, migration, and capillary tube formation of human umbilical vein
endothelial cells. J Biol Chem. 2004;279(18):18600-7.
94. Nakashio A, Fujita N, Tsuruo T. Topotecan inhibits VEGF- and bFGF-induced vascular
endothelial cell migration via downregulation of the PI3K-Akt signaling pathway. Int J Cancer.
2002;98(1):36-41.
95. Soldi R, Mitola S, Strasly M, Defilippi P, Tarone G, Bussolino F. Role of alphavbeta3
integrin in the activation of vascular endothelial growth factor receptor-2. EMBO J.
1999;18(4):882-92.
96. Dong F, Tian H, Yan S, Li L, Dong X, Wang F, et al. Dihydroartemisinin inhibits
endothelial cell proliferation through the suppression of the ERK signaling pathway. Int J Mol
Med. 2015;35(5):1381-7.

149
97. Hoeflich KP, O'Brien C, Boyd Z, Cavet G, Guerrero S, Jung K, et al. In vivo antitumor
activity of MEK and phosphatidylinositol 3-kinase inhibitors in basal-like breast cancer models.
Clin Cancer Res. 2009;15(14):4649-64.
98. Teven CM, Farina EM, Rivas J, Reid RR. Fibroblast growth factor (FGF) signaling in
development and skeletal diseases. Genes Dis. 2014;1(2):199-213.
99. Clark JA, Black AR, Leontieva OV, Frey MR, Pysz MA, Kunneva L, et al. Involvement
of the ERK signaling cascade in protein kinase C-mediated cell cycle arrest in intestinal epithelial
cells. J Biol Chem. 2004;279(10):9233-47.
100. Thomas SM, DeMarco M, D'Arcangelo G, Halegoua S, Brugge JS. Ras is essential for
nerve growth factor- and phorbol ester-induced tyrosine phosphorylation of MAP kinases. Cell.
1992;68(6):1031-40.
101. Kawakami Y, Kitaura J, Yao L, McHenry RW, Newton AC, Kang S, et al. A Ras activation
pathway dependent on Syk phosphorylation of protein kinase C. Proc Natl Acad Sci U S A.
2003;100(16):9470-5.
102. Ueda Y, Hirai S, Osada S, Suzuki A, Mizuno K, Ohno S. Protein kinase C activates the
MEK-ERK pathway in a manner independent of Ras and dependent on Raf. J Biol Chem.
1996;271(38):23512-9.
103. Ho PD, Zechner DK, He H, Dillmann WH, Glembotski CC, McDonough PM. The Raf-
MEK-ERK cascade represents a common pathway for alteration of intracellular calcium by Ras
and protein kinase C in cardiac myocytes. J Biol Chem. 1998;273(34):21730-5.
104. Corbit KC, Trakul N, Eves EM, Diaz B, Marshall M, Rosner MR. Activation of Raf-1
signaling by protein kinase C through a mechanism involving Raf kinase inhibitory protein. J Biol
Chem. 2003;278(15):13061-8.
105. Wen-Sheng W. Protein kinase C alpha trigger Ras and Raf-independent MEK/ERK
activation for TPA-induced growth inhibition of human hepatoma cell HepG2. Cancer Lett.
2006;239(1):27-35.
106. Koch S, Claesson-Welsh L. Signal transduction by vascular endothelial growth factor
receptors. Cold Spring Harb Perspect Med. 2012;2(7):a006502.
107. Mac Gabhann F, Popel AS. Model of competitive binding of vascular endothelial growth
factor and placental growth factor to VEGF receptors on endothelial cells. Am J Physiol Heart
Circ Physiol. 2004;286(1):H153-64.

150
108. Weddell JC, Imoukhuede PI. Quantitative characterization of cellular membrane-receptor
heterogeneity through statistical and computational modeling. PLoS One. 2014;9(5):e97271.
109. Zimmermann S, Moelling K. Phosphorylation and regulation of Raf by Akt (protein kinase
B). Science. 1999;286(5445):1741-4.
110. Moelling K, Schad K, Bosse M, Zimmermann S, Schweneker M. Regulation of Raf-Akt
Cross-talk. J Biol Chem. 2002;277(34):31099-106.
111. Chiu D, Ma K, Scott A, Duronio V. Acute activation of Erk1/Erk2 and protein kinase B/akt
proceed by independent pathways in multiple cell types. FEBS J. 2005;272(17):4372-84.
112. C VS, Babar SM, Song EJ, Oh E, Yoo YS. Kinetic analysis of the MAPK and PI3K/Akt
signaling pathways. Mol Cells. 2008;25(3):397-406.
113. Adlung L, Kar S, Wagner MC, She B, Chakraborty S, Bao J, et al. Protein abundance of
AKT and ERK pathway components governs cell type-specific regulation of proliferation. Mol
Syst Biol. 2017;13(1):904.
114. Ammann KR, DeCook KJ, Tran PL, Merkle VM, Wong PK, Slepian MJ. Collective cell
migration of smooth muscle and endothelial cells: impact of injury versus non-injury stimuli. J
Biol Eng. 2015;9:19.
115. Braun A, Caesar NM, Dang K, Myers KA. High-resolution Time-lapse Imaging and
Automated Analysis of Microtubule Dynamics in Living Human Umbilical Vein Endothelial
Cells. J Vis Exp. 2016(114).
116. Koon YL, Zhang S, Rahmat MB, Koh CG, Chiam KH. Enhanced Delta-Notch Lateral
Inhibition Model Incorporating Intracellular Notch Heterogeneity and Tension-Dependent Rate of
Delta-Notch Binding that Reproduces Sprouting Angiogenesis Patterns. Sci Rep. 2018;8(1):9519.
117. Walpole J, Mac Gabhann F, Peirce SM, Chappell JC. Agent-based computational model
of retinal angiogenesis simulates microvascular network morphology as a function of pericyte
coverage. Microcirculation. 2017;24(8).
118. Blanco R, Gerhardt H. VEGF and Notch in tip and stalk cell selection. Cold Spring Harb
Perspect Med. 2013;3(1):a006569.
119. Heiss M, Hellström M, Kalén M, May T, Weber H, Hecker M, et al. Endothelial cell
spheroids as a versatile tool to study angiogenesis in vitro. FASEB J. 2015;29(7):3076-84.

151
120. Liebler SS, Feldner A, Adam MG, Korff T, Augustin HG, Fischer A. No evidence for a
functional role of bi-directional Notch signaling during angiogenesis. PLoS One.
2012;7(12):e53074.
121. West DC, Rees CG, Duchesne L, Patey SJ, Terry CJ, Turnbull JE, et al. Interactions of
multiple heparin binding growth factors with neuropilin-1 and potentiation of the activity of
fibroblast growth factor-2. J Biol Chem. 2005;280(14):13457-64.
122. Jih YJ, Lien WH, Tsai WC, Yang GW, Li C, Wu LW. Distinct regulation of genes by
bFGF and VEGF-A in endothelial cells. Angiogenesis. 2001;4(4):313-21.
123. Brütsch R, Liebler SS, Wüstehube J, Bartol A, Herberich SE, Adam MG, et al. Integrin
cytoplasmic domain-associated protein-1 attenuates sprouting angiogenesis. Circ Res.
2010;107(5):592-601.
124. Herberich SE, Klose R, Moll I, Yang WJ, Wüstehube-Lausch J, Fischer A. ANKS1B
Interacts with the Cerebral Cavernous Malformation Protein-1 and Controls Endothelial
Permeability but Not Sprouting Angiogenesis. PLoS One. 2015;10(12):e0145304.
125. Nacak TG, Alajati A, Leptien K, Fulda C, Weber H, Miki T, et al. The BTB-Kelch protein
KLEIP controls endothelial migration and sprouting angiogenesis. Circ Res. 2007;100(8):1155-
63.
126. Bazzazi H, Popel AS. Computational investigation of sphingosine kinase 1 (SphK1) and
calcium dependent ERK1/2 activation downstream of VEGFR2 in endothelial cells. PLoS Comput
Biol. 2017;13(2):e1005332.
127. Zhang Y, Kontos CD, Annex BH, Popel AS. Angiopoietin-Tie Signaling Pathway in
Endothelial Cells: A Computational Model. iScience. 2019;20:497-511.
128. Eklund L, Saharinen P. Angiopoietin signaling in the vasculature. Exp Cell Res.
2013;319(9):1271-80.
129. Loughna S, Sato TN. Angiopoietin and Tie signaling pathways in vascular development.
Matrix Biol. 2001;20(5-6):319-25.
130. Saharinen P, Eklund L, Alitalo K. Therapeutic targeting of the angiopoietin-TIE pathway.
Nat Rev Drug Discov. 2017;16(9):635-61.
131. Wu Q, Finley SD. Mathematical Model Predicts Effective Strategies to Inhibit VEGF-
eNOS Signaling. J Clin Med. 2020;9(5).

152
132. Bazzazi H, Zhang Y, Jafarnejad M, Isenberg JS, Annex BH, Popel AS. Computer
Simulation of TSP1 Inhibition of VEGF-Akt-eNOS: An Angiogenesis Triple Threat. Front
Physiol. 2018;9:644.
133. Fukumura D, Gohongi T, Kadambi A, Izumi Y, Ang J, Yun CO, et al. Predominant role of
endothelial nitric oxide synthase in vascular endothelial growth factor-induced angiogenesis and
vascular permeability. Proc Natl Acad Sci U S A. 2001;98(5):2604-9.
134. Duda DG, Fukumura D, Jain RK. Role of eNOS in neovascularization: NO for endothelial
progenitor cells. Trends Mol Med. 2004;10(4):143-5.
135. Vega R, Carretero M, Travasso RDM, Bonilla LL. Notch signaling and taxis mechanisms
regulate early stage angiogenesis: A mathematical and computational model. PLoS Comput Biol.
2020;16(1):e1006919.
136. Kühn C, Checa S. Computational Modeling to Quantify the Contributions of VEGFR1,
VEGFR2, and Lateral Inhibition in Sprouting Angiogenesis. Front Physiol. 2019;10:288.
137. Weinstein N, Mendoza L, Gitler I, Klapp J. A Network Model to Explore the Effect of the
Micro-environment on Endothelial Cell Behavior during Angiogenesis. Front Physiol. 2017;8:960.
138. Nacev BA, Liu JO. Synergistic inhibition of endothelial cell proliferation, tube formation,
and sprouting by cyclosporin A and itraconazole. PLoS One. 2011;6(9):e24793.
139. Liao H, He H, Chen Y, Zeng F, Huang J, Wu L. Effects of long-term serial cell passaging
on cell spreading, migration, and cell-surface ultrastructures of cultured vascular endothelial cells.
Cytotechnology. 2014;66(2):229-38.
140. Roy M, Finley SD. Metabolic reprogramming dynamics in tumor spheroids: Insights from
a multicellular, multiscale model. PLoS Comput Biol. 2019;15(6):e1007053.
141. Letort G, Montagud A, Stoll G, Heiland R, Barillot E, Macklin P, et al. PhysiBoSS: a multi-
scale agent-based modelling framework integrating physical dimension and cell signalling.
Bioinformatics. 2019;35(7):1188-96.
142. Athale C, Mansury Y, Deisboeck TS. Simulating the impact of a molecular 'decision-
process' on cellular phenotype and multicellular patterns in brain tumors. J Theor Biol.
2005;233(4):469-81.
143. Zhang L, Athale CA, Deisboeck TS. Development of a three-dimensional multiscale agent-
based tumor model: simulating gene-protein interaction profiles, cell phenotypes and multicellular
patterns in brain cancer. J Theor Biol. 2007;244(1):96-107.

153
144. Page DJ, Clarkin CE, Mani R, Khan NA, Dawson JI, Evans ND. Injectable nanoclay gels
for angiogenesis. Acta Biomater. 2019.
145. Wolfe A, O'Clair B, Groppi VE, McEwen DP. Pharmacologic characterization of a kinetic
in vitro human co-culture angiogenesis model using clinically relevant compounds. J Biomol
Screen. 2013;18(10):1234-45.
146. Dellian M, Witwer BP, Salehi HA, Yuan F, Jain RK. Quantitation and physiological
characterization of angiogenic vessels in mice: effect of basic fibroblast growth factor, vascular
endothelial growth factor/vascular permeability factor, and host microenvironment. Am J Pathol.
1996;149(1):59-71.
147. Bazzazi H, Isenberg JS, Popel AS. Inhibition of VEGFR2 Activation and Its Downstream
Signaling to ERK1/2 and Calcium by Thrombospondin-1 (TSP1):. Front Physiol. 2017;8:48.
148. Mack JJ, Iruela-Arispe ML. NOTCH regulation of the endothelial cell phenotype. Curr
Opin Hematol. 2018;25(3):212-8.
149. Meyer RD, Mohammadi M, Rahimi N. A single amino acid substitution in the activation
loop defines the decoy characteristic of VEGFR-1/FLT-1. J Biol Chem. 2006;281(2):867-75.

Abstract (if available)

Abstract Angiogenesis, the formation of new blood capillaries from pre-existing blood vessels, plays an important role in the survival of tissues, because vessels carry blood throughout the body to provide oxygen and nutrients required by the resident cells. Angiogenesis has received great attention, as angiogenic-based strategies are beneficial in many contexts, including in both tissue engineering and cancer treatment. Specifically, promoting new blood vessel formation to ensure adequate vasculature is critical for the survival of artificial tissues; on the other hand, inhibiting angiogenesis is an important strategy for cancer treatment, as vessels not only deliver nutrients to support cancer growth, they also provide routes for metastasis. However, the outcomes of these pro- and anti-angiogenic therapies are not all effective. Lack of therapeutic response, low efficacy, or drug resistance have been observed. One major obstacle for angiogenesis-based therapies is that the angiogenic response is driven by signals integrated from multiple relevant signaling pathways. Targeting one pathway could be insufficient, as alternative pathways may compensate, diminishing the overall effect of the treatment strategy. For example, there are multiple proteins that promote signaling needed for angiogenesis. Fibroblast growth factor (FGF) and vascular endothelial growth factor (VEGF), which are two major pro-angiogenic factors, mediate mitogen-activated protein kinase (MAPK) and phosphatidylinositol 3-kinase/protein kinase B (PI3K/Akt) pathways. Signaling through these pathways leads to activated extracellular regulated kinase (ERK) and Akt, which are important in angiogenic responses such as endothelial cell survival, proliferation, and migration. However, there is a limited understanding of how these promoters combine together to stimulate angiogenesis. Therefore, understanding FGF- and VEGF-mediated angiogenic signaling systematically can be beneficial in improving current angiogenic strategies. ? We seek to achieve optimal angiogenic outcomes and inform the development of effective angiogenic therapies by constructing experimentally validated mathematical models. In this work, I present a series of mechanistic models to address questions regarding the interactions of FGF- and VEGF-mediated angiogenic signaling in endothelial cells, the effects of the crosstalk in angiogenic cellular responses mediated by these factors, and potential targets for angiogenic-based therapies.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Understanding anti-angiogenic signaling and treatment for cancer through mechanistic modeling

PDF

Mechanistic modeling of angiogenic factors network and cancer therapy

PDF

Exploring the role of vegfr transcripts during vascular development in the zebrafish embryo

PDF

Mechanistic model of chimeric antigen receptor T cell activation

PDF

Insights from a computational approach to natural killer cell activation

PDF

Targeting vessel maturation: an anti-angiogenesis based cancer therapy

PDF

Understanding the impact of cell-to-cell variability on intracellular signaling in CAR cells through mathematical modeling

PDF

Membrane-bound regulation of hematopoietic stem cells

PDF

Metabolic profiling of single hematopoietic stem cells for developing novel ex vivo culture strategies

PDF

Metabolomic and proteomic approaches to understanding senescence and aging in mammary epithelial cells

PDF

Importance of alpha-catenin and cell-cell interactions in hair follicle stem cells homeostasis

PDF

Chitosan nanoparticle mediated oral drug delivery of kidney targeting micelles for polycystic kidney disease

PDF

Physiologically-based modeling to predict monoclonal antibody pharmacokinetics in humans from physiochemical properties

PDF

Proliferation and maturation events in second heart field cells during cardiovascular development activated by the Delta like ligand-4 mediated notch signaling

PDF

System biology approaches to cancer and diabetes metabolism

PDF

Improving antitumor efficacy of chimeric antigen receptor-engineered immune cell therapy with synthetic biology and combination therapy approaches

PDF

Nonlinear dynamical modeling of single neurons and its application to analysis of long-term potentiation (LTP)

PDF

Elements of photoreceptor homeostasis: investigating phenotypic manifestations and susceptibility to photoreceptor degeneration in genetic knockout models for retinal disease

PDF

Kaposi's sarcoma associated herpes-virus induces cellular proxi expression to modulate host gene expression that benefits viral infection and oncogenesis

PDF

Modulation of transcription and receptor function with synthetic small molecules and multi-finctional integrin ligands

Asset Metadata

Creator Song, Min (author)

Core Title Mechanistic investigation of pro-angiogenic signaling and endothelial sprouting mediated by FGF and VEGF

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Biomedical Engineering

Degree Conferral Date 2021-08

Publication Date 07/22/2021

Defense Date 04/19/2021

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag angiogenesis,cell signaling,computational modeling,endothelial cell,OAI-PMH Harvest,sprouting

Format application/pdf (imt)

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Finley, Stacey (committee chair), D’Argenio, David (committee member), Fraser, Scott (committee member), Newton, Paul (committee member), Shen, Keyue (committee member)

Creator Email m.song4746@gmail.com

Permanent Link (DOI) https://doi.org/10.25549/usctheses-oUC15614230

Unique identifier UC15614230

Legacy Identifier etd-SongMin-9818

Document Type Dissertation

Format application/pdf (imt)

Rights Song, Min

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the author, as the original true and official version of the work, but does not grant the reader permission to use the work if the desired use is covered by copyright. It is the author, as rights holder, who must provide use permission if such use is covered by copyright. The original signature page accompanying the original submission of the work to the USC Libraries is retained by the USC Libraries and a copy of it may be obtained by authorized requesters contacting the repository e-mail address given.

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA

Repository Email cisadmin@lib.usc.edu