VEGF and soluble VEGF receptor-1 (sFlt-1) distributions in peripheral arterial disease: an in silico model

Florence T. H. Wu, Marianne O. Stefanini, Feilim Mac Gabhann, Christopher D. Kontos, Brian H. Annex, Aleksander S. Popel


Vascular endothelial growth factor (VEGF) is a key regulator of angiogenesis, the growth of new capillaries from existing microvasculature. In peripheral arterial disease (PAD), lower extremity muscle ischemia develops downstream of atherosclerotic obstruction. A working hypothesis proposed that the maladaptive overexpression of soluble VEGF receptor 1 (sVEGFR1) in ischemic muscle tissues, and its subsequent antagonism of VEGF bioactivity, may contribute to the deficient angiogenic response in PAD, as well as the limited success of therapeutic angiogenesis strategies where exogenous VEGF genes/proteins are delivered. The objectives of this study were to develop a computational framework for simulating the systemic distributions of VEGF and sVEGFR1 (e.g., intramuscular vs. circulating, free vs. complexed) as observed in human PAD patients and to serve as a platform for the systematic optimization of diagnostic tools and therapeutic strategies. A three-compartment model was constructed, dividing the human body into the ischemic calf muscle, blood, and the rest of the body, connected through macromolecular biotransport processes. Detailed molecular interactions between VEGF, sVEGFR1, endothelial surface receptors (VEGFR1, VEGFR2, NRP1), and interstitial matrix sites were modeled. Our simulation results did not support a simultaneous decrease in plasma sVEGFR1 during PAD-associated elevations in plasma VEGF reported in literature. Furthermore, despite the overexpression in sVEGFR1, our PAD control demonstrated increased proangiogenic signaling complex formation, relative to our previous healthy control, due to sizeable upregulations in VEGFR2 and VEGF expression, thus leaving open the possibility that impaired angiogenesis in PAD may be rooted in signaling pathway disruptions downstream of ligand-receptor binding.

  • therapeutic angiogenesis
  • molecular systems biology
  • computational model
  • mathematical model

vegf, a family of cytokines that include VEGF-A, is central to the regulation of competing stimulatory and inhibitory biochemical reactions that keep angiogenesis, the growth of new capillaries from pre-existent microvasculature, at a critical homeostasis in normal physiology (11, 43, 83). There exist many endogenous splice variants of human VEGF-A; among the predominant proangiogenic isoforms are VEGF121 and VEGF165 (28, 43, 54). Once secreted into the extracellular matrix, VEGF121 is generally considered a freely diffusible isoform, whereas VEGF165 can be reversibly sequestered through its heparin-binding domain at interstitial and cell surface proteoglycans in significant quantities. On the endothelial cell surface, the specific receptor tyrosine kinases VEGF receptor 2 (VEGFR2; also known as Flk-1 in mice or KDR in humans) and VEGF receptor 1 (VEGFR1; also known as Flt-1) are thought to be the major transducers of stimulatory and modulatory angiogenic signals, respectively, upon binding and activation by VEGF ligands. Signaling VEGF-VEGFR complexes can additionally be stabilized by the nonsignaling transmembrane coreceptor neuropilin-1 (NRP1), which couples directly with VEGFR1 and indirectly with VEGFR2 through nonoverlapping binding sites on VEGF165 (42).

Alternative splicing of the VEGFR1 gene also produces soluble forms of VEGFR1 (sVEGFR1 or sFlt-1) that are truncated just before the transmembrane region and thus lack the intracellular signaling tyrosine kinase domains (34, 66). The retained extracellular domains confer ligand-binding affinity and receptor dimerization capacity to sVEGFR1; hence sVEGFR1 is thought to be a potent antagonist of VEGF bioactivity through two mechanisms: direct reversible sequestration of VEGF ligands or dominant-negative heterodimerization with surface VEGFR monomers (3, 12, 33, 37). A detailed review of the physiological and pathological functions of sVEGFR1 was recently published (79).

Therapeutic Angiogenesis for PAD

Peripheral arterial disease (PAD) is a chronic complication of systemic atherosclerosis that typically results in muscle ischemia downstream to occluded lower extremity arteries (30). PAD affects 20% of the population over 55 years of age (52) and is strongly associated with diabetes mellitus (DM); 20–30% of diabetic patients have PAD, and 20–30% of PAD patients have diabetes (48, 71). Patients of PAD have either intermittent claudication (leg pain or cramping with walking that is relieved with rest) or the more severe form of critical limb ischemia (where patients have leg pain at rest, ischemic ulceration, and gangrene) (59). The inadequacy of traditional treatments targeting the primary atherosclerotic obstruction (e.g., pharmacotherapy, surgical revascularization, endovascular interventions) has motivated research on proangiogenic therapies aimed at directly improving collateral perfusion at downstream ischemic tissues (41, 59). However, the numerous human clinical trials since 1996 for therapeutic angiogenesis in PAD, via intra-arterial, intravenous, or intramuscular delivery of the proteins or genes of angiogenic growth factors [VEGF121, VEGF165, and fibroblast growth factor-2 (FGF-2)] have not been able to reproduce the consistent efficacy observed in animal studies (15, 83). One possible reason for the clinical disappointments is that PAD patients may suffer from an inability to respond to changes in ligand concentration due to impaired receptor signaling and/or altered expression of antagonists, rather than the lack of angiogenic growth factors themselves (15). Another reason supposed that the pharmacokinetics of administrations were not optimal for localizing proangiogenic responses within ischemic tissue and that systemic elevation of angiogenic growth factors in blood could have a countereffect in further promoting angiogenesis at the vasa vasorum feeding the growth of primary atherosclerotic plaques (7, 15, 53). Both cases describe conditions of chronic pervasive vascular disease that are difficult to mimic accurately in current animal models of PAD, such as surgically induced acute hindlimb ischemia in diabetic mice (15, 29). Hence there is tremendous potential for computational models, featuring detailed ligand-receptor binding interaction networks, whole body multiscale compartmentalization, and pharmacokinetics, in overcoming the limitations of animal models, as a platform for preclinical research into therapeutic angiogenesis for human PAD.

VEGF and sVEGFR1 Distributions in PAD

It is still not fully understood how dysfunctionalities in the VEGF ligand and receptor system contribute to the impaired angiogenic response in PAD. There has been general consensus that plasma concentration of VEGF is elevated in PAD patients relative to healthy subjects (6, 7, 45); a recent study further correlated the magnitude of increase in plasma VEGF165 with PAD disease severity (25). On the other hand, reports of PAD-associated changes in plasma levels of sVEGFR1 were less consistent; older studies have documented significant decreases (6, 7), whereas newer studies did not observe significant changes (25, 45). Hence it remains uncertain whether circulating sVEGFR1 levels have prognostic value for angiogenic status in PAD.

Even less is known about the intramuscular state of VEGF system regulation in PAD. Currently, there are no noninvasive experimental techniques for readily measuring the in vivo interstitial concentrations of VEGF and sVEGFR1, or the endothelial surface density of signaling VEGF-VEGFR complexes, within the ischemic muscles of human PAD patients. Studies on mouse models of PAD have suggested that the maladaptive overexpression of sVEGFR1 and VEGFR1 in ischemic muscles may account for the blunted angiogenic response in PAD (29).

Computational Modeling: Objectives and Achievements

We have previously modeled VEGF ligand interactions with endothelial receptors in the context of PAD in single-tissue (rat extensor digitorum longis skeletal muscle) computational models (36). In recent model extensions, we developed whole body multicompartmental models for predicting the systemic distributions of VEGF and sVEGFR1 in healthy human subjects between three major compartments: the calf muscle, blood, and the rest of the body (70, 80, 81). Model features included detailed biomolecular interactions between VEGF121, VEGF165, VEGFR1, VEGFR2, NRP1, sVEGFR1, and interstitial matrix binding sites, as well as intercompartmental macromolecular transport through vascular permeability, lymphatic drainage, and direct plasma clearance. For completeness, the model is presented in the supplemental material (supplemental material can be found with the online version of this article).

