|Home | About | Journals | Submit | Contact Us | Français|
Most physiological processes are subjected to molecular regulation by growth factors, which are secreted proteins that activate chemical signal transduction pathways through binding of specific cell-surface receptors. One particular growth factor system involved in the in vivo regulation of blood vessel growth is called the vascular endothelial growth factor (VEGF) system. Computational and numerical techniques are well-suited to handle the molecular complexity (the number of binding partners involved, including ligands, receptors, and inert binding sites) and multi-scale nature (intra-tissue vs. inter-tissue transport and local vs. systemic effects within an organism) involved in modeling growth factor system interactions and effects. This paper introduces a variety of in silico models that seek to recapitulate different aspects of VEGF system biology at various spatial and temporal scales: molecular-level kinetic models focus on VEGF ligand-receptor interactions at and near the endothelial cell surface; meso-scale single-tissue 3D models can simulate the effects of multi-cellular tissue architecture on the spatial variation in VEGF ligand production and receptor activation; compartmental modeling allows efficient prediction of average interstitial VEGF concentrations and cell-surface VEGF signaling intensities across multiple large tissue volumes, permitting the investigation of whole-body inter-tissue transport (e.g., vascular permeability and lymphatic drainage). The given examples will demonstrate the utility of computational models in aiding both basic science and clinical research on VEGF systems biology.
At the molecular and cellular levels, growth factors are extracellularly secreted polypeptides which, upon binding to specific cell-surface target receptors, trigger intracellular signal transduction pathways that regulate cell proliferation, differentiation and survival (Lodish et al., 2004). At the tissue and organ levels, growth factors are responsible for orchestrating many physiological processes in complex multicellular organisms.
Of utmost importance in human physiology and pathology is the process known as angiogenesis – the growth of new capillaries or microvessels from pre-existing blood vasculature – which critically supports organogenesis during embryonic development (Haigh, 2008); physiological growth and repair in adult tissues (such as in wound healing (Bao et al., 2009), muscular adaptation to exercise (Brown and Hudlicka, 2003), or endometrial regeneration (Girling and Rogers, 2005)); as well as the malignant growth of tumor tissues (Kerbel, 2008).
Sprouting angiogenesis is a well-coordinated and complex cascade of molecular, cellular, and tissue-level events (Qutub et al., 2009): First, tissue ischemia is converted into a chemical cue for angiogenesis, as hypoxic cells (e.g., cancer cells in tumor tissue or myocytes in ischemic muscle) transcribe and secrete growth factors in response to hypoxia-inducible factor 1 (HIF1) activation. The growth factor ligands then diffuse throughout the extracellular fluid, where some become sequestered to matrix proteoglycans and others bind cell-surface receptors on the capillary endothelium. Cell surface receptor-bound ligands initiate vessel sprouting by turning quiescent endothelial cells into migratory tip cells. Extracellular matrix-bound and freely diffusing ligands form chemotatic gradients that guide the migration of tip cell filopodia in the capillary sprout. In the next two sections, we further explore the multi-scale nature and molecular complexity of angiogenic growth factor interactions.
Many growth factor systems are involved in angiogenic regulation, including the vascular endothelial growth factor (VEGF) system of at least 5 ligands (VEGF-A, PlGF, VEGF-B, VEGF-C, VEGF-D) and 3 receptors (VEGFR1, VEGFR2, VEGFR3) (Roy et al., 2006; Mac Gabhann and Popel, 2008); the fibroblast growth factor (FGF) system of at least 18 ligands (FGF1 to FGF10 and FGF16 to FGF23) and 4 receptors (FGFR1 to FGFR4) (Beenken and Mohammadi, 2009); the angiopoietin (Ang) system of at least 4 ligands (ANG1 to ANG4) and 2 receptors (TIE1 and TIE2) (Augustin et al., 2009); the platelet-derived growth factor (PDGF) system of at least 4 ligands (PDGF-A to PDGF-D) and 2 receptors (PDGFR-α and PDGFR-β) (Andrae et al., 2008); and the insulin-like growth factor (IGF) system of at least 2 ligands (IGF1 and IGF2) and 2 receptors (IGF1R and IGF2R) (Mazitschek and Giannis, 2004; Pollak, 2008).
There are organizational similarities between these growth factor systems – all of the above receptors except for the IGF2R are transmembrane receptor tyrosine kinases (RTKs) that are activated by ligand-induced dimerization and transphosphorylation of tyrosine residues (Gschwind et al., 2004); co-receptors (e.g., neuropilin-1 (NRP1) for VEGFRs and syndecan-4 for FGFRs) and endothelial integrins (e.g., αvβ5 for VEGFRs and αvβ3 for FGFRs) often modulate receptor signaling (Simons, 2004); heparan sulfate proteoglycans (HSPGs) are involved in the extracellular matrix sequestration of ligands in the VEGF, PDGF, FGF systems (Roy et al., 2006; Andrae et al., 2008; Beenken and Mohammadi, 2009).
Intracellularly, details are also emerging on the convergent and integrative crosstalk between and within growth factor systems in angiogenic signaling. Among the overlap in their transcriptional profiles downstream of RTK activation, all aforementioned growth factor systems can activate the canonical Ras-MAPK signaling pathway (Simons, 2004). VEGF and FGF2 were observed to induce the expression of each other (Simons, 2004). Within the VEGF system, the existence of heterodimeric ligands (e.g., VEGF-A/PlGF and VEGF-A/VEGF-B) and heterodimeric receptors (e.g., VEGFR1/VEGFR2) is expected to introduce new signal transduction pathways in addition to those downstream of classic homodimeric ligand-receptor activation (Mac Gabhann and Popel, 2008; Cao, 2009). While VEGFR1 is mainly a negative regulator of angiogenesis, its possible pro-angiogenic roles are suspected to involve the intermolecular transphosphorylation of VEGFR2 by PlGF-activated VEGFR1 (Autiero et al., 2003).
In this paper, we will introduce computational frameworks that are well suited for the quantitative modeling of highly complex molecular interaction networks such as that of the VEGF system in angiogenesis. While our examples focus on the VEGF system, the mathematical frameworks are generally adaptable for any of the organizationally similar growth factor systems introduced above, with the potential of further integration between the VEGF, FGF, PDGF, IGF, Ang-Tie system modules themselves.
In understanding the biology of angiogenic growth factors, it is of equal importance to identify the key molecular players and to distinguish where in the body the molecular interactions take place. The spatial range of activity can vary between growth factors (or between isoforms of the same growth factor) depending on their propensity for intra-tissue and inter-tissue transport.
The intra-tissue transport of a secreted growth factor – i.e., its diffusive and convective transport within the extracellular matrix (ECM) – is dependent on the growth factor’s molecular size, the pore sizes of the ECM, and its chemical affinity for ECM proteoglycans. For instance, heparin-binding affinity of VEGF, which determines the extent of its sequestration by ECM heparan sulfate proteoglycans, is traditionally thought to be encoded in the VEGF-A gene on exons 6 and 7 (Harper and Bates, 2008). Hence, the pro-angiogenic VEGF121 and anti-angiogenic VEGF121b splice isoforms which skip exons 6 and 7 are mostly freely diffusible once secreted into the extracellular fluid; while the higher molecular-weight isoforms VEGF145(b), VEGF148, VEGF165(b), VEGF183(b), VEGF189(b), and VEGF206 have progressively higher heparin-binding affinity due to their inclusion of increasingly greater portions of exons 6 and 7 (Harper and Bates, 2008). Yet these higher molecular-weight isoforms in their matrix-bound state can be subjected to proteolytic cleavage by plasmin or matrix metalloproteinases (MMPs), which releases active fragments of 110 to 113 amino acids in length and with similar angiogenic properties as VEGF121 (Ferrara and Davis-Smyth, 1997; Lee et al., 2005; Qutub et al., 2009).
On a larger scale, the inter-tissue transport of a growth factor is first affected by its rate of entry into or exit from the blood or lymphatic vasculature, i.e., its permeability through the blood capillary endothelium or its lymphatic drainage rate from interstitial spaces. Once in the bloodstream or lymphatic fluid, the inter-tissue transport of a growth factor may be further facilitated by specific carrier proteins or circulating cells. For VEGF, its potential carriers in blood include soluble forms of its normal receptors (sVEGFR1, sVEGFR2, sNRP1) (Gagnon et al., 2000; Ebos et al., 2004; Sela et al., 2008), plasma fibronectin (Wijelath et al., 2006), as well as platelets (Verheul et al., 2007).
All together, these transport properties influence the distance over which a growth factor can signal. Autocrine signaling occurs when growth factors act upon the same cells that produced them; juxtacrine growth factors act upon adjacent cells after secretion; paracrine growth factors diffuse through the extracellular fluid and target cells within the same tissue but of a different cell type than that which secreted them; whereas growth factors that are transported through the bloodstream to distant target tissues act in an endocrine manner (Lauffenburger and Linderman, 1993; Lodish et al., 2004).
Angiogenic VEGF signaling occurs predominantly in a paracrine fashion: epithelial cells in fenestrated organs (e.g. glomerular podocytes), mesenchymal cells (e.g., skeletal myocytes in ischemic muscles), vascular stromal cells (e.g., pericytes and smooth muscle cells) and hypoxic tumor cells are all known to secrete VEGF, which then diffuses to and activates the VEGF receptors on nearby endothelial cell surfaces (Maharaj and D'Amore, 2007; Kerbel, 2008). However, there is also evidence of autocrine VEGF signaling loops in VEGF-producing endothelial cells (Martin et al., 2009). Furthermore, autocrine VEGF signaling involving intracellular VEGF receptors (“intracrine signaling”) has been documented in breast carcinoma cells (Lee et al., 2007) and haematopoietic stem cells (Gerber et al., 2002); although in these contexts, VEGF functions as a cell survival signal rather than a pro-angiogenic factor.
While there have not been formally established specific endocrine functions of VEGF in normal physiology, aberrantly high circulating levels of VEGF may have deleterious systemic effects. In clinical trials administering VEGF via intravascular infusion to stimulate therapeutic angiogenesis for ischemic muscle diseases, unintended side effects of the high systemic VEGF concentrations such as hypotension and macular edema have been transiently and sporadically observed (Collinson and Donnelly, 2004). The VEGF concentrations in the plasma of cancer patients are also known to be several-fold higher than healthy baseline levels, although it is uncertain whether the tumors are themselves the source of the elevated circulating VEGF, or whether conversely the elevated circulating VEGF triggered the malignant growth of tumors (Kut et al., 2007; Stefanini et al., 2008). Therefore, a complete understanding of VEGF biology – including its pathogenic role in cancer and its therapeutic potential in ischemic diseases – necessitates an appreciation of the dynamic distribution of VEGF in the human body (Kut et al., 2007).
The computational models presented in this paper are complementary, capturing the biology of VEGF interactions at different length scales: the molecular-level kinetic models can predict the local intensity of VEGF-VEGFR complex formation on endothelial cell surfaces as a marker of angiogenic activation or vessel sprouting initiation; the meso-scale single-tissue 3D models involving multiple cell types can predict the spatial gradients of VEGF in the extracellular space and simulate paracrine signaling distributions that guide capillary sprout migration; and the multi-tissue compartmental models can be used to simulate whole-body VEGF distributions and to investigate the possibility of endocrine VEGF effects.
Figure 1 summarizes the multi-scale nature and complexity of the molecular interactions involved in the systems biology of the VEGF ligand-receptor system. In the following sections, we introduce models for investigating emergent behavior at progressively higher spatial scales: from molecular (subcellular) level models, to meso-scale (intra-tissue) models, to whole-body (inter-tissue) models. The chosen examples also illustrate the versatility of computational modeling for investigating basic science questions (e.g., simulating the molecular mechanisms underlying the PlGF-VEGF synergy and NRP1-VEGFR2 synergy) and assisting in the design of translational medicine (e.g., comparing the therapeutic efficacy of cell-based versus protein delivery of VEGF; optimizing the dose of anti-VEGF therapy).
In our first two examples, models were developed to investigate specific molecular functions and interactions of key players in the VEGF ligand-receptor system – PlGF and NRP1 – by recapitulating in vitro experiments. The spatial scope of these models focused on molecular behaviors near the endothelial cell surface, including extracellular ligand diffusion and cell-surface ligand-receptor binding.
Mathematical theory and formulations for kinetic modeling of cell surface ligand-receptor binding and cell-surface receptor/ligand trafficking have been presented in classical texts (Lauffenburger and Linderman, 1993). The standard description for the binding kinetics of ligand L to receptor R to form complex C involves characterization of the complex association and dissociation rate constants kon and koff:
Endocytotic internalization of free receptors and complexes are generally characterized by first-order rate constants:
Free receptor insertion rates are typically introduced through zero-order source terms; in the following models, they are chosen to maintain a steady total population (free and bound) of receptors in the absence of added ligand.
Our first case study sought to decipher the molecular mechanisms behind PlGF’s observed ability to augment the angiogenic response to VEGF-A in in vitro assays for endothelial cell survival, proliferation and migration. Details and full references can be found in (Mac Gabhann and Popel, 2004). These two members of the VEGF family have different receptor-binding properties: VEGF-A (hereinafter referred to as simply “VEGF”) binds with both VEGFR1 and VEGFR2; while PlGF only binds VEGFR1. Two proposed mechanisms for the PlGF-VEGF synergy were: (a) “ligand shifting”, where PlGF displaces VEGF from VEGFR1, effectively freeing more VEGF to bind the more pro-angiogenic VEGFR2; and (b) PlGF-VEGFR1 signaling, where PlGF activation of VEGFR1 may transduce qualitatively different (pro-angiogenic) signals than that from VEGF activation (generally inhibitory of angiogenic signaling). An in silico model was thus constructed to quantify these mechanistic contributions to the VEGF-PlGF synergy.
The in silico model formulation mimicked in vitro assay geometry and conditions as illustrated in Fig. 2A. At the bottom of a cell culture well, a confluent layer of endothelial cells expressing receptors VEGFR1 and VEGFR2 on the surface (z=0) was exposed to the fluid media. Into the cell culture media (from z=0 to z=h), ligands were administered at time zero – either VEGF alone (“PlGF-” case) or VEGF and PlGF (“PlGF+” case) – to assess the synergistic effects of PlGF. Mathematically, each molecular species was represented by either a volumetric concentration (V for VEGF, P for PlGF) or cell-surface concentration (R1 for VEGFR1, R2 for VEGFR2, VR1 for VEGF·VEGFR1 complex, PR1 for PlGF·VEGFR1 complex, VR2 for VEGF·VEGFR2 complex), using the continuum approach. The initial value problem was essentially a single spatial dimension problem (z-direction) as molecular concentrations are assumed uniform in the plane parallel to the endothelial cell surface.
Coupled diffusion and reaction equations (Eqns. 3 and 4 respectively) were used to describe the time evolution of extracellular ligand transport and cell-surface molecular binding interactions. D represented ligand diffusivity (cm2/s); sR is the insertion rate or receptor R (mol/cm2/s); kint is the receptor or complex internalization rate (s−1); kon and koff are the rate constants of complex association (M−1s−1) and dissociation (s−1).
Boundary conditions were given by Eqn. 5. qV is the endothelial secretion rate of VEGF (mol/cm2/s).
The initial conditions describing exogenous ligand administration and total receptor densities were specified by Eqn. 6.
The predicted final concentration of cell-surface ligated VEGF·VEGFR2 complexes served as a surrogate marker for achieved angiogenic response.
Representative values used for the model parameters were based on experimental literature as detailed in (Mac Gabhann and Popel, 2004). Numerical solution of the coupled nonlinear differential equations (Eqns. 3–6) was achieved using an iterative implicit finite-difference scheme. Sample results are shown in Fig. 2B. Crucially, the simulations predicted that PlGF addition only increased VEGF·VEGFR2 complex formation up to 5% at the peak percentile change between PlGF+ and PlGF− cases; i.e., the “ligand shifting effect” is expected to be minimal. On the other hand, the transient increase in total ligated VEGFR1 complexes was as high as 43%, during which, the magnitude increase in PlGF·VEGFR1 complexes exceeded that of the decrease in VEGF·VEGFR1 complexes. This would suggest that “VEGFR1 signaling” plays a more prominent role in the observed PlGF-VEGF synergy: that PlGF critically alters the VEGFR1 signaling profile in both the absolute quantity of signaling VEGFR1 complexes and the signaling quality of VEGFR1 complexes (elevated pro-angiogenic PlGF·VEGFR1 signaling and reduced modulatory VEGF·VEGFR1 signaling). Experimental support for these computational predictions have been found (Autiero et al., 2003), as intermolecular crosstalk was reported to occur downstream of PlGF-VEGFR1 binding, leading to the transphosphorylation of VEGFR2 and amplification of pro-angiogenic VEGF·VEGFR2 signaling.
Our second case study investigated another molecular player in the systems biology of VEGF, neuropilin-1 (NRP1). As illustrated in Fig. 2C, NRP1-binding affinity is conferred to the higher molecular-weight isoforms of VEGF such as VEGF165 predominantly through transcription of exon 7; the VEGF121 isoform, which lacks exon 7, is generally considered to have negligible affinity for NRP1. Because the NRP1-binding domains of VEGF165 do not overlap with its VEGFR-binding domains, VEGF165 can act as a bridge in the formation of a ternary complex: VEGF165·VEGFR2·NRP1. A reduced interaction network model involving VEGF121, VEGF165, VEGFR2 and NRP1 (Fig. 2C) was thus constructed to quantify the role of VEGF165-bridged VEGFR2-NRP1 coupling in generating VEGF isoform-specific angiogenic responses. Details and full references can be found in (Mac Gabhann and Popel, 2005).
Geometrically, the experimental setup again involved a confluent layer of endothelial cells on the bottom of a cell culture well, as in Fig. 2A. As before, the initial value problem in one spatial dimension was formulated as a system of coupled diffusion and reaction equations (Eqns. 7 and 8–13 respectively). A new parameter, kC, represents the coupling rate between VEGFR2 and NRP1 via VEGF165.
Boundary conditions were given by Eqn. system 14.
The predicted final VEGF·VEGFR2 concentration again served as a marker for the strength of pro-angiogenic signal transduction.
The experimental and theoretical derivations of model parameter values are detailed in (Mac Gabhann and Popel, 2005). Numerical solution of the coupled nonlinear differential equations (Eqns. 7–14) was achieved using an iterative implicit finite-difference scheme. Sample results are shown in Fig. 2D. The in silico modeling of stepwise reaction kinetics (Fig. 2C) allowed the prediction of differential anti-angiogenic efficacies from therapeutic interference of two distinct aspects of NRP1 function: VEGF binding and VEGFR2 coupling. Simulated blockade of VEGF165 binding to NRP1 (blocking reactions ‘Rx1’ and ‘Rx2’ in Fig. 2C) resulted in the convergence of VEGF165 response to that of VEGF121 in terms of VEGFR2 activation; however, simulated blockade of NRP1-VEGFR2 coupling (blocking reactions ‘Rx1’ and ‘Rx3’ in Fig. 2C) converted NRP1 into a VEGF165 sink (through intact reaction ‘R2’ in Fig. 2C), further reducing VEGF165 response to below that of VEGF121 (Fig. 2D).
The study of VEGF binding to receptors on cells in vitro, and the validation of the VEGF kinetic interaction network between multiple ligands and multiple receptors, leads us to ask the question: how does this network behave in vivo? In sections 4 and 5 we will discuss the transport of VEGF between tissues and around the body, but here we will focus first on the behavior of VEGF in a local volume of tissue. This multicellular milieu requires significant additions to our model in order to accurately simulate the local transport of VEGF, including diffusion of VEGF ligands over significant distances, extracellular matrix sequestration and variable production rates of VEGF throughout the tissue. We place these all in an anatomically-based 2D or 3D multicellular tissue geometry.
The models can predict the creation of interstitial VEGF gradients due to the non-uniform nature of the tissue anatomy. This is of particular interest because VEGF is believed to be a chemotactic guiding agent for blood vessels, but also because local variability in VEGF concentration can lead to local variation in VEGF receptor ligation and signaling, allowing for focal activation of endothelial cells. The model framework can be adapted to most tissues; here we present a case with parameters specifically selected to represent a skeletal muscle experiencing ischemia (specifically, rat extensor digitorum longus, or EDL, for a rodent model of hindlimb artery ligation), and we describe how to computationally test several therapeutic interventions including gene therapy and exercise.
A cross-section (for 2D) or a volume of tissue (for 3D) is reconstructed from histological and other microanatomical information (Figure 3A–C). The major relevant features of the tissue are the blood vessels, the parenchymal cells (from here on, we will assume these are skeletal myocytes, i.e. long multinucleated cells) and the interstitial space between them. From a computational modeling point of view, the tissue comprised of volumes and surfaces, defined as those portions of the tissue where molecules can move in all directions (volumes) and those portions where the movement of molecules is restricted to a plane (e.g. receptors inserted in cell membrane can move only laterally). There are three major volumes of the tissue for our purposes: the vascular space (i.e. inside the blood vessels, determined by the density of blood vessels and their diameters); the intracellular space (whether inside of parenchymal cells or endothelial cells); and the interstitial space between cells, which is itself divided into three volumes (each of which is not contiguous), based on the density of the fibrous matrix present: the extracellular matrix, and basement membrane regions surrounding the endothelial cells and the myocytes.
There are two major surfaces; again, these are not contiguous. First, the combined cell surfaces of the skeletal myocytes, which are assumed to be cylindrical (diameter 37.5 µm, consistent with rat histology), and arranged in a regular hexagonal grid formation, accounting for almost 80% of total tissue cross-sectional area. VEGF is secreted from the myocytes’ surfaces. Second, the surface of the endothelial cells that make up the blood vessels, specifically, the abluminal surface (the luminal surface faces the blood stream and we neglect it for now). Again, the blood vessels are assumed to be cylindrical, and although most (but not all) are parallel to the muscle fibers, they do not occupy every possible position between fibers, but instead have a stochastic, nonuniform arrangement (based on experimentally measured capillary-to-fiber ratios, capillary-to-fiber distances and histology) occupying 2.5% of total tissue volume (leaving ~18% as interstitial space). On this endothelial surface, VEGF receptors are expressed. Thus, VEGF must diffuse from the myocyte surface where it is secreted, through basement membranes and extracellular matrix, to the endothelial surface where it ligates its cognate receptors.
To model tissues at the meso-scale, we use the above microanatomical information as an input to a set of integrated models of blood flow, oxygen transport, and VEGF transport (Fig. 3D).
Blood flow and hematocrit value calculations are based on Pries et al.’s two-phase continuum model (Pries and Secomb, 2005), and reduces to a system of nonlinear algebraic equations (two per vessel) that are solved iteratively (Ji et al., 2006). The Fahraeus-Lindqvist effect and non-uniform hematocrit distribution at vascular bifurcations are included in the blood flow model. Higher blood flow rates are used for exercising conditions, to represent the increased perfusion (and enhanced oxygen delivery) to exercising muscles. In addition, exercise-trained rats have higher average capillary blood velocity.
Oxygen transport in the tissue in detailed in (Ji et al., 2006) and (Goldman and Popel, 2000). Oxygen arrives in the tissue via the blood vessel network, and the partial pressure of oxygen in the vessels, Pb, is described by
where vb is mean blood velocity; αb is oxygen solubility in blood; HD is the discharge hematocrit; are the oxygen binding capacity and oxygen saturation of the red blood cell; ξ is the longitudinal position in the vessel; R is vessel radius; Jwall is the oxygen flux across the vessel walls (i.e. into the tissue). Oxygen diffuses across the endothelial cells, and freely throughout the tissue (both interstitial and intracellular). Within the cells, it may also be consumed by binding to Myoglobin (Mb). The local partial pressure of oxygen in the tissue, P, is described by
where DO2 and DMb are the diffusivities of oxygen and myoglobin; α is the oxygen solubility in tissue; CMbbind and SMb are the oxygen binding capacity and oxygen saturation of myoglobin; and M(P) represents Michaelis-Menten kinetic consumption of oxygen.
The VEGF ligands, VEGF120 (rodent form of VEGF121) and VEGF164 (rodent form of VEGF165), can both diffuse through the interstitium following secretion, however the longer isoform also binds to glycoproteins in the extracellular matrix, becoming reversibly sequestered. The equations are thus identical to those of section 2, with the addition of binding and unbinding terms:
where Ci and Cj are concentrations of two interstitial molecules, i and j. In the rat EDL model, Ci=[V120] or [V164]; Cj=[GAG]. Concentration of proteins in the thin endothelial or myocyte basement membranes is given by an equation of the form:
where si is the secretion rate from the cell (typically from myocytes); Rm and Rim are the concentrations of receptor m and of the i-m complex on the cell surface (typically on endothelial cells); Jout is the Fickian diffusive flux from BM to ECM of VEGF; and dBM is the basement membrane thickness.
The ligand-receptor interactions that take place are precisely those that were outlined in section 2, and that will be used in Sections 4 and 5: VEGF120 and VEGF164 bind to VEGFR1 and VEGFR2, while only the longer isoform binds Neuropilin-1 and extracellular matrix. The general form of the receptor and receptor complex equations is therefore:
where sm and kint,m are the membrane insertion rate and internalization rate of receptor m; kcouple,m,n and kdissoc,mn are the kinetic rate of binding and unbinding of two surface receptors m and n to each other. Note in particular that the concentration of the ligand (Ci) in each case is the concentration in the basement membrane region closest to the receptor. Thus, the receptor occupancy varies from cell to cell across the capillary network. Examples of specific individual equations can be found in section 2.
The production and secretion of VEGF has been observed to be inducible by hypoxia (Forsythe et al., 1996; Qutub and Popel, 2006). Here, we use an empirical relationship (Mac Gabhann et al., 2007b) for the increase in the baseline secretion rate of VEGF (S0) based on the observed upregulation of VEGF mRNA and protein during hypoxia in cells and tissues (Jiang et al., 1996; Tang et al., 2004):
Intracellular VEGF is not included in these simulations; that includes both post-internalization VEGF and pre-secretion VEGF. In addition, we neglect the intravasation of VEGF into the bloodstream, either by endothelial cell secretion or through paracellular routes, e.g. permeability. Lymphatic transport of VEGF is also neglected. These additional transport routes could be accommodated in the above model structure with the addition of new surfaces or terms. Although endothelial VEGF production and parenchymal VEGFR expression have been observed in recent years (Lee et al., 2007; Bogaert et al., 2009), these are not included as part of these simulations; there is no technical obstacle to doing so.
It is important to note that the spatial averages of VEGF concentrations at the endothelial cell surface and of VEGFR activation in the meso-scale models match well with the values in single-compartment models (section 4) that do not include diffusion or VEGF gradients. Thus, it may be possible to calculate the average receptor activation using less computationally intensive compartment models, and use the meso-scale models to estimate the spatial gradients.
In order to improve the perfusion and healing of ischemic muscle tissue with impaired angiogenic response, several therapies have been suggested, typically involving the delivery of VEGF (one or more isoforms) to the muscle.
The first of these, gene therapy, increases the VEGF secretion by adding additional VEGF-encoding genes to the cells that are transfected. By transfecting multiple copies, or by judicious choice of VEGF promoters and enhancers in the new construct, significant increases in VEGF secretion can be obtained. We have modeled both uniform upregulation of VEGF (increasing VEGF secretion at every myocyte surface point in the model) and stochastic upregulation, in which each cell has a randomly increased VEGF production within a certain range (using the myonuclear density, we know the size of the myocyte surface that is under the control of each nucleus; thus, we can assign a random number to each region, that stays constant through the simulation) (Mac Gabhann et al., 2007a). These increases in VEGF production result in increased VEGFR2 activation, however the VEGF gradients are not significantly increased (Fig. 3E,F); in this case blood vessels might be induced to sprout, but have no directional cues. Further simulations restricting the VEGF transfection to a specific region of the muscle demonstrates increased VEGFR2 activation coupled with very high VEGF gradients towards the transfected tissue, but only in a narrow region between transfected and nontransfected tissue (Mac Gabhann et al., 2007a). This suggests that VEGF gene delivery needs to be effectively localized with a high degree of spatial accuracy to allow the gradients of VEGF to bring the new vessels to the affected volume.
Another route to bringing more VEGF to the tissue, and one which may allow for more spatial specificity, is the delivery of VEGF-overexpressing cells, e.g. myoblasts that will effectively integrated into the existing muscle and produce excess VEGF locally. To simulate this, we select specific myocytes in the model to overexpress VEGF, and distribute these distantly or close together (Mac Gabhann et al., 2006; Mac Gabhann et al., 2007a). That is, since the secretion rate of VEGF can have a different value for every spatial location on the myocytes surface in our model, we can upregulate VEGF in a specific subset of these cells.
For this therapy, we observe in the simulations both increased VEGFR2 binding and increased VEGF gradients (Fig. 3E,F), but only within approximately one to two myocyte diameters from the new VEGF overexpressing cells (Mac Gabhann et al., 2006; Mac Gabhann et al., 2007a). In addition, cells close together synergize while distant ones do not. In this way, we can see that a small number of cells, or cells distributed too broadly, would have a low probability of attracting perfusion from a neighboring region; however, a large mass of cells, at the right location, could serve as a local chemoattractant.
The results described in sections 3.2 and 3.3, for therapies reliant on VEGF upregulation alone, mirror the outcome of the several clinical trials of VEGF isoforms in humans for coronary artery disease (CAD) or peripheral artery disease (PAD); these trials have not had the success that was expected of them. Instead, the standard of care for PAD continues to be exercise, and it is this therapy that we consider next.
Exercise training in rats has been shown to not only restore the ability of hypoxic, ischemic tissue to upregulate VEGF following injury, but also increases the expression levels of the VEGF receptors (Lloyd et al., 2003). Thus, we used our model to simulate the exercise-dependent upregulation of both the ligands and the receptors, using experimentally measured increases (Ji et al., 2007; Mac Gabhann et al., 2007a; Mac Gabhann et al., 2007b).
In this case, we increase the secretion rate of VEGF isoforms from each point on the myocyte surface, during exercise; in addition, we increase the insertion rate of the VEGF receptors at every point on the endothelial cell surface at all times (as a result of exercise training). The results of these simulations are quite different from those before: first, during exercise, both the VEGFR2 activation and the VEGF gradients are increased, not just locally but across the upregulated tissue (Ji et al., 2007; Mac Gabhann et al., 2007a); second, during rest periods, while VEGF upregulation ceases and the occupancy of VEGFR2 returns to lower levels, the high VEGF gradients are maintained (Fig. 3E,F). This suggests that the activation step for attracting new blood vessels may be during a smaller window of time, while the guidance of the new vessel to its destination can take place continuously.
This observation that our current best strategy for PAD, exercise, increases both ligand expression and receptor activation, leaves us with the possibility of developing combined ligand-receptor therapy (especially for those who cannot exercise).
From section 3, we saw that the investigative use of VEGF models for therapeutic applications require larger spatial scale modeling of in vivo growth factor transport and effects at the tissue level. However, the full spatial resolution afforded by 3D modeling becomes computationally intensive when the volume of interest grows from tissue regions to whole tissues and organs. Here we demonstrate the strategy of compartmental modeling – where tissue fluid volumes (e.g., interstitial fluid volume) are approximated as well-mixed compartments of uniform protein concentrations, based on the assumption that diffusion occurs on faster timescales than that of molecular binding kinetics (Damkohler number <1 (Mac Gabhann and Popel, 2006)). While compartmental models cannot predict diffusion-limited VEGF gradients, they allow the prediction of average interstitial VEGF levels and average cell-surface VEGF-VEGFR binding within single tissues (section 4) or in multiple interconnected tissues (section 5).
In converting interstitial volumetric protein concentrations to their appropriate per-tissue-volume basis, it is useful to consider extravascular regions as porous media (Truskey et al., 2004) and to use the following standard definitions. ε is the porosity, defined as the fractional void space within the total tissue volume. Φ is the partition coefficient, representing the fractional interstitial fluid volume that is available or accessible to macromolecules, i.e., excluding all isolated or impenetrable pores. Together, they define the available volume fraction, KAV, which is used to convert interstitial concentrations from per-interstitial-fluid-volume (M) to per-tissue-volume (mol/(cm3 tissue)) basis:
Endothelial cell-surface receptor and complex densities also have to be converted from their per-surface-area units (moles per cm2 cell surface area) to per-tissue-volume (moles per cm3 tissue) basis for consistency within the mathematical equations given below. For this purpose, other geometrical attributes (e.g., microvessel density and diameters, surface area per endothelial cell, and endothelial surface area per tissue volume) have to be quantified to reflect the particular tissue architecture modeled, as documented in (Mac Gabhann and Popel, 2006) for breast tumor tissue geometry.
In this example of a single-tissue compartment model, three strategies of inhibiting different aspects of NRP1 functionality were simulated to predict their relative in vivo anti-angiogenic efficacy as cancer treatments in breast tumor tissues. The first strategy of inhibiting NRP1 expression (Fig. 4A), which could be achieved clinically through siRNA methods, was simulated by reducing the NRP1 insertion rate, sN. The second strategy of impeding VEGF binding to NRP1 (Fig. 4B), which could be implemented clinically through a PlGF fragment (P2Δ) that binds NRP1 exclusively, was modeled with additional equations representing the competition between VEGF and P2Δ for NRP1 (Eqn. 22) and corresponding modifications to the original NRP1 equation. The third strategy of blocking NRP1-VEGFR coupling while allowing NRP1-VEGF binding (Fig. 4C), which had been done experimentally using a NRP1 antibody, was modeled with additional equations for the antibody interactions (Eqn. 23) and corresponding modifications to the original VEGF165, NRP1, and VEGF165·NRP1 equations. Detailed references and full equations can be found in (Mac Gabhann and Popel, 2006).
Numerical solution of the full system of coupled nonlinear ordinary differential equations was achieved using a Runge-Kutta integration scheme (Mac Gabhann and Popel, 2006). Sample results are shown in Fig. 4D. The results predict that the third strategy of blocking NRP1-VEGFR coupling will induce the most effective and sustained inhibition of VEGFR2 ligation and signaling, with the differential advantage being most discernible in tumor types where endothelial expression of VEGFR1 is low. This case study demonstrates that model incorporation of molecularly detailed interaction networks allow the sophisticated optimization of therapeutic strategies, from the fine-tuning of therapeutic agent dosing to the customization of targeted molecular mechanisms according to the diseased tissue type with its characteristic receptor expression levels.
In this section we introduce models that compartmentalize the whole body into multiple tissue compartments. The two examples given below each have three compartments (Fig. 5): the blood; the tissue of interest (breast tumor tissue in section 5.2 and calf muscle in section 5.3); and the rest of the body, which we call the “normal compartment”. Future model extensions can further subdivide the normal compartment into individual organ compartments. In this paper we demonstrate that the modeling of two major inter-compartment transport processes – vascular permeability and lymphatic drainage – is essential in predicting whole-body pharmacokinetics of macromolecules.
Unlike small molecules that can easily traverse microvessel endothelia through interendothelial cleft junctions (Fu and Shen, 2003), the transendothelial permeation of macromolecules such as peptide growth factors occurs through caveolar structures and vesiculo-vacuolar organelles at much slower rates that are dependent on protein size and VEGF-induction (VEGF is also known as VPF, vascular permeability factor) (Garlick and Renkin, 1970; Feng et al., 2002; Fu and Shen, 2003). The equations relevant to protein transport between the available fluid volumes of the tissue compartment and the blood compartment via vascular permeability are:
where represents the microvascular permeability rate for the protein of interest going from the tissue (T) to the blood (B) compartment (cm/s) across STB (cm2), the total endothelial surface area (i.e. the tissue-blood interface).
Microvascular permeability is modeled as a passive bidirectional transport process, represented on Figure 5 by a double-arrow, with identical intravasation and extravasation rates, . The in vivo value of kp for a given protein is rarely found in the literature. Instead, we extrapolate from calibration curves (Garlick and Renkin, 1970) correlating permeability rates with protein size (Stokes-Einstein radius of the molecule). For instance, the calculated Stokes-Einstein radius (ae) of the 45kDa globular VEGF protein is 30.2 Å according to 0.483 × (Molecular weight in Da)0.386 as given in (Venturoli and Rippe, 2005). Our extrapolation methods yield a baseline permeability rate of 4.3 × 10−8 cm/s for VEGF (Stefanini et al., 2008). The microvascular permeability for VEGF in tumor tissue (for use in modeling breast tumor as the tissue of interest in section 5.2) can be estimated from corresponding values for similar-sized molecules, e.g., ovalbumin (45 kDa; ae = 30.8 Å) was measured to permeate at about 5.77 × 10−7 cm/s (Yuan et al., 1995). Details can be found in (Stefanini et al., 2008).
Lymphatic drainage is a major route by which interstitial proteins are transported to the blood because size-dependent transendothelial permeability restricts their intravasation into blood capillaries. In contrast, there is no macromolecular impedance in the filling of the initial lymphatics, hence the protein concentrations drained through the lymphatics are assumed to be continuous with interstitial concentrations at the lymphatic entrance. Mathematically, we describe unidirectional lymphatic drainage as:
where kL is the lymphatic drainage rate in cm3/s. Detailed derivation and parameterization can be found in (Wu et al., 2009b).
The pharmacokinetics of anti-VEGF therapy can be studied via whole-body compartmental modeling based on the detailed biochemical reactions between VEGF ligands and their receptors described above. Specifically, we can simulate and study the VEGF isoform specificity of anti-VEGF agents, ligand-agent binding configurations (e.g., whether such binding is monomeric or multimeric), agent biodistribution (if the anti-VEGF agent is confined to one or several compartments), as well as various therapeutic regimen designs (varying dosage, frequency of administration, the site of injection).
This example models bevacizumab, a humanized monoclonal antibody to VEGF, the characteristics and properties of which have been reported (150 kDa; Kd of 1.8 nM in (Presta et al., 1997); half-life = 21 days in (Gordon et al., 2001)). The diseased tissue of interest is a tumor. In the absence of the anti-VEGF agent, a total of 40 ordinary differential equations describes the compartmental system (19 ordinary differential equations (ODEs) for each tissue; 2 ODEs for the blood compartment). When the anti-VEGF agent is added and confined to the blood compartment (non-extravasating agent), 3 more equations are added to the model (blood compartment), representing the chemical interactions of the anti-VEGF with VEGF121 and VEGF165, as well as free anti-VEGF.
where KAV,blood,i is the available volume fraction of the molecule i (A=anti-VEGF, V=VEGF, VA=complex VEGF/anti-VEGF). Note that UAV = KAV U. The first term represents the concentration of anti-VEGF drug injected into the patient bloodstream. The second term represents the clearance of the anti-VEGF agent by the organs (for example, kidneys or liver). The equations related to the complex form of the anti-VEGF drug are of the form:
Figure 6 illustrates the transient dynamics resulting from the intravenous injection of the VEGF-antibody; in these simulations a tumor of 2 cm diameter is considered. In the single-dose treatment (10 mg/kg), intravenous injection leads to a rapid decrease in the free VEGF concentration (Fig. 6A). That level returns to baseline after about 3 to 4 weeks. For daily smaller doses of treatment or “metronomic” treatment (1 mg/kg for 10 days), a new lower pseudo-steady state for the plasma VEGF level emerges during the duration of treatment (Fig. 6B). Following treatment cessation (10 days), the plasma VEGF level returns to that before treatment after about 3 weeks. The metronomic injection also delays the peak of maximum formation of VEGF-antiVEGF complex as compared to a single-dose treatment (Figure 6C–D).
Equations (26–27) describe an intravenous injection of anti-VEGF antibodies that would be confined to the plasma. If the anti-VEGF agent is injected intravenously and can extravasate, terms of the form of Eqn. (24) are added to equations (26–27) and additional equations corresponding to the free and bound antibody concentrations are needed for each tissue compartment in the systems of ODEs. Interestingly, the addition of extravasation to the model drastically changes the form of the response. Transiently after injection, the free VEGF level in plasma drops drastically (data not shown). This is due to the binding of the antibody to the free VEGF present in the plasma. Following this drop there is a several-fold increase of free VEGF concentration in plasma. The apparent “rebound” effect is due to the amount of drug delivered and its extravasation. Briefly, while some antibodies bind to the free VEGF in plasma, another portion extravasates and binds to the free VEGF present in the available interstitial fluid volume in the tissues (healthy tissue and breast tumor). Although some of the formed complexes subsequently dissociate within the same tissue, significant quantities are brought into the bloodstream (via microvascular permeability and lymphatic drainage) where the complex dissociates, leading to more free VEGF in the blood compartment. Such a counter-intuitive increase in serum VEGF following intravenous administration of anti-VEGF agents has been observed in experiments (Gordon et al., 2001; Willett et al., 2005; Segerstrom et al., 2006) and our model is, to our knowledge, the first to explain this phenomenon by an intrinsic mechanism of inter-tissue transport.
This example illustrates how computational models can provide useful insights that are not easily accessible by in vitro or in vivo experiments. This model can also be extended to examine the effects of drug treatment via alternate routes of anti-VEGF administration (e.g. intramuscular injection).
In this last example of a multi-tissue compartmental model, we investigated the molecular mechanisms by which sVEGFR1, a truncated soluble variant of the endothelial cell-surface VEGFR1, inhibits VEGF signaling. The two prevailing postulated mechanisms are: direct VEGF ligand sequestration, reducing the VEGF available for VEGFR activation (Fig. 7A middle); and heterodimerization with cell-surface VEGFR monomers, rendering the receptor dimer non-functional as trans-phosphorylation of paired intracellular domains of full-length VEGFRs is necessary for activating signal transduction (Fig. 7A bottom). The model as described in detail in (Wu et al., 2009a) simulated the first mechanism, to assess the anti-angiogenic potential of sVEGFR1’s ligand trapping capacity alone.
sVEGFR1 in its monomeric (110 kDa) and dimeric (220 kDa) forms are about 2–5 times larger than VEGF; thus free sVEGFR1 and the sVEGFR1·VEGF complex have lower vascular permeability rates than VEGF, while sharing the same lymphatic drainage rates as VEGF. Full mathematical equations describing these transport properties, along with the sVEGFR1-VEGF binding and sVEGFR1-NRP1 coupling interactions, can be found in (Wu et al., 2009a). Sample equations for the concentrations of free sVEGFR1 and sVEGFR1·VEGF165 in tissue j are given here:
In this model, calf muscle tissue was chosen as the “tissue of interest” compartment (Fig. 5), in order to investigate the effects of endogenous sVEGFR1 produced from the calf muscle on local (calf compartment) and global (normal compartment) VEGF signaling complex formation. Such predictions were expected to provide insight on whether pathological upregulation in sVEGFR1 expression in the calf may contribute to the dampened VEGF response observed in ischemic calf muscles in peripheral arterial disease. Therapeutic intravascular (IV) delivery of exogenous sVEGFR1 was also simulated to assess its efficacy in lowering systemic levels of VEGF (Wu et al., 2009a).
While, intuitively, intravascular sVEGFR1·VEGF complex formation following simulated IV infusion of sVEGFR1 would be expected to lead to a sustained reduction of plasma free VEGF, in fact sample results in Fig. 7B show that the permeability rates and lymphatic drainage rates of free sVEGFR1 and sVEGFR1·VEGF are predicted to critically determine whether this reduction will take place or not. In other words, these simulations suggest that prior to the clinical translation of administering exogenous sVEGFR1 to lower systemic VEGF levels in anti-angiogenic therapy, extensive experimental research is needed to exclude the computationally predicted possibilities where sVEGFR1 delivery counter-intuitively elevates plasma free VEGF. As in section 5.2, this computational study demonstrates that the inter-tissue transport properties of proteins significantly affect their whole-body pharmacokinetic effects.
In this paper, we summarized several computational models for investigating different spatial and temporal aspects of VEGF systems biology. In section 2, we described molecular network models that simulated in vitro endothelial cell-surface interaction experiments to investigate the roles of the PlGF ligand and NRP1 co-receptor within the VEGF system. In section 3, we presented meso-scale models for investigating the effects of in vivo tissue architecture on VEGF ligand and receptor interactions, and for predicting the intramuscular response and relative therapeutic efficacy of various modalities of pro-angiogenic therapies (gene vs. cell vs. exercise therapy) for ischemic muscle diseases. In sections 4 and 5, molecular-detailed compartmental modeling was introduced as a method for efficient prediction of average molecular concentrations within tissue subcompartments (e.g., interstitial or plasma VEGF concentrations; intramuscular cell-surface density of signaling complexes) and investigation of inter-compartment transport processes (e.g., microvascular permeability and lymphatic drainage). We described several compartmental model studies predicting in vivo effects of intra-tissue trafficking and whole-body pharmacokinetics on angiogenic response to treatments using NRP1, anti-VEGF and sVEGFR1 as therapeutic targets/agents. While the VEGF system models presented in this paper were limited to the representation of two ligand isoforms (VEGF121 and VEGF165) and three receptors (VEGFR1, VEGFR2, NRP1), these model frameworks can be readily extended to include other VEGF ligand isoforms and receptors (e.g.,VEGFR3, NRP2). Furthermore, similar computational modeling techniques are applicable and have contributed to the study of other growth factor systems, including that of FGF (Forsten et al., 2000; Filion and Popel, 2004; Filion and Popel, 2005) and EGF (Wiley et al., 2003).
Research described in this paper was supported by NIH grants R01 HL079653, R33 HL0877351, and R01 CA138264.