A Hybrid Model of the Role of VEGF Binding in Endothelial Cell Migration and Capillary Formation

Vascular endothelial growth factor (VEGF) is the most studied family of soluble, secreted mediators of endothelial cell migration, survival, and proliferation. VEGF exerts its function by binding to specific tyrosine kinase receptors on the cell surface and transducing the effect through downstream signaling. In order to study the influence of VEGF binding on endothelial cell motion, we develop a hybrid model of VEGF-induced angiogenesis, based on the theory of reinforced random walks. The model includes the chemotactic response of endothelial cells to angiogenic factors bound to cell-surface receptors, rather than approximating this as a function of extracellular chemical concentrations. This allows us to capture biologically observed phenomena such as activation and polarization of endothelial cells in response to VEGF gradients across their lengths, as opposed to extracellular gradients throughout the tissue. We also propose a novel and more biologically reasonable functional form for the chemotactic sensitivity of endothelial cells, which is also governed by activated cell-surface receptors. This model is able to predict the threshold level of VEGF required to activate a cell to move in a directed fashion as well as an optimal VEGF concentration for motion. Model validation is achieved by comparison of simulation results directly with experimental data.

1. Introduction

Motility – random, directed, and collective – is a fundamental property of cells. Coordinated cellular motion leads to all physiological tissue patterns, a consequence of integration across multiple temporal and spatial scales. However, when this integration is aberrant, the properties that emerge lead to a critical bifurcation point in cancer progression: angiogenesis. Angiogenesis, the formation of new blood vessels from pre-existing ones, provides the necessary blood supply for the growth and nourishment of solid tumors beyond a few millimeters in diameter (Hanahan and Weinberg, 2000; Augustin et al., 2009). Tumor angiogenesis is associated with an extremely complex, yet well-ordered series of events at the center of which is the enhanced replicative potential and motility of endothelial cells (ECs) that line the inner surface of blood vessels (Folkman, 1985; Hanahan and Weinberg, 2000).

The multistep process associated with successful angiogenesis can be summarized as EC degradation of the adjacent basement membrane, migration (sprouting), proliferation, alignment, tube formation, branching that increases near the tumor leading to a brush-border, anastomosis (fusion of vessels), synthesis of new basement membrane, recruitment of parenchymal cells, network remodeling, and a return to quiescence (Folkman, 1985; Yancopoulos et al., 2000; Conway et al., 2001; Augustin et al., 2009). Precise coordination and integration of molecular, cellular, and tissue level interactions is required for angiogenesis to be successful from initiation to stabilization of a functional vascular plexus.

Under conditions of hypoxia, tumor cells induce angiogenesis by releasing a wide variety of polypeptide angiogenic factors that stimulate EC activation, survival, proliferation, migration, and maturation. Members of the vascular endothelial growth factor (VEGF) family have been identified as the predominant amongst these angiogenic factors that regulate EC phenotype (Yancopoulos et al., 2000; Jain, 2002; Ferrara, 2004; Hicklin and Ellis, 2005). VEGF has been implicated across a range of human cancer and preclinical studies have shown that VEGF stimulates survival of existing vessels, promotes new vessel growth, and contributes to vascular abnormalities such as tortuousness and hyperpermeability. The angiogenic effects of the VEGF pathway are primarily initiated through the interaction of VEGFA and its natural, endothelial cell specific receptor, VEGFR2, which is up-regulated during angiogenesis (Neufeld et al., 1999; McMahon, 2000; Conway et al., 2001). Dimerization and activation of VEGFR2 results in mitogenic, chemotactic, and prosurvival signals (Nor et al., 1999; Ferrara et al., 2003; Ferrara, 2004), which help to determine endothelial cell phenotype.

The various steps of the angiogenic cascade require endothelial cells to take on spatio-temporally varying phenotypes; that is, at any given time and at any specific spatial location within a developing sprout, ECs can have a proliferative, migratory, or quiescent phenotype. For example, tip cells are highly migratory and lead the extending sprout through the extracellular matrix (ECM), whereas stalk cells, which form the vessel lumen and recruit support cells, can be either proliferative or quiescent. It has been shown that endothelial cells compete for the tip cell position through relative levels of VEGF-receptors (Jakobsson et al., 2010). While much is known about the sequential morphogenetic processes required for angiogenesis and the growth factors that drive it, far less is known about how cellular and molecular mechanisms are coordinated to control cell motility decisions and phenotype choices. In order to advance the understanding and manipulating of the processes that occur during angiogenesis, it is critical to understand how individual cells interpret the biochemical signals that come from their unique microenvironment.

For decades, mathematical models have been employed to help address some of the pressing questions associated with tumor angiogenesis. As discussed in detail in Jackson and Zheng (2010), Zheng et al. (2013), existing models of tumor-induced angiogenesis can be characterized as continuous approaches (Balding and McElwain, 1985; Byrne and Chaplain, 1995, 1996; Anderson and Chaplain, 1998a,b; Holmes and Sleeman, 2000; Levine et al., 2001; Arakelyan et al., 2002; Sleeman and Wallis, 2002; Manoussaki, 2003; Plank and Sleeman, 2003, 2004; Plank et al., 2004; Levine and Nilsen-Hamilton, 2006; Schugart et al., 2008; Billy et al., 2009; Xue et al., 2009; Travasso et al., 2011), wherein cells are assumed to have a continuous distribution; discrete or hybrid models (Stokes and Lauffenburger, 1991; Anderson and Chaplain, 1998b; Tong and Yuan, 2001; Plank and Sleeman, 2003, 2004; Sun et al., 2005; Bartha and Rieger, 2006; Gevertz and Torquato, 2006; Frieboes et al., 2007; Milde et al., 2008; Capasso and Morale, 2009; Owen et al., 2009; Perfahl et al., 2011), wherein cells are modeled as individual agents and diffusible chemicals are modeled as a continuum; and cell-based formulations (Peirce et al., 2004; Bauer et al., 2007; Bentley et al., 2009; Qutub and Popel, 2009; Wcislo et al., 2009; Jackson and Zheng, 2010; Liu et al., 2011) wherein explicit incorporation of different properties of individual cells allows collective behavior of cell clusters to be predicted from the behavior and interactions of individual cells. Reviews of these models that appeared in or before 2009 can be found in Mantzaris et al. (2004), Peirce (2008), Qutub et al. (2009). However, these models suffer from the following limitations. Continuum descriptions of biological motion such as chemotaxis are derived by averaging quantities such as cellular and vascular densities, and therefore apply to the macroscopic behavior of a large number of cells. However, the initial stages of new capillary development requires only a small number of cells in a highly discrete arrangement, which is better described by treating cells as individual agents. Further, even when a hybrid or cell-based approach has been adopted, endothelial cell movement, and/or microvessel formation speed and direction is typically assumed to depend on extracellular chemokine concentrations, whereas it is known that cells integrate the chemical signal via receptors on their surfaces in order to make behavioral decisions (Nor et al., 1999; Ferrara et al., 2003; Ferrara, 2004).

In order to study the influence of VEGF binding on EC motion, we develop here a hybrid model of VEGF-induced angiogenesis that is based on the theory of reinforced random walks. We will include in our model, the chemotactic response of endothelial cells to angiogenic factors bound to cell-surface receptors, rather than approximating this as a function of extracellular chemical concentrations. This will allow us to capture biologically observed phenomena such as the activation and polarization of endothelial cells in response to VEGF gradients across their lengths. We will also propose a novel and more biologically reasonable functional form for the chemotactic sensitivity of cells, which is also governed by activated cell-surface receptors.

In the sections that follow, we will first develop a model to describe the motion of a single EC that has a tip cell phenotype. The motion of an EC across a 1 mm² domain will be simulated, and average time taken by the cell to reach the tumor (VEGF source) computed as a function of source strength. We will then extend this model to study the early stages of blood vessel formation and sprout extension. In particular, we will capture the following growth pattern of developing sprouts, typically observed in experiments. Sprouts arising from parent vessels are observed to grow parallel to each other initially (Paweletz and Knierim, 1989), with anastomoses between tip cells and stalk cells, or between two tip cells observed to occur a certain distance into the stroma. As the developing vessels near the source of chemoattractants, new sprouts emerge in a process called sprout branching, which increases with proximity to the source. This has been described as a “brush-border” effect (Muthukkaruppan et al., 1982; Sholley et al., 1984). We will also investigate the effect on vascular development of the source strength of VEGF. Model validation will follow from a qualitative comparison of simulations with experimental data on neovascularization of the rat cornea taken from Sholley et al. (1984).

2. Materials and Methods

2.1. Model Development: Single Cell Motion under the Influence of VEGF

A vital characteristic of all cells is their ability to sense their environment and respond to it, such as motion toward or away from an external, chemical stimulus. The response of endothelial cells to a chemokine like VEGF involves the following two major steps that a mathematical description of this process needs to account for: (i) detection of the signal (via gradient of bound VEGF to cell-surface receptors) and (ii) transduction of the external signal into an internal signal that controls the pattern of movement (Mantzaris et al., 2004). The theory of reinforced random walks, where a master equation governing cell movement is derived directly from the governing biology, as opposed to discretizing a continuous equation of macroscopic motion, provides a natural framework for modeling the movement of individual endothelial cells that initiate vascular sprout development. We remark that for simplicity, the effects of the extracellular matrix on EC motion are not explicitly considered at this time so that the cell is assumed to migrate on a homogeneous and isotropic medium. Further, on the time-scale of interest, cell proliferation, and death are assumed to be negligible.