The purpose of this study was to develop PAD counterparts to elements found in our previous healthy human model, through reparameterization of the calf compartment to reflect the pathological characteristics of PAD-affected calf muscles. We first investigated whether the concentrations of circulating VEGF and sVEGFR1 clinically observed in PAD patients could be reproduced in silico by simulating the pathological levels of intramuscular protein expression extrapolated from mouse models of PAD. At the PAD control point, we then examined whether overexpression of sVEGFR1 diminished formation of the VEGF-VEGFR signaling complexes in the calf or, alternately, whether the impaired angiogenic response in chronic human PAD is more likely rooted in pathological disruptions of signaling pathways downstream of ligand-receptor binding. Furthermore, we demonstrated the clinically translatable value of our computational model of VEGF and sVEGFR1 systems biology in the context of PAD via two examples: 1) the simulation of system responses to upregulated calf production of VEGF and sVEGFR1 to facilitate the optimization of proangiogenic therapies and 2) integrative molecular distribution analyses to assist in the diagnostic identification of underlying molecular disturbances to the VEGF system within PAD-affected ischemic tissues.


Compartmental Model Formulation

To model the systemic distributions of VEGF and sVEGFR1 associated with PAD patients, a physiologically based pharmacokinetic model was constructed that divided the human body into three compartments: 1) the calf compartment, where the classical symptom of PAD, intermittent claudication, is most commonly localized (55); 2) the normal (rest of the body) compartment, consisting of all other solid tissues; and 3) the blood compartment, which connects the other two compartments. Further spatial distinctions within the tissue compartments into parenchymal space, interstitial space, and capillary space were necessary to distinguish tissue-specific pools of VEGF and sVEGFR1 into their free versus matrix-bound versus receptor-bound subpopulations. As illustrated in Fig. 1A, histological data from human skeletal muscle cross sections were used to characterize the relevant volumes and surface areas of these spatial subdivisions for the solid tissue compartments. The parameters characterizing healthy tissue are included in the supplemental material Tables S. 1–S. 5; the diseased calf compartment was parameterized based on gastrocnemius histology of PAD patients as described below in Tables 19.

Fig. 1.

Schematic of multitissue model of VEGF and soluble VEGF receptor-1 (sVEGFR1) distributions. A: adapting the whole body compartmental model for a peripheral arterial disease (PAD) patient. Calf muscles are common symptomatic sites of impaired ischemic muscle angiogenesis in PAD. With the assumption that systemic manifestations of PAD were negligible, the normal (rest of the body) compartment of the multitissue model was kept at its previous healthy control parameterization. In contrast, the calf compartment was reparameterized to encapsulate pathological adaptations to PAD, including histological changes to ischemic skeletal muscle tissues, as illustrated here. BM, basement membrane; EBM, endothelial BM; MBM, myocyte BM; ECM, extracellular matrix; EC, endothelial cell; RBC, red blood cell. Illustrations adapted from muscle man series from Andreas Vesalius, De Humani Corpis Fabrica, 1543, courtesy of the National Library of Medicine; histological micrographs from Baum et al. J Vasc Res 44: 202–213, 2007. Note, BM thicknesses versus molecule sizes are not drawn to scale. B: mass flows through three compartment model. VEGF and sVEGFR1 were secreted from parenchymal and endothelial cell sources, respectively. Both were internalized upon binding with abluminal endothelial surface receptors. All soluble species were subjected to lymphatic drainage from interstitial space into the blood, bidirectional permeability through the endothelia, and direct clearance from the blood. C: interactions between VEGF121 (yellow), VEGF165 (blue), sVEGFR1 (orange), interstitial matrix binding sites (purple), and the abluminal endothelial cell surface receptors VEGFR1 (red), VEGFR2 (blue), and NRP1 (green). The VEGF-trapping interactions of sVEGFR1 could directly affect free VEGF concentrations (T1a or T2), NRP1 availability (T1b or T3), and interstitial matrix-bound reservoirs of VEGF (D1a,b). This implementation neglected sVEGFR1 heterodimerization with surface VEGFRs.

Within the two tissue compartments, VEGF was secreted from the myocytes in a 1:10 VEGF121:VEGF165 ratio [roughly corresponding with the mRNA expression ratio of freely diffusible isoforms (VEGF120) vs. heparin-binding isoforms (VEGF164 + VEGF188) in mice (54)], whereas sVEGFR1 was secreted abluminally from the endothelial cells into the interstitial space (Fig. 1B). Potential contributions to total body sVEGFR1 production from direct luminal endothelial secretion or proteolytic shedding of abluminal surface-bound VEGFR1 (10) were neglected in this implementation. Tissue concentrations of VEGF and sVEGFR1 were then distributed, through competing binding interactions summarized in Fig. 1C, over three populations: 1) diffusing in free or complexed (VEGF-sVEGFR1) form within the available interstitial fluid; 2) sequestered at heparan sulfate glycosaminoglycan binding-sites within the distinctly composed subregions of the interstitial space, i.e., endothelial basement membrane (EBM), extracellular matrix (ECM), and parenchymal basement membrane (PBM); or 3) bound to abluminal endothelial cell surface receptors. Since the current implementation neglected heterodimerization of sVEGFR1 with surface VEGFRs, internalization of sVEGFR1 from the endothelial surface could only occur through binding to surface NRP1s, whereas VEGF could be internalized through the six signaling VEGFR complexes in addition to uncoupled NRP1s (Fig. 1, B and C). Additionally, the soluble species, free VEGF, free sVEGFR1, and VEGF-sVEGFR1 complex, were subjected to intercompartmental transport processes: 1) lymphatic drainage from the interstitial space into plasma, 2) bidirectional vascular permeability between interstitial fluid and plasma, and 3) direct clearance from plasma (Fig. 1B).

In this spatially averaged model, gradients of soluble species were neglected within the interstitial fluid and plasma [an assumption justified in (23)], as was spatial variability in matrix- or receptor binding. As in our previous models (44, 70, 81), VEGF receptors surface receptors (VEGFR1 and VEGFR2) and the diffusible sVEGFR1 were present only in preassociated homodimeric forms. In particular, sVEGFR1 heterodimerization with surface VEGFR1 or VEGFR2 was neglected in this study to specifically investigate the VEGF-trapping effects of sVEGFR1 (Fig. 1C, interactions T1–3) on the VEGF system.

Mathematical Formulation and Notations

This human PAD model was based on the mathematical formulation used for our previous compartmental models (70, 80, 81; see supplemental material for complete mathematical equations); the solution of a system of 59 coupled nonlinear ordinary differential equations predicts changes in molecular concentrations over time as a function of binding and transport kinetic constants. In the equations, all blood and tissue concentrations of molecular species X are expressed in the notations of [X]j (mol/cm3 of tissue j; j = N for Normal or D for PAD calf) and [X]B (mol/cm3 of blood), respectively. In the discussion of results below, separate notations [X]IS,j (mol/l of available interstitial fluid in tissue j) and [X]pl (mol/l of plasma) explicitly refer to the concentrations of soluble species X (e.g., V for VEGF or sR1 for sVEGFR1) after unit conversion to their relevant fluid volumes, given in molar (M) units.

Model Parameterization for a Sedentary PAD Patient

This section describes the parameterization of PAD-specific compartmental geometry, kinetic rates, biotransport rates, and protein expression levels necessary to solve the model system of ordinary differential equations. All comparisons with the healthy model refer our previous parameterization for a healthy human subject (81), the values of which are summarized in the supplemental Tables S. 1–S. 5.



There had been increasing evidence that PAD is more than an impairment of blood flow or transport characteristics, since angioplasty at the atherosclerotic site alone did not improve walking performance as well as exercise therapy did, presumably via systemic metabolic adaptations (60). Tissue structural changes accompanying PAD, which may be associated with documented abnormalities in cellular metabolism (31), were modeled as follows for the calf compartment of a sedentary PAD patient.

