This article appeared in a journal published by Elsevier. The attached copy is furnished to the author for internal non-commercial research and education use, including for instruction at the authors institution and sharing with colleagues. Other uses, including reproduction and distribution, or selling or licensing copies, or posting to personal, institutional or third party websites are prohibited. In most cases authors are permitted to post their version of the article (e.g. in Word or Tex form) to their personal website or institutional repository. Authors requiring further information regarding Elsevier’s archiving and manuscript policies are encouraged to visit: http://www.elsevier.com/copyright

Author's personal copy ARTICLE IN PRESS

Journal of Theoretical Biology 252 (2008) 593–607 www.elsevier.com/locate/yjtbi

Spatial distribution of VEGF isoforms and chemotactic signals in the vicinity of a tumor Alexander R. Smalla,b,, Adrian Neaguc,d, Franck Amyota, Dan Sacketta, Victor Chernomordika, Amir Gandjbakhchea a

Laboratory of Integrative and Medical Biophysics, National Institute of Child Health and Human Development, National Institutes of Health, Bethesda, MD 20892, USA b Department of Physics, California State Polytechnic University, Pomona, CA 91768, USA c Department of Physics and Astronomy, University of Missouri-Columbia, Columbia, MO 65211, USA d Department of Biophysics and Medical Informatics, Victor Babes University of Medicine and Pharmacy, 1900 Timisoara, Romania Received 18 May 2007; received in revised form 23 January 2008; accepted 7 February 2008 Available online 16 February 2008

Abstract We propose a mathematical model that describes the formation of gradients of different isoforms of vascular endothelial growth factor (VEGF). VEGF is crucial in the process of tumor-induced angiogenesis, and recent experiments strongly suggest that the molecule is most potent when bound to the extracellular matrix (ECM). Using a system of reaction–diffusion equations, we study diffusion of VEGF, binding of VEGF to the ECM, and cleavage of VEGF from the ECM by matrix metalloproteases (MMPs). We find that spontaneous gradients of matrix-bound VEGF are possible for an isoform that binds weakly to the ECM (i.e. VEGF165), but cleavage by MMPs is required to form long-range gradients of isoforms that bind rapidly to the ECM (i.e. VEGF189). We also find that gradient strengths and ranges are regulated by MMPs. Finally, we find that VEGF molecules cleaved from the ECM may be distributed in patterns that are not conducive to chemotactic migration toward a tumor, depending on the spatial distribution of MMP molecules. Our model elegantly explains a number of in vivo observations concerning the significance of different VEGF isoforms, points to VEGF165 as an especially significant therapeutic target and indicator of a tumor’s angiogenic potential, and enables predictions that are subject to testing with in vitro experiments. r 2008 Elsevier Ltd. All rights reserved. PACS: 87.18.Bb; 82.39.k; 87.15.Vv; 87.17.Jj Keywords: Tumor-induced angiogenesis; Vascular endothelial growth factor (VEGF); Matrix metalloproteases (MMPs); Reaction–diffusion equations; Chemotaxis

1. Introduction Tumor-induced angiogenesis, in which a tumor recruits nearby capillaries to feed it by growing new sprouts, is a crucial process in the development and metastasis of tumors (Folkman, 1971, 2003). When tumors become too large to satisfy their metabolic requirements with nutrients from preCorresponding author at: Department of Physics, California State Polytechnic University, Pomona, CA 91768, USA. Tel.: +1 909 869 5202. E-mail addresses: [email protected] (A.R. Small), [email protected] (A. Gandjbakhche).

0022-5193/$ - see front matter r 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.jtbi.2008.02.009

existing capillaries, they release a variety of pro-angiogenic molecules, the most important of which is known as vascular endothelial growth factor A (VEGFA, commonly shortened to VEGF). VEGF diffuses through the extracellular matrix (ECM) surrounding the tumor and stimulates the proliferation and chemotactic migration of endothelial cells from nearby capillaries (Ng et al., 2006). The new vasculature enables the tumor to take in more nutrients, continue its growth, and metastasize via the bloodstream. In an effort to both prevent metastasis and starve the tumor, angiogenesis (and VEGF in particular) has been identified as a target for cancer therapies (Ferrara et al., 2004).

Author's personal copy ARTICLE IN PRESS 594

A.R. Small et al. / Journal of Theoretical Biology 252 (2008) 593–607

However, there are open questions concerning the roles of different forms of VEGF, including questions regarding the gradients that guide chemotactic migration of endothelial cells (Fleury et al., 2006), and the resulting vascular morphology (Lee et al., 2005; Grunstein et al., 2000). Alternative splicing leads to several different VEGF isoforms, sharing a common receptor-binding portion but differing in a matrix-binding portion (Ferrara et al., 2003; Park et al., 1993; Houck et al., 1992). The most commonly studied isoforms, denoted by the number of amino acids in their sequence, are VEGF121, VEGF165, and VEGF189. VEGF121 diffuses freely without binding to the ECM, VEGF189 binds strongly to ECM fibers, and VEGF165 is an intermediate case that binds weakly, and is generally considered the most potent form of VEGF (Grunstein et al., 2000). Matrix-bound VEGF can be cleaved from the ECM by a variety of matrix metalloproteases (MMP), in particular matrix metalloprotease 9 (MMP9) (Bergers et al., 2000; Sternlicht and Werb, 2001; Deryugina and Quigley, 2006). Cleavage from the ECM releases a 110 or 113 amino acid fragment, which retains the receptorbinding portion. The secreted, matrix-bound, and cleaved VEGF can all bind to the receptor VEGFR2 (known to guide chemotactic migration of endothelial cells, Ferrara et al., 2003; Ng et al., 2006) and stimulate the production of new vasculature (Lee et al., 2005). The basic process that we examine in this work is illustrated in Fig. 1. An aggregate of hypoxic tumor cells secretes VEGF, which diffuses through the ECM. Upon encountering an ECM fiber, a VEGF molecule can bind to the fiber. Bound molecules can subsequently unbind spontaneously (since binding is a reversible process) and continue diffusing until the next binding event. Meanwhile,

Fig. 1. Schematic of the interaction of VEGF, the ECM, and MMPs: (1) A tumor secretes VEGF (open circles), which diffuses through the matrix. (2) VEGF binds to ECM fibers (lines). (3) MMPs (dark circles, produced near capillary) cleave VEGF from the ECM, releasing a cleaved VEGF molecule (triangle). (4) Endothelial cells migrate from the capillary to the tumor, in response to VEGF gradients.

MMP molecules are secreted primarily from cells adjacent to pre-existing vasculature (we review the experimental evidence for this statement below), and diffuse through the ECM while cleaving VEGF molecules from ECM fibers. Proteolytic cleavage (unlike spontaneous unbinding) releases a smaller molecule, which also diffuses through the ECM but cannot bind to ECM fibers. Finally, in response to one or more forms of VEGF (uncleaved-diffusing, matrix-bound, or cleaved-diffusing), endothelial cells in the capillary proliferate and migrate toward the tumor. The question we address is: Which of the ECM-associated (i.e. binding to or cleaved from the ECM) forms of VEGF are distributed in a manner that can guide chemotactic migration toward the tumor? (Note that we do not consider VEGF121. Although VEGF121 plays an important role in tumor-induced angiogenesis, its spatial distribution is a straightforward diffusion problem, and so a detailed theoretical treatment is less necessary. Also, there is a significant body of experiments regarding the ECM-associated isoforms and their effects, and so the ECM-associated isoforms are a theoretically and biologically significant topic on their own. In addition, although stromal cells and neutrophils can also produce significant amounts of VEGF, Fukumura et al., 1998, we focus on tumor-derived VEGF, in keeping with the key experimental studies that we seek to model in this work.) Our efforts to model this process are motivated by a number of findings regarding the potency of matrix-bound VEGF (Lee et al., 2005), the role of MMP9 in triggering the ‘‘angiogenic switch’’ (Hanahan and Folkman, 1996; Bergers et al., 2000), and the different vascular morphologies produced in the vicinity of tumors expressing different VEGF isoforms (Lee et al., 2005; Grunstein et al., 2000). While tumor vasculature is in general of an inefficient and leaky nature, the morphology does depend to a noticeable degree on the effects of different isoforms. Mouse tumors expressing the rapidly binding isoform VEGF188 (mouse VEGF has one less amino acid than human VEGF) are surrounded by a very dense disorganized vasculature. On the other hand, tumors expressing VEGF164 are surrounded by a comparatively sparser but somewhat more organized vasculature, in which new vessels proceed from the parent vessel to the tumor in a less tortuous manner (Grunstein et al., 2000). Interestingly, expression of VEGF120 by tumor cells often leads to the dilation of existing vessels rather than growth of new vessels into the ECM and tumor. Similar effects were observed when Lee et al. (2005) generated tumor cells that express VEGF113 (normally produced by proteolytic cleavage of matrix-bound VEGF) and found that VEGF113 (and presumably VEGF110) primarily stimulate the dilation of existing vessels rather than the formation of new vessels. The importance of binding to the ECM was illustrated when Lee et al. (2005) created an MMP-resistant form of VEGF that remains bound to the ECM. When tumors expressing this isoform