We begin our model development by first simulating a single EC moving under the influence of VEGF. The EC is interpreted as a (biased) random walker that adapts its motility decisions under the influence of activated VEGF-receptors on its surface. We consider a 2-dimensional spatial domain, with a tumor located at x = 1 serving as the source of VEGF, and a parent vessel located at x = 0 providing individual ECs to begin sprout development, as shown in Figure 1A. The tumor secretes VEGF under the condition of hypoxia, which diffuses toward the parent vessel. VEGF is taken up by cells lining the parent vessel, transforming them into sprout tip cells. These migrate up its chemical gradient, pulling behind them the developing capillaries. The principal dynamics that we wish to capture with the model are the binding and uptake of VEGF by the sprout tip cells, the subsequent activation of cell-surface receptors, and the chemotactic response of the cells to this stimulation. To our knowledge, this level of molecular detail has not been implemented previously in a model of tumor-induced angiogenesis.

FIGURE 1

Figure 1. (A) Geometry of the model domain. A source of VEGF (e.g., a tumor) is situated at x = 1. VEGF diffuses toward a parent vessel located at x = 0, and is taken up by endothelial cells lining it. The activated cells migrate up gradients of VEGF, elongating behind them capillaries. (B–D) Motion of a cell on a 2-d lattice in response to VEGF stimulus. The cell begins at position (n, m) in (B). VEGF molecules are shown in red. The number of activated VEGF-receptors is greatest on the cell surface at lattice site (n + 1/2, m), and the probability of motion in this direction is the greatest. Consequently, the cell is likely to move from site (n, m) to (n + 2, m), as shown in (C,D).

The cell is located initially at spatial position x = 0, y = 0.5, and will move in response to a local, cellular gradient of VEGF, which has its source at x = 1. A schematic of this process is shown in Figures 1B–D. Following Plank et al. (Plank and Sleeman, 2003; Plank et al., 2003, 2004), we base our spatial discretization on purely biological considerations. As per the approach developed in Plank and Sleeman (2003), Plank et al. (2003), Plank et al. (2004), Othmer and Stevens (1997), the following master equation is used to describe a biased random walk (in two dimensions) of the endothelial cell, moving under the influence of VEGF in its local environment:

$\begin{array}{ll}\frac{\partial p_{n,m}}{\partial t}=& {\widehat{\mathcal{T}}}_{n-1,m}^{H+}\left(W\right)p_{n-1,m}+{\widehat{\mathcal{T}}}_{n+1,m}^{H-}\left(W\right)p_{n+1,m}\\ & +{\widehat{\mathcal{T}}}_{n,m-1}^{V+}\left(W\right)p_{n,m-1}+{\widehat{\mathcal{T}}}_{n,m+1}^{V-}\left(W\right)p_{n,m+1}\\ & -\left({\widehat{\mathcal{T}}}_{n,m}^{H+}\left(W\right)+{\widehat{\mathcal{T}}}_{n,m}^{H-}\left(W\right)+{\widehat{\mathcal{T}}}_{n,m}^{V+}\left(W\right)\right.\\ & \left.+{\widehat{\mathcal{T}}}_{n,m}^{V-}\left(W\right)\right)p_{n,m}.\end{array}$ (1)

Here, p_n,m(t) describes the probability that a cell is at site (n, m), at time t. ${\widehat{\mathcal{T}}}_{n,m}^{H\pm }\left(\cdot \right),$ and ${\widehat{\mathcal{T}}}_{n,m}^{V\pm }\left(\cdot \right)$ are the transition probabilities per unit time for a one step horizontal jump to (n ± 1, m), or a one step vertical jump to (n, m ± 1) respectively. The vector W gives the concentration of the chemoattractant C, at the lattice sites. In order for the master equation to translate to the standard diffusion-chemotaxis equation for cell movement in the continuum limit, it is assumed that the dependence of transition rates at lattice site (n, m) is localized to chemoattractant concentration at sites (n ± 1/2, m) and (n, m ± 1/2). This is reasonable, since we may think of a cell present at lattice site (n, m), with its boundaries extending to half the lattice length. The cell can therefore sense the chemical concentrations at these half-lattice sites, and make a decision where to move, as illustrated in Figures 1B–D. Under these assumptions, W = (…, C_−n−1/2,m, C_−n,m, C_−n+1/2,m, C_−n+1,m, …).

The mean waiting time at the (n, m)th site is given by $1∕\left({\widehat{\mathcal{T}}}_{n,m}^{H+}\left(W\right)+{\widehat{\mathcal{T}}}_{n,m}^{H-}\left(W\right)+{\widehat{\mathcal{T}}}_{n,m}^{V+}\left(W\right)+{\widehat{\mathcal{T}}}_{n,m}^{V-}\left(W\right)\right).$ We make the assumption that the decision of where to move in space is independent of the decision when to move in time. Mathematically, this is equivalent to setting

${\widehat{\mathcal{T}}}_{n,m}^{H+}\left(W\right)+{\widehat{\mathcal{T}}}_{n,m}^{H-}\left(W\right)+{\widehat{\mathcal{T}}}_{n,m}^{V+}\left(W\right)+{\widehat{\mathcal{T}}}_{n,m}^{V-}\left(W\right)=k.$ (2)

That is to say, the cell makes a decision to move (or to stay still) after a constant amount of time, k. One way to satisfy these assumptions is the following choice of transition probabilities, as made by Othmer and Stevens (1997)):

$\begin{array}{c}{\mathcal{T}}_{n,m}^{H\pm }=\frac{1}{k}\frac{\tau \left(C_{n\pm 1∕2,m}\right)}{\begin{array}{c}\tau \left(C_{n+1∕2,m}\right)+\tau \left(C_{n-1∕2,m}\right)\hfill \\ +\tau \left(C_{n,m+1∕2}\right)+\tau \left(C_{n,m-1∕2}\right)\hfill \end{array}},\\ {\mathcal{T}}_{n,m}^{V\pm }=\frac{1}{k}\frac{\tau \left(C_{n,m\pm 1∕2}\right)}{\begin{array}{c}\tau \left(C_{n+1∕2,m}\right)+\tau \left(C_{n-1∕2,m}\right)\hfill \\ +\tau \left(C_{n,m+1∕2}\right)+\tau \left(C_{n,m-1∕2}\right)\hfill \end{array}},\end{array}$ (3)

for some function τ(C) of the chemoattractant. The choice of the functional form for τ(C) is based on the particular form of the chemotactic sensitivity desired, and is explained in the following section. A grid of mesh size h is chosen, thereby fixing x = nh. Passing to the continuum limit h → 0 and 1/4k → ∞ such that h²/4k = D_p, where D_p is the diffusion coefficient of ECs, Othmer and Stevens (1997) show that the master equation (1) translates to the familiar diffusion-chemotaxis equation (4), for cell motion.

$\frac{\partial p}{\partial t}=D_p\Delta p-\nabla \cdot \left(p\chi \left(C\right)\nabla C\right),$ (4)

where the chemotactic sensitivity χ(C) = D_p(ln τ(C))’. To get a completely discretized model of the motion of the cell, the time derivative of p in equation (1) is approximated by a simple forward difference scheme, with k as the time step, given by

$k=\frac{h^2}{4D_p}.$ (5)

A diagrammatic representation of the motion of the cell is shown in Figure 1. The cell starts out at time t at the lattice site (n, m) (Figure 1B). Endothelial cells are large enough to detect gradients of chemoattractants across their length, which is typically 20 μm (Vadapalli et al., 2000). In contrast to existing models of cellular chemotaxis, in which cells typically respond to free, extracellular chemokine concentrations, or their gradients in the surrounding tissue, the model developed here will capture the response of ECs to VEGF that is bound to cell-surface receptors. VEGF-VEGFR2 binding is known to be the signal that initiates endothelial motility, therefore incorporating this molecular response is important for a realistic description of cell motion.

The cell detects bio-available VEGF by taking it up at the half-lattice sites (Figure 1B). The model will thus make the crucial distinction between VEGF that is free to bind to the cell, versus VEGF that might be sequestered in the underlying extracellular matrix, unavailable to the cell. Based on the numbers of activated receptors at its four sides, the cell becomes polarized and attains a bias in a particular direction. It correspondingly elongates in this direction (Figure 1C). Finally, the rear of the cell detaches from the underlying matrix and contracts, and the cell has now moved to the site (n + 1, m) (Figure 1D). We remark that it has been observed experimentally that ECs may respond to chemoattractant concentration differences of as small as 2% across their length, frequently at concentrations at which molecular fluctuations are significant (Mantzaris et al., 2004). In our model, fluctuations of the order of 100 molecules of VEGF per cell are significant enough to alter its polarization, and hence its direction of motion.

2.2. A Novel Chemotactic Sensitivity Function