An experimentally measured 40% increase in capillary-to-fiber ratio (CF) from healthy gave a CF of 1.63 for PAD, whereas the capillary density for PAD was taken directly as measured at 286 capillaries/mm2 tissue (50). Increased capillarization, relative to healthy, might indicate an adaptive angiogenic response to muscle ischemia, albeit insufficient to alleviate symptoms of PAD. However, a reduction in capillary lumen size has been observed in PAD, from 18% in mild to 43% in severe intermittent claudication (46). With the assumption of an average reduction of 36% from healthy (46), the lumen cross-sectional area (CSA) became 7.91 μm2 in PAD, whereas the endothelial portion of the capillary CSA remained unchanged from its healthy value at 11.6 μm2 (75, 81).

PAD-associated muscle atrophy has been reported at a 15% reduction in leg muscle mass (1), as well as individual muscle fiber CSA (FCSA) reductions averaged at 17% (1, 50, 63), relative to healthy; thus we assumed a compartment volume of 738 cm3 (representative of the combined volume of unilateral gastrocnemius and soleus muscles under PAD conditions) and FCSA of 3,464 μm2. In addition, the reduced muscle fiber density of 175 fibers/mm2 tissue versus healthy might be indicative of muscle cell death consistent with observed PAD-associated muscle denervation (8, 63).

All together, the following changes in volume fractions were reflected in our parameterization of the PAD-affected calf muscle histology compared with healthy: a decreased muscle fiber space of 0.608 cm3/cm3 tissue, similar capillary and vascular spaces of 0.006 and 0.002 cm3/cm3 tissue, respectively, along with a more than doubled interstitial space of 0.387 cm3/cm3 tissue, as summarized in Table 1 and Fig. 1A.

View this table:
Table 1.

Geometric parameterization for the PAD calf compartment

Within the interstitial space (IS), although PBM thickness was relatively unchanged, EBM thickness increased drastically depending on concurrent diabetic status (69), age (78), and PAD severity, up to 1,650 nm in late-stage PAD requiring amputation (5). The available fluid fractions of IS subvolumes (Table 1) were calculated using the same solid fractions and partition coefficients as before (70, 81).


By confining the symptomatic muscle ischemia and PAD-associated histological changes to the calf compartment of this PAD model (i.e., assuming negligible effects of PAD on skeletal muscle far away from the calf), the normal compartment was assigned the same set of geometric parameters previously established for the healthy human model (81) (Tables S. 1–S. 5). Similarly, the blood compartment was maintained at a volume of 5 l, 60% of which was plasma.

Binding and coupling kinetics.

All binding and coupling kinetic rates were based on the same in vitro values used in our previous healthy compartmental model (81) (Table S. 4); tissue values for the calf compartment, converted based on PAD calf geometry, are tabulated in Table 2. Full sensitivity analysis on the effects of binding affinities and kinetic constants on steady-state protein distributions has been performed and published previously in the context of health (81).

View this table:
Table 2.

Kinetic parameters of the VEGF system (control values)



The molecular size-dependent basal vascular permeability rates used in this PAD model were the same as those used in the healthy model (81) (Table S. 2). In light of evidence suggestive of decreased muscle perfusion in PAD patients, including measurements of decreased red blood cell velocity (17), increased time to peak intensity of bolus-injected ultrasound contrast agents (19), decreased calf muscle perfusion at peak exercise (35), decreased skin perfusion pressure (68), and increased percentage of no-flow capillaries (9), we assume the effective vascular permeability rate for the PAD calf at control (supine position) to be at 50% of healthy basal level by adjusting the endothelial surface area recruitment factor to γ = 0.5 (Table 3). However, the loss of the venoarterial reflex (VAR) had been documented in PAD patients, such that the lower-than-normal perfusion of the supine foot about doubled (i.e., increased, rather than decreased as should happen with a normal VAR) to a higher-than-normal perfusion in the dependent position (14). Hence we note the possibility that the effective vascular permeability rates could reach 100% of healthy levels (γ = 1) with PAD patients in dependent positions (sitting, standing, etc.).

View this table:
Table 3.

Transport parameters: surface area recruitment factor (γ) for effective vascular permeability rates for PAD calf


The nighttime (asleep; supine position) lymph flow rate for the calf compartment, adjusted for the mass of the PAD calf, was 0.0004 cm3/s (Table 4). With the maintenance of the relationship that awakeness alone could increase the average basal lymph flow rates to about fivefold of night values (57), 0.0022 cm3/s, the basal daytime (awake; supine posture) rate, was taken as the control lymphatic drainage rate for the PAD calf in simulations (Table 4). We did not explore higher muscle activity-dependent lymphatic flow rates in the PAD model, as we did in our study of healthy subjects (81), with the assumption that PAD patients with symptomatic intermittent claudication or critical limb ischemia may experience reduced mobility relative to healthy subjects. The higher exercise-dependent ranges of lymphatic drainage rates should be investigated in future comprehensive simulations of exercise therapy, along with other parameter changes in response to exercise therapy (e.g., altered ligand/receptor protein expression, higher effective permeability due to exercise recruitment of permeable surface area, histological adaptations to long-term exercise, etc.).

View this table:
Table 4.

Transport parameters: lymphatic flow rates


The plasma clearance rates for soluble species used in this PAD model were unchanged from the healthy model (81) (Table S. 2).

We have previously investigated the steady-state and dynamic system sensitivity to all transport parameters (vascular permeability, lymphatic drainage, plasma clearance) in the context of health (80, 81); similar behaviors are expected for the PAD model.

Initial concentrations/protein expression levels.


As in the healthy model (81), we estimated the densities of interstitial matrix binding sites for VEGF and sVEGFR1 based on only heparan sulfate proteoglycan (HSPG) contributions. For the ECM, the same estimated molar density of VEGF-binding sites from the healthy model was converted using the new ECM volume fraction of the PAD calf compartment (Table 5). For the BMs, we modeled the HSPG content of the PAD calf compartment based on measurements from Engelbreth-Holm-Swarm (EHS) sarcomas in diabetic (DM) mice, where the decreased HSPG content and synthesis were postulated as trigger events for compensatory increase in collagen synthesis and BM thickening in DM (64). Justification for using the EHS and DM models were as follows: 1) EHS mouse tumors produce a homogenous interstitial matrix of basement membrane material, making EHS interstitial matrix a typical in vitro model for physiological BMs, which are usually too thin for content analysis (38, 64) and 2) a significant proportion of PAD patients are diagnosed with concomitant DM, and the pathological changes to BM geometry are similar in PAD and DM, as described in previous sections. Using the conversion from HSPG content into VEGF-binding site densities (81), we derived a 5-μM density of VEGF-binding sites for the BMs of the PAD calf (Table 5), at four times lower than that of the healthy calf.

View this table:
Table 5.

Initial conditions for interstitial VEGF/sVEGFR1-binding site densities, based on matrix HSPGs


The previous healthy receptor densities (70, 81), 10,000 VEGFR1/EC, 10,000 VEGFR2/EC, and 10,000 NRP1/EC, were maintained for the normal compartment in the present PAD model. To the best of our knowledge, pathological changes to VEGFR protein expression on endothelial cell surfaces have not been extensively quantified in human PAD subjects. Thus in the PAD calf compartment, we approximated pathological receptor densities by extrapolating from animal models of PAD, specifically, based on protein expression changes observed in the tibialis anterior muscle between normal (NC for normal chow or WT for wild-type) versus postischemic diabetic (DM10dI for 10-day postischemia diabetic or db14dI for 14-day postischemia diabetic) mice (29, 65), as summarized in Table 6.

View this table:
Table 6.

Initial conditions for surface receptor densities