Author's personal copy ARTICLE IN PRESS A.R. Small et al. / Journal of Theoretical Biology 252 (2008) 593–607

were implanted in mice, the resulting vasculature was very dense and disorganized as compared with the vasculature around tumors expressing VEGF164. Interestingly, however, tumors expressing MMP-resistant VEGF experienced the most rapid growth. Taken together, these findings strongly suggest that ECM-bound VEGF plays a significant role in endothelial cell migration in the tumor microenvironment. Also, the fact that both VEGF164 and VEGF188 are at least as potent as VEGF120 or VEGF113 in stimulating cell migration suggests that the unique effects of VEGF164 cannot be entirely explained by its ability to bind to neuropilin on the cell surface (which is known to have a chemotactic effect) (Soker et al., 1998). We will therefore focus much of our attention on the spatial distribution of matrix-bound VEGF. Also, the dependence of morphology on the VEGF isoform secreted by the tumor suggests that VEGF164 is distributed with a long-range gradient (providing strong directional cues to migrating cells), while VEGF188 and the MMP-resistant form are distributed in a pattern with a locally uniform concentration profile. The presence of gradients throughout the microenvironment should cause highly directional sprout growth, while the absence of gradients should lead to disorganized sprout growth analogous to a random walk. We will show that the different rates of binding to the ECM for different VEGF isoforms do indeed lead to these different spatial distributions, suggesting a chemotactic explanation for the influence of different VEGF isoforms on vascular morphology. In regard to proteolytic cleavage of VEGF, the release of MMP9 is often associated with the onset of angiogenesis (Bergers et al., 2000), and it has been hypothesized that a key function of MMP9 is to cleave VEGF from the ECM and make it available to endothelial cells in a cleaved form (Bergers et al., 2000; Sounni et al., 2003; Deryugina and Quigley, 2006). Cleaved VEGF may very well have a mitogenic effect, but experiments with MMP-resistant VEGF by Lee et al. show that proteolytic cleavage of VEGF is not necessary for the generation of new vasculature around a xenografted tumor. Also, there is considerable experimental evidence that during the ‘‘angiogenic switch’’ proteolytic cleavage of VEGF from the ECM primarily occurs near the parent vasculature rather than near the tumor surface. Mouse models for a variety of cancers show that when MMP9 production triggers the onset of angiogenesis, the protease is released not by tumor cells but rather by cells adjacent to the parent vasculature, often associated with an inflammatory response (Bergers et al., 2000; Hiratsuka et al., 2002; Nozawa et al., 2006; Jodele et al., 2005; Huang et al., 2002; Kaplan et al., 2005; Coussens et al., 2000). Tumor cells can indeed produce MMP9, but MMP9 production by tumor cells typically occurs during the invasive stage of tumor growth (Deryugina and Quigley, 2006; Kondraganti et al., 2000; Lakka et al., 2002), after the tumor has already recruited a vascular network (Folkman, 2003).

595

Regarding the other vertebrate MMPs associated with tumor growth, Lee et al. showed that only three (MMP3, MMP7, and MMP19) efficiently release matrix-bound VEGF164 from the ECM, and two others (MMP1 and MMP16) also release VEGF but do so far less effectively. MMP1 and MMP3 seem to follow a similar pattern of release, in which they are not produced by tumor cells until the invasive stage (Deryugina and Quigley, 2006; Van Themsche et al., 2004; McCawley et al., 2004). MMP16 is membrane-bound and hence cannot diffuse through the ECM (Sounni et al., 2003), making it less relevant to the phenomena we describe here. MMP19 is apparently not ubiquitous in tumor-induced angiogenesis, and can act to limit tumor growth (Jost et al., 2006; Bister et al., 2004; Blackhall et al., 2005), so it is not clear how to model its effects. Only MMP7 is associated with early tumor growth, particularly intestinal polyps and breast tumors (Wilson et al., 1997; Hulboy et al., 2004; Ii et al., 2006). However, it is not clear that significant MMP7 production by tumor cells is ubiquitous in the early stages of angiogenesis. It therefore seems reasonable to restrict our attention to MMP9, on the grounds of its significance as well as similarities between its pattern of activity and that of other MMPs. To model these phenomena quantitatively, we develop here a system of reaction–diffusion equations. Our model has four inputs: the rate of VEGF production by the tumor (which can be highly variable, and can be regulated by drugs), the concentration distribution of MMPs (which can be regulated by various inhibitors, Sounni et al., 2003; Deryugina and Quigley, 2006), and the rates of binding to and spontaneous unbinding from the ECM (which vary among VEGF isoforms). We work in the steady state, since the characteristic time scale for diffusion in our system is less than an hour, while most observations of tumorinduced angiogenesis occur on time scales of days or weeks (Bergers et al., 2000; Lee et al., 2005; RodriguezManzaneque et al., 2001; Muthukkaruppan et al., 1982). Using our model, we show that MMPs are required to sustain a gradient of matrix-bound VEGF189 in the steady state, while spontaneous unbinding from the ECM can produce gradients of VEGF165 without proteolytic cleavage. We also find that VEGF165, which binds to the ECM more slowly than VEGF189, is distributed over a longer distance from the tumor, with a correspondingly longerrange gradient. This is consistent with observations on the different vascular morphologies produced in response to VEGF165 and VEGF189 (Grunstein et al., 2000). Finally, our calculations show that gradients of cleaved VEGF are typically either flat or point away from the tumor, depending on the level of MMP production and number of cells secreting MMPs. 2. Theory We model the concentration of four species: soluble secreted VEGF (either VEGF165 or VEGF189), matrix-bound

Author's personal copy ARTICLE IN PRESS 596

A.R. Small et al. / Journal of Theoretical Biology 252 (2008) 593–607

VEGF, cleaved VEGF, and MMPs. We denote their respective concentrations as C s , C b , C c , and C m . For simplicity, we assume that all molecules (except matrixbound VEGF) have similar diffusion coefficients D, a practice used elsewhere in the literature (Mac Gabhann and Popel, 2005). Researchers modeling related phenomena have used estimates of D ranging from 105 to 104 mm2 =s (Mac Gabhann and Popel, 2005; Anderson and Chaplain, 1998; Fleury et al., 2006). For the tumor geometry described here, this gives a characteristic time t ¼ R2 =D in the range of 10 min to 2 h for R  0:25 mm or less. Even if D is reduced by another order of magnitude, the time scale is still shorter than the typical time scale in experimental studies of blood vessel growth (days or weeks). Because VEGF expression by the tumor precedes MMP9 activity (Hiratsuka et al., 2002; Bergers et al., 2000) and continues during at least the early stages of blood vessel growth (while the tumor is hypoxic), the short time scale for VEGF diffusion justifies an assumption that the VEGF concentration profile quickly reaches a steady state in the early stages of the angiogenic process. We will therefore restrict our attention to the steady state. Moreover, even if the level of VEGF production by the tumor changes over time, as long as it changes on a time scale longer than t a quasi-steady approximation remains valid. Note that we assume here a tumor radius of 0.25 mm, and that the distance from the tumor surface to the nearest capillary is also 0.25 mm. This radius is an upper bound, chosen to place limits on the relevant time scales. In experiments with xenotransplanted tumor cells, this may be a very realistic assumption. In other in vivo experiments, the distance from an aggregate of hypoxic cancer cells to the nearest capillary is probably closer to 0.1 mm (Carmeliet and Jain, 2000), but we can rescale all distances, times, and rates appropriately (multiply all length scales by a factor a, all time scales by a2 , and all rates by a2 ) and get the same qualitative results. 2.1. Governing equations The time-dependence of C s is controlled by three processes: diffusion through the ECM, binding to the ECM, and spontaneous unbinding from the ECM. The rate of binding to the ECM is proportional to a rate constant kon , the concentration of soluble molecules C s , and the area available for binding to the ECM. The available area must decrease as C b increases, since the ECM offers a large but finite number of sites for VEGF to bind to. The finite number of binding sites corresponds to a maximum value of C b that we will call C sat . To incorporate this effect, we assume that the second term is proportional to C sat  C b . The saturation value C sat can presumably vary considerably, depending on matrix composition (which in turn depends on the particular tissue, and the levels of fibroblast and MMP activity, both of which can vary considerably during tumor growth).

Meanwhile, the rate of spontaneous unbinding is proportional to the concentration of bound VEGF molecules and a rate constant koff . We will discuss below the relationship between kon and koff for different VEGF isoforms. (Note that for simplicity we assume the tumor secretes only a single VEGF isoform. This condition has been realized in several important experimental studies, including Grunstein et al., 2000 and Lee et al., 2005.) These effects lead to the following form for the time derivative of C s : q C s ¼ Dr2 C s  kon C s ðC sat  C b Þ þ koff C b qt

(1)

The concentration of bound VEGF is also governed by three processes: binding of soluble VEGF to the ECM, spontaneous unbinding of VEGF from the ECM, and proteolytic cleavage by MMP molecules. The rates of binding and unbinding are just the negative of the second and third terms in Eq. (1). The rate of cleavage by MMPs should be proportional to C b and the MMP concentration C m multiplied by a rate constant kc . We get the equation q C b ¼ kon C s ðC sat  C b Þ  koff C b  kc C m C b qt

(2)