An important difference that sets this model apart from those preceding it, is the choice of the chemotactic sensitivity function χ(C). Various choices have been proposed thus far in the modeling literature for χ(C), for a review of the most commonly used functional forms (see Ford and Lauffenburger, 1991). The simplest choice is to assume that the chemotactic sensitivity is constant, χ(C) = χ₀ (Keller and Segel, 1971a,b). However, this implies that the chemotactic sensitivity is unchanging in the presence of the chemoattractant, and does not account for the desensitization of cells which has been experimentally observed to occur in regions of high chemokine concentrations (Kuppuswamy and Pike, 1989; Wang et al., 2000; Kurt et al., 2001). To overcome this, Lapidus and Schiller (1976), and later Murray (2003) used the functional form χ(C) = χ₀/(K + C)², also known as the receptor-kinetic law. This has the advantage that it is able to account for the desensitization of receptors when c is large. Yet another popular phenomenological choice is χ(C) = χ₀/(K + C), where K is the dissociation constant of the chemokine binding to the receptors (Balding and McElwain, 1985; Anderson and Chaplain, 1998b; Plank et al., 2003).

Although the choices mentioned above have been widely used in the angiogenesis literature, there are a few biological issues that these choices do not address. Firstly, they indicate that when no chemokine is present at a site, the chemotactic sensitivity is the greatest, and that the sensitivity decreases as chemokine concentrations increase. However, there is experimental evidence which points toward the existence of a minimum threshold level of chemical stimulus required for the cell-surface receptors to become activated, and for the cell to start moving in a directed fashion (Favier et al., 2006; Liu et al., 2006). This threshold has been incorporated in the cell-based model of tumor angiogenesis by Bauer et al. (2007). Secondly, the above functions do not account for the fact that the amount of chemokine required to desensitize cells depends on the concentration of cells present at the lattice site. For instance, while 10 fg of VEGF is enough to desensitize a single EC, it is not enough for 10 cells. Finally, for an external chemical signal to elicit a chemotactic response from a cell, it needs to be detected by the cell, and transduced into an internal signal controlling cell motion. Neutrophils have been shown to sense chemical gradients of 1% across their lengths, under optimal conditions (Wang et al., 2004; Levine et al., 2006), while this number can be as low as 0.1% for axons (Wang et al., 2004). In general, eukaryotic cells are reported to be able to polarize and migrate in a directed fashion in alignment with chemical gradients of about 2% across their lengths (Franz et al., 2002). It is therefore biologically more reasonable to assume that the chemotactic response of cells is dependent on the gradients of activated receptor complexes formed on the cell surface when the chemokine binds to its receptors, rather than gradients of free chemokine concentration throughout the tissue.

To address these concerns, we propose that the chemotactic sensitivity function should in fact be a function of the activated receptor concentration, A. In this case, equation (4) for the motion of a cell in 2 dimensions transforms to the following:

$\frac{\partial p}{\partial t}=D_p\Delta p-\nabla \cdot \left(\chi \left(A\right)p\nabla A\right).$ (6)

Correspondingly, the transition probabilities in equation (1) will now be functions of concentration of activated receptor complexes on the cell surface a and not extracellular VEGF. That is,

$\begin{array}{c}{\mathcal{T}}_{n,m}^{H\pm }=\frac{1}{k}\frac{\tau \left(A_{n\pm 1∕2,m}\right)}{\begin{array}{c}\tau \left(A_{n+1∕2,m}\right)+\tau \left(A_{n-1∕2,m}\right)\hfill \\ +\tau \left(A_{n,m+1∕2}\right)+\tau \left(A_{n,m-1∕2}\right)\hfill \end{array}},\\ {\mathcal{T}}_{n,m}^{V\pm }=\frac{1}{k}\frac{\tau \left(A_{n,m\pm 1∕2}\right)}{\begin{array}{c}\tau \left(A_{n+1∕2,m}\right)+\tau \left(A_{n-1∕2,m}\right)\hfill \\ +\tau \left(A_{n,m+1∕2}\right)+\tau \left(A_{n,m-1∕2}\right)\hfill \end{array}}.\end{array}$ (7)

We have to add equations for the binding of chemokine to their cell-surface receptors, which will need to be solved wherever a cell is present (see Section 2.3). Biologically, the chemotactic sensitivity χ(A)▽A can be interpreted by breaking it down as follows: a velocity χ(A) imparted to the cell due to the presence of bound chemokine on its surface, and a gradient ▽A which governs the direction of motion. This gradient simply means that the cell is able to sense the amount of chemokine bound to its various faces, and is correspondingly able to align itself for motion in this direction. Therefore, a is in fact taken to be the amount of activated receptors per cell face. We choose a velocity function that satisfies the requirements that there can be no chemotaxis in the absence of a signal, and that the cell gets desensitized in the presence of excess signal. One such functional form is:

$\chi \left(A\right)={\chi }_{0}A{e}^{-A∕K}.$ (8)

The maximum of this function occurs at A = K, while its maximum value is given by χ₀ Ke⁻¹. In order for this choice to be consistent with the discrete formulation, the function τ(A), from equation (3) must be taken as follows;

$\tau \left(A\right)=exp\left[\frac{{\chi }_{0}K}{D_p}\left(K-\left(K+A\right)e^{-A∕K}\right)\right].$ (9)

The parameters χ₀ and K are unknown in our model formulation, and would ideally be determined from experimental observations. K specifies the fractional occupancy of the receptors on the cell surface at which its chemotactic response is the greatest, while χ₀ determines the maximum value of this response. Here, values of these parameters are chosen to produce biologically realistic simulation results. Figure 2B plots the chemotactic sensitivity equations (8) as a function of the fraction of activated receptors on a cell face, for a particular choice of K and χ₀. We can see that at zero fractional activation, the cell remains inactive. The sensitivity peaks at 5% fractional activation of receptors, and decays thereafter. Also shown for comparison are the receptor-kinetic law, and constant chemotactic sensitivity.

FIGURE 2

Figure 2. (A) Typical profile of unbound VEGF, the source of which is located at x = 1. (B) Various choices for the chemotactic sensitivity of an endothelial cell to VEGF bound to its surface receptors, as a function of the fraction of activated VEGFR2 per cell face. (C–F) Average migration times (in hours) for a single cell to travel across a 1 mm × 1 mm domain as a function of increasing the maximum free VEGF concentration at a lattice site, for various choices of the chemotactic sensitivity function. Cells are assumed to respond to activated VEGFR2 on their surfaces, with chemotactic sensitivity taken as: (C) as proposed in equation (8); (D) cells are assumed to respond to free VEGF with chemotactic sensitivity as defined in equation (8) with activated receptor concentration A replaced by free VEGF concentration C; (E) receptor-kinetic law, for which τ(A) = e^{χ0A/(KDp(K+A))}, χ₀ = 0.4416 (pg/mm³) mm²/h, K = 2 pg/mm³; and (F) constant, χ₀ = 0.0046 mm²/(pg/mm³)/h.

2.3. VEGF-VEGFR2 Binding Dynamics

We now describe the equations governing the rates of change of the concentrations of free VEGF (C), free VEGFR2 (R), VEGF-VERFR2 monomers (M), and activated VEGF-VEGFR2 dimer complexes (A). Beginning with free VEGF, we assume that the processes of diffusion and natural decay dominate the dynamics, which are represented by the reaction-diffusion equation:

$\frac{\partial C}{\partial t}=D_c\Delta C-{\alpha }_{c}C-f\left(p\right)C.$ (10)

Here, D_c is the diffusion coefficient of VEGF, and α_c is its rate of decay in tissue. The uptake of VEGF by the migrating EC has also been accounted for via the term f(p)C, which is derived in the following discussion. As in Anderson and Chaplain (1998b), a line source of tumor cells is assumed at x = 1 that produces VEGF at a constant rate, so that C(1, y, t) = C₀. At each of the remaining domain boundaries x = 0 and y = 0, 1, a no-flux condition is imposed on free VEGF. Figure 2A shows the distribution of VEGF, expressed in non-dimensional terms, as determined by equation (10) across the domain.

At each half-lattice site surrounding a site occupied by an EC, free extracellular VEGF (C) binds to free cell-surface receptors, VEGFR2 (R), to form activated dimerized receptor complexes (A). Following Jain et al. (2008), we assume ligand-induced dimerization to be the dominant mechanism by which VEGF activates VEGFR2, as represented by the following chemical reactions:

$\begin{array}{ccc}\hfill C\hfill & \hfill +\hfill & \hfill R\underset{k_{r1}}{\overset{k_{f1}}{\rightleftharpoons }}M\hfill \\ \hfill M\hfill & \hfill +\hfill & \hfill R\underset{k_{r2}}{\overset{k_{f2}}{\rightleftharpoons }}A\hfill \\ \hfill \hfill & \hfill \hfill & \hfill A\stackrel{k_p}{\to }2R\hfill \\ \hfill \hfill \end{array}$

the rates of forward reactions are indicated above the reaction arrows, while those of reverse reactions are indicated below the reaction arrow.

This reaction diagram can be converted to the following system of differential equations using principles of mass balance:

$\begin{array}{lll}\frac{dC}{dt}& =-2{\eta }_1{k}_{f1}CR+{\eta }_2{k}_{r1}M,& \end{array}$ (11)

$\begin{array}{lll}\frac{dR}{dt}& =-2k_{f1}CR+{\eta }_3{k}_{r1}M-k_{f2}MR+2{\eta }_4{k}_{r2}A+2{\eta }_4{k}_{p}A,& \end{array}$ (12)