Similarly, pathological changes to interstitial free VEGF and sVEGFR1 concentrations have not been quantified for human PAD subjects. We thus estimated the target ranges for pathological interstitial concentrations by extrapolating from animal models of PAD, specifically based on protein changes observed in tibialis anterior muscle homogenates between normal (NC) mice versus 3- or 10-day postligation diabetic (DM3/10dI) mice (29), as summarized in Table 7. The limitations of these estimations were noted: 1) homogenate measurements could not discriminate between interstitial and intracellular proteins, 2) ligated diabetic mice were imperfect models for human PAD subjects, and 3) chronic concentration changes might not have been evident at 10 days postligation. Thus, in the case that plasma and interstitial concentrations could not be simultaneously modeled, we opted to fit targeted ranges for plasma concentrations rather than the interstitial targets.

View this table:
Table 7.

Steady-state interstitial concentrations of VEGF and sVEGFR1


In PAD patients, given that the majority of literature data supported several-fold increases in plasma VEGF concentration and several-fold decreases in plasma sVEGFR1 concentration compared with healthy subjects (6, 7), we targeted a range of 3 to 6 pM of VEGF and ∼30 pM of sVEGFR1 in the blood for PAD simulations (Table 8), i.e., at ∼3 times higher and lower than our previous healthy targets, respectively.

View this table:
Table 8.

Steady-state plasma concentrations of VEGF and sVEGFR1

Numerical Solutions and Computational Simulations

All simulations and data plots presented in this article were performed on the numerical computing platform, Matlab 7.3.0 R2006b (MathWorks, Natick, MA) and run on a desktop PC. The system of coupled nonlinear ordinary differential equations was solved as an initial value problem to steady state using ode15s solver routine of Matlab, with a relative error tolerance set at 10−6 (0.0001% accuracy) and an initial step size of 10−4.

The first part of our simulation experiments involved a search for the secretion rates of VEGF and sVEGFR1 from the calf compartment that could reproduce a PAD control patient's expected ranges of interstitial and plasma concentrations of VEGF and sVEGFR1 (Tables 7 and 8), with all other parameters set at values characterizing a supine PAD patient (Tables 16). The basal molecular distribution and flow profiles at the established PAD control were then compared with those of our healthy control results (summarized in supplemental material Tables S. 7 and S. 8); their expected influence on simulated system responses to variation in secretion rates were noted. Finally, an etiological study was done to explore how individual parameter changes (e.g., in geometry vs. secretion rates) from healthy to PAD state contributed to the overall pathological changes in VEGF and sVEGFR1 distributions. From this physiologically based sensitivity analysis, the predicted correlation patterns between the responses of plasma free VEGF, plasma free sVEGFR1, and VEGF signaling potentials in the calf allowed the design of better differential diagnostic tools to pinpoint the underlying etiological defect within the VEGF and sVEGFR1 system (e.g., pathological receptor expression vs. pathological secretion rates) in PAD patients.


Targeting PAD-associated VEGF and sVEGFR1 Concentrations as Functions of Their Secretion Rates from the PAD Calf Compartment

To simulate a PAD patient, we kept the normal compartment parameters identical to the healthy subject's parameter set but changed the relevant calf parameters, specifically, the geometric volume distributions, permeability rates, basement membrane matrix site densities, and receptor densities, to their PAD-associated values summarized in Tables 1, 3, 5, and 6. Although the VEGF and sVEGFR1 secretion rates from the normal tissue were kept at the control values (Table S. 6) from our previous model of a healthy subject, we attempted to target the expected plasma and calf interstitium concentrations of VEGF and sVEGFR1 for PAD patients (Tables 7 and 8) by varying VEGF and sVEGFR1 secretion rates from the calf compartment alone (qTotalVEGF,Calf and qsR1,Calf), as illustrated in Fig. 2.

Fig. 2.

Targeting calf secretion rates of VEGF (qTotalVEGF,Calf) and sVEGFR1 (qsR1,Calf) to match the basal profile of PAD patient. AF: concentrations of free VEGF (AC) and free sVEGFR1 (DF) in the normal interstitial fluid (A and D), plasma (B and E), and PAD calf interstital fluid (C and F). Orange spheres: although calf secretion rates were kept at their healthy control values, interstitial concentrations of VEGF and sVEGFR1 were found to be much reduced and elevated with respect to their healthy control levels, due to PAD-associated changes to the geometric, protein expression, and transport parameters of the calf compartment. Yellow spheres: qTotalVEGF,Calf and qsR1,Calf were increased to reach the upper bounds of the PAD target ranges for the calf interstitial VEGF and sVEGFR1 concentrations. Purple spheres: PAD-target concentration of 4.5 pM VEGF in plasma was reached by further synergistic increases to both qTotalVEGF,Calf and qsR1,Calf. This point was henceforth defined as the PAD control point; we noted that the simulated deviations from targeted PAD ranges in plasma sVEGFR1, calf interstitial VEGF, and normal interstitial sVEGFR1 concentrations could not be reconciled. Green spheres: even imposing a compensatory decrease in sVEGFR1 secretion rate from the normal compartment (qsR1,Normal = 0) was not able to sufficiently dampen the fluctuation in the interstitial sVEGFR1 concentration of the normal compartment caused by pathological qsR1,Calf. Orange/black arrows indicate inverse/direct relations between sVEGFR1 concentrations and VEGF secretion rate. Blue/black arrows indicate inverse/direct relations between VEGF concentrations and sVEGFR1 secretion rate. Dashed arrows indicate relations that differed from the healthy model. +Control; max, min, and mean, PAD target ranges; EC, endothelial cell; V121, VEGF121; V165, VEGF165; R1, VEGFR1; sR1, sVEGFR1; R2, VEGFR2.

Starting off with qTotalVEGF,Calf and qsR1,Calf at their healthy control values (Fig. 2, orange spheres), the plasma concentrations ([V]pl = 1.5 pM; [sR1]pl = 100 pM) and the interstitial concentrations in the normal compartment ([V]IS,Normal = 10 pM; [sR1]IS, Normal = 36 pM) remained consistent with those predicted in our previous healthy control model (81) (Table S. 7), whereas the calf interstitial concentrations deviated considerably (i.e., [V]IS,Calf drastically reduced to 0.8 pM; [sR1]IS,Calf elevated to 43 pM) as a result of the PAD-associated parameterization for the calf compartment. These deviations in [V]IS,Calf and [sR1]IS,Calf (Fig. 2, C + F) were respectively in opposite direction and in lesser magnitude than expected from the experimental changes observed in VEGF and sVEGFR1 expressions in postischemic diabetic mice models of PAD (Table 7) (29). The lack of change in [V]pl and [sR1]pl also did not correspond to the clinically observed changes in PAD patients, i.e., increased [V]pl and decreased [sR1]pl (Table 8) (6, 7, 39, 45).

Thus our first step was to increase protein secretion rates in the calf (both qTotalVEGF,Calf and qsR1,Calf) in attempt to reach the upper bound of the targeted PAD ranges of interstitial VEGF and sVEGFR1 concentrations in the calf (Table 7); [V]IS,Calf = 63 pM and [sR1]IS,Calf = 505 pM were successfully reached at secretion rates of qTotalVEGF,Calf = 3.74 molecules·MD−1·s−1 and qsR1,Calf = 0.705 molecules·EC−1·s−1 (Fig. 2, C + F, yellow spheres). However, the concurrent changes to their plasma concentrations did not match those clinically observed in PAD (Table 8), i.e., [V]pl only elevated slightly from healthy control to 2.4 pM, falling short of the PAD-targeted lower bound of 3 pM, whereas [sR1]pl increased to 179 pM rather than decreased (Fig. 2, B + E).