If we set the time derivative in Eq. (2) to zero, and solve for C b , we get the steady-state result Cb ¼

kon C s C sat kon C s þ koff þ kc C m

(3)

A steady-state equation for C s can be obtained by setting the time derivative in Eq. (1) to zero and using Eq. (3) Dr2 C s ¼

kon C s kc C m C sat kon C s þ koff þ kc C m

(4)

Solving Eq. (4) is our principal task. We impose two boundary conditions. The first is that C s goes to zero at large distances from the tumor. While the background VEGF concentration is typically non-zero in vivo, we are only interested in phenomena close to the tumor, where the VEGF concentration will be significantly higher than the background. We can therefore choose any value we like for the VEGF concentration at sufficiently distant boundary, as long as it is consistent with a concentration profile that decreases over short distances from the tumor. For convenience, we will assume that C s ¼ 0 at a distant boundary 10 tumor radii from the center of the computational window. (Note that absorption by the distant boundary is not the only loss mechanism governing C s , as binding to the ECM and proteolytic conversion to a cleaved form is also included in our model.) The second boundary condition is that there is a steady flux of VEGF from the surface of the tumor. The flux of VEGF per unit area and unit time is given by DrC s in accordance with Fick’s Law. Often the more easily measured parameter, however, is the total rate of VEGF

Author's personal copy ARTICLE IN PRESS A.R. Small et al. / Journal of Theoretical Biology 252 (2008) 593–607

Fig. 2. Geometry of our model: size and position of tumor and cell producing MMP, as well as our cylindrical coordinate system. The system is azimuthally symmetric about the z-axis.

expression by the tumor (flux integrated over the surface of the tumor), which we will denote E. We will assume a spherical tumor of radius R ¼ 0:25 mm (Fig. 2). The boundary condition can then be expressed as E Dn^  rC s jsurface ¼ (5) 4pR2 where n^ is a unit vector normal to the tumor surface. As discussed above, the geometry and length scales depicted in Fig. 2 are probably most directly relevant to experiments with xenotransplanted tumor cells. However, even inside a tumor containing capillaries there will typically be spatial separation of the hypoxic tumor cells and the capillaries from which new sprouts grow. The typical separation in that case will be of order 0.1–0.2 mm, the diffusion distance of oxygen in vivo (Carmeliet and Jain, 2000). 2.2. Relationship between kon and koff In the dilute limit (i.e. C s ; C b 5C sat ), and in the absence of MMPs, the thermodynamic equilibrium concentrations of soluble and bound VEGF are dictated solely by energetic considerations rather than proteolytic cleavage or competition for binding sites. In this limit, Eq. (3) reduces to C b kon C sat ¼ Cs koff

(6)

If the system is in thermodynamic equilibrium we can apply the Boltzmann relation to C b =C s and obtain a relationship between kon and koff : kon ¼

koff C b koff ¼ expðDG=kB TÞ C sat C s C sat

(7)

where kB T is the temperature multiplied by Boltzmann’s constant and DG is the free energy change when a molecule of VEGF unbinds from the ECM. Because binding is energetically favorable, DG40, and so koff okon C sat . Also, as the binding energy increases, kon should increase (since the probability of sticking to the ECM increases) and koff should decrease (since molecules remain stuck to the ECM for longer periods of time).

597

It is difficult to be sure what the most relevant measurements are for kon , koff , C sat , and DG, as they depend on factors such as pH (Goerges and Nugent, 2004) and matrix composition, both of which can vary among tumors, and during the different stages of tumor growth. However, it is generally reported that the bound and soluble forms of VEGF165 coexist under physiologically relevant conditions, while VEGF189 exists primarily in the matrix-bound form under physiologically relevant conditions (Ferrara et al., 2003; Ng et al., 2006). In addition, the robustness of binding to the ECM suggests that DG must be at least several times kB T. In calculations for VEGF165 we will therefore vary koff between 0:01kon C sat and 0:1kon C sat . On the other hand, in calculations for VEGF189 we will assume that koff is effectively zero, since the overwhelming prevalence of bound VEGF189 relative to the soluble form indicates that koff 5kon C sat for VEGF189. 2.3. Gradients of matrix-bound VEGF To obtain the spatial distribution of ECM-bound VEGF, we must solve Eq. (4) for C s and then find C b from the steady-state condition (Eq. (3)). However, we can obtain a significant result without solving those equations. If we take the gradient of Eq. (3), we get that rC b ¼ kon C sat

ðkoff þ kc C m ÞrC s  C s kc rC m ðkon C s þ koff þ kc C m Þ2

(8)

This equation has three implications. The first is that a non-negligible steady-state gradient of matrix-bound VEGF requires either significant proteolytic cleavage or non-negligible unbinding (i.e. non-negligible value for koff ). This implies that a gradient of VEGF189 cannot be sustained without the presence of MMP molecules, but the significant rate of spontaneous unbinding of VEGF165 does enable the formation of gradients. This is consistent with the finding that vascular networks around tumors expressing VEGF165 tend to be somewhat less disorganized than vascular networks around tumors expressing VEGF189, which suggests that VEGF165 may produce stronger chemotactic signals (Grunstein et al., 2000; Lee et al., 2005). The second implication is that when the gradient of the MMP concentration points opposite to the gradient of the soluble VEGF concentration, the gradient of matrix-bound VEGF will be stronger. This suggests a possible further role for MMPs in the initiation of the angiogenic switch, by regulating gradients of matrix-bound VEGF. Finally, at high levels of MMP secretion, the magnitude of the gradient will be inversely proportional to the rate at which MMP molecules are produced. Although MMP molecules are needed to sustain a gradient, when the rate of cleavage by MMP molecules is fast relative to binding (kc C m Xkon ) the effect of adding more MMP molecules is to deplete the ECM-bound VEGF, reducing the magnitude of the gradient.

Author's personal copy ARTICLE IN PRESS 598

A.R. Small et al. / Journal of Theoretical Biology 252 (2008) 593–607

2.4. MMP concentration profile and model geometry In addition to specifying the tumor geometry, we must also specify the distribution of MMP molecules. As discussed above, the MMPs relevant to the onset of tumor-induced angiogenesis are typically secreted by discrete cells located adjacent to pre-existing vasculature. If MMP sources are approximately uniformly distributed between the tumor surface and the pre-existing vessels, then the MMP concentration profile will be peaked near the outermost MMP sources (pre-existing vessels). To see this, consider Poisson’s equation (steady-state diffusion) with a source term that is isotropic in a spherical shell, where the outer radius of the shell is the characteristic distance from the center of the tumor to a blood vessel. The solution inside the shell will have a constant term and a quadratic term. It is thus reasonable to assume that proteolytic cleavage will be most significant in the vicinity of MMP-secreting cells near the parent capillaries. Rather than specifying the number of MMP-secreting cells and the position of each cell, when interested in the spatial location of proteolytic cleavage we will restrict our attention to phenomena in the vicinity of a single MMPsecreting cell (typically located approximately one tumor radius away from the tumor surface). Focusing our attention on the vicinity of a single MMP-secreting cell is motivated by the fact that proteolytic cleavage is localized near the parent vasculature, which is the same place where cell migration commences (and hence the most significant region of interest when studying the ‘‘angiogenic switch’’). Also, the considerations outlined above suggest that the effects of other MMP sources can be averaged and treated as a background added to the local source. Our simplification yields the geometry shown in Fig. 2. The associated MMP concentration profile is given by 0 1 f B C C m ¼ C m;0 @1  f þ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiA 2 ðz  2RÞ þ r2 þ r2cell

(9)

The parameter C m;0 is related to the overall rate of MMP secretion in our system. The parameter f, restricted to the range 0pfp1, determines the background level of MMP production relative to MMP production by a single cell. The case f ¼ 1 corresponds to a single MMP-secreting cell (and zero background, or a sparse distribution of MMPsecreting cells), while the case f ¼ 0 corresponds to a uniform MMP concentration profile (densely distributed MMP-secreting cells). Finally, rcell , which is the radius of an endothelial cell ( 5 mm, or 0.02 in our units), is inserted to avoid singularities in our numerical calculations. Although this parameter is somewhat artificial, it reflects the fact that on short length scales the MMP concentration does not diverge, but rather approaches a finite value at the surface of a cell secreting MMP molecules. The value of rcell does not qualitatively affect results when considered on