$\begin{array}{lll}\frac{dM}{dt}& =2{\eta }_5{k}_{f1}CR+k_{r1}M-{\eta }_5{k}_{f2}MR+2{\eta }_6{k}_{r2}A,& \end{array}$ (13)

$\begin{array}{lll}\frac{dA}{dt}& ={\eta }_7{k}_{f2}MR-2k_{r2}A-k_{p}A.& \end{array}$ (14)

The multiplicative factor 2 in some of the equations accounts for the possibility that there may be two ways for that product to form. The constants η_i represent the ratios of weights of different molecules and have been introduced to express chemical concentrations in units of pg/mm³. The values of these constants are given in Table 1 and have been estimated from Ferrara et al. (2003), Stewart et al. (2003).

TABLE 1

Table 1. Parameter values relating to the molecular weights of VEGF and VEGFR2.

Since EC migration and sprout elongation occurs on a time-scale of several hours to days, and the biochemical reactions equations (11–14) occur on a time-scale of several scones to minutes, we assume that the VEGF-receptor complex concentrations M and A are at quasi steady state. This is equivalent to setting the left hand sides of equations (13) and (14) to zero, and solving for M and A. Further, by conservation of total receptor number, we have R + η₃ M + 2η₄ A = R_fN, where R_f is the total number of receptors per EC face and N is the number of ECs (N = 1 in the case of single cell migration, and N = the total number of tip cells in the case of capillary formation). We therefore deduce that at quasi steady state,

$A=\frac{\begin{array}{c}-2\alpha \delta -\gamma +\theta \beta \hfill \\ +\sqrt{{\left(2\alpha \delta +\gamma -\theta \beta \right)}^2+4\alpha \beta \left({\eta }_4+\delta \right)\left(\theta -\alpha \right)}\hfill \end{array}}{2\beta \left({\eta }_4+\delta \right)},$ (15)

where,

$\begin{array}{llll}\alpha & =\frac{k_{r1}{R}_f}{2k_{f1}C+k_{r1}},\beta =\frac{{\eta }_4\left(k_p-2k_{r1}\right)}{2k_{f1}C+k_{r1}},& & \\ \gamma & =\frac{2k_{r2}+k_p}{{\eta }_5{\eta }_7{k}_{f2}}\delta ={\eta }_4+\beta ,\theta =R_{f}N.& \end{array}$ (16)

In all that follows, equation (15) will be used to estimate the concentration of activated VEGF-VEGFR2 dimers in the domain. Adding equations (12) and (14) and substituting in equation (11) gives the following equation for the uptake of VEGF by ECs:

$\frac{dC}{dt}=-{\eta }_0{k}_{p}A.$ (17)

Therefore, in equation (10) the cellular uptake function f(p)C = −η₀ k_pAI(p), where I(p) is an indicator function that has a value of 1 at half-lattice sites where EC boundaries are present and is zero otherwise. Observing that η₀ k_p ≪ α_C (see parameter values in Table 2), we make a final simplifying assumption that due to the constant production and rapid diffusion of extracellular VEGF, cellular uptake will not significantly effect its concentration, that is, f(p)C is neglected. Thus, the equation governing free VEGF dynamics is taken to be

$\frac{\partial C}{\partial t}=D_c\Delta C-{\alpha }_{c}C.$ (18)

TABLE 2

Table 2. List of parameter values for single cell motion.

2.3.1. Summary of model equations

The principle variables in our model are: p(n, m, t), the probability that a cell occupies lattice site (n, m) at time t; C(x, y, t), the concentration of free VEGF at position (x, y) and at time t in pg per lattice site volume, where each lattice site has a height of 1 mm and a base equal to the surface area of a cell; R(i, j, t), the concentration of free VEGFR2 at half-lattice sites (i, j) and at time t in pg per lattice site volume; M(i, j, t), the concentration of VEGF-VEGFR2 monomers at half-lattice sites (i, j) and at time t in pg per lattice site volume; and A(i, j, t), the concentration of free activated VEGFR2-VEGF-VEGFR2 dimers at half-lattice sites (i, j) and at time t in pg per lattice site volume. We remark that R, M, and A can only take positive values at neighboring half-lattice sites where a cell is present, and are 0 otherwise. Equations (1), (7), and (9) that describe the biased random walk of a cell under the influence of activated VEGF-receptors have already been discussed. The following conditions are imposed on the transition probabilities for cell motion as described by equation (1) to ensure that no cell exits the domain:

${\mathcal{T}}_{1,m}^{H-}\left(\cdot \right)={\mathcal{T}}_{N_s+1,m}^{H+}\left(\cdot \right)={\mathcal{T}}_{n,1}^{V-}\left(\cdot \right)={\mathcal{T}}_{n,N_s+1}^{V+}\left(\cdot \right)=0.$ (19)

Here, N_s = 1/h, h being the lattice size so that 1 ≤ n, m ≤ N_s + 1.

2.4. Parameter Estimation

