The splice isoforms of vascular endothelial growth A (VEGF) each have different affinities for the extracellular matrix (ECM) and the coreceptor NRP1, which leads to distinct vascular phenotypes in model systems expressing only a single VEGF isoform. ECM-immobilized VEGF can bind to and activate VEGF receptor 2 (VEGFR2) directly, with a different pattern of site-specific phosphorylation than diffusible VEGF. To date, the way in which ECM binding alters the distribution of isoforms of VEGF and of the related placental growth factor (PlGF) in the body and resulting angiogenic signaling is not well-understood. Here, we extend our previous validated cell-level computational model of VEGFR2 ligation, intracellular trafficking, and site-specific phosphorylation, which captured differences in signaling by soluble and immobilized VEGF, to a multi-scale whole-body framework. This computational systems pharmacology model captures the ability of the ECM to regulate isoform-specific growth factor distribution distinctly for VEGF and PlGF, and to buffer free VEGF and PlGF levels in tissue. We show that binding of immobilized growth factor to VEGF receptors, both on endothelial cells and soluble VEGFR1, is likely important to signaling in vivo. Additionally, our model predicts that VEGF isoform-specific properties lead to distinct profiles of VEGFR1 and VEGFR2 binding and VEGFR2 site-specific phosphorylation in vivo, mediated by Neuropilin-1. These predicted signaling changes mirror those observed in murine systems expressing single VEGF isoforms. Simulations predict that, contrary to the ‘ligand-shifting hypothesis, ’ VEGF and PlGF do not compete for receptor binding at physiological concentrations, though PlGF is predicted to slightly increase VEGFR2 phosphorylation when over-expressed by 10-fold. These results are critical to design of appropriate therapeutic strategies to control VEGF availability and signaling in regenerative medicine applications.
The primary objectives of this study were: (1) to predict the distribution of VEGF and PlGF within the body, (2) to understand the effect of VEGF and PlGF on the balance of VEGFR1 and VEGFR2 ligation and VEGFR2 phosphorylation, (3) to quantify the effect of matrix-bound VEGF & PlGF binding to endothelial and soluble receptors on VEGFR signaling, and (4) to study the impact of changes in VEGF & PlGF isoform expression on absolute and relative VEGFR1 & VEGFR2 activation and site-specific phosphorylation of VEGFR2, as a result of isoform-specific matrix- and NRP1-binding properties, all within the context of a healthy human body. The computational systems pharmacology model developed in this study is based on previously-developed computational models of VEGF distribution and receptor binding in vivo. These models have included VEGF165, VEGF121, VEGFR1, VEGFR2, soluble VEGFR1 (sR1), NRP1, and sites in the interstitial matrix to which some growth factors and sR1 can bind [56–58]. The distribution of these proteins and their complexes has been examined in tissues of therapeutic interest (healthy or PAD calf [58], or tumor [59,60]), the blood, and non-diseased tissue (main body mass) [56,57], in humans or mice [61,62], incorporating transport between these compartments via vascular permeability and lymphatic drainage of tissues, and clearance of proteins from the plasma. By including multiple tissue compartments, we can compare quantities in a tissue of interest to those in the bulk of body tissue. In the present study, we greatly expand upon previous models to further capture the complexity of VEGF distribution and VEGF receptor activation in the body. For the first time, we include two isoforms of placental growth factor (PlGF1 & PlGF2), and the VEGF isoform VEGF189. Additionally, we account for binding of matrix-immobilized ligands in the endothelial basement membrane (EBM) to cell-surface receptors (VEGFR1 & VEGFR2), binding of immobilized ligands throughout the interstitial space to soluble sR1, and the ability of sR1, when sequestered in the interstitial matrix, to bind some VEGF isoforms. To capture these effects, we simulate receptor trafficking and VEGFR2 tyrosine site-specific phosphorylation following ligand binding or unbinding explicitly, implementing the reactions in a previously-developed in vitro computational model that captures differences in VEGFR2 phosphorylation following stimulation with soluble or matrix-bound VEGF165 [20]. Finally, we leverage recent measurements to update endothelial cell surface receptor densities [63].
To capture the pharmacokinetics of VEGF, PlGF, and sR1 distribution in the human body, we divide the body into three compartments: a healthy calf muscle (gastrocnemius + soleus muscles), blood, and the main body mass (the rest of the tissues), approximated with the properties of skeletal muscle (Fig 1A). Transport between compartments occurs via bi-directional vascular permeability and lymphatic drainage of tissues (into the blood), while growth factors and sR1 are cleared from the blood (via the liver and kidneys), using rates previously determined (S10 Table). Each tissue compartment includes physiological proportions of interstitial space, extracellular matrix (ECM), endothelium, other parenchymal cells, and basement membranes for both the endothelium and parenchyma (endothelial- EBM, and parenchymal- PBM). Within each tissue, we incorporate molecularly-detailed pharmacodynamics, including secretion into the interstitial space of VEGF and PlGF by parenchymal cells and sR1 by endothelial cells. In the interstitium, these diffusible proteins can then bind to heparan sulfate proteoglycans (HSPGs) in the ECM and basement membranes (see Fig 1B, S11 Table), bind to receptors on endothelial cells (ECs), or be removed from the compartment via physiological transport processes (Fig 1A). VEGF and PlGF isoforms have different affinities for matrix sites and for the coreceptor NRP1, which are included (Table 1), to account for isoform-specific ligand distribution and receptor activation. On the surface of and within endothelial cells, we simulate binding of sR1 to NRP1, binding of PlGF to VEGFR1 and/or NRP1, and binding of VEGF to VEGFR1, VEGFR2, and/or NRP1, based on the binding properties of each protein (summarized in Tables 1–3 and Fig 1E). Endothelial cell surface receptors are continually produced, internalized, recycled, and degraded, with trafficking rates that depend on ligation status and complex formation with NRP1 (Fig 1C). We include detailed VEGFR2 trafficking based on a previous in vitro computational model (S8 Table). Surface receptor production rates were tuned to match experimental measurements of cell surface receptor levels in human umbilical vein endothelial cells (Table 4). We also explicitly include phosphorylation and site-specific dephosphorylation of VEGFR2 (Fig 1D), which is dependent on receptor trafficking, with higher net activation at Y1214 than Y1175 on the cell surface, and higher Y1175 phosphorylation in early (Rab4/5) endosomes (S9 Table), as a result of differential dephosphorylation of Y1175 and Y1214 on the cell surface and in early endosomes [20]. This allows us to study phosphorylation explicitly, instead of using receptor occupancy as a surrogate, and to look at relative activation of downstream signaling pathways leading to proliferation (pY1175 via ERK1/2) and migration (pY1214 via p38). Due to the spatially-averaged nature of this model, gradients and heterogeneity in growth factor, soluble receptor, and cell surface receptor patterning are neglected. Instead, we examine the tissue-averaged behavior within the context of the human body. We neglect secretion of sR1 directly into the bloodstream, receptors present on the luminal side of ECs, and degradation of growth factors by proteases. All parameters are based on or fit to experimental data, either newly here or previously for other computational models. By building on previous modeling efforts, we have built more molecular detail into our models, while adding only a modest number of new parameters (indicated in bold in Tables 1–4 and S7 Table). To simulate the time-course of each molecular species in each tissue and the blood, this model includes 635 nonlinear ordinary differential equations that are solved simultaneously. The model equations can be found in S1 Equations. The full set of differential equations was solved in Fortran using the Livermore Solver for Ordinary Differential Equations with Automatic method switching for stiff and nonstiff problems (LSODA), on a laptop PC, with a relative error tolerance of 10−6.
The ligand secretion and receptor production rates necessary to hit baseline (healthy) targets had to be fit simultaneously, due to the highly non-linear nature of the system. At our baseline steady-state, the VEGF production rate is 0. 2830 molecules/myonuclear domain/s, the PlGF production rate is 0. 0146 molecules/myonuclear domain/s, and the sR1 production rate is 0. 0893 molecules/EC/s (see Table 4). The VEGF and sR1 production rates here are higher than previous estimates. This is unsurprising, given the changes in receptor levels, trafficking, and growth factor isoforms. Surprisingly, the PlGF production rate is lower than that for VEGF, despite a higher target plasma level (see Flux Analysis section for the mechanism by which this occurs). To illustrate the nonlinearity of our model, we perturbed each ligand secretion and receptor production rate slightly (2%), and examined changes in plasma ligand and tissue receptor levels. As shown in Fig 2A, plasma VEGF and tissue VEGFR2 are highly sensitive to changes in either VEGF secretion or VEGFR2 production in the main body mass, with changes of 11–25% per percent change in input. As VEGF levels increase, more VEGFR2 becomes occupied, internalized, and degraded, reducing VEGFR2 levels and decreasing VEGF consumption (Fig 2B and S1 Fig). Similarly, as VEGFR2 production increases, more VEGF is bound to VEGFR2, internalized, and degraded, reducing VEGF levels and thus increasing EC surface VEGFR2. This super-sensitivity was not present in previous models, where surface VEGFR2 levels were fixed (see S1 Fig). This new, emergent result suggests that, lacking upregulation of VEGFR2 in response to VEGF, VEGFR2 levels would be highly sensitive to even small fluctuations in local VEGF concentration (Fig 2), highlighting the importance of dynamic adjustments to ligand and receptor expression in vivo. In the calf muscle, perturbing VEGFR2 production has a large impact on EC surface VEGFR2, but little effect on plasma VEGF, due to the smaller size of the compartment. Changes in receptor production in one tissue compartment have little effect on receptor levels in the other tissue compartment. In this model, we assume the same rates for ligand production in both the healthy calf muscle and the main body mass. As such, perturbing the VEGF secretion rate (in both compartments) alters the receptor levels in both tissues (Fig 2). Due to differences in the geometric parameterizations of the calf and other tissues (S7 Table), using the same ligand secretion rates results in different interstitial VEGF, sR1, and PlGF levels (Fig 3D). We focus primarily on quantities measured in the “Main Body Mass” compartment, which, due to its larger size, represents the primary determinant of plasma VEGF, sR1, and PlGF levels. After establishing the secretion and production rates required to achieve basal targets, we next examined the steady-state distribution of VEGF, PlGF, and sR1. Having examined the distribution of VEGF, PlGF, and sR1, we next zoomed in to examine the effect of these proteins and their distributions on the binding and activation of endothelial VEGFR1 and VEGFR2 within healthy tissue. It is clear that the different proteins—ligands, soluble receptors, and co-receptors—regulating VEGFR1 and VEGFR2 activation do not act in isolation. Changes to any single feature affect the total multi-factor system in a way that is difficult to predict without the use of a computational model. Here, we perturb several interactions that are of interest therapeutically, and/or are included in this model for the first time. The most convincing evidence to date of differential signaling by VEGF isoforms is the distinct vascular phenotypes of mice or human tumors (implanted in mice) expressing only single isoforms of VEGF, with VEGF121-only tissues producing high diameter, sparsely branched networks, VEGF165-only tissue a relatively normal phenotype, and VEGF189-only tissues networks of thin, highly branched vessels. Endothelial cells isolated from these single isoform-expressing mice also display distinct signaling and behavior in cell culture [81]. It is assumed that similar regulation occurs in humans. To better understand VEGF isoform-specific signaling in the context of the human, as well as to qualitatively validate our model, we simulated expression of a single VEGF isoform in the human body. While no significant changes in VEGFR1 or VEGFR2 mRNA were observed in the muscle of mice expressing only VEGF120 [82] (equivalent to human VEGF121), we re-fit our model for each case, in order to maintain target ligand and receptor levels (S13 Table). The need for these changes in receptor production and ligand secretion rates may be a result of differences between humans and mice, or underlying compensation mechanisms and physiological changes in the engineered mice [82] not included in this model. Consistent with observations in mice, ligand distribution and VEGFR2 activation are more similar to wild type (baseline) in the VEGF165-only than the VEGF121-only or VEGF189-only cases (Fig 8A and 8B). Similar to the baseline case (Fig 5), where all three isoforms are expressed, with single VEGF isoform expression the ratio of migratory to proliferative signaling downstream of VEGFR2 (pY1214/pY1175) is predicted to increase with isoform length, paralleling the observed phenotypes (Fig 8C). The model’s ability to capture this trend provides qualitative validation of our isoform-specific signaling predictions in vivo. Interestingly, the model also predicts other changes, in free VEGF levels in tissue interstitium (Fig 8A) and in relative activation of VEGFR1 and VEGFR2 (Fig 8B and 8D, S1 File).
Our model predicts that, based on their binding properties and in vivo concentrations, PlGF and VEGF have distinct distributions within the body. PlGF2, binding to the ECM more strongly than VEGF, is bound to interstitial matrix sites at very high levels (~1 nM in tissue: soluble + ECM-bound + EC-bound predicted, Fig 3C), forming a large reservoir available for proteolytic release. Despite high tissue PlGF levels, our simulations predict that only about 30% of ligated EC surface VEGFR1 is bound to PlGF. As a result, while most VEGF removal from tissue is predicted to occur via binding to endothelial receptors, only 25% of PlGF was predicted to bind to and be subsequently degraded by endothelial cells. PlGF also binds VEGFR1 on other cells, e. g. monocytes and macrophages, that are implicated in arteriogenesis [26,83]. We found that removing PlGF or increasing PlGF secretion has only a modest effect on predicted VEGFR2 phosphorylation, while substantially altering VEGFR1 activation (Fig 6A). This result suggests that observed physiological PlGF-dependent pro-angiogenic effects are likely mediated directly by VEGFR1, either on ECs or other cells, and not via changes in VEGFR2 signaling, contrary to the ‘ligand-shifting hypothesis’. This result implicates VEGFR1 in the impaired angiogenic responses to ischemia, wound healing, and cancer [21] observed in mice lacking PlGF. It also implicates VEGFR1 in diseases where PlGF levels are known to change or to be predictive of prognosis, e. g. pre-eclampsia [42] and breast cancer [84]. The pro-angiogenic effects of PlGF likely also rely on its ability to up-regulate other growth factors, including VEGF, FGF2, and PDGF [85,86]. This result is not inconsistent with recent work by the Alitalo group showing that therapeutic over-expression of VEGFB (which like PlGF binds only VEGFR1) in mice improves metabolic health even following endothelial Flt1 gene deletion, and inhibits doxorubicin-induced cardiotoxicity [54,87]. Competition between ligands is concentration-dependent, and in these studies, VEGFB protein levels were elevated 20-fold or more in serum, heart, liver, and white adipose tissue. Our model predicts that competition is not a driver of PlGF signaling in physiological conditions, but does not preclude the existence of competition following supraphysiologic therapy. Indeed, at 10-fold PlGF over-expression, outside of the concentration range likely to be observed in untreated healthy or diseased tissue [42], the model does begin to predict an effect on VEGFR2 signaling. Both the ECM and sR1 regulate tissue levels of free interstitial VEGF and PlGF, the amount of growth factor available to bind ECs, and the steady-state distribution of ligand throughout the body (Fig 3). The model predicts that sR1 modulates the magnitude of EC receptor ligation, potentially also altering the balance of signaling via VEGFR1 vs. VEGFR2 (Fig 6G). This is of therapeutic interest because ratios of VEGF or PlGF to sR1 levels in plasma are increasingly of interest as a biomarker (e. g. in pre-eclampsia) [70], and sR1 levels increase in diabetic mice following hindlimb ischemia [88]. Including binding of immobilized ligands to sR1, and binding of immobilized sR1 to VEGF121 and PlGF1, increases total extracellular VEGF and PlGF stored in tissue (Fig 7). While there is not yet evidence to prove the existence of such complexes, the heparin- and ligand-binding sites on sR1 are distinct, as are the heparin- and receptor-binding domains on VEGF and PlGF, and therefore these complexes are likely. Unlike matrix-ligand-sR1 complexes, VEGF immobilized to both surfaces and ECM proteins has been shown to bind and activate VEGFR2 in vitro, preferentially increasing VEGFR2 activation of tyrosine Y1214, upstream of p38 phosphorylation and migratory cell behavior, demonstrating an important role for physical immobilization of VEGF in signal regulation in vitro [18,19,89]. However, whether VEGFR2 ligation by immobilized VEGF would occur to any notable extent in vivo, and what the physiological impact on EC receptor signaling would be, have been unknown. Here, we saw that including these reactions increased EC receptor ligation and altered VEGFR2 signaling (Fig 7). While the number of available sites in the EBM is not well-established, our model suggests that these M-L-R complexes may make up a small but significant portion of ligated EC receptors (Fig 4D). To improve our estimates of the extent of EC receptor ligand by EBM-bound growth factor, it is necessary to obtain better estimates of heparin-binding sites in basement membranes. Interestingly, the fraction of ligated VEGFR1 bound to immobilized ligand was predicted to be higher than that for VEGFR2, owing largely to the strong M-PlGF2 affinity (Fig 6F). To date, the impact of VEGFR1 ligation by immobilized ligand has not been studied. However, as these are largely PlGF2-VEGFR1 complexes (Fig 6F), EBM binding site density may shift relative ligation of VEGFR1 by VEGF versus PlGF, which is known to alter VEGFR1-mediated signaling [44]. Spatial patterning of receptor ligation by soluble and immobilized ligand is also likely to be important, but cannot be examined with this model. Additionally, the potential roles for HSPGs and NRP1 expressed on other cells engaging with VEGFR2 in trans [90,91] are of interest for future study. We were interested in differences in signaling between VEGF isoforms upon binding to VEGFR1 and VEGFR2. Explicitly simulating VEGFR2 trafficking and site-specific phosphorylation, placed in the context of physiological geometry and transport processes, allowed us to predict isoform-specific VEGFR2 signaling in vivo (Fig 5). Immobilization in the matrix alters VEGF distribution and the resulting signaling, while NRP1 alters VEGF-receptor binding and trafficking. By including these isoform-specific properties, the model predicts that VEGF121 induces a shift in VEGFR2 distribution towards early signaling endosomes, decreasing the signaling ratio pY1214/pY1175, and shifting the net cellular signaling towards proliferation. Conversely, a larger portion of VEGFR2 bound to VEGF189 was localized on the EC surface at steady-state, increasing pY1214/pY1175, and shifting the balance towards pro-migratory signaling (Fig 5C). This isoform-specific patterning in VEGFR2 signaling was seen in both the baseline case (Fig 5C), with all three VEGF isoforms present, and in cases where only single isoforms of VEGF were expressed (Fig 8C). This is key validation, as our simulated signaling predictions in humans match the observed vascular phenotypes in mice or tumors expressing single VEGF isoforms (Fig 9B). Interestingly, in the single isoform cases, change in relative activation of VEGFR1 and VEGFR2 were also predicted (Fig 8B), which may contribute to these phenotypes [92,93]. This is in line with another interesting model prediction; while all VEGF isoforms can bind to both VEGFR1 and VEGFR2, physiologically it appears that VEGF165 and VEGF189 bind almost exclusively to VEGFR2, while VEGF121 comprises a large portion of the ligand on VEGFR1, and also binds VEGFR2 to an extent (Fig 4D). This segregation of ligands suggests that, while ligand levels are limiting for receptor binding, VEGFR1 and VEGFR2 don’t directly compete for VEGF in vivo, instead binding to largely distinct subsets of ligands dictated primarily by isoform-specific NRP1-binding properties (Fig 9C). The relative levels of VEGF isoforms are not yet extensively-characterized, but they are known to vary by tissue and to change in disease [69,82,94,95]. As such, this model can be used to understand splicing-induced tissue- and disease-specific changes in VEGF receptor signaling. Our model is built upon experimental data and a validated model of VEGFR2 signaling in vitro, and provides new insight into distribution of and signaling by VEGF and PlGF isoforms in vivo. However, when interpreting the results, it is important to acknowledge mismatch between model predictions and experimental measurements, which may result from limitations of our modeling approach, uncertainly in interpretation of experimental measures, and/or missing understanding of underlying biological mechanism. Similar to previous models, our predicted interstitial VEGF concentrations when fitting the model to measured plasma VEGF levels are higher than those measured in tissues using microdialysis. This discrepancy could be due to: difficulty in obtaining accurate measurements for high molecular weight proteins using microdialysis; production of VEGF by blood sources (e. g. PBMCs, platelets) or specific organs (e. g. highly fenestrated tissue), reducing the requisite VEGF production by skeletal muscle; or degradation of VEGF by tissue-resident proteases and/or other cell types expressing VEGF receptors (modeled in [96,97]). Inclusion of proteases in the model would reduce immobilized growth factor stores at steady state. Additionally, as in previous models, the predicted fraction of plasma sR1 bound to ligand was higher than the experimentally-measured fraction. There are other soluble receptors that may be important to consider and are not included here. There may also be limitations with the experimental method that make these in vivo measurements inaccurate. To quantify the importance of some difficult-to-measure parameters, as well as reactions included in this model for the first time (some of which have not been explicitly demonstrated experimentally), we analyzed the sensitivity of many new or poorly characterized parameters (see S3 Fig and S1 File). In order to achieve simulation at the whole body scale, compartment models neglect spatial effects, instead predicting only average values for tissue. The interstitial space of the tissue, the cell surface of endothelial cells and the cell surface of myocytes are still independent entities in this case and each is treated as well-mixed. Detailed study of gradients in interstitial space and along cell surfaces, which are difficult to measure in vivo but are likely key to angiogenic signaling, requires development of detailed 2- and 3-dimensional models of tissue and experimental set-ups, calibrated to match predicted average concentrations from compartment models [98–101] such as the one presented here. Much work remains to fully understand the role of spatial gradients of VEGF distribution and receptor activation in health, disease, and response to therapy. This model integrates detailed regulation of VEGF and PlGF distribution and binding to EC VEGFR1 and VEGFR2 by sR1, the ECM, and NRP1 into a multi-scale pharmacokinetic/pharmacodynamic (PK/PD) framework. The resulting model predicts that all of these features interact, and contribute to regulation of tissue-level VEGF family signaling. While many model predictions are difficult to validate in vivo, the mechanisms included were first modeled using detailed in vitro measurements, and validated in many cases on the cellular level, before being put in a physiological context using an existing PK/PD framework. By progressively adding complexity, we can study the impact of each contribution, and compare simulation results to quantities that are measurable and to observable phenotypes, such as the vascular morphologies in mice expressing single isoforms of VEGF. By the same turn, this model provides a window into details of growth factor distribution and signaling that are essentially impossible to measure (especially on the protein level), though in many cases implicated in disease-related impairment in angiogenic response, or targeted by potential therapies. The lack of approved pro-angiogenic therapies to date makes it clear that a better understanding of the molecular mechanisms driving disease is critical to identify more effective drug targets, optimize drug properties (e. g. affinity), and avoid off-target effects leading to toxicity and drug failure [55]. This work can be extended to disease applications with changes in VEGF splicing, and to compare results in humans versus mice, to aid in translation of therapeutics targeting the VEGF system and to further validate the model against data obtained in mice.
Angiogenesis, the growth of new blood vessels from the existing vasculature, is critical for maintenance of health and response to injury. In ischemic disease, this process is impaired, but therapies targeting a key family of proteins, the vascular endothelial growth factors (VEGF), have failed to translate clinically. This suggests a need for deeper understanding of the complex regulation underlying angiogenic signaling. Here, we translate a previously developed and validated model of VEGF family signaling into a human, whole-body framework. The different splice isoforms of VEGF and the related PlGF proteins have different affinities for the extracellular matrix (ECM) and the co-receptor Neuropilin-1. Using our model, we examine the effect of these different binding properties on the distribution of each isoform in tissue, and subsequent receptor signaling. The model predicts isoform-specific receptor activation that is consistent with observed vascular phenotypes in mice expressing a single VEGF isoform; non-ECM-binding isoforms lead to signaling that promotes cell proliferation, while strong ECM-binding promotes migratory signaling and increased vessel branching. This understanding is critical for design of biomaterials that manipulate VEGF-ECM binding to control growth factor delivery, and for understanding of splicing-induced changes in VEGF family signaling in different tissues and in disease.
cell binding
phosphorylation
cell physiology
medicine and health sciences
placental growth factor
body fluids
vegf signaling
endothelial cells
epithelial cells
endocrine physiology
growth factors
cellular structures and organelles
animal cells
proteins
extracellular matrix
endocrinology
biological tissue
blood plasma
biochemistry
signal transduction
blood
cell biology
post-translational modification
anatomy
physiology
epithelium
biology and life sciences
cellular types
cell signaling
A computational analysis of in vivo VEGFR activation by multiple co-expressed ligands