To further elevate plasma VEGF ([V]pl), without relying solely on increasing calf production of VEGF (qTotalVEGF,Calf), which would soon exponentially increase [V]IS,Calf to beyond physiological bounds, production of sVEGFR1 (qsR1,Calf) was simultaneously increased to contain [V]pl at the PAD target of 4.5 pM. This point, henceforth defined as our PAD control point, was achieved at secretion rates of qTotalVEGF,Calf = 6.6 molecules·MD−1·s−1 and qsR1,Calf = 2 molecules·EC−1·s−1 (Table 9; Fig. 2, purple spheres). At this PAD control point, [V]IS,Normal, [V]pl, and [sR1]IS,Calf were successfully targeted (Fig. 2, A, B, and F), although [V]IS,Calf, [sR1]IS,Normal, and [sR1]pl were not (Fig. 2, CE), for reasons as follows. First, [sR1]pl was 317 pM at the PAD control point, about 3× higher than that at our previous healthy control (Table S. 7), because [sR1]pl exhibited a strong directly proportional relationship with qsR1,Calf in the range we considered, such that an upregulated calf expression of sVEGFR1, which was necessary to assist in the elevation of [V]pl, inevitably led to a drastically increased plasma concentration of sVEGFR1 (Fig. 2E). Second, although [sR1]IS,Calf (at 506 pM) was able to stay within targeted PAD range, we noted that a higher than expected elevation in [V]IS,Calf (to 242 pM; Fig. 2C) was necessary to achieve the PAD target [V]pl of 4.5 pM (Fig. 2B), which we accepted as part of the uncertainty inherent with modeling human PAD patients using interstitial data extrapolated from mice models. Additionally, [V]IS,Normal (at 10 pM) remained relatively insulated from the changes in the calf secretion rates (Fig. 2A), whereas [sR1]IS,Normal responded strongly to increased qsR1,Calf with its own elevation (to 81 pM; Fig. 2D). In other words, localized pathological changes to the secretion rate of sVEGFR1 (qsR1,Calf) affected not only its local tissue concentration ([sR1]IS,Calf) but also caused its global tissue concentration ([sR1]IS,Normal) to deviate significantly from its healthy control level (36 pM, Table S. 7 and Fig. 2D). Even when a compensatory shut-down in global production rate of sVEGFR1 (qsR1,Normal = 0 molecules·EC−1·s−1) was introduced (Fig. 2, green spheres), both [sR1]pl (at 226 pM, Fig. 2D) and [sR1]IS,Normal (at 45 pM, Fig. 2E) were still elevated relative to healthy control (Table S. 7).

View this table:
Table 9.

VEGF- and sVEGFR1 secretion rates (fitted kinetic parameters)

Altered Basal Molecular Distribution Profiles for a PAD Patient

The molecular distributions established at the PAD control point were compared with those of the healthy control; see PAD control versus healthy control in supplemental material Fig. S. 3 for details.

VEGF distribution.

Although free VEGF concentrations in plasma and calf interstitial fluid were higher in the PAD control, their VEGF165 isoform fractions were actually lower than that in the healthy control (Fig. S. 3A). Among the total VEGF in plasma, the sVEGFR1-bound fraction increased from 77% in healthy to 91% in PAD (Fig. S. 3B), much higher than that experimentally measured by Belgore et al. (6) in the plasma of PAD patients (∼0.6% mol fraction, based on 403 pg/ml free VEGF and 8 pg/ml complexed sVEGFR1). The total amount of VEGF in muscle did not change significantly for the normal compartment but was drastically higher in the PAD calf than the healthy calf as a result of higher availability of reservoirs (sVEGFR1, matrix, surface receptors; Fig. S. 3, C and D). On top of the increase in absolute concentrations of total VEGF, its distribution between sVEGFR1, matrix sites, and surface receptors was also changed from healthy to PAD control. In particular, total VEGF121 in the PAD calf had a notably increased sVEGFR1-bound fraction (51% vs. 5%) and lowered surface receptor-bound fraction (46% vs. 90%), among the latter, VEGF121·VEGFR2 had taken over as the predominant surface receptor-bound species instead of VEGF121·VEGFR1·NRP1. On the other hand, total VEGF165 in the PAD calf had a greatly increased sVEGFR1-bound fraction (23% vs. 1.5%), a lower matrix-bound fraction (52% vs. 72%), and a similar surface receptor-bound population (23% vs. 25%). Among the latter, VEGF165·VEGFR2 had taken over as the predominant surface receptor-bound species instead of VEGFR2·VEGF165·NRP1. Overall, in the context of PAD, sVEGFR1 became a much more significant volumetric reservoir of VEGF, whereas the role in VEGF·VEGFR complex formation of NRP1 was reduced in favor of unassisted direct binding to VEGFR2.

sVEGFR1 distribution.

Total sVEGFR1 only tripled in plasma and in normal interstitium from healthy to PAD control but increased ∼28-fold in the calf interstitium (Fig. S. 3B). VEGF occupancies of total sVEGFR1 in PAD were also considerably higher at 13% versus 5% in plasma and 17% versus 0.6% in calf muscle but relatively unchanged in normal muscle (<0.5%). Again, the predicted VEGF occupancy of total plasma sVEGFR1 was higher than that experimentally measured in the plasma of PAD patients (<2% mol fraction, based on 8 ng/ml free sVEGFR1 vs. 18 pg/ml complexed sVEGFR1) (6). The increased prominence of sVEGFR1-VEGF complexes at the PAD control was reflective of the fact that when PAD-associated upregulation of abluminal receptors strongly sequestered VEGF interstitially, sVEGFR1 became a necessary vehicle in transporting VEGF into the plasma to sustain any concurrent elevation in plasma VEGF.

Occupancy of matrix sites.

The fractional occupancies of total matrix sites remained <0.4% in the normal compartment's interstitium but increased about 15-fold to ∼3% in PAD calf interstitium (uniformly across the EBM, ECM, and PBM; Fig. S. 3, EG). In other words, at steady state, ∼97% of the matrix sites in the PAD calf at control were still not utilized in the sequestration of VEGF or sVEGFR1, despite being exposed to high concentrations of free VEGF (242 pM) and high free sVEGFR1 (506 pM). However, we could not rule out the possibility that total matrix site occupancy could be higher than 3%, had we allowed matrix sites to directly buffer sVEGFR1·VEGF complexes (currently at 3,610 pM without matrix sequestration) in addition to the complex constituents.

Occupancy of VEGFRs and NRP1.

Relative to the healthy calf, the PAD calf at control was predicted to have dramatically increased VEGF-ligated non-NRP1-coupled fractions of VEGFR1 (VEGF·VEGFR1: 83% vs. 6%) and VEGFR2 (VEGF·VEGFR2: 56% vs. 3%). This resulted in lower NRP1-coupled VEGF-bound fractions of VEGFR1 (VEGF121·VEGFR1·NRP1) and VEGFR2 (VEGFR2·VEGF165·NRP1) despite almost all (97%) NRP1s being coupled to signaling VEGFRs (Fig. S. 3, HJ). The overall signaling profile of the calf compartment at PAD control was still proangiogenic, but by a wider margin than that in the healthy control; specifically, there were ∼591 pmoles/cm3 tissue more VEGF-bound VEGFR2s than VEGF-bound VEGFR1s at PAD control versus 9 pmoles/cm3 tissue more at healthy control (Fig. S. 3K). As was the case in the healthy model, VEGF occupancies of VEGFRs were directly proportional to interstitial free VEGF concentrations; hence qTotalVEGF,Calf was greatly effective in altering the signaling potentials of the PAD calf, whereas qsR1,Calf was not (compare Fig. S. 2K vs. S. 1K).

Altered Basal Concentration Gradients and Flow Profiles of PAD model

At the PAD control point, the concentration gradients reversed directions in two areas (Fig. 3, dashed arrows) compared with the healthy control flow profile (Table S. 8): 1) free sVEGFR1 in the PAD calf interstitium became higher than that in plasma and 2) plasma sVEGFR1-VEGF complex in a PAD patient model became higher than that in the normal interstitium. Hence, aside from the altered magnitudes of macromolecular flows, in particular those leaving and entering the calf interstitium, the biggest changes in the basal flow profile for a PAD patient relative to a healthy subject were predicted to be the reversed directions of net permeability of free sVEGFR1 across the calf endothelium and of the sVEGFR1-VEGF complex across the normal endothelium.