A list of parameter values and sources is given in Table 2. The random motility coefficient of endothelial cells has been estimated to lie within the range 7.2 × 10⁻⁴–7.2 × 10⁻³ mm²/h (Anderson and Chaplain, 1998b). Consequently, intermediate value of 1.44 × 10⁻⁴ mm²/h is assumed. The rates k_r1 and k_f1 are chosen to ensure that the equilibrium disassociation constant k_D = k_r1/k_f1 has a value of 30.375 pg/mm³ (Wang et al., 2002). VEGF binding is known to induce receptor aggregation; therefore, as in Jain et al. (2008) we assume that the rate of formation of a dimerized VEGF-VEGFR2 complex is greater than the rate of formation of a monomer VEGF-VEGFR2 (that is, k_f2 ( k_f1). Further, because the dimerized complex A is the signaling form of VEGFR2, it is reasonable to assume that A is more stable than the monomer complex M, that is, k_r2 = k_r1. The size h of the lattice on which the cell moves is taken to be 20 μm, since typical microvascular endothelial cell volume is about 400 μm (Vadapalli et al., 2000), while its thickness is about 1 μm (Levine et al., 2002). Finally, the parameters χ₀ and K relating to chemotactic sensitivity are chosen to reproduce cell motion and capillary formation profiles that are biologically realistic.

2.5. Method of Simulation

The time interval over which the movement of the cell is simulated, is divided into subintervals of length k, given by the mean waiting time of the cell at any lattice site. The cell moves on a lattice of step-size h. Activated VEGFR2 concentrations are calculated at half-lattice sites, neighboring a site where the cell is currently situated. The method of simulation of cell movement is based on that described in Plank et al. (2003). Briefly, at each time step, the movement of the cell is simulated according to the master equation (1), with the probabilities of moving up, down, left, and right calculated according to equation (3). Equation (9) quantifies the dependence of the transition probabilities on the levels of activated VEGFR2 on each cell face as given by equation (15). The full interval [0, 1] is divided into five subintervals, each of length proportional to the probabilities of moving or staying still. A random number q lying within this interval is generated, and depending on the sub-interval in which it lies, the cell either executes a motion in the corresponding direction or stays stationary. Thus, the cell moves left if $q\in \left[0,{\mathcal{T}}_{n,m}^{H-}\right),$ moves right if $q\in \left[,{\mathcal{T}}_{n,m}^{H-}+{\mathcal{T}}_{n,m}^{H+}\right),$ and so on.

The results of the single cell motion model are discussed in Section 3.1.

2.6. Adaptation of Single Cell Motion Model to Simulate Capillary Formation

To simulate capillary formation in response to a VEGF stimulus from a tumor source, we modify the single cell motion model described above as follows. As mentioned earlier, we motivate our model of capillary formation by the experiments of Sholley et al. (1984) wherein inflammatory neovascularization of the rat cornea was induced by cauterization using silver nitrate and levels of EC proliferation and degree of vascular profusion measured periodically. From these experiments, the average rate of sprout extension into the cornea is estimated to be 0.26 mm/day or 0.78 mm in 3 days. In our model, we do not account for vessel maturation; a process that typically occurs after 3 days of vessel formation. Consequently, we simulate vessel growth for lengths ≤0.78 mm. With this constraint, a parent vessel, from which sprout tips will migrate toward the tumor, is assumed at x = 0.22 mm. As in the single cell model, a line of tumor cells is assumed at x = 1, providing a constant source of VEGF. For ease of computations, the domain size is reduced to 0.5 mm in the y-direction, and initially 4 sprouts are assumed to have formed along the parent vessel at y = 0.1, 0.2, 0.3, and 0.4 mm.

It is known that specialized ECs situated at the tips of the sprouts, called tip cells, are activated by, and respond to VEGF, by chemotactic migration (Hangai et al., 2002; Gerhardt et al., 2003). We therefore keep track of these leading cells in our simulations. As a tip cells moves, it pulls behind it a developing vessel. Hence, receptors on its tail are made unavailable for binding VEGF at any given time. This eliminates the possibility for the tip cell to back-track. By keeping track of all the lattice sites a tip cell visits, we know the location of the newly formed vessel behind it.

The processes of branch formation and anastomoses formation of loops by capillary sprouts are also included explicitly in our model. At each time step, as the tip cells migrate under the influence of VEGF, probabilities of motion to adjacent lattice sites are calculated. Anastomoses between the tip cell and a sprout may occur if a sprout is present at a site which the tip cell wants to move to. We assume that the probability of tip cell loss as a result of such an event is 1%. Likewise, as in Anderson and Chaplain (1998b), it is assumed that if another tip cell is encountered at a site, only one of these cells continues to grow (with a probability 99%), while the rest of the time, a loop is created with the loss of both cells.

Sholley et al. (1984) have demonstrated that sprout extension cannot occur in the absence of mitosis. While we do not explicitly model cell division, the dependence of capillary extension on it is accounted for in the processes of capillary elongation and branch formation as follows. The proliferation of cells is known to be regulated by total concentration of activated cell-surface receptors (Gerhardt et al., 2003). Thus, in our model, the tip cell integrates the total VEGF bound to it and sprout extension via tip cell motion, and branch formation is only possible if there are enough activated VEGFR2 on its surface. The effect of proliferation on tip cell motion is simulated by introducing a scaling factor of P_m(A_t) that multiplies the movement probabilities of each cell, where A_t is the total concentration of activated VEGFR2 on the cell. P_m is assumed to be a positive, increasing, and saturating function of A_t, with a saturating value of 1. Thus, for small values of A_t, the probability of capillary extension will be ∼ 0 due to an insufficient proliferation stimulus. Here, we take $P_m(A_t)=1/(1+{\mu }_m{e}^{-A_t})$ which is plotted as a function of the fraction of total activated VEGFR2 per tip cell in Figure 5D.

We further assume that the generation of new sprouts occurs only from existing sprout tips. This is in keeping with the fact that there is a region of proliferating cells just behind the tip cell (Sholley et al., 1984), which could give rise to new branches. As in the case for P_m, the branching probability P_b is also taken to be an increasing and saturating function of A_t. This will result in the creation of the brush-border effect. Similar rules for branching have been applied previously by Anderson and Chaplain (1998b). Here, we take P_b(A_t) = 1/(μ_b1 + e^{−μb2(At − A0)}) which is plotted as a function of the fraction of total activated VEGFR2 per tip cell in Figure 5D.

2.6.1. Parameter estimation for capillary formation model

A list of parameter values that are different or new in the capillary formation model is given in Table 3. For consistency with the single cell model, we keep the time step-size k unchanged at 0.07 h. Further, the diffusion rate of a tip cell, say D_t, has been estimated to be much smaller than the diffusion rate D_p of an individual EC (Anderson and Chaplain, 1998b; Levine et al., 2001). Therefore, the lattice size h for the capillary model needs to be altered accordingly. Equation (5) is used to estimate $h=\sqrt{4k{D}_t}\approx 0.001\text{mm}\text{.}$ Finally, the parameters μ_m, μ_b1, μ_b2, and A₀ relating to the movement probability and branching probability are chosen to reproduce capillary formation profiles that are biologically realistic.

TABLE 3

Table 3. List of parameter values for capillary formation.

2.6.2. Simulation methodology for capillary formation model

The simulation methodology is similar to that of single cell motion described in section 2.5, with the additional computation of accounting for branching and anastomoses for each tip cell, and at each time step. Briefly, in addition to generating a random number q, which is used to determine the direction of tip cell motion, two further random numbers are generated (q_a and q_b) by uniformly sampling the interval [0, 1]. We use q_a to determine whether or not anastomoses occurs, and q_b is used to determine whether a new branch forms, in accordance with the rules described above.

The results of the capillary formation simulation are discussed in Section 3.2.

3. Results

3.1. Single Cell Motion

Simulations of the system governing a single endothelial cell migrating up a gradient of VEGF, as described by equations (1), (7), (9), and (15), were run in two dimensions, with unbound VEGF profile described by equation (18). The average time in hours it takes for the cell to travel across the domain is plotted in Figures 2C–F, as a function of C₀, the maximum free VEGF concentration at a lattice site, for various possible choices of the chemotactic sensitivity function χ(·). Standard deviations and average times are computed over 500 runs of the model.

Figures 2C,E,F depict the cases when χ(·) is assumed to be a function of activated VEGFR2 on the cell surface. When χ(·) is as defined in equation (8), the model captures the existence of a minimum level of VEGF stimulus required for directed cell motion, as well as desensitization of VEGFR2 at high VEGF concentrations (see Figure 2C). As C₀ increases from 0.002 pg/mm³, the average EC migration time is observed to first decrease and then increase, attaining a minimum of 8.23 h at C₀ = 0.015 pg/mm³. A typical cellular trajectory is plotted in Figure 3C for this optimal value of C₀, and the corresponding movement probabilities at any lattice site are plotted in Figure 3D. Note that since the VEGF profile is invariant along the y-direction, the movement probabilities are also invariant along this axis – they only vary as x varies. The probabilities show a large bias toward stepping to the right, while steps to the left are very unlikely to occur. This is because: (i) the chemokine gradient across the cell length has an average value of 1.32%, over the entire domain, which lies within the reported value of 1–2% at which eukaryotic cells become polarized; and (ii) the fraction of activated receptors on any cell face is sufficiently large, with an average value of 9%, over the entire domain.

FIGURE 3

Figure 3. (A,C,E) Typical trajectories of a cell migrating across a 2-dimensional domain under the influence of VEGF. Here, A represents the concentration of activated VEGFR2 per cell face and A_max represents the maximum value A can take so that A/A_max is the fraction of activated VEGFR2 per cell face expressed here as a percentage, and <A_grad> represents the gradient of A across a cell length, averaged over the entire domain. (B,D,F) Corresponding movement probabilities for various values of maximum VEGF concentration.

As C₀ is decreased below 0.015, the average migration time is predicted to increase exponentially. For instance, when C₀ = 0.002 pg/mm³, the average migration time is predicted to be 37.46 h, and the cell exhibits a high degree of randomness in its motion, as evident from a typical cellular trajectory shown in Figure 3A. The corresponding movement probabilities at any lattice site plotted in Figure 3B show that a definite bias is apparent for motion to the right only close to x = 1. This is because the fraction of activated VEGFR2 on any cell face is very low, with a maximum of < 2%, even though the chemokine gradient across the cell length has an average value of 1.46%. Thus, the model is able to account for the fact that if chemokine concentrations are too low, cell-surface receptors do not achieve a sufficient degree of activation.

As C₀ is increased beyond its optimal value of 0.015–0.08 pg/mm³, the model replicates the desensitization effect which has been observed to occur when receptors are over-exposed to chemokines. It now takes the cell an average of 35.77 h to migrate across the domain. From Figure 3E, we observe that typical cell trajectories exhibit a large degree of random motion. Now, activated receptor gradients across the cell have an average value of only 1%. Further, the fraction of activated receptors that vary between 18 and 34% across the domain so that the negative exponential in equation (8) dominates resulting in a very slight bias of movement to the right (see Figure 3F).

For comparison, we also consider the cases where χ(A) = χ_0/(K + A)² or the receptor-kinetic law and when χ(A) = constant = χ₀. As can be seen from Figure 2E, while the receptor-kinetic law captures the desensitization of VEGFR2 at high concentrations of VEGF, the cell still displays a high degree of directed motion for very low values of C₀. For instance, when C₀ = 0.002 pg/mm³ the average migration time is as low as 16.08 h as compared to 37.46 h in the earlier case. In contrast, the existence of a minimum activation threshold for VEGF is predicted by assuming χ(A) = χ₀, as evident from Figure 2F. However, this model is unable to capture receptor desensitization at high values of C₀, and in fact, the average migration time is predicted to decrease monotonically with C₀.

Finally, for illustration purposes, we also consider the case when χ(·) has the same qualitative properties as in equation (8), but the cell now responds to free VEGF rather than activated VEGFR2, that is, χ(C) = χ_0Ce^−C/K (see Figure 2D). While the graph is qualitatively similar to Figure 2C, the fastest migration of the EC across the domain is occurs at C₀ = 7.5 pg/mm³. This is biologically implausible since for such high receptor activation levels, the fraction of activated VEGFR2 on any cell face > 0.97 throughout the domain, and the cell should be completely desensitized to the chemical gradient around it.

3.1.1. Effect of receptor expression level on cell migration

An important parameter in our simulations of EC migration is R_T, the expression level of VEGFR2 per cell. This is known to be highly variable across cell lines, and it is even possible to find different values for R_T for the same cell line. We therefore conduct a sensitivity analysis on the migration times of an EC across the domain as R_T is varied, the results of which are graphed in Figure 4. For the baseline simulations discussed earlier, a value of R_T = 230,000 receptors per cell or 0.08 pg/cell (Stewart et al., 2003; Mac Gabhann and Popel, 2004) was used (see Figures 2 and 3). We now simulate the effect on cell migration of increasing R_T from a minimum of 46,000 to a maximum of 1,115,000 receptors per cell for various values of C₀, the maximum free VEGF concentration at a lattice site. For each of the cases when R_T = 115,000 (Figure 4B), R_T = 230,000 (Figure 2C), R_T = 460,000 (Figure 4C), and R_T = 1,150,000 (Figure 4D), the average migration time of the EC is predicted to first decrease and then increase, as C₀ is increases. Thus, for a large range of values of R_T, the model captures the existence of an activation threshold of VEGFR2, and their desensitization when exposed to high VEGF concentrations. However, when R_T is very low (46,000/cell, Figure 4A), receptor desensitization is not predicted. This is possibly due to a high value of the parameter K, which is held fixed in all our simulations. As can be seen from equation (8), K determines the concentration of activated VEGFR2 per cell face at which the chemotactic sensitivity χ(A) is maximum.

FIGURE 4

Figure 4. (A–D) Average migration times (in hours) for a single cell to travel across a 1 mm × 1 mm domain as a function of increasing the maximum free VEGF concentration at a lattice site, for different values of R_T, the total number of VEGFR2 per cell. Baseline simulations correspond to R_T = 230,000/cell, and are shown in Figure 2C. (E) Minimum migration time for a single cell to cross the domain as a function of increasing R_T, solid triangle corresponds to baseline simulations. (F) The maximum free VEGF concentration at a lattice site (C₀) at which cell migration is fastest, as a function of increasing R_T, solid triangle corresponds to baseline simulations.

Next, as can be seen from a plot of fastest migration times versus receptor expression in Figure 4E, the EC migrates more rapidly across the domain as R_T increases. The fastest migration time is predicted to be 8.25 h, for R_T ≥ 460,000/cell. Interestingly, the maximum free VEGF concentration at which EC migration is fastest decreases with increasing R_T (Figure 4F). Thus receptor over-expression is predicted to lower the activation threshold for ECs, possibly because gradients of activated VEGFR2 become more pronounced across the cell.

3.2. Capillary Formation

Simulations of the system governing capillary network formation under the influence of VEGF, described in section 2.6, were run in two dimensions. Averages and standard deviations of all observed quantities were calculated from 100 runs of the model.

The first case considered is when the sprout tips move across the domain in the least amount of time. This occurs when the maximum concentration of unbound VEGF C₀ = 0.015 pg/mm³, as deduced from the single cell simulations. The results from a typical simulation are shown in Figure 5A. We begin with 4 initially formed sprout tips at x = 0.22 mm. As the tip cells migrate across the domain, they lay down behind them capillary sprouts. As the vascular network penetrates deeper into the stroma, branching is observed to occur leading to the brush-border effect. The model predicts that it takes on average 1170 ± 27 steps or 3.38 ± 0.08 days for the vasculature to reach the tumor source at x = 1 mm. Our model is validated by the experiments in Sholley et al. (1984) where the vascular sprouts traveled the same distance in 3 days. Further validation follows by observing that the vascular networks generated by our model are qualitatively similar to those observed by Sholley et al. (1984).

FIGURE 5

Figure 5. (A–C) Typical vascular networks formed by 4 initial sprouts located along x = 0.22 at positions y = 0.1, 0.2, 0.3, 0.4; x being plotted along the abscissa and y along the ordinate – migrating across a 2-dimensional domain under the influence of VEGF for various values of C₀, the maximum VEGF concentration per lattice site. (A) Optimal VEGF concentration, C₀ = 0.015 pg/lattice volume. The bias of movement is overwhelmingly in the forward direction. Branching and anastomoses are observed to occur as the vasculature penetrates deeper into the stroma. The resulting networks are qualitatively similar to those observed experimentally in Sholley et al. (1984). (B) C₀ = 0.005 pg/lattice volume. The amount of VEGF is too low to induce proliferation or polarization of the tip cell, leading to a poorly developed and stunted vasculature that does not reach the VEGF source within the time frame of simulations (3 days). (C) C₀ = 0.030 pg/lattice volume. Due to a high VEGF concentration, over-stimulation of endothelial cells occurs, and extensive branching, anastomoses and lateral movement of the tip cell is observed. Due to excessive lateral movement, the vasculature that does not reach the VEGF source within the time frame of simulations (3 days). (D) Assumed branching probability P_b of the migrating tip cell (black curve), and assumed extension probability P_m of the capillary (blue curve), expressed as functions of total fraction of activated VEGFR2 per cell.

Next, the effects of low (0.005 pg/mm³) and high (0.030 pg/mm³) maximum VEGF concentrations on vascular formation are investigated. As remarked earlier, we do not model vessel maturation, which is typically observed after ∼3 days of vessel formation. Model simulations are run for a maximum of 1170 time steps and the average degree of vascular penetration into the stroma, along with the fraction of sprouts that remain viable (that is, have at least one active tip cell) at the end of this time is computed. When C₀ = 0.005 pg/mm³, the average lengths of sprouts formed is predicted to be 0.2 ± 0.1 mm, with only 40% of the initial sprouts still viable after 1170 time steps. Sprouts that remain viable after 1170 time steps extend a greater distance (0.4 ± 0.03 mm) into the stroma. However, these display virtually no branching, with the average number of branches per sprout only 1.1 ± 0.4. Figure 5B shows the results of a typical simulation. As can be seen, there has been no branching and all but the first sprout have anastomosed with themselves to form closed loops. This is due to an insufficient bias to move forwards, coupled with a low value of the scaling factor P_m (see section 2.6).

Finally, when C₀ = 0.030 pg/mm³, the average lengths of sprouts formed after 1170 time steps is predicted to be 0.6 ± 0.1 mm. As can be seen from a typical simulation shown in Figure 5C, extensive branching and anastomoses are observed. Given the density of vessel branches, it is reasonable to expect that several of these may fuse into one another resulting in thicker and more dilated vessels, which is a morphology consistent with vascular hyperplasia, as seen in Lee et al. (2005). The higher VEGF concentration implies that the vessels have a weaker bias for forward motion, and lateral movement of vessels as well as movement against the gradient of VEGF are observed to occur. These phenomenon have been observed in vivo, and have been numerically simulated previously (Anderson and Chaplain, 1998b; Plank and Sleeman, 2004; Sun et al., 2005; Zheng et al., 2013).

4. Conclusion

We have developed a hybrid model of cellular chemotaxis and capillary formation under the influence of VEGF. The migrating cell, whether by itself or as the tip cell “pulling” behind it a developing sprout, was treated as an agent. Its movement was simulated stochastically with movement probabilities based on the theory of biased random walks. On the other hand, due to its fast diffusion coefficient, VEGF dynamics were governed by a continuum reaction-diffusion equation. Using this approach, we first simulated the motion of a single cell on a two-dimensional grid, following the gradient of VEGF laid down by a constant source. Next, our model was adapted to simulate the formation of new vessels from pre-formed sprouts along a parent vessel, also under the influence of a constant source of VEGF, such as a tumor. Events such as branching and anastomoses, which are observed to occur in vivo, were incorporated explicitly in the model. The rate of vessel formation closely matched that observed experimentally (Sholley et al., 1984) under an optimal VEGF concentration. Additionally, as the forming vessels neared the VEGF source, a brush-border effect due to increased branching was predicted, thus proving both quantitative and qualitative validation of our approach. Using this framework, we also tested the effects of excessive as well as low levels of VEGF signaling on vascular development. Insufficient chemotactic and mitotic cues from VEGF resulted in stunted and solitary vessels, while an over-stimulation induced a high degree of branching and lateral movement.

An important difference that sets our model apart from similar hybrid models of chemotaxis is the inclusion of a molecular level detail of interaction between VEGF and its cell-surface receptor VEGFR2, the activation of which polarizes the cell and induces directed motion. This has been observed experimentally as well – endothelial cells respond to gradients of chemokines across their lengths, rather than to free chemokine concentrations. These gradients have been shown to be between 1 and 2%, which was seen in the numerical simulations as well, thus validating our model. Crucially, a chemotaxis sensitivity function was proposed that incorporated biological detail hitherto ignored by commonly used sensitivity functions currently. The model could thus capture realistic dynamics, such as the requirement of a minimum activation level of cell-surface receptors and receptor desensitization in high concentrations of VEGF.

Angiogenesis, both physiological and pathological, is a highly complex process, and understanding its mechanisms can lead to significant breakthroughs in the treatment of diseases such as cancer that depend on it. To this end, it is vital that modeling efforts keep up with current advances in experimentation. Our model provides such a framework, in which it is easy to build in biochemical and biomechanical forces guiding vessel formation. In fact, a number of highly detailed and complex hybrid models of vascular tumor growth have recently been proposed (Frieboes et al., 2007; Owen et al., 2009; Perfahl et al., 2011) and a significant strength of our model is that it can easily be incorporated into these. The inclusion of greater biological detail would only increase confidence in the predictive power of such models.

In addition, a number of refinements of the model proposed here are under active consideration. For instance, EC response to cell-surface bound VEGF has already been explicitly included. However, for the ease of computation, certain simplifying assumptions were made. Most notably, activated VEGFR2 were assumed to be in quasi steady state. Further, only the tip cell was tracked, while VEGF uptake by stalk cells was ignored. Cell death was also omitted, while the processes of cell proliferation, branching, and anastomoses were included phenomenologically. We plan to extend this model by relaxing some of these assumptions. Lattice-based models of angiogenesis face the criticism that the capillary networks generated by them are artificial to a certain extent, as they are forced to follow the lattice used to discretize the model. A first step would therefore be to develop a lattice-free version of our model of capillary formation, in which the ECs move without geometric constraints. Such models have been applied to capillary formation previously (Plank and Sleeman, 2004; Frieboes et al., 2007).

Other model refinements include incorporation of the relation between extra cellular matrix or ECM and vascular morphology. ECs require the ECM to gain traction in order to move. To facilitate their migration, ECs also secrete proteolytic enzymes such as matrix metalloproteinases (MMPs), that degrade collagen and elastin and clear a path for the ECs to follow. As ECs interact with the matrix, they also cause the release of matrix bound angiogenic factors such as VEGF, which are then available to induce further pro-angiogenic activity (Mantzaris et al., 2004). Further, pericytes, macrophages, and angiopoietins are also important determinants of developing vascular morphology and maturation (Levine et al., 2000; Plank and Sleeman, 2003), and need to be considered explicitly. The framework presented here is highly flexible, and would allow for the inclusion of the above processes, grounding it further in biology, and enhancing its usefulness as a tool to understanding the process of angiogenesis.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

References

Anderson, A. R. A., and Chaplain, M. A. J. (1998a). A mathematical model for capillary network formation in the absence of endothelial cell proliferation. Appl. Math. Lett. 11, 109–114.

CrossRef Full Text

Anderson, A. R. A., and Chaplain, M. A. J. (1998b). Continuous and discrete mathematical models of tumor-induced angiogenesis. Bull. Math. Biol. 60, 857–900.

CrossRef Full Text

Arakelyan, L., Vainstein, V., and Agur, Z. (2002). A computer algorithm describing the process of vessel formation and maturation, and its use for predicting the effects of anti-angiogenic and anti-maturation therapy on vascular tumor growth. Angiogenesis 5, 203–214.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Augustin, H. G., Koh, G. Y., Thurston, G., and Alitalo, K. (2009). Control of vascular morphogenesis and homeostasis through the angiopoietin-tie system. Nat. Rev. Mol. Cell Biol. 10, 165–177.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Balding, D., and McElwain, D. L. S. (1985). A mathematical model of tumor-induced capillary growth. J. Theor. Biol. 114, 53–73.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Bartha, K., and Rieger, H. (2006). Vascular network remodeling via vessel cooption, regression and growth in tumors. J. Theor. Biol. 21, 903–918.

CrossRef Full Text

Bauer, A., Jackson, T., and Jiang, Y. (2007). A cell-based model exhibiting branching and anastomosis during tumor-induced angiogenesis. Biophys. J. 92, 3105.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Bentley, K., Mariggi, G., Gerhardt, H., and Bates, P. A. (2009). Tipping the balance: robustness of tip cell selection, migration and fusion in angiogenesis. PLoS Comput. Biol. 5:e1000549. doi:10.1371/journal.pcbi.1000549

CrossRef Full Text

Billy, F., Ribba, B., Saut, O., Morre-Trouilhet, H., Colin, T., Bresch, D., et al. (2009). A pharmacologically based multiscale mathematical model of angiogenesis and its use in investigating the efficacy of a new cancer treatment strategy. J. Theor. Biol. 260, 545–562.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Byrne, H. M., and Chaplain, M. A. J. (1995). Mathematical models for tumour angiogenesis: numerical simulations and nonlinear wave solutions. Bull. Math. Biol. 57, 461–486.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Byrne, H. M., and Chaplain, M. A. J. (1996). Explicit solutions of a simplified model of capillary sprout growth during tumor angiogenesis. Appl. Math. Lett. 9, 69–74.

CrossRef Full Text

Capasso, V., and Morale, D. (2009). Stochastic modelling of tumour-induced angiogenesis. J. Math. Biol. 58, 219–233.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Conway, E. M., Collen, D., and Carmeliet, P. (2001). Molecular mechanisms of blood vessel growth. Cardiovasc. Res. 49, 507–521.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Favier, B., Alam, A., Barron, P., Bonnin, J., Laboudie, P., Fons, P., et al. (2006). Neuropilin-2 interacts with VEGFR-2 and VEGFR-3 and promotes human endothelial cell survival and migration. Blood 108, 1243–1250.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Ferrara, N. (2004). Vascular endothelial growth factor: basic science and clinical progress. Endocr. Rev. 25, 581–611.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Ferrara, N., Gerber, H. P., and LeCouter, J. (2003). The biology of VEGF and its receptors. Nat. Med. 9, 669–676.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Folkman, J. (1985). Tumor angiogenesis. Adv. Cancer Res. 43, 175–203.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Ford, R. M., and Lauffenburger, D. A. (1991). Analysis of chemotactic bacterial distributions in population migration assays using a mathematical model applicable to steep or shallow attractant gradients. Bull. Math. Biol. 53, 721–749.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Franz, C. M., Jones, G. E., and Ridley, A. J. (2002). Cell migration in development and disease. Dev. Cell 2, 153–158.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Frieboes, H. B., Lowengrub, J. S., Wise, S., Zheng, X., Macklin, P., Bearer, E. L., et al. (2007). Computer simulation of glioma growth and morphology. Neuroimage 37, 59–70.

CrossRef Full Text

Gerhardt, H., Golding, M., Fruttiger, M., Lundkvist, A., Ruhrberg, C., Abramsson, A., et al. (2003). VEGF guides angiogenic sprouting utilizing endothelial tip cell filopodia. J. Cell Biol. 161, 2263–1167.

CrossRef Full Text

Gevertz, J. L., and Torquato, S. (2006). Modeling the effects of vasculature evolution on early brain tumor growth. J. Theor. Biol. 243, 517–531.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Hanahan, D., and Weinberg, R. A. (2000). The hallmarks of cancer. Cell 100, 57–70.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Hangai, M., Kitaya, N., Xu, J., Chan, C. K., Kim, J. J., Werb, Z., et al. (2002). Matrix metalloproteinase-9-dependent exposure of a cryptic migratory control site in collagen is required before retinal angiogenesis. Am. J. Pathol. 161, 1429–1437.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Hicklin, D. J., and Ellis, L. M. (2005). Role of the vascular endothelial growth factor pathway in tumor growth and angiogenesis. J. Clin. Oncol. 23, 1011–1027.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Holmes, M. J., and Sleeman, B. D. (2000). A mathematical model of tumour angiogenesis incorporating cellular traction and viscoelastic effects. J. Theor. Biol. 202, 95–112.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Jackson, T. L., and Zheng, X. (2010). A multiscale model of cell elongation, proliferation and maturation during angiogenesis. Bull. Math. Biol. 72, 830–868.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Jain, H. V., Nor, J. E., and Jackson, T. L. (2008). Modelling the VEGF-Bcl-2-CXCL8 pathway in intra-tumoral angiogenesis. Bull. Math. Biol. 70, 89–117.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Jain, R. K. (2002). Tumor angiogenesis and accessibility: role of vascular endothelial growth factor. Semin. Oncol. 29(Suppl. 16), 3–9.

CrossRef Full Text

Jakobsson, L., Franco, C. A., Bentley, K., Collins, R. T., Ponsioen, B., Aspalter, I. M., et al. (2010). Endothelial cells dynamically compete for the tip cell position during angiogenic sprouting. Nat. Cell Biol. 12, 943–953.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Keller, E. F., and Segel, L. A. (1971a). Models of chemotaxis. J. Theor. Biol. 30, 225–234.

CrossRef Full Text

Keller, E. F., and Segel, L. A. (1971b). Traveling bands of chemotactic bacteria: a theoretical analysis. J. Theor. Biol. 30, 235–248.

CrossRef Full Text

Kuppuswamy, D., and Pike, L. J. (1989). Ligand-induced desensitization of 125I-epidermal growth factor internalization. J. Biol. Chem. 264, 3357–3363.

Pubmed Abstract | Pubmed Full Text

Kurt, R. A., Baher, A., Wisner, K. P., Tackitt, S., and Urba, W. J. (2001). Chemokine receptor desensitization in tumor-bearing mice. Cell. Immunol. 207, 81–88.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Lapidus, I. R., and Schiller, R. (1976). A model for the chemotactic response of a bacterial population. Biophys. J. 16, 779–789.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Lee, S., Jilani, S. M., Nikolova, G. V., Carpizo, D., and Iruela-Arispe, M. L. (2005). Processing of VEGF-A by matrix metalloproteinases regulates bioavailability and vascular patterning in tumors. J. Cell Biol. 169, 681–691.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Levine, H., Kessler, D. A., and Rappel, W. J. (2006). Directional sensing in eukaryotic chemotaxis: a balanced inactivation model. Proc. Natl. Acad. Sci. U.S.A. 103, 9761–9766.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Levine, H. A., and Nilsen-Hamilton, M. (2006). “Angiogenesis – a biochemial/mathematical perspective,” in Tutorials in Mathematical Biosciences III, Number 1872 in Lecture Notes in Mathematics, Chap. 2, ed. A. Friedman (Heidelberg: Springer), 65.

Levine, H. A., Pamuk, S., Sleeman, B. D., and Nilsen-Hamilton, M. (2001). Mathematical modeling of capillary formation and development in tumor angiogenesis: penetration into the stroma. Bull. Math. Biol. 63, 801–863.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Levine, H. A., Sleeman, B. D., and Nilsen-Hamilton, M. (2000). A mathematical model for the roles of pericytes and macrophages in the initiation of angiogenesis. I. The role of protease inhibitors in preventing angiogenesis. Math. Biosci. 168, 77–115.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Levine, H. A., Tucker, A. L., and Nilsen-Hamilton, M. (2002). A mathematical model for the role of cell signal transduction in the initiation and inhibition of angiogenesis. Growth Factors 20, 155–175.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Liu, G., Qutub, A., Vempati, P., Mac Gabhann, F., and Popel, A. (2011). Module-based multiscale simulation of angiogenesis in skeletal muscle. Theor. Biol. Med. Model. 8, 6.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Liu, Y., McCarthy, J., and Ladisch, S. (2006). Membrane ganglioside enrichment lowers the threshold for vascular endothelial cell angiogenic signaling. Cancer Res. 66, 10408–10414.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Mac Gabhann, F., and Popel, A. S. (2004). Model of competitive binding of vascular endothelial growth factor and placental growth factor to VEGF receptors on endothelial cells. Am. J. Physiol. Heart Circ. Physiol. 286, H153–1564.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Mac Gabhann, F., and 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, H2851–60.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Manoussaki, D. (2003). A mechanochemical model of angiogenesis and vasculogenesis. Esaim Math. Model. Numer. Anal. 37, 581–599.

CrossRef Full Text

Mantzaris, N., Webb, S., and Othmer, H. G. (2004). Mathematical modeling of tumor-induced angiogenesis. J. Math. Biol. 49, 111–187.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

McMahon, G. (2000). VEGF receptor signaling in tumor angiogenesis. Onoclogist 5(Suppl. 1), 3–10.

CrossRef Full Text

Milde, F., Bergdorf, M., and Koumoutsakos, P. (2008). A hybrid model for three-dimensional simulations of sprouting angiogenesis. Biophys. J. 95, 3146–3160.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Murray, J. D. (2003). Mathematical Biology II: Spatial Models and Biomedical Applications. New York: Springer-Verlag.

Muthukkaruppan, V. R., Kubai, L., and Auerbach, R. (1982). Tumor-induced neovascularization in the mouse eye. J. Natl. Cancer Inst. 69, 699–708.

Pubmed Abstract | Pubmed Full Text

Neufeld, G., Cohen, T., Gengrinovitch, S., and Poltorak, Z. (1999). Vascular endothelial growth factor (VEGF) and its receptors. FASEB J. 13, 9–22.

Pubmed Abstract | Pubmed Full Text

Nor, J. E., Christensen, J., Mooney, D. J., and Polverini, P. J. (1999). Vascular endothelial growth factor (VEGF)-mediated angiogenesis is associated with enhanced endothelial cell survival and induction of Bcl-2 expression. Am. J. Pathol. 154, 375–384.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Othmer, H. G., and Stevens, A. (1997). Aggregation, blowup and collapse: the ABC’s of generalized taxis in reinforced random walks. SIAM J. Appl. Math. 57, 1044–1081.

CrossRef Full Text

Owen, M. R., Alarcon, T., Maini, P. K., and Byrne, H. M. (2009). Angiogenesis and vascular remodelling in normal and cancerous tissues. J. Math. Biol. 58, 689–721.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Paweletz, N., and Knierim, M. (1989). Tumor-related angiogenesis. Crit. Rev. Oncol. Hematol. 9, 197–242.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Peirce, S. M. (2008). Computational and mathematical modeling of angiogenesis. Microcirculation 15, 739–751.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Peirce, S. M., Van Gieson, E. J., and Skalak, T. C. (2004). Multicellular simulation predicts microvascular patterning and in silico tissue assembly. FASEB J. 18, 731–733.

Pubmed Abstract | Pubmed Full Text

Perfahl, H., Byrne, H. M., Chen, T., Estrella, V., Alarcon, T., Lapin, A., et al. (2011). Multiscale modelling of vascular tumour growth in 3D: the roles of domain size and boundary conditions. PLoS ONE 6:e14790. doi:10.1371/journal.pone.0014790

CrossRef Full Text

Plank, M. J., and Sleeman, B. D. (2003). A reinforced random walk model of tumor angiogenesis and anti-angiogenesis strategies. IMA J. Math. Appl. Med. Biol. 20, 135–181.

CrossRef Full Text

Plank, M. J., and Sleeman, B. D. (2004). Lattice and non-lattice models of tumour angiogenesis. Bull. Math. Biol. 66, 1785–1819.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Plank, M. J., Sleeman, B. D., and Jones, P. F. (2003). A mathematical model of an in vitro experiment to investigate endothelial cell migration. J. Theor. Med. 4, 251–207.

CrossRef Full Text

Plank, M. J., Sleeman, B. D., and Jones, P. F. (2004). A mathematical model of tumour angiogenesis, regulated by vascular endothelial growth factor and the angiopoietins. J. Theor. Biol. 229, 435–454.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Qutub, A. A., Mac Gabhann, F., Karagiannis, E. D., Vempati, P., and Popel, A. S. (2009). Multiscale models of angiogenesis. IEEE Eng. Med. Biol. Mag. 28, 14–31.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Qutub, A. A., and Popel, A. S. (2009). Elongation, proliferation & migration differentiate endothelial cell phenotypes and determine capillary sprouting. BMC Sys. Biol. 3:13. doi:10.1186/1752-0509-3-13

CrossRef Full Text

Schugart, R. C., Friedman, A., Zhao, R., and Sen, C. K. (2008). Wound angiogenesis as a function of tissue oxygen tension: a mathematical model. Proc. Natl. Acad. Sci. U.S.A. 105, 2628–2633.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Serini, G., Ambrosi, D., Giraudo, E., Gamba, A., Preziosi, L., and Bussolino, F. (2003). Modeling the early stages of vascular network assembly. EMBO J. 22, 1771–1779.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Sholley, M. M., Ferguson, G. P., Seibel, H. R., Montour, J. L., and Wilson, J. D. (1984). Mechanisms of neovascularization. vascular sprouting can occur without proliferation of endothelial cells. Lab. Invest. 51, 624–634.

Pubmed Abstract | Pubmed Full Text

Sleeman, B. D., and Wallis, I. P. (2002). Tumour induced angiogenesis as a reinforced random walk: modeling capillary network formation without endothelial cell proliferation. Math. Comput. Model. 36, 339–358.

CrossRef Full Text

Stewart, M., Turley, H., Cook, N., Pezzella, F., Pillai, G., Ogilvie, D., et al. (2003). The angiogenic receptor KDR is widely distributed in human tissues and tumours and relocates intracellularly on phosphorylation. Histopathology 43, 33–39.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Stokes, C. L., and Lauffenburger, D. A. (1991). Analysis of the roles of microvessel endothelial cell random mobility and chemotaxis in angiogenesis. J. Theor. Biol. 152, 377–403.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Sun, S., Wheeler, M. F., Obeyesekere, M., and Patrick, C. (2005). A deterministic model of growth factor-induced angiogenesis. Bull. Math. Biol. 67, 313–337.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Tong, S., and Yuan, F. (2001). Numerical simulations of angiogenesis in the cornea. Microvasc. Res. 61, 14–27.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Travasso, R. D. M., Corvera Poire, E., Castro, M., Rodrguez-Manzaneque, J. C., and Hernandez-Machado, A. (2011). Tumor angiogenesis and vascular patterning: a mathematical model. PLoS ONE 6:e19989. doi:10.1371/journal.pone.0019989

CrossRef Full Text

Vadapalli, A., Pittman, R. N., and Popel, A. S. (2000). Estimating oxygen transport resistance of the microvascular wall. Am. J. Physiol. Heart Cric. Physiol. 279, 657–671.

Wang, D., Donner, D. B., and Warren, R. S. (2000). Homeostatic modulation of cell surface KDR and Flt1 expression and expression of the vascular endothelial cell growth factor (VEGF) receptor mRNAs by VEGF. J. Biol. Chem. 275, 15905–15911.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Wang, D., Lehman, R. E., Donner, D. B., Matli, M. R., Warren, R. S., and Welton, M. L. (2002). Expression and endocytosis of VEGF and its receptors in human colonic vascular endothelial cells. Am. J. Physiol. Gastrointest. Liver Physiol. 282, G1088–G1096.

Pubmed Abstract | Pubmed Full Text

Wang, S. J., Saadi, W., Lin, F., Minh-Canh Nguyen, C., and Li Jeon, N. (2004). Differential effects of EGF gradient profiles on MDA-MB-231 breast cancer cell chemotaxis. Exp. Cell Res. 300, 180–189.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Wcislo, R., Dzwinel, W., Yuen, D., and Dudek, A. (2009). A 3-D model of tumor progression based on complex automata driven by particle dynamics. J. Mol. Model. 15, 1517–1539.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Xue, C., Friedman, A., and Sen, C. K. (2009). A mathematical model of ischemic cutaneous wounds. Proc. Natl. Acad. Sci. U.S.A. 106, 16782–16787.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Yancopoulos, G. D., Davis, S., Gale, N. W., Rudge, J. S., Wiegand, S. J., and Holash, J. (2000). Vascular-specific growth factors and blood vessel formation. Nature 407, 242–248.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Zheng, X., Koh, G. Y., and Jackson, T. L. (2013). A continuous model of angiogenesis: initiation, extension, and maturation of new blood vessels modulated by vascular endothelial growth factor angiopoietins, platelet-derived growth factor-b, and pericytes. Discrete Continuous Dyn. Syst. Ser. B, 18, 1109–1154.

CrossRef Full Text