a coarse-grained scale larger than rcell . We are primarily interested in vascular network morphology, and hence VEGF concentration gradients on length scales comparable to the distances that new blood vessel sprouts must traverse to reach the tumor. The use of Eq. (9) involves an explicit decision to neglect MMP7 production by tumors, as well as a decision to neglect endogenous tissue inhibitors of metalloproteases (TIMPs). These decisions were made primarily for simplicity. The neglect of MMP7 production by tumor cells is justified by the fact, discussed above, that most studies implicate MMP9 (released by cells adjacent to the parent vasculature) in the angiogenic switch, and only a handful of studies (mostly in particular organs) find a role for MMP7 production in the early stages of tumor growth. In the discussion of our results we will comment on how certain findings would change in the presence of MMP7 production by tumors during the angiogenic switch. We will also discuss the effects of neglecting endogenous TIMPs, and argue that TIMP activity would, in many cases, support the phenomena that we describe. Also, a uniform MMP concentration profile amounts to dealing with the heterogeneity of the tumor microenvironment via averaging. The tumor microenvironment in vivo is highly heterogeneous, with a variety of MMP sources, an inhomogeneous ECM, MMP inhibitors (Deryugina and Quigley, 2006; Sounni et al., 2003; Sternlicht and Werb, 2001), and secretion of variable levels of MMPs by tumor cells. In that complex environment, the most meaningful experimental work would involve average measurements of [VEGF] as a function of distance from the tumor. Finally, it is important to solve this problem in 3 dimensions rather than 2. 2D mathematical models are often used to study tumor-induced angiogenesis (Anderson and Chaplain, 1998; Levine et al., 2001; Sun et al., 2005), but such models are inappropriate if two oppositely oriented gradients may fully or partially cancel out (i.e. Eq. (8)). Spatial dimensionality has a significant effect on gradients around an isotropic source of a diffusing species, with the gradient decreasing as one over the distance from the origin in 2D and one over the square of the distance from the origin in 3D. It is crucial to use a 3D model (or, in our case, a 2D model with azimuthal symmetry) to accurately calculate gradients. 2.5. Concentration profile of cleaved VEGF After solving Eq. (4) for C s and using the solution to obtain C b from Eq. (3), we can calculate the concentration profile of cleaved VEGF molecules by solving the steadystate diffusion equation with a source term Dr2 C c ¼ kc C m C b

(10)

The concentration C c a position x can be obtained by convolving the source term with Green’s function 1=jx  x0 j. When we do the integration in cylindrical coordinates, and integrate over the azimuthal angle to

Author's personal copy ARTICLE IN PRESS A.R. Small et al. / Journal of Theoretical Biology 252 (2008) 593–607

account for symmetry, we get Z 1Z 1 kc C m ðr0 ; z0 ÞC b ðr0 ; z0 Þ C c ðr; zÞ / 1

599

For the boundary condition (Eq. (5)), we introduce a dimensionless counterpart to the rate of VEGF expression E: Q  Et=ðC sat R3 Þ. This gives the boundary condition

0

r0  qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðz  z0 Þ2 þ ðr þ r0 Þ2 þ r2cell 0

n^  rC s jsurface ¼ 1

0

4rr B C K @qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiA dr0 dz0 2 2 ðr þ r0 Þ þ ðz  z0 Þ þ r2cell

(11)

where KðmÞ is the complete elliptic integral of the first kind with parameter m (frequently called the square of the modulus) (Abramowitz et al., 1964), the factor r2cell is again included in the denominator to avoid singularities in a physically reasonable manner, and the proportionality symbol / is used because we have omitted factors of 4p, D, etc. 2.6. Non-dimensionalization

Q 4p

(15)

where the gradient is taken with respect to the dimensionless coordinates. Finally, non-dimensionalization of Eq. (11) is straightforward, since scaling factors such as 4p, D, etc. were neglected: Z 1Z 1 V c ðr; z0 Þ / Mðr0 ; z00 ÞV b ðr0 ; z00 Þ 1

0

r0  qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðz0  z00 Þ2 þ ðr þ r0 Þ2 þ ðrcell =RÞ2 0

1

0

4r  r B C K @qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiA dr0 dz00 2 2 2 0 0 00 ðr þ r Þ þ ðz  z Þ þ ðrcell =RÞ (16)

Our work is easier if we can eliminate some of the parameters from our model, and cast our equations in dimensionless form (Murray, 2003). To work in dimensionless time and length units, we first introduce the dimensionless variables r  r=R, ðx0 ; y0 ; z0 Þ  ðx=R; y= R; z=RÞ, and t0  t=t ¼ tD=R2 . All lengths will hereafter be given in units of the tumor radius R, and all times and rates will hereafter be given in terms of the characteristic diffusion time t ¼ D=R2 . We will also use dimensionless VEGF concentrations V s  C s =C sat , V b  C b =C sat , and V s  C c =C sat . All VEGF concentrations will hereafter be given relative to the ECM saturation concentration C sat . For the rate constants kon and koff we introduce the dimensionless counterparts k1  kon C sat t and k2  koff t. For the MMP concentration, because C m is always premultiplied by the MMP cleavage rate constant kc here, we will use the dimensionless MMP activity rate M  kc C m t. This gives the following dimensionless form for the MMP concentration profile: 0 1 f B C M ¼ kc C m ¼ M 0 @1  f þ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiA 2 2 ðz0  2Þ þ r2 þ ðrcell =RÞ (12) where M 0 ¼ kc C m;0 t. The governing equation (4) becomes r2 V s ¼

k1 V s M k1 V s þ k2 þ M

(13)

where the gradient is taken with respect to the dimensionless coordinates. Eq. (3) for the bound VEGF concentration becomes Vb ¼

k1 V s k1 V s þ k2 þ M

(14)

The dimensionless quantities introduced here are summarized in the Appendix. 2.7. Parameter values After non-dimensionalization our model has five adjustable parameters: k1 , k2 , Q, M 0 , and f. The parameters Q, M 0 , and f, reflecting the levels and patterns of VEGF and MMP production, can vary considerably when comparing different tumors, and different stages of tumor growth. The parameters k1 and k2 depend on the VEGF isoform under study, and also on matrix composition and pH (Goerges and Nugent, 2004). However, we can estimate physiologically relevant and qualitatively interesting ranges for these parameters. We can estimate Q, the flux of VEGF from the tumor, using the generally accepted observation that most of the VEGF present in vivo is matrix-bound (i.e. V s o1 at the surface of the tumor). A spherical source of diffusing molecules releases molecules at a rate proportional to the surface concentration multiplied by 4pRD (Berg, 1993). In dimensionless units, we get the result that Qo1=4p. We thus considered values of Q between 102 and 1. In regard to k1 , we are aware of experiments indicating that the half-life of soluble VEGF must be no longer than a few hours (Park et al., 1993). This imposes a lower bound of k1 4101 . For an upper bound, examining Eq. (13) shows that the length scale over which V s decays will be on 1=2 the order of k1 (if M4k1 V s ) or M 1=2 (if Mok1 V s ). In order to obtain a significant gradient of matrix-bound VEGF near a capillary, that length scale must be larger than  13 (in units of tumor radius R), so k1 or M must be less than or equal to 10. Our upper bound for k1 is presumably valid for VEGF165, which binds to the ECM only weakly, but may not hold for VEGF189.

Author's personal copy ARTICLE IN PRESS 600

A.R. Small et al. / Journal of Theoretical Biology 252 (2008) 593–607

We can see from Eq. (14) that in order to reduce V b significantly below the saturation level (without bringing it to nearly zero) and achieve a gradient, k2 þ M (and hence M 0 , k2 , or both) and k1 V s must be of the same order. We therefore typically varied the ratio M 0 =ðk1 Q=4pÞ between 102 and 102 . Also, as discussed above, we will use the parameter range k1 4k2 40:1k1 when studying VEGF165, and k2  0 for VEGF189. Finally, in regard to the parameter f, we will use f ¼ 0 (spatially uniform MMP concentration profile) for convenience when studying the effect of varying k1 (rate of binding to the ECM). To study the distribution of cleaved VEGF in the vicinity of an MMP source, we will use nonzero f, and frequently consider the value f ¼ 0:2. This value is consistent with the images of sparse MMPsecreting cells presented by Bergers et al. (2000). 3. Methods

4. Numerical results Our numerical investigation focuses on two questions. The first is how the rate of binding to the ECM affects the distribution of matrix-bound VEGF. This question is motivated by the observation that different VEGF isoforms give rise to different vascular morphologies, suggesting different directional cues for endothelial cell migration. The second question that we ask is how the cleaved VEGF is distributed in the vicinity of a cell secreting MMP molecules. It has been suggested that the role of MMP molecules in tumor-induced angiogenesis is to cleave matrix-bound VEGF from the ECM, making it available to endothelial cells (Deryugina and Quigley, 2006; Sounni et al., 2003). We ask whether the cleaved VEGF is distributed in a manner that can provide directional cues to guide endothelial cell migration. 4.1. Varying VEGF release rate We begin by considering the effect of the VEGF release rate Q, which sets our boundary condition at the tumor surface. We show in Fig. 4 plots of V s versus distance from the center of the tumor for k1 ¼ 1, k2 ¼ 0:1, M ¼ 102 ( k1 Q=4p for Q ¼ 101 ), and Q between 102 and 100 . We see that V s is roughly proportional to the VEGF release rate Q in the range of values under consideration, in 10−1 10−2 soluble VEGF conc.

For the case of a uniform (or average) MMP concentration profile (f ¼ 0) we solved Eq. (13) in spherical coordinates using Matlab 7:0:4. We used the function ‘‘bvp4c’’, which is a solver for non-linear boundary value problems. At r ¼ 10 we imposed the boundary condition V s ¼ 0. We used a grid of approximately 100 points in each calculation. For the case of a non-uniform MMP concentration profile (fa0) we used Femlab Version 2.3. We used the 2D axisymmetric mode (cylindrical coordinates) of the Chemical Engineering module, due to the azimuthal symmetry of the 3D problem depicted in Fig. 2. We used the stationary non-linear solver for Eq. (13). Our mesh, depicted in Fig. 3, has 2987 triangular elements, the smallest of which are  20 mm on a side in the vicinity of the MMP-secreting cell. (Note that this is larger than a typical endothelial cell, but we are calculating concentrations and gradients over the length scales of vascular networks.) Due to azimuthal symmetry we only need to solve the problem for rX0. Along the outer circular boundary in Fig. 3 we imposed the boundary condition V s ¼ 0, and at