Fig. 3.

PAD patient: basal steady-state flow profiles of free VEGF (left), sVEGFR1-VEGF complex (middle), free sVEGFR1 (right). Solid-colored arrows represent intracompartmental flows (i.e., secretion, internalization) and intercompartmental flows (i.e., net vascular permeability, lymph flow, plasma clearance), with their relative magnitudes reflected by arrow widths. Dashed arrows indicate flows with directions reversed compared with the basal profile of the healthy model. Graded-colored arrows between columns indicate mass transfer flows between species (i.e., net association of free VEGF and free sVEGFR1 to form sVEGFR1-VEGF complexes or net dissociation of the complex back into its constituents).

The basal flow profile influences how the VEGF/sVEGFR1 system responded to perturbations in their secretion rates and transport parameters. In comparing the PAD system's responses to increasing qTotalVEGF,Calf and qsR1,Calf (Fig. 2) with the responses of the healthy system, two differences were noted: 1) at high qTotalVEGF,Calf increasing qsR1,Calf slightly increased, rather than decreased, free VEGF in the normal interstitium (Fig. 2A, black dashed arrow) and 2) at high qsR1,Calf increasing qTotalVEGF,Calf beyond about 2 molecules·MD−1·s−1 started to slightly decrease, rather than increase, plasma free sVEGFR1 (Fig. 2E, orange dashed arrow). Both cases could potentially have resulted from the influence of the PAD-associated net extravasation of sVEGFR1·VEGF complex into the normal interstitium. In the former case, near the PAD control point, the tendency for interstitial free sVEGFR1 to associate with free VEGF might have lessened due to the influx of complexes from plasma, thus reducing the effect of local complex formation in lowering interstitial free VEGF as seen in the healthy model. In the latter case, complex extravasation into the normal compartment presented an alternative process to local complex dissociation for plasma to equilibrate the extra complexes intravasating from the calf when qTotalVEGF,Calf was increased; thus instead of responding with a gain in free VEGF from complex dissociation as in the healthy model, the PAD model instead responded with an overall loss of free VEGF to normal tissue.

Relative Contributions of Individual Parameter Changes to the Overall Pathological Molecular Distribution Profiles

Multiple pathological characteristics associated with PAD, e.g., changes in histological geometry, ligand or receptor protein expression levels, and basal vascular permeability rates, were collectively modeled at the PAD control point as presented above. To examine the individual contributions of these etiological factors in transforming the VEGF and sVEGFR1 distribution profiles from healthy control to PAD control state, a series of steady-state physiologically based sensitivity analyses were performed (Fig. 4 and supplemental Fig. S. 3). In this series of simulations labeled PAD + {Healthy XCalf}, all parameters were set at PAD control values except for one calf attribute X that was set at its healthy control value.

Fig. 4.

Contributions of individual parameter changes to PAD-associated free VEGF, free sVEGFR1, and signaling profiles: a diagnostic application. A, C, and E: free VEGF. B, D, and F: free sVEGFR1. GJ: signaling profiles. AJ: simulations between PAD control (ctrl) and healthy control had all parameters set at PAD control values except for calf attribute X set at its healthy control value. Sensitivity of a particular distribution to attribute X was quantified by how different/similar the steady state value of run PAD + {Healthy XCalf} was to the PAD control/healthy control run. For instance, the PAD-associated increase (from healthy) in plasma free VEGF was mostly attributed to pathological secretion rates of VEGF and sVEGFR1, while restrained by pathological changes to histological geometry and VEGFR upregulation. Differential sensitivities were observed between biological fluid compartments (e.g., plasma vs. calf interstitium) or between molecular species (e.g., VEGF vs. sVEGFR1 vs. signaling VEGFRs). K: a clinical application of results from AJ, using combinatorial analysis of multiple markers for etiological differential diagnosis at the molecular level. Kp, vascular permeability rate; [BM], density of matrix binding sites for VEGF or sVEGFR1; [R1], [R2], and [NRP1], surface receptor densities; qsR1 and qVEGF, secretion rates of sVEGFR1 and VEGF; Geom, geometric parameters. Percentages in brackets are the fractional occupancies of total VEGFR1 or total VEGFR2. See supplemental Fig. S. 3 for the complete set of molecular distribution plots.

Within each molecular distribution plot, the X attribute of the PAD + {Healthy XCalf} simulation that produced results most similar to the healthy control case or most different from the PAD control case was the most significant contributing factor to the overall effect of PAD on that distribution. For instance, pathological changes in VEGF and sVEGFR1 secretion rates were the predominant contributing factors to the PAD-associated increase in [V]pl (Fig. 4A). Additionally, some factors produced opposing effects, suggestive of compensatory feedback mechanisms, e.g., pathological VEGFR upregulation and histological geometry changes actually lowered [V]pl, returning it closer to the healthy control value (Fig. 4A). Such intradistribution analyses of the direction and extent to which individual pathological factors contributed to our PAD control profiles could provide insight into patient-specific diagnoses, e.g., a PAD patient's [V]pl being much lower than our PAD control's 4.5 pM could indicate a smaller upregulation in VEGF/sVEGFR1 secretion or a greater overexpression of VEGFR than that represented in our PAD control model.

Moreover, interdistribution analyses showed that sensitivities to individual pathological factors were not uniform across biological fluid compartments or between molecular species. Two examples of nonuniform concentrational changes for the same molecular species across fluid compartments include: 1) [V]IS,Calf was much less sensitive than [V]pl to PAD-associated sVEGFR1 secretion rates and was in fact affected in the opposite direction (Fig. 4, A and E) and 2) the effects of PAD-associated geometry, vascular permeability rates, and VEGF secretion rates on [sR1]IS,Calf were not similarly reflected in [sR1]pl (Fig. 4, B and F). Examples of nonuniform changes in the same fluid compartment of different molecular species include unlike [V]pl, which decreased with PAD-associated VEGFR expression levels, both [sR1]pl and the calf density of VEGF-VEGFR complexes (or signaling potential) increased as a result (Fig. 4, A, B, I, and J). This implied that the use of either free VEGF or free sVEGFR1 in the plasma alone as a diagnostic marker for PAD may not always accurately reflect the sensitivity of surface signaling complexes, i.e., the two circulating angiogenic markers may not fully encapsulate the interstitial angiogenic status in all scenarios.

However, the combinatorial use of plasma free VEGF and sVEGFR1 markers may still have integrative value for differential diagnosis. For instance, Fig. 4 showed that by exploiting the differential sensitivities of VEGF and sVEGFR1, the predominant pathological factor may be identified in individual patients: 1) a drop in [V]pl paired with a rise in [sR1]pl could uniquely indicate increased VEGFR expression, 2) a rise in [V]pl but not in [sR1]pl could suggest increased VEGF secretion, whereas 3) a rise in both [V]pl and [sR1]pl could signify increased sVEGFR1 secretion. Yet another alternative might be to monitor the free versus bound ratios of VEGF or sVEGFR1 as indication of changes in VEGFR expression: the higher the VEGFR overexpression, the lower the complexed fraction plasma sVEGFR1 or the lower the free fraction of plasma VEGF (supplemental Fig. S. 3).


Targeting PAD Control Point: Insights on Clinical Discrepancies and Model Limitations

In developing this novel computational model of the VEGF ligand and receptor distributions in human PAD patients, we began by reparameterizing the calf compartment of our previous healthy human model to reflect characteristics of critical limb ischemia and assuming that the normal compartment, as representative of the rest of the body, was unaffected by disease. More precisely, in the normal compartment, the geometric parameters and protein expression levels were kept at previous healthy control levels, whereas in the calf compartment, the surface receptor densities were first adjusted according to experimental data from mice models of PAD, then soluble protein secretion rates were tuned in an attempt to fit the PAD-associated target ranges of VEGF and sVEGFR1 concentrations in the calf interstitium (extrapolated from mice data) and plasma (human clinical data) simultaneously. However, several predicted trends prevented the simultaneous achievement of all targeted concentrations.