r ¼ 0 we imposed no-flux boundary conditions, reflecting the symmetry of our problem. After solving Eq. (13) in Femlab, the solution for V s was exported to Matlab and interpolated onto a square mesh with the function ‘‘postinterp’’ for more convenient postprocessing. During post-processing we calculated V b from Eq. (14) and V c from Eq. (16).

10−3 10−4 10−5

Q = 100 Q = 10−1 Q = 10−2

10−6 1

2

3

4

5

6

7

8

9

10

r

Fig. 3. Mesh used for finite element simulations. Note that the mesh is finer in the vicinity of the MMP source.

Fig. 4. V s versus r for different values of Q. k1 ¼ 1, k2 ¼ 0:1, and M ¼ 102 . In all cases, the MMP concentration profile is spatially uniform. r ¼ 1 corresponds to the surface of the tumor, and V s is scaled to the saturation concentration of the ECM (C sat ).

Author's personal copy ARTICLE IN PRESS A.R. Small et al. / Journal of Theoretical Biology 252 (2008) 593–607

are shown in Fig. 5. Although larger values of k1 can give much steeper gradients at some particular distance r (if M is low enough), the gradients have a consistently longer range for smaller values of k1 . The steeper gradient at a particular distance corresponds to a V b profile that is mostly flat except for a steep transition region. This is consistent with a scenario where isoforms that bind rapidly to the ECM diffuse only until they encounter a vacant binding site. Increasing M has two effects. The first is that it attenuates the gradient of V b at long distances, irrespective of the value of k1 . The second effect is that at higher values of M there is no steep transition where the gradient spikes. Both effects are due to the action of MMP molecules cleaving VEGF from the ECM, depleting the supply of matrix-bound VEGF. However, increasing M does not change our qualitative finding that only slowly binding VEGF isoforms are distributed through the ECM in a manner that gives rise to long-range concentration gradients of matrix-bound VEGF. We also show concentration profiles of soluble VEGF in Fig. 6 for the same values of k1 , k2 , and M used in Fig. 5. Similar results hold, with weaker binding producing stronger gradients over longer distances, irrespective of

spite of the non-linearity of Eq. (13). We are primarily interested in V b , which is a function of the product k1 V s , so for our purposes varying k1 and varying Q will have the same effect on V b . We will therefore work with Q ¼ 101 in the remainder of this paper. 4.2. Comparing VEGF165 with VEGF189 We wish to compare the distribution of matrix-bound VEGF165 with matrix-bound VEGF189 in the presence of a tumor. These isoforms are distinguished by their affinities for binding to the ECM, and hence the values of k1 and k2 , the rate constants for binding to and unbinding from the ECM. As discussed above, the specific numbers will depend on factors such as ECM composition and pH, but commonly accepted experimental observations imply that the value of k1 for VEGF165 must be significantly smaller than the value for VEGF189, that k2 must be essentially zero for VEGF189, and that k2 must be non-negligible in comparison with k1 for VEGF165. We examined the spatial profile of V b and its gradient, varying k1 between 1 and 100. We considered a range of values for M, ranging from 103 ( 0:1k1 Q=4p for k1 ¼ 1) to 101 ( 10k1 Q=4p for k1 ¼ 1). Representative results

M = 10−3

1

0.7

k1 = 100, k2 = 0.01

0.6

k1 = 100, k2 = 0.1

k1 = 102, k2 = 0

0.4

k1 = 101, k2 = 0.01

bound VEGF conc.

bound VEGF conc.

0.45

k1 = 102, k2 = 0

0.8

0.5 0.4 0.3

k1 = 101, k2 = 0.01

0.35

k1 = 100, k2 = 0.01

0.3

k1 = 100, k2 = 0.1

0.25 0.2 0.15

0.2

0.1

0.1

0.05 0

0 1

1.5

2

2.5

3 r

3.5

4

4.5

1

5

M = 10−3

101

k1 = 101, k2 = 0.01 k1 = 100, k2 = 0.01

100

k1 = 100, k2 = 0.1

10−1 10−2 10−3 10−4 1

1.5

2

2.5

3 r

3.5

4

4.5

1.5

2

2.5

5

3 r

3.5

4

4.5

5

3.5

4

4.5

5

M = 10−1

101

k1 = 102, k2 = 0

gradient of bound VEGF conc.

gradient of bound VEGF conc.

M = 10−1

0.5

0.9

601

k1 = 102, k2 = 0

100

k1 = 101, k2 = 0.01 k1 = 100, k2 = 0.01

10−1

k1 = 100, k2 = 0.1

10−2 10−3 10−4

1

1.5

2

2.5

3 r

Fig. 5. (a,c) V b for various values of k1 , k2 , M. (b,d) gradient of V b for various values of k1 , M. In all cases the MMP concentration profile is spatially uniform. r ¼ 1 corresponds to the surface of the tumor.

Author's personal copy ARTICLE IN PRESS A.R. Small et al. / Journal of Theoretical Biology 252 (2008) 593–607

602

7

x 10−3

M = 10

−3

6

5

2

k1 = 10 , k2 = 0 1

k1 = 10 , k2 = 0.01

4

k1 = 100, k2 = 0.01 0

3

k1 = 10 , k2 = 0.1

2

2

4

k1 = 10 , k2 = 0 1 k1 = 10 , k2 = 0.01

3

k1 = 10 , k2 = 0.01

0 0

k1 = 10 , k2 = 0.1

2 1

1 0

0 1

1.5

2

2.5

3 r

3.5

4

4.5

5

1

−3

10−2

M = 10

10−3

2

k1 = 10 , k2 = 0 1 k1 = 10 , k2 = 0.01

10−4

0

k1 = 10 , k2 = 0.01 0 k1 = 10 , k2 = 0.1

1

1.5

2

2.5

3 r

1.5

2

2.5

3.5

4

4.5

5

3 r

M = 10

10−2 gradient of soluble VEGF conc.

gradient of soluble VEGF conc.

−1

M = 10

5 soluble VEGF conc.

soluble VEGF conc.

6

x 10−3

3.5

4

4.5

5

−1 2

k1 = 10 , k2 = 0 1

k1 = 10 , k2 = 0.01 k1 = 100, k2 = 0.01

10−3

0

k1 = 10 , k2 = 0.1

10−4

1

1.5

2

2.5

3 r

3.5

4

4.5

5

Fig. 6. (a,c) V s for various values of k1 , k2 , M. (b,d) gradient of V s for various values of k1 , M. In all cases the MMP concentration profile is spatially uniform. r ¼ 1 corresponds to the surface of the tumor.

the value of M. Importantly, decreasing k2 can attenuate V s for all r (by decreasing the amount of VEGF released from the ECM) but does not significantly change the length scale over which V s decays. Taken together, these results suggest that VEGF165 can guide chemotactic migration of endothelial cells toward tumors, but VEGF189 cannot. VEGF189 is mostly located near the tumor, so that the concentration gradient is appreciable in only a narrow region. On the other hand, molecules of VEGF165 bind to the ECM more slowly, and can travel longer distances before binding. This gives rise to directional cues that can guide endothelial cell migration. A simple check of our theory, which would not require in vivo measurements of VEGF concentration profiles, would be to measure the rates of binding to and cleavage from a gel containing realistic concentrations of VEGF, MMPs, and heparan sulfate proteoglycans. Our explanation for the different vascular morphologies induced by VEGF165 and VEGF189 (short-range versus long-range gradients of matrix-bound VEGF) requires that for VEGF165, k1 is of order 1 in our dimensionless units. To translate this result to a form that can be compared with experiment, the lifetime ðkon C sat Þ1 for a VEGF165 molecule in an

unsaturated ECM must be of order t1 ¼ D=R2 , if our explanation for the different effects of VEGF165 and VEGF189 is correct. Our model also assumes measurable production of cleaved VEGF, implying that MX0:01k1 . Physically, this implies that the rate at which MMPs cleave VEGF from the ECM must be at least  1% of the rate at which VEGF binds to the ECM. Also, the results obtained make it easy to qualitatively understand what would happen if we generalized our mathematical model to the case of a tumor secreting multiple isoforms simultaneously: The ECM will be saturated with matrix-bound VEGF, primarily matrixbound VEGF189, out to some distance rsat . Although VEGF165 can also bind to the ECM in this region, the concentration profile will be dominated by VEGF189, as a VEGF189 molecule encountering a binding site is far more likely to bind than a VEGF165 molecule. This distance can be estimated by setting the rate of VEGF189 production by the tumor (Q189 ) equal to the rate at which a given volume of matrix-bound VEGF189 is proteolytically cleaved. In this saturated region, the concentration profile of matrix-bound VEGF will be flat, and hence the gradient will be negligible. Calculations to show these mathematical predictions explicitly are underway.