First, under no circumstances could simulated changes in protein secretion rates from the calf compartment predict concurrently increased VEGF and decreased sVEGFR1 levels in plasma, as clinically observed in PAD patients relative to healthy subjects (6, 7). As shown in Fig. 2, increasing secretion of sVEGFR1 in the calf (qsR1,Calf) generally elevated plasma concentrations of both VEGF and sVEGFR1 ([V]pl and [sR1]pl), whereas increasing VEGF secretion by the calf (qTotalVEGF,Calf) effectively raised [V]pl but had negligible effect on [sR1]pl. The present study did not rule out the possibility that concomitant alterations in protein secretion rates from the normal (rest of the body) compartment could contribute to the set of plasma concentrations associated with PAD; this was not explored for the lack of physiological evidence of nonlocalized (distal to ischemic sites) or systemic disturbances to the VEGF system as a result of PAD. Additionally, there had been more recent clinical studies reporting plasma concentrations inconsistent with the studies on which our target ranges were based; in fact, Fig. 2 would give more support to trends observed by Makin et al. (45) and Findley et al. (25), both reporting a much lesser increase in [V]pl and an unchanged [sR1]pl in PAD patients relative to healthy subjects. These disparities among clinical studies remain to be resolved to establish the true plasma target ranges for PAD simulations.

To sustain an elevated plasma VEGF level at the mean of the PAD target range (4.5 pM), the interstitial VEGF concentration in the calf compartment (242 pM) had to exceed its PAD target range; even with the assistance of increased secretion of sVEGFR1 in the calf (qsR1,Calf) in reaching a higher plasma VEGF concentration ([V]pl), an amplification of VEGF secretion by the calf (qTotalVEGF,Calf) was unavoidable, ultimately resulting in a calf interstitial sVEGFR1 concentration ([sR1]IS,Calf) at four times the upper bound of its targeted range. However, we noted the compromised stringency of our PAD calf interstitial target ranges due to their extrapolation from mice models of PAD with surgically induced acute hindlimb ischemia. Until experimental techniques are available to quantitate skeletal muscle interstitial concentrations in human PAD patients, it remains uncertain whether interstitial VEGF is as high as predicted for ischemic muscle in chronic PAD. Nonetheless, although 242 pM of interstitial VEGF would appear to be an extreme increase relative to healthy control levels (10 pM; Table S. 7), it does not appear to be physiologically impossible: VEGF receptor occupancies were high (84% VEGFR1; 69% VEGFR2) but not completely saturated at the PAD control point (Fig. 4, S. 1 and S. 2). Additionally, in vitro studies have shown that it takes more than 40 ng/ml (870 pM) of VEGF to saturate receptor binding on microvascular endothelial cells, albeit at apparently higher total receptor densities (13).

Although the decision to fix the protein secretion rates of the normal compartment at healthy control rates was intended to keep the interstitial levels of VEGF and sVEGFR1 close to healthy levels in noncalf tissues, interstitial sVEGFR1 concentration in the normal compartment ([sR1]IS,Normal) inevitably increased away from its healthy benchmark to 81 pM at the PAD control point, due to its sensitivity to the calf secretion rate of sVEGFR1 (qsR1,Calf). Although the mass balance of sVEGFR1 in and out of the tissue compartments (in particular the relative magnitudes of the permeability vs. lymphatic vs. secretion vs. internalization flows) remains to be validated in vivo, current simulations suggest that systemic sVEGFR1 concentrations (e.g., in the body compartment) may be affected by localized changes to production rates (e.g., from the calf compartment).

Altered Basal Molecular Distributions: Implications for Impaired Angiogenesis in PAD

Several key differences in the basal molecular distribution profiles at the PAD control point were noted in comparison with our previous healthy control distributions (81).

Among them, simulations predicted that in PAD, the sVEGFR1-VEGF complex would constitute higher percentages of both total VEGF (91% vs. 77% in plasma) and total sVEGFR1 (13% vs. 5% in plasma), partly because complex formation provided a synergistic route for shuttling interstitial VEGF into the blood, and consequentially was a crucial means of elevating plasma VEGF to PAD target range, aside from independent upregulation in VEGF production. As in our previous healthy model, predictions were much higher than experimentally measured by Belgore et al. (6), suggesting that either significant in vivo binding partners of VEGF [e.g., sVEGFR2 (20), plasma fibronectin (77)] and sVEGFR1 [e.g., PlGF (72)] were critically neglected in our computational model, or alternatively, that the VEGF-sVEGFR1 complexed fractions were underestimated experimentally, with the counterimplication that measurements of free VEGF should not be assumed to approximate circulating levels of total VEGF (6). Both possibilities indicate that further experimental investigation is needed to assess the full effects of circulating carriers, including sVEGFR1, in the physiological regulation of VEGF bioavailability, especially in the context of VEGF-dependent diseases.

As a direct consequence of the sizeable upregulations in VEGFR2 expression (extrapolated from mice data) and VEGF production (tuned to fit the PAD target concentration of plasma VEGF) in the calf compartment, the PAD calf was predicted to have an increased margin in its overall proangiogenic signaling potential [approximated by the quantity of VEGF-activated VEGFR2 complexes minus VEGF-activated VEGFR1 complexes, at 0.6 in PAD vs. 0.009 pmoles/(cm3 of calf tissue) in healthy]. This current prediction in which the angiogenic compensatory response was intact up to the level of proangiogenic signaling complex formation, even in the presence of overexpressed sVEGFR1, points to a need to explore alternative hypotheses. One alternative hypothesis suggests that impaired angiogenesis may be rooted in signaling pathway disruptions downstream of ligand-receptor binding in PAD. Early evidence showed decreased downstream signaling through AKT and increased phosphoinositide (PI) phosphatase activity in a mouse model of PAD, suggesting that PI3-kinase/AKT pathway-dependant VEGF resistance may be involved in the impairment of ischemia-induced angiogenesis (29). A second hypothesis is that sVEGFR1 may interfere with angiogenic signaling through dominant-negative heterodimerization with cell surface VEGFR1 and VEGFR2 monomers, rather than through VEGF ligand sequestration (3, 12, 33, 37). Experimental data on sVEGFR1-VEGFR heterodimerization are presently lacking, and theoretical considerations involving significant extensions of our previous computational models of cell surface VEGF receptor dimerization may be warranted.

PAD-specific System Responses to Upregulated Protein Expressions: Implications for Therapeutic Angiogenesis

Comparing the basal flow profiles for the soluble species at the PAD control point (Fig. 3) with healthy flow profiles [Fig. 4 in (81)], directionally reversed net flows of free sVEGFR1 across the PAD calf endothelium and of the sVEGFR1-VEGF complex across the normal endothelium were observed. We postulated that these flow reversals played a role in the altered steady-state system responses to perturbed protein secretion rates in certain cases (Fig. 2). At high calf secretion rate of VEGF, increasing calf secretion rate of sVEGFR1 began to increase interstitial VEGF globally (normal compartment); at high calf secretion rate of sVEGFR1, increasing calf secretion rate of VEGF began to decrease circulating sVEGFR1 (blood compartment). In general, most of the net flow directions in this PAD model were conserved relative to the healthy model, such that we expect the preservation of the mechanism whereby sVEGFR1 shuttles interstitial VEGF into the plasma (80).

Experimental validation of the concentration gradients, flow profiles, and shuttling mechanism in the context of PAD would be critical for accurate prediction of the systemic responses to therapeutic angiogenesis strategies including the intramuscular delivery of VEGF (similar to increasing qTotalVEGF,Calf) and the intramuscular inhibition of sVEGFR1 expression (similar to decreasing qsR1,Calf). For instance, the molecular interplay between VEGF and sVEGFR1, which was accountable for the proposed shuttling mechanism, may carry significant implications for the optimization of proangiogenic therapies; Fig. 2 predicted that the intramuscular delivery of VEGF to the ischemic calf would be less effective when calf production of sVEGFR1 is high, i.e., the slope of increase in attainable VEGF concentration is less steep at the target site (calf interstitium) but more steep in nontargeted sites (plasma and normal interstitium), which may inadvertently lead to systemic side effects (e.g., luminal exposure to plasma VEGF may exacerbate the angiogenic growth of primary atherosclerotic plaques; abluminal exposure to globally elevated interstitial VEGF may incite angiogenic imbalances in noncalf body tissues).

Systems Biology Aids Design of Integrative and Patient-specific Diagnostic Tools for Disease

Through an integrative examination of how individual parameter changes associated with PAD pathology differentially affected each molecular distribution, we proposed that surveillance of just the circulating (plasma) concentration of either free VEGF alone or free sVEGFR1 alone would not provide a reliable assessment of the intramuscular signaling potential, e.g., an increased signaling potential may be reflected in a simultaneous increase in [V]pl at times, as in the case of increasing calf production of VEGF (run X = qVEGF in Fig. 4), yet accompanied by a concurrent decrease in [V]pl in other times, as in the case of increasing VEGFR expression (run X = [R1],[R2] in Fig. 4). It appeared that no combination of markers, other than the direct measurement of both calf interstitial VEGF concentration and total surface VEGFR densities, could allow the reliable estimation of intramuscular signaling potential (Fig. S. 3), for which there is no currently available noninvasive technology for measuring in situ under clinical settings. However, we also proposed that the combinatorial analysis of multiple markers showed promise for patient-specific etiological differential diagnosis at the molecular level, e.g., the temporal tracking of fluctuations in the two circulating markers, [V]pl and [sR1]pl, could differentiate between altered VEGFR expression, VEGF secretion, or sVEGFR1 secretion as the main molecular disturbance to the VEGF system within PAD-affected muscle tissues (Fig. 4K); a survey of the free versus bound fractions of plasma VEGF or sVEGFR1 could further corroborate the identification of altered VEGFR expression. Subsequent computational simulation of the underlying pathological molecular adaptations, as diagnosed from circulating angiogenic markers, can then potentially predict the intramuscular signaling potential and assist the prognostic evaluation of disease.

Other Limitations

Several simplifications made in deriving the current computational model that can potentially affect its clinical translatability need to be further considered, including the exclusion of nonangiogenic or nonendothelial functions of VEGF; the omission of complex interactions of VEGF system molecules with other angiostatic factors, receptor antagonists, and upstream transcriptional regulators; and the extrapolation of human parameter values from animal experiments.

Physiological functions of VEGF independent of angiogenesis and potentially crucial in PAD pathogenesis, which are not current modeled, might account for some discrepancy between our computational predictions and clinical measurements in PAD. The current model of VEGF interactions in PAD can be extended to include inflammatory pathways in addition to angiogenic pathyways. PAD is typically a manifestation of atherosclerosis, where lower extremity arterial occlusions by atherosclerotic plaques lead to downstream skeletal muscle ischemia (49, 73). Highlighting the role of vascular inflammation in the atherosclerosis of PAD, studies have correlated serum elevations of an inflammatory biomarker, C-reactive protein, with disease severity (49, 73). Although the current model focuses on the angiogenic role of VEGF in downstream ischemic muscle in PAD, it neglects the potential contribution of VEGF to the vascular inflammation that drives upstream atherosclerotic plaque formation in PAD (7, 15, 53). Reports have also indicated a role for VEGF in mediating smooth muscle cell proliferation, another critical process in atherosclerotic pathogenesis, through induction of sympathetic innervation (16, 47) or autocrine loops (at >0.8–1 ng/ml or >17–22 pM) (58). The various pathways through which VEGF can be involved in PAD pathogenesis (angiogenic vs. inflammatory; endothelial vs. nonendothelial) underscore the difficulty of modeling such a fundamentally heterogeneous and polygenic disease.

In addition to the VEGF system of ligands and receptors, a multitude of other endogenous molecules play a role in the physiological complexity of angiogenic regulation, including angiostatic factors [e.g., angiostatin and matrix metalloproteinases (18)], receptor antagonists [e.g., TIMP-3, a VEGR2 antagonist (82)], and upstream transcriptional regulators [e.g., HIF-1 (67)]. Some of these systems have previously been independently modeled in our laboratory (62, 74), and ongoing efforts seek to integrate these modules into multiscale models of angiogenesis (61).

The current PAD model assumed pathological levels of cell surface receptor protein concentrations that were extrapolated from experimental data on mouse models of PAD, due to a lack of in vivo data from human PAD patients. This could potentially introduce inaccuracies in the current model parameterization, since surgically induced acute hindlimb ischemia in diabetic mice may not be able to fully recapitulate the chronic intracellular maladaptations that are symptomatic of human PAD (15, 29). This limitation of data source should be addressed at the earliest availability of in vivo human data of protein expression levels in PAD.


In this study, we present a novel computational framework for simulating VEGF and sVEGFR1 in PAD patients. The model was constructed based on experimentally derived parameter values where available and physiologically based assumptions where experimental measurements were lacking. The apparent discrepancies between predicted phenomena from our current model and experimental data are reflective of an incomplete experimental understanding of the system and point to areas where the missing information lies and further physiological experiments are needed, as summarized in Table 10. In summary, more comprehensive clinical data from human PAD patients, of both circulating and interstitial levels of various angiogenic markers, are still needed to validate the simulated molecular distribution and flow profiles of our PAD control model. Furthermore, the alternate role of sVEGFR1 in dominant-negative with surface VEGFRs needs to be investigated to conclusively rule out sVEGFR1 contribution to the impaired angiogenesis of PAD. Nevertheless, this study was a proof-of-concept for the value of an in silico model of whole body VEGF systems biology as a platform for conducting personalized medicine and for evaluating treatment strategies.

View this table:
Table 10.

Comparing computational predictions with experimental data


This work was supported by the National Heart, Lung, and Blood Institute (NHLBI) Grants R33-HL087351 and R01-HL079653, a Natural Sciences and Engineering Research Council of Canada postgraduate PGS-M scholarship (to F. T. H. Wu), and NHLBI Grant R00-HL-093219 (to F. Mac Gabhann).


No conflicts of interest, financial or otherwise, are declared by the author(s).


We thank Dr. Amina A. Qutub, David Noren, Prakash Vempati, and all members of the Popel Lab for helpful discussions.


  1. 1.
  2. 2.
  3. 3.
  4. 4.
  5. 5.
  6. 6.
  7. 7.
  8. 8.
  9. 9.
  10. 10.
  11. 11.
  12. 12.
  13. 13.
  14. 14.
  15. 15.
  16. 16.
  17. 17.
  18. 18.
  19. 19.
  20. 20.
  21. 21.
  22. 22.
  23. 23.
  24. 24.
  25. 25.
  26. 26.
  27. 27.
  28. 28.
  29. 29.
  30. 30.
  31. 31.
  32. 32.
  33. 33.
  34. 34.
  35. 35.
  36. 36.
  37. 37.
  38. 38.
  39. 39.
  40. 40.
  41. 41.
  42. 42.
  43. 43.
  44. 44.
  45. 45.
  46. 46.
  47. 47.
  48. 48.
  49. 49.
  50. 50.
  51. 51.
  52. 52.
  53. 53.
  54. 54.
  55. 55.
  56. 56.
  57. 57.
  58. 58.
  59. 59.
  60. 60.
  61. 61.
  62. 62.
  63. 63.
  64. 64.
  65. 65.
  66. 66.
  67. 67.
  68. 68.
  69. 69.
  70. 70.
  71. 71.
  72. 72.
  73. 73.
  74. 74.
  75. 75.
  76. 76.
  77. 77.
  78. 78.
  79. 79.
  80. 80.
  81. 81.
  82. 82.
  83. 83.
View Abstract