Author's personal copy ARTICLE IN PRESS A.R. Small et al. / Journal of Theoretical Biology 252 (2008) 593–607

Beyond rsat , the concentration profile will be dominated by VEGF165, as the supply of available VEGF189 has been exhausted. Since concentration profiles of matrix-bound VEGF165 tend to decay slowly, there will be a large region over which the concentration profile attenuates, and hence a large region in which the concentration gradient is nonnegligible. It is in this region, where the matrix-bound VEGF is primarily the VEGF165 isoform, that migrating endothelial cells will encounter the strongest chemotactic signals. A more detailed and quantitative study of these phenomena, and implications for developmental processes, is underway.

4.3. Cleaved VEGF We also calculated the spatial distribution of soluble, matrix-bound, and cleaved VEGF molecules in the vicinity of a cell secreting MMP molecules. We held Q constant at 0:1 as before, and fixed k2 at a value 0:01 (in keeping with our finding that varying k2 does not have a significant qualitative effect on the spatial profile of V s . We varied k1 between 1 and 100. When we vary the parameters M 0 (proportional to the rate of MMP secretion), and f (from Eq. (12)) we find three distinct qualitative regimes of behavior.

603

When MMP-secreting cells are sparse and proteolytic cleavage is tightly localized (f  1), the gradient of cleaved VEGF can point toward or away from the tumor, depending on the rate of cleavage by MMPs relative to the rate at which VEGF binds to the ECM. For low values of the ratio M 0 =ðk1 Q=4pÞ, cleavage of VEGF from the ECM is localized near the cell secreting MMP molecules. This leads to a concentration profile of cleaved VEGF that is strongest in the vicinity of the MMP-secreting cell. On the other hand, for larger values of the ratio M 0 =ðk1 Q=4pÞ, cleavage of VEGF from the ECM is rapid everywhere in the vicinity of the tumor. The strongest concentration of cleaved VEGF is then found near the surface of the tumor, where VEGF is released. We illustrate these findings in Fig. 7 for k1 ¼ 1 and 10. The value of M 0 =ðk1 Q=4pÞ at which the maximum concentration of cleaved VEGF shifts from the MMPsecreting cell to the tumor surface depends on k1 . Increasing k1 shifts the concentration profile of soluble (secreted) VEGF (and hence matrix-bound VEGF) closer to the tumor. This, in turn, localizes binding of VEGF to the ECM (and hence cleavage of VEGF from the ECM) near the tumor. Still, regardless of the value of k1 , the same trend is observed: If MMP-secreting cells are sparsely distributed, then at sufficiently low levels of MMP secretion, the cleaved species VEGF is localized near an

Fig. 7. Concentration profile of cleaved VEGF in the vicinity of the tumor (white semicircle) and a single MMP-secreting cell (black circle). All concentrations are in arbitrary units. (a) k1 ¼ 1, M 0 ¼ 3  104 , (b) k1 ¼ 1, M 0 ¼ 5  103 , (c) k1 ¼ 10, M 0 ¼ 103 , (d) k1 ¼ 10, M 0 ¼ 102 .

Author's personal copy ARTICLE IN PRESS 604

A.R. Small et al. / Journal of Theoretical Biology 252 (2008) 593–607

MMP-secreting cell. Increasing the level of MMP secretion tends to localize cleaved VEGF near the surface of the tumor. It is possible that the localization of cleavage near a tumor surface may give rise to vascular morphologies that differ according to whether the tip of the sprout is near or far from a tumor surface. However, to properly treat this it will be necessary to develop a model that combines our VEGF transport calculations with cell migration, to fully explore issues of vascular morphology. In the case where the local concentration of MMP is not dominated by a single MMP-producing cell (fo1 in Eq. (12)), the concentration gradients of cleaved VEGF no longer point toward the MMP-producing cell. Instead, as we show in Fig. 8, cleaved VEGF is concentrated near the surface of the tumor (Fig. 8(a)). A closer view (Fig. 8(b)) shows that the concentration profile is not quite uniform around the surface of the tumor. Near an MMP-secreting cell the concentration profile is somewhat flatter than elsewhere. This is consistent with our result in Eq. (8), where oppositely oriented gradients of soluble VEGF and MMP cancel out. Interestingly, we found that the concentration profile of cleaved VEGF in general decays quite slowly compared with the concentration profiles of soluble (secreted) and matrix-bound VEGF. One way to compare the nearly flat gradient of cleaved VEGF with the gradients of bound VEGF and the diffusing VEGF species secreted by the tumor is to compute the ratio of VEGF concentration at the tumor surface to VEGF concentration at the MMP-secreting cell. We do this in Table 1 for different levels of MMP secretion, i.e. varying M 0 in Eq. (12). (For all these calculations, k1 ¼ 1, k2 ¼ 0:01, and f ¼ 0:2.) At every level of MMP secretion, the concentration profile of cleaved VEGF is nearly flat between the tumor and the MMP-secreting cell, increasing only very slowly (but monotonically) from the MMP source to the tumor surface. The results are similar if we rescale k1 and M 0 proportionally. Overall, our calculations show that gradients of cleaved VEGF are significantly weaker (by at least an order of magnitude)

than gradients of soluble (secreted) and matrix-bound VEGF. Taken together, our results suggest that gradients of cleaved VEGF cannot account for chemotactic migration of endothelial cells toward tumors. The gradients are at best quite weak, and frequently point in the wrong direction for chemotactic migration toward the tumor. Also, even if endothelial cells do not distinguish cleaved VEGF from other soluble forms of VEGF (e.g. secreted VEGF165 or VEGF189), and respond only to the total concentration profile of soluble VEGF (irrespective of isoform), the addition of a flat or nearly flat cleaved VEGF concentration profile will contribute little to the total gradient of soluble VEGF, and will diminish the relative gradient of soluble VEGF. We therefore conclude that cleaved VEGF, acting either alone or in concert with other soluble VEGF isoforms, is not distributed in a manner that is conducive to endothelial cell migration toward tumors. 4.4. Location of MMP-secreting cells We also varied the location of the cell secreting MMP molecules. Varying the location of the MMP-secreting cell

Table 1 Ratio of VEGF concentration at tumor surface to VEGF concentration at location of MMP-secreting cell M0

105 104 103 102 101 100

Concentration ratio Soluble (165 or 189)

Bound

Cleaved

2.25 2.30 2.64 3.97 5.16 5.40

1.74 1.80 2.45 7.60 16.92 19.39

1.06 1.06 1.08 1.17 1.25 1.26

Calculations are performed for the case k1 ¼ 1 and f ¼ 0:2.

Fig. 8. Concentration profile of cleaved VEGF in a system with a tumor (white circle), an MMP-secreting cell (black circle), and a background concentration of MMP (f ¼ 0:2 in Eq. (12)). All concentrations are in arbitrary units. Calculations performed for k1 ¼ 1, M 0 ¼ 104 . (a) Concentration profile around tumor. (b) Close-up near MMP-secreting cell.

Author's personal copy ARTICLE IN PRESS A.R. Small et al. / Journal of Theoretical Biology 252 (2008) 593–607

does not appreciably alter the MMP concentration profile in the vicinity of the tumor surface, as the MMP concentration decays only as the inverse power of the distance. Small displacements of the MMP-secreting cell also fail to change the qualitative behavior described above: At low MMP concentrations (slow cleavage rates), production of cleaved VEGF is localized at the cell secreting MMP molecules, and at high MMP concentrations (fast cleavage), production of cleaved VEGF is localized at the tumor surface. However, moving the MMP-secreting cell sufficiently far from the tumor surface actually localizes production of cleaved VEGF at the tumor surface. Far away from the tumor (a few tumor radii for the typical parameters considered in this paper) the concentration of ECM-bound VEGF is low, and so no significant production of cleaved VEGF can occur there. The only significant production of cleaved VEGF occurs near the tumor surface. These results show that the most significant MMP sources, for the purpose of studying the distribution of cleaved VEGF, are those located near but not necessarily at the tumor surface.

5. Discussion Our results demonstrate that the spatial distribution of VEGF molecules is highly dependent on the rate at which they bind to the ECM, the rate at which they are cleaved from the ECM by MMPs, and the spatial distribution of cells secreting MMPs. VEGF isoforms that bind to the ECM rapidly (e.g. VEGF189) tend to cluster near the tumor, in both their soluble and matrix-bound forms. VEGF isoforms that bind to the ECM more slowly (e.g. VEGF165) are distributed over greater distances, with significant gradients over correspondingly longer distances. Molecules of cleaved VEGF, produced by MMPs, can either cluster near the cells secreting MMPs or spread out with a nearly uniform concentration profile in the vicinity of the tumor, depending on the MMP concentration profile. Together, our calculations suggest that among the ECMassociated isoforms of VEGF, only VEGF165 is consistently distributed in a manner that is conducive to chemotactic migration of endothelial cells over long distances. Our calculations also suggest that the spatial distribution of cleaved VEGF cannot account for chemotactic migration of endothelial cells toward tumors. This is consistent with literature reports concerning the different effects of the various VEGF isoforms (Grunstein et al., 2000; Lee et al., 2005). Our calculations do not address VEGF121, but as discussed above, the ECM binding forms of VEGF give rise to distinct vascular morphologies, and are at least as potent as VEGF121 in stimulating endothelial cell migration. The logical conclusion is that VEGF165 is more effective than any other ECM-associated VEGF isoform at guiding chemotactic migration toward tumors, due in large part to its spatial distribution.

605

Full experimental tests of our calculations and conclusions would require direct measurements of the concentration profiles of the different forms of VEGF in vivo, which would be very technically challenging. However, measurements of the rate at which different VEGF isoforms bind to the ECM, and the rate of cleavage by MMPs, could be done with tissue samples. Such data would help determine whether the parameter values used and the qualitative behaviors described are indeed relevant to real biological systems. We have not considered the effects of TIMPs (Sternlicht and Werb, 2001) in our model. Qualitatively, since TIMPs inhibit the activity of MMPs, the main effect of TIMPs would be to localize MMP activity. This would result in highly localized production of cleaved VEGF even in cases where the cells producing MMPs are densely distributed. More localized production of cleaved VEGF would imply concentration profiles that guide endothelial cells toward MMP sources rather than the tumor surface. This analysis suggests that the neglect of TIMP activity does not significantly change the applicability of our simulations. We are planning to analyze in more detail the implications of TIMPs for our model. Another obvious extension of this work would be to model the effects of MMP production by the tumor during the initial stages of angiogenesis, as there is evidence that some tumors produce MMP7 during the early stages of their growth. From Eq. (8), it is clear that parallel gradients of MMPs and soluble VEGF would lead to weaker gradients of bound VEGF (and hence flatter concentration profiles). Also, localizing proteolytic cleavage to the tumor surface would lead to greater production of cleaved VEGF production near the tumor surface, and hence cleaved VEGF gradients that can indeed guide endothelial cell migration toward the tumor. A more complete quantitative study of these issues is being planned, to contrast the effects of MMP7 release from the tumor with the effects of MMP9 release near the vasculature. A different extension of this work would be to use the calculated VEGF concentration profiles in simulations of endothelial cell migration that include both random and directed motility (Anderson and Chaplain, 1998; Levine et al., 2001; Sun et al., 2005). It would be especially interesting to see if such a simulation, with gradients of chemoattractants as predicted by our model, could reproduce the different vascular morphologies produced in response to VEGF165 and VEGF189. It may also be that phenomena in the vicinity of a cell secreting MMPs can explain the formation of loops in tumor-induced angiogenesis. Going beyond tumor-induced angiogenesis, extensions of our model may be of use in tissue engineering and wound healing. The incorporation of VEGF and proteoglycans in collagen matrices is an area of current interest (Yao et al., 2006). Accurate models of VEGF transport and binding could help in the design of surgical products. Finally, our model of VEGF gradient formation may be useful in attempts to model vascular growth in the

Author's personal copy ARTICLE IN PRESS 606

A.R. Small et al. / Journal of Theoretical Biology 252 (2008) 593–607

embryo, as it is known that expression of different VEGF isoforms during embryonic development can give rise to very different (even pathological) vascular morphologies (Ruhrberg et al., 2002). In conclusion, we have a model for the interaction of VEGF molecules with the ECM and MMPs. Our model predicts very different spatial distributions for the various forms of VEGF. We hypothesize that a number of experimental observations concerning tumor-induced angiogenesis can be explained in part by the way that different forms of VEGF are distributed to guide endothelial cell migration. We also hypothesize, based on our model, that the most important matrix-binding VEGF isoform for chemotactic migration of endothelial cells is VEGF165. The assumptions underlying our model can be tested with straightforward in vitro experiments. Our model may also be useful in studying angiogenesis in wound healing, tissue engineering, and vascular development in the embryo. Acknowledgments This work was supported in part by the intramural program of the National Institute of Child Health and Human Development. Adrian Neagu acknowledges the support of National Science Foundation grant NSF-FIBR0526854. We dedicate this paper to the memory of science educator Don Herbert, who in his role as Mr. Wizard encouraged numerous young people to study science, including the first author of this paper. Appendix A. Definitions of dimensionless quantities In the list below we summarize the key dimensionless quantities introduced in Section 2.6. Note that these definitions involve the quantities C sat (saturation concentration of matrix-bound VEGF), R (tumor radius), and t (diffusion time scale  R2 =D). (1) r ¼ r=R: Distance from z-axis in cylindrical coordinates, normalized to tumor radius R. (2) V s ¼ C s =C sat : Concentration of VEGF in solution, normalized to C sat , the concentration of matrix-bound VEGF in a saturated ECM. (3) V b ¼ C b =C sat : Concentration of VEGF bound to the ECM, normalized to the saturation concentration. (4) V c ¼ C c =C sat : Concentration of cleaved VEGF in solution, normalized to the concentration of VEGF in a saturated ECM. (5) M ¼ kc C m t: Local rate of MMP activity, normalized to the characteristic diffusion time t. When M ¼ 1, the time constant for proteolytic cleavage of VEGF from the ECM is 1=t. (6) M 0 ¼ kc C m;0 t: Overall rate of MMP activity, normalized to the characteristic diffusion time t. (7) k1 ¼ kon C sat t: Rate constant for VEGF to bind to the ECM, normalized to the saturation concentration C sat

and characteristic diffusion time t. When k1 ¼ 1 and the ECM is empty of VEGF (V b ¼ 1), the lifetime of a VEGF molecule (uncleaved) in solution is t. (8) k2 ¼ koff t: Rate at which VEGF un-binds from the ECM, normalized to the characteristic diffusion time t. When k2 ¼ 1, the lifetime of an ECM-bound VEGF molecule is k2 . (9) Q ¼ Et=ðR3 C sat Þ: Rate of VEGF production by the tumor, normalized to the characteristic diffusion time t and the amount of VEGF present when a volume R3 is saturated to a concentration C sat . When Q ¼ 1 the tumor produces enough VEGF to saturate a volume R3 to a concentration C sat in a time t. It is worth noting that the relationship between the rate constants kon and koff given in Eq. (7) leads to the following relationship for their dimensionless counterparts: kon C sat k1 ¼ ¼ expðDG=kB TÞ. koff k2

(A.1)

Because DG40, it follows that k1 4k2 , just as kon 4koff . References Abramowitz, M., Stegun, I.A., U. S. N. B. of Standards, 1964. Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables. Applied Mathematics Series, vol. 55. National Bureau of Standards, Washington, DC. Anderson, A.R., Chaplain, M.A., 1998. Continuous and discrete mathematical models of tumor-induced angiogenesis. Bull. Math. Biol. 60 (5), 857–899, 0092-8240 Journal Article. Berg, H.C., 1993. Random Walks in Biology, expanded edition. Princeton University Press, Princeton. Bergers, G., Brekken, R., McMahon, G., Vu, T.H., Itoh, T., Tamaki, K., Tanzawa, K., Thorpe, P., Itohara, S., Werb, Z., Hanahan, D., 2000. Matrix metalloproteinase-9 triggers the angiogenic switch during carcinogenesis. Nat. Cell Biol. 2 (10), 737–744. Bister, V.O., Salmela, M.T., Karjalainen-Lindsberg, M.L., Uria, J., Lohi, J., Puolakkainen, P., Lopez-Otin, C., Saarialho-Kere, U., 2004. Differential expression of three matrix metalloproteinases, MMP-19, MMP-26, and MMP-28, in normal and inflamed intestine and colon cancer. Dig. Dis. Sci. 49 (4), 653–661. Blackhall, F., Pintilie, M., Liu, N., Pennington, C., Zhu, C., Shepherd, F., Edwards, D., Tsao, M., 2005. Expression profiling of matrix metalloproteinases (MMPs) and inhibitors reveals a novel prognostic role for MMP-19 in patients with non-small cell lung cancer. Lung Cancer 49, S43–S44. Carmeliet, P., Jain, R.K., 2000. Angiogenesis in cancer and other diseases. Nature 407 (6801), 249–257. Coussens, L.M., Tinkle, C.L., Hanahan, D., Werb, Z., 2000. MMP-9 supplied by bone marrow-derived cells contributes to skin carcinogenesis. Cell 103 (3), 481–490. Deryugina, E.I., Quigley, J.P., 2006. Matrix metalloproteinases and tumor metastasis. Cancer Metastasis Rev. 25 (1), 9–34. Ferrara, N., Gerber, H.P., LeCouter, J., 2003. The biology of VEGF and its receptors. Nat. Med. 9 (6), 669–676. Ferrara, N., Hillan, K.J., Gerber, H.P., Novotny, W., 2004. Discovery and development of bevacizumab, an anti-VEGF antibody for treating cancer. Nat. Rev. Drug Discovery 3 (5), 391–400. Fleury, M.E., Boardman, K.C., Swartz, M.A., 2006. Autologous morphogen gradients by subtle interstitial flow and matrix interactions. Biophys. J. 91 (1), 113–121.

Author's personal copy ARTICLE IN PRESS A.R. Small et al. / Journal of Theoretical Biology 252 (2008) 593–607 Folkman, J., 1971. Tumor angiogenesis: therapeutic implications. N. Engl. J. Med. 285 (21), 1182–1186. Folkman, J., 2003. Fundamental concepts of the angiogenic process. Curr. Mol. Med. 3 (7), 643–651. Fukumura, D., Xavier, R., Sugiura, T., Chen, Y., Park, E.C., Lu, N., Selig, M., Nielsen, G., Taksir, T., Jain, R.K., Seed, B., 1998. Tumor induction of VEGF promoter activity in stromal cells. Cell 94 (6), 715–725. Goerges, A.L., Nugent, M.A., 2004. pH regulates vascular endothelial growth factor binding to fibronectin—a mechanism for control of extracellular matrix storage and release. J. Biol. Chem. 279 (3), 2307–2315. Grunstein, J., Masbad, J.J., Hickey, R., Giordano, F., Johnson, R.S., 2000. Isoforms of vascular endothelial growth factor act in a coordinate fashion to recruit and expand tumor vasculature. Mol. Cell. Biol. 20 (19), 7282–7291. Hanahan, D., Folkman, J., 1996. Patterns and emerging mechanisms of the angiogenic switch during tumorigenesis. Cell 86 (3), 353–364. Hiratsuka, S., Nakamura, K., Iwai, S., Murakami, M., Itoh, T., Kijima, H., Shipley, J.M., Senior, R.M., Shibuya, M., 2002. MMP9 induction by vascular endothelial growth factor receptor-1 is involved in lung-specific metastasis. Cancer Cell 2 (4), 289–300. Houck, K.A., Leung, D.W., Rowland, A.M., Winer, J., Ferrara, N., 1992. Dual regulation of vascular endothelial growth-factor bioavailability by genetic and proteolytic mechanisms. J. Biol. Chem. 267 (36), 26031–26037. Huang, S.Y., Van Arsdall, M., Tedjarati, S., McCarty, M., Wu, W.J., Langley, R., Fidler, I.J., 2002. Contributions of stromal metalloproteinase-9 to angiogenesis and growth of human ovarian carcinoma in mice. J. Natl. Cancer Inst. 94 (15), 1134–1142. Hulboy, D.L., Gautam, S., Fingleton, B., Matrisian, L.M., 2004. The influence of matrix metalloproteinase-7 on early mammary tumorigenesis in the multiple intestinal neoplasia mouse. Oncol. Rep. 12 (1), 13–17. Ii, M., Yamamoto, H., Adachi, Y., Maruyama, Y., Shinomura, Y., 2006. Role of matrix metalloproteinase-7 (matrilysin) in human cancer invasion, apoptosis, growth, and angiogenesis. Exp. Biol. Med. 231 (1), 20–27. Jodele, S., Chantrain, C.F., Blavier, L., Lutzko, C., Crooks, G.M., Shimada, H., Coussens, L.M., DeClerck, Y.A., 2005. The contribution of bone marrow-derived cells to the tumor vasculature in neuroblastoma is matrix metalloproteinase-9 dependent. Cancer Res. 65 (8), 3200–3208. Jost, M., Folgueras, A.R., Frerart, F., Pendas, A.M., Blacher, S., Houard, X., Berndt, S., Munaut, C., Cataldo, D., Alvarez, J., Melen-Lamalle, L., Foidart, J.M., Lopez-Otin, C., Noel, A., 2006. Earlier onset of tumoral angiogenesis in matrix metalloproteinase-19deficient mice. Cancer Res. 66 (10), 5234–5241. Kaplan, R.N., Riba, R.D., Zacharoulis, S., Bramley, A.H., Vincent, L., Costa, C., MacDonald, D.D., Jin, D.K., Shido, K., Kerns, S.A., Zhu, Z.P., Hicklin, D., Wu, Y., Port, J.L., Altorki, N., Port, E.R., Ruggero, D., Shmelkov, S.V., Jensen, K.K., Rafii, S., Lyden, D., 2005. VEGFR1-positive haematopoietic bone marrow progenitors initiate the pre-metastatic niche. Nature 438 (7069), 820–827. Kondraganti, S., Mohanam, S., Chintala, S.K., Kin, Y., Jasti, S.L., Nirmala, C., Lakka, S.S., Adachi, Y., Kyritsis, A.P., Ali-Osman, F., Sawaya, R., Fuller, G.N., Rao, J.S., 2000. Selective suppression of matrix metalloproteinase-9 in human glioblastoma cells by antisense gene transfer impairs glioblastoma cell invasion. Cancer Res. 60 (24), 6851–6855. Lakka, S.S., Rajan, M., Gondi, C., Yanamandra, N., Chandrasekar, N., Jasti, S.L., Adachi, Y., Siddique, K., Gujrati, M., Olivero, W., Dinh, D.H., Kouraklis, G., Kyritsis, A.P., Rao, J.S., 2002. Adenovirus-

607

mediated expression of antisense MMP-9 in glioma cells inhibits tumor growth and invasion. Oncogene 21 (52), 8011–8019. Lee, S., Jilani, S.M., Nikolova, G.V., Carpizo, D., Iruela-Arispe, M.L., 2005. Processing of VEGF-a by matrix metalloproteinases regulates bioavailability and vascular patterning in tumors. J. Cell Biol. 169 (4), 681–691. Levine, H.A., Pamuk, S., Sleeman, B.D., Nilsen-Hamilton, M., 2001. Mathematical modeling of capillary formation and development in tumor angiogenesis: penetration into the stroma. Bull. Math. Biol. 63 (5), 801–863. Mac Gabhann, F., Popel, A.S., 2005. Differential binding of VEGF isoforms to VEGF receptor 2 in the presence of neuropilin-1: a computational model. Am. J. Physiol.—Heart Circ. Physiol. 288 (6), H2851–H2860. McCawley, L.J., Crawford, H.C., King, L.E., Mudgett, J., Matrisian, L.M., 2004. A protective role for matrix metalloproteinase-3 in squamous cell carcinoma. Cancer Res. 64 (19), 6965–6972. Murray, J.D., 2003. Mathematical Biology. Springer, Berlin. Muthukkaruppan, V.R., Kubai, L., Auerbach, R., 1982. Tumor-induced neovascularization in the mouse eye. J. Natl. Cancer Inst. 69 (3), 699–708. Ng, Y.S., Krilleke, D., Shima, D.T., 2006. VEGF function in vascular pathogenesis. Exp. Cell Res. 312 (5), 527–537. Nozawa, H., Chiu, C., Hanahan, D., 2006. Infiltrating neutrophils mediate the initial angiogenic switch in a mouse model of multistage carcinogenesis. Proc. Natl. Acad. Sci. USA 103 (33), 12493–12498. Park, J.E., Keller, G.A., Ferrara, N., 1993. Vascular endothelial growthfactor (VEGF) isoforms—differential deposition into the subepithelial extracellular-matrix and bioactivity of extracellular matrix-bound VEGF. Mol. Biol. Cell 4 (12), 1317–1326. Rodriguez-Manzaneque, J.C., Lane, T.F., Ortega, M.A., Hynes, R.O., Lawler, J., Iruela-Arispe, M.L., 2001. Thrombospondin-1 suppresses spontaneous tumor growth and inhibits activation of matrix metalloproteinase-9 and mobilization of vascular endothelial growth factor. Proc. Natl. Acad. Sci. USA 98 (22), 12485–12490. Ruhrberg, C., Gerhardt, H., Golding, M., Watson, R., Ioannidou, S., Fujisawa, H., Betsholtz, C., Shima, D.T., 2002. Spatially restricted patterning cues provided by heparin-binding VEGF-a control blood vessel branching morphogenesis. Genes Dev. 16 (20), 2684–2698. Soker, S., Takashima, S., Miao, H.Q., Neufeld, G., Klagsbrun, M., 1998. Neuropilin-1 is expressed by endothelial and tumor cells as an isoformspecific receptor for vascular endothelial growth factor. Cell 92 (6), 735–745. Sounni, N.E., Janssen, M., Foidart, J.M., Noel, A., 2003. Membrane type1 matrix metalloproteinase and TIMP-2 in tumor angiogenesis. Matrix Biol. 22 (1), 55–61. Sternlicht, M.D., Werb, Z., 2001. How matrix metalloproteinases regulate cell behavior. Annu. Rev. Cell Dev. Biol. 17, 463–516. Sun, S.Y., Wheeler, M.F., Obeyesekere, M., Patrick, C., 2005. Multiscale angiogenesis modeling using mixed finite element methods. Multiscale Modeling Simulation 4 (4), 1137–1167. Van Themsche, C., Potworowski, E.F., St-Pierre, Y., 2004. Stromelysin-1 (MMP-3) is inducible in T lymphoma cells and accelerates the growth of lymphoid tumors in vivo. Biochem. Biophys. Res. Commun. 315 (4), 884–891. Wilson, C., Heppner, K., Labosky, P., Hogan, B., Matrisian, L., 1997. Intestinal tumorigenesis is suppressed in mice lacking themetalloproteinase-matrilysin 10.1073/pnas.94.4.1402. Proc. Natl. Acad. Sci. 94 (4), 1402–1407. Yao, C., Roderfeld, M., Rath, T., Roeb, E., Bernhagen, J., Steffens, G., 2006. The impact of proteinase-induced matrix degradation on the release of VEGF from heparinized collagen matrices. Biomaterials 27 (8), 1608–1616.