Mathematical models of the VEGF receptor and its role in cancer therapy

Tomás Alarcón, Karen M Page


We present an analysis of a stochastic model of the vascular endothelial growth factor (VEGF) receptor. This analysis addresses the contribution of ligand-binding-induced oligomerization, activation of src-homology 2 domain-carrying kinases and receptor internalization in the overall behaviour of the VEGF/VEGF receptor (VEGFR) system. The analysis is based upon a generalization of a Wentzel–Kramers–Brillouin (WKB) approximation of the solution of the corresponding master equation. We predict that tumour-mediated overexpression of VEGFRs in the endothelial cells (ECs) of tumour-engulfed vessels leads to an increased sensitivity of the ECs to low concentrations of VEGF, thus endowing the tumour with increased resistance to anti-angiogenic treatment.


1. Introduction

Angiogenesis is the process whereby new blood vessels are generated from the existing vasculature in response to substances secreted and released by the surrounding tissues. These substances are special types of cytokines called growth factors (GFs). Endothelial cells (ECs) possess surface receptors specific for each of these GFs. There are many such GFs and cell-surface receptors involved in angiogenesis, but there is a particularly important one, vascular endothelial growth factor (VEGF), as the VEGF receptor (VEGFR) is expressed only by ECs. Angiogenesis can occur in a variety of biological settings, both normal and pathological, ranging from wound healing to cancer.

The onset of angiogenesis is controlled by the so-called angiogenic switch. The usual picture of the angiogenic switch is a scale measuring the levels of angiogenic factors and anti-angiogenic substances. When the former are found in excess of the latter, angiogenesis is triggered (Berger & Benjamin 2003).

Here, we argue that this image of the angiogenic switch might be incomplete. The VEGFR is a receptor tyrosine kinase (RTK). Within their cytoplasmic domains, RTKs have regions which, upon phosphorylation, exhibit tyrosine kinase activity. Activation of these regions, however, occurs only upon receptor oligomerization (Helmreich 2001; Alberts et al. 2002). Most GF molecules are multivalent ligands, i.e. one molecule of GF has more than one receptor-binding domain and, therefore, it can engage in binding with as many receptors as it possesses binding regions. Mathematical models of multivalent ligand/multivalent receptor systems have been formulated and analysed (see Lauffenburger & Linderman (1993), Posner et al. (1995), Woolf & Linderman (2004) and references therein). One remarkable property exhibited by these models concerns the behaviour of the response curve (which, roughly speaking, represents the probability of a cell within a population responding to a given concentration of ligand). Whereas the response curve for receptors that do not depend on receptor oligomerization for activation (e.g. G-protein-linked receptors) is monotonic, saturating for high ligand concentration, the response curve for multivalent ligand/multivalent receptors is bell-shaped: cellular responses are inactivated at high concentrations of ligand. With the help of the models presented here, we aim to discuss the implications of this property of RTKs in relation to the onset of angiogenesis and anti-angiogenic therapy.

In spite of the initial enthusiasm raised by anti-angiogenic therapy, the actual results obtained on patients in clinical practice have been poor and its impact on the life expectancy of cancer patients has been very disappointing, in particular, when the anti-angiogenic drugs were used alone (see the review by Jain (2005) and references therein). This lack of results, especially in contrast with the success registered on laboratory animals, has been puzzling. A commonly adopted explanation for such a failure is that, whereas anti-angiogenic drugs can kill many cancer cells, they do not eradicate the tumour completely and the remaining tumour cells will eventually trigger angiogenesis anew (Hampton 2005). One of our aims is to use our models to try to produce plausible explanations of this failure.

Tumour vasculature, whether tumour vessels are the product of tumour-induced angiogenesis or native vessels of the host which have been engulfed by the growing tumour mass, presents many structural abnormalities in comparison to its normal counterpart (Jain 2005). An example of such deregulation is an overexpression of the VEGF surface receptor (Ferrara 2002; Cross et al. 2003). Further evidence for this can be found in experiments carried out on retinal microvascular ECs under stimulation with oestrogen (Suzuma et al. 1999). Oestrogen is known to promote proliferation of some types of breast cancer cells (Amlal et al. 2006) and, therefore, it is plausible that the same mechanism upregulates VEGFR. Moreover, recent experiments by Zhang et al. (2005) show that the platelet-derived growth factor (PDGF) receptor, a system analogous to the VEGFR in every significant aspect, is upregulated in ECs of hepatocellular carcinoma. In this paper, we aim to study the effect of overexpression of surface VEGFR on anti-angiogenic therapy. We consider two possible mechanisms for overexpression of surface VEGFR, namely increased rate of VEGFR synthesis (Suzuma et al. 1999; Zhang et al. 2005) and downregulation of receptor endocytosis (Polo et al. 2004).

The models presented here are formulated in terms of Markov processes and analysed by means of a Wenzel–Kramers–Brillouin (WKB) approximation of the master equation (Kitahara 1973; Kubo et al. 1973). Our models include ligand/receptor binding, ligand-induced receptor dimerization, receptor internalization and binding of enzymes carrying src-homology 2 (SH2) domains (e.g. members of the Src tyrosine kinase family) to activated (dimerized) receptors. This last process constitutes the earliest event in RTK-activation-induced signalling. Our analysis allows us to discern the contribution of each of these processes to the overall behaviour of the VEGFR system as well as to assess the plausible roles of each of them in resistance to anti-angiogenic therapy in solid tumours.

There are several recent studies on models of GF/RTK ligation dynamics. Mac Gabhann & Popel (2004) study a model of competitive binding of VEGF and placental growth factor (PlGF) to VEGFR. There is some evidence of synergy (i.e. enhancement of cell response) between PlGF and VEGF in pathological situations. The mathematical models presented by Mac Gabhann & Popel (2004) help to elucidate the mechanisms of this synergy. Their models take into account receptor internalization but do not account for VEGFR dimerization. Mac Gabhann & Popel (2005a) have studied the system VEGF/VEGF receptor 2 (VEGFR-2)/neuropilin-1 (NRP-1). Both VEGFR-2 and NRP-1 are found on the surface of ECs; they do not interact directly, but can be cross-linked by a VEGF isoform which has binding sites for both. This model considers cross-linking between VEGFR-2 and NRP-1, but does not account for either receptor internalization or VEGFR dimerization.

A model of the PDGF/PDGF receptor (PDGFR) has been proposed by Park et al. (2003). This model incorporates some early events in the signalling cascade triggered by PDGF/PDGFR binding, including phosphoinositide 3-kinase-dependent activation of Akt. The authors also incorporate an alternative model for receptor dimerization in which dimerization is mediated by receptor domains that are active or exposed only when the receptors are bound, forming a sort of ‘pre-dimer’. Ligand and receptor are supposed to associate and dissociate rapidly. The dissociation of one of the ligands from its receptor within a pre-dimer leads to the formation of a stable dimerized complex. The model presented by Park et al. (2003) exhibits good agreement with experimental data.

Mac Gabhann & Popel (2005b) have proposed a stochastic analysis of VEGF binding to cell-surface receptors. In physiological circumstances, VEGF is usually found in very low concentrations, typically of the order of the picomolar. These concentrations imply less than one ligand molecule in each cubic micrometer of fluid. Such low concentrations lead them to consider the validity of the excess of ligand assumption (Sulzer et al. 1996) and the effects of the fluctuations in cellular response, especially in a scenario in which response is threshold triggered. They find agreement between stochastic and deterministic models in the range of VEGF concentrations handled in in vitro experiments (of the order of the nanomolar), but argue that in in vivo situations the effects of fluctuations might be more important.

This paper is organized as follows. Section 2 is devoted to giving details of our model formulation and a brief summary of the necessary biological background. In §3, the stochastic models formulated in §2 are analysed by means of an asymptotic analysis (a generalization to arbitrary dimension of the work by Kubo et al. (1973)). This analysis produces a set of Ordinary Differential Equations (ODEs) for the first and second moments which are then solved numerically. Section 3 also contains details of estimation of parameter values. In §4, we present numerical simulations of the response of our models to anti-angiogenic therapy, assuming both a physiological and a pathological scenario, the latter one characterized by overexpression of surface receptors by inhibition of endocytosis. Another possible source of overexpression of receptors is upregulation of receptor synthesis. This situation is analysed in §5. Section 6 presents an analysis of the fluctuations. Finally, in §7 we summarize and discuss our results.

2. Biological background and model formulation

Here, we briefly summarize the biological background necessary to understand how our models are set up and then we discuss our model formulation.

2.1 Biological background

VEGF denotes a large family of dimeric glycoproteins which consists of five mammalian and one virus-encoded members. VEGF-A was the first member to be discovered and has been shown to be involved in a large number of processes, with both physiological and pathological functions. VEGF-A, in turn, is expressed as four isoforms of different lengths. The shortest of them (VEGF-A121, 121 amino acids long) differs from the other three in its lack of ability to bind to the extracellular matrix and, therefore, it diffuses freely (Hicklin & Ellis 2005).

Regarding the VEGFRs, there are three different types VEGFR-1, -2 and -31. ECs in tumour blood vessels express mostly VEGFR-2, although VEGFR-1 and -3 might also be expressed. In physiological conditions, the vascular endothelium expresses VEGFR-1 and -2, whereas the lymphatic endothelium expresses VEGFR-2 and -3 (Cross et al. 2003). Of the two receptors expressed on ECs, only VEGFR-2 seems to contribute to intracellular signalling, with the function of VEGFR-1 most probably being sequestering (excess) VEGF (Cross et al. 2003).

In order to keep our model as simple as possible and stay focused on the study of how ligand/receptor-binding dynamics affect the early events of the VEGF-binding-induced signalling cascade, we concentrate on the effects of diffusible forms VEGF-A and their binding to VEGFR-2. This particular system appears to make a major contribution to tumour-induced angiogenesis. Thus, hereafter, for simplicity in the notation, the system VEGF-A/VEGFR-2 will be referred simply as VEGF/VEGFR.

In the case of the VEGF/VEGFR system, the ligand (VEGF molecule) is bivalent and the receptor (VEGFR) is monovalent, meaning that one VEGF molecule binds two VEGFRs, while each VEGFR can bind a single VEGF molecule.

This property provides a mechanism for RTK activation, which is elusive from a purely structural perspective: receptors are oligomerized (in the particular case of the VEGFR, dimerized) upon ligand binding. The receptors within the oligomer are brought into close proximity which leads to receptor cross-phosphorylation (Alberts et al. 2002). Cross-phosphorylation yields attachment of phosphate to the tyrosine kinase domains within the cytoplasmic tail of the RTKs, providing high-affinity docking sites for selected substrates to bind. These substrates, usually members of the Src family of tyrosine kinases, carry the SH2 domain, which has high specificity for the phosphorylated domains within the RTKs, and are themselves tyrosine kinases activated by binding to phosphorylated RTKs. These are the earlier events in the signalling cascade triggered by GF/RTK binding. Activated SH2-carrying kinases relay the signal on to other tyrosine kinases, which lead to activation of the corresponding pathways and the alteration of cell behaviour.

Each VEGFR has two kinase domains. We consider that each of these has only one tyrosine residue that is cross-phosphorylated under ligand-induced dimerization, thus providing four high-affinity docking sites for SH2 domains (Cross et al. 2003). Actually, dimerized receptors exhibit more than four possible docking sites (six or more according to Cross et al. 2003). We have made this approximation in order to keep the model as simple as possible. Below, we comment on the effects of this approximation.

According to Sawyer (1998), there are basically two types of SH2-bearing tyrosine kinase: those carrying only one SH2 domain, hereafter to be referred to as SH2 monomers, and those carrying two SH2 domains (e.g. ZAP70 or PI3K). In this paper, only the former ones are considered.

2.2 Model formulation

The stochastic models we analyse in this paper are specified in terms of three quantities, namely the state vector X whose components are the number of molecules of each of the species involved, the probability per unit time corresponding to each of the reactions involved in the process being modelled, Wi, and the corresponding vector ri. The components of ri are the increments in the number of molecules when the ith reaction occurs. To summarize, the occurrence of the ith reaction induces the change in the state vector XX+ri and occurs with probability proportional to Wi. The system is then described by the probability density of the system being in state X at time t, Ψ(X,t), whose dynamics is given by the master equationEmbedded Image(2.1)

Next, we present a description of the three models to be analysed.

2.3 Receptor-binding model

We use a version of the stochastic model for multivalent ligand-induced receptor oligomerization developed by Alarcón & Page (2006). Here, the model corresponds to a bivalent ligand and a univalent receptor, which corresponds to the case of the VEGF/VEGFR system. The stochastic dynamics of this model is summarized in table 1, where the precise form of the transition rates for the different events involved in the ligand–receptor-binding model are given, and depicted in figure 1, where the different reactions involved in the ligand/receptor binding are represented schematically. In figure 1, U is the number of unbound receptors, B is the number of bound receptors and X is the number of dimers (Embedded Image, where Embedded Image is the number of surface receptors). In table 1, uU/N, bB/N and xX/N.2 Here, N=NS+NR is the total number of molecules in our simulation. NS is the number of SH2-carrying enzymes (table 2). L is the concentration (in moles per litre) of free ligand, which is assumed to be constant, i.e. ligand is supplied at a rate that matches its rate of binding to the surface of the cells.

Figure 1

Schematic of our receptor-binding model (table 1) including ligand-induced dimerization (receptor activation). The black rectangles represent the cytoplasmic kinase domains of the receptors. The asterisks denote that the dimerized receptors have undergone cross-phosphorylation and hence provide high-affinity docking sites for SH2 domain-carrying kinases. See text for a definition of Graphic.

View this table:
Table 1

Reaction probability per unit time, WiW(X,ri,t), i=1, … , 4, for the four elementary reaction steps involved in our model. The initials ‘p.u.t.’ stand for per unit time. See table 2 for a summary of parameter values.

View this table:
Table 2

Parameter values for the VEGFR models.

The actual model (i.e. transition rates for each reaction and the corresponding vectors ri) used for receptor binding is summarized in table 1 and figure 1. This model will be hereafter referred to as Model 1. The transition rate corresponding to reaction r3 (Embedded Image in figure 1) needs further clarification (Alarcón & Page 2006). The transition rate for this reaction, which corresponds to the formation of a dimer, is obtained as the product of two factors: the rate of binding between an unbound receptor and a ligand–receptor heterodimer and the probability of finding another receptor within a characteristic distance Δ of the ligand–receptor heterodimer. The latter is given by πΔ2ρ, with Embedded Image the surface density of receptors on the cell surface and R the average radius of an EC.

2.4 Src-homology 2 binding to dimerized receptors

We consider that each VEGFR provides two high-affinity docking sites for tyrosine kinases carrying SH2 domains upon ligand-induced dimerization, thus providing four high-affinity docking sites for SH2 domains.

In figure 2 and table 3, SFNSS stands for the number of free SH2 domains, i.e. those that are not bound to dimerized receptor. S is the number of bound SH2 domains, X is the total number of dimers, X1 is the number of dimers bound to a single SH2 domain, X2 is the number of dimers bound to two SH2 domains, X3 is the number of dimers bound to three SH2 domains, X4 is the number of dimers bound to four SH2 domains and X0 is defined as the number of ‘free’ (i.e. not bound to SH2) receptor dimers. Taking these definitions into account, it is straightforward to see that Embedded Image and Embedded Image.

Figure 2

Schematic of our receptor-binding model (table 3). The black rectangles represent the cytoplasmic kinase domains of the receptors. The asterisks denote that the dimerized receptors have undergone cross-phosphorylation and hence provide high-affinity docking sites for SH2 domain-carrying kinases. See text for a definition of Graphic, i=1, … , 4.

View this table:
Table 3

Reaction probability per unit time, Graphic, i=1, …, 12 and Graphic. This second reaction scheme corresponds to the interaction of dimerized (phosphorylated) receptors with proteins carrying the SH2 domain, e.g. the Src family of tyrosine kinases, which are instrumental in relaying the signal initiated by receptor binding. N refers to the number of receptors plus the number of proteins carrying a SH2 domain and NA is Avogadro's number. We will assume that each pair of dimerized receptors carries a pair of phosphorylated tyrosines to which two SH2 domains can bind. The initials ‘p.u.t.’ stand for per unit time. See table 2 for a summary of parameter values.

The rate constants Embedded Image, i=1, … ,4 that appear in figure 2 need further clarification. According to table 3, Embedded Image, where Embedded Image is the binding rate of a SH2 domain to a phosphotyrosine residue on a receptor dimer (table 2). Since this constant is given in M−1s−1, we need to use the factor VcNA, where Vc is the volume of an EC and NA is Avogadro's number, to convert to appropriate units. The factor 4−(i−1) corresponds to the number of free docking sites that are left on the receptor dimer. We note that we assume that only receptor dimers free of SH2 domains can dissociate.

The model summarized in table 3 will be hereafter referred to as Model 2.

2.5 Endocytosis of surface receptors

We introduce a third element in our model, namely, receptor internalization (Teis & Huber 2003). RTKs undergo clearance from the surface upon ligand-mediated activation (i.e. dimerization) in a very efficient manner. Inactivated receptors (i.e. unbound and non-dimerized bound receptors) also undergo internalization.

Receptor endocytosis is a complex process involving a sophisticated network of protein interactions of which not all the details are known. As we intend to produce a simple model of receptor endocytosis capturing its essential features, we only give a very brief and incomplete summary of the biology of receptor internalization. The reader is referred to Helmreich (2001) and Teis & Huber (2003) for more detailed reviews.

Both inactive (i.e. non-dimerized) RTKs and RTK dimers undergo internalization by essentially the same mechanism.3 The first step is the formation of a structure called a clathrin-coated pit around the RTK dimer. This pit eventually pinches off the membrane forming a vesicle containing the RTK dimer. Once these vesicles are formed, the RTKs enter the so-called early endosome. Although the mechanism of early RTK internalization is the same for both non-dimerized and dimerized receptors, the rate of endocytosis of the latter is much in excess of the former. This indicates that a protein network regulating endocytosis is upregulated upon RTK dimerization (Teis & Huber 2003).

After entering the early endosome, dimer and non-dimer RTKs follow different pathways. Non-dimer RTKs are rapidly recycled to the membrane (we will assume that they do so as unbound receptors). However, most of the dimerized RTKs are transported to the so-called late endosomes, from where they pass into the lysosomes, where they undergo degradation (Teis & Huber 2003).

As the total number of cell-surface receptors seems to stay constant, parallel to endocytosis, RTK production must be sustained by the cell at some given rate (Lauffenburger & Linderman 1993).

Further to these assumptions, we make the hypothesis that somewhere down the endocytic pathway, the SH2-carrying enzymes potentially attached to RTK dimers are detached and released into the cytoplasm. This assumption allows us to ensure that the total number of SH2 domains stays constant. A further simplifying assumption will be that the rate of internalization and degradation for dimerized receptors is independent of the number of SH2 domains bound to their active sites. We also assume that unbound RTKs and non-dimerized ligand/receptor complexes are internalized and recycled back to the membrane at the same rate.

Our stochastic model is inspired by a model by Lauffenburger & Linderman (1993). Our model, however, contains some simplifications with respect to Lauffenburger & Linderman (1993). We will assume that all the internalized RTK dimers go to degradation in the lysosomes (without considering the two intermediate compartments described above) and that only unbound RTKs and non-dimerized ligand/receptor complexes undergo recycling. We further assume that none of these pass into the lysosome and that they are not degraded. Moreover, we assume that the rate at which receptors are synthesized is such that receptor degradation is exactly balanced by receptor synthesis. This assumption is introduced in order to have a constant number of ‘particles’ in our model. This assumption could be relaxed by considering NR as a random variable included in the model rather than as a model parameter.

In table 4, the variables bearing the superindex ‘i’ are the internalized counterparts of the surface variables, which bear no index. The physical meaning of the ‘surface’ variables is the same as in Model 2. For example, X1 is the number of surface dimers bound to a single SH2 domain, whereas Embedded Image is the number of internalized dimers bound to a single SH2 domain. The total number of bound SH2 dimers is now Embedded Image.

View this table:
Table 4

Reaction probability per unit time, Graphic, i=1, … , 12 and Graphic. N refers to the number of receptors plus the number of proteins carrying a SH2 domain and NA is Avogadro's number. We assume that each receptor dimer carries four phosphorylated tyrosines to which two SH2 domains can bind. See table 2 for a summary of parameter values.

The model summarized in table 4 will be hereafter referred to as Model 3.

3. Model analysis: WKB approximation

The methodology we use to analyse the models presented in §2, originally proposed within the field of chemical physics, is due to Kubo et al. (1973) and Van Kampen (1992).4 This technique essentially consists of extending the form of the equilibrium probability density to a non-equilibrium setting. In thermodynamic systems, the equilibrium probability density is given byEmbedded Image(3.1)where C is the normalization constant, X is a set of extensive variables, whose values determine the state of the system, and the function Φe(X) has the properties of a thermodynamic potential, i.e. is a homogeneous function,Embedded Image(3.2)where N is the size of the system which, for example, can correspond to the number of particles. From these two equations, we can see that Ψe(X) is the probability density for fluctuations of the macroscopic extensive variables X with respect to the equilibrium state, as Ψe(X) is a Boltzmann-like function, i.e. the exponential of a homogeneous function which plays the role of thermodynamic potential.

Kubo et al. (1973) have proved that, under the appropriate scaling substitution, the time-dependent solution of the Master Equation (ME) can be approximated by a function of the same form as its equilibrium solution (equation (3.1)), namely the exponential of a homogeneous function, which we call S, of X,Embedded Image(3.3)

Intuitively, the accuracy of this approximation can be assessed in terms of the comparison between the characteristic time-scale associated with the disturbance that drives the system out of equilibrium and the time-scale of the phenomena occurring locally in the system. If the former is much shorter than the latter (which usually happens in weak noise limits), the system may be considered locally in equilibrium. From a more rigorous point of view, Van Kampen (1992) has shown that the extensive variables exhibit the following asymptotic behaviour:Embedded Image(3.4)where 〈X(t)〉 is the solution of the macroscopic equations, i.e. the average value of X and ξ is a Gaussian random variable ∼N(0, 1). This equation reveals that out of equilibrium the fluctuations of an extensive Markovian variable around its mean value depend on the volume in the same way as in equilibrium. Thus, this fact justifies the substitution of equation (3.3).

Let us consider a system whose stochastic dynamics is described by the ME,Embedded Image(3.5)

Kubo et al. (1973) have shown that the transition rates W(X,r,t) must be homogeneous functions of X to obtain a solution of the ME of the form of equation (3.3),Embedded Image(3.6)Accordingly, the probability of a given reaction to occur within an infinitesimal interval of time is proportional to the size of the system and is determined only by the state of the system, represented by the set of intensive variables x. The definitionEmbedded Image(3.7)together with equation (3.6), enables us to write the ME (3.5) in WKB formEmbedded Image(3.8)where we have used that Embedded Image is the generator of the translations in the space of states of the system.

For arbitrary n, the scaling substitution for the cumulants of the probability distribution (for example, (q1)i=〈xi〉; Embedded Image and so on) Embedded Image yields a consistent expansion leading to balanced equations for the cumulants qn(t). Eventually, this leads to the following equation for the leading-order term for q1(t):Embedded Image(3.9)where w(v,r,t) is the Fourier transform of a(x,r,t) and the quantity Embedded Image(q11,t) is defined byEmbedded Image(3.10)

Note that this is the result predicted by the Law of Mass Action. Likewise, the equation for the cumulants of order 2 is given byEmbedded Image(3.11)where Embedded Imageij≡(q21)ij and the first term on the right-hand side has been symmetrized (q21 is a symmetrical matrix). Equations (3.9) and (3.11) are our final result and constitute the generalization to arbitrary dimension of the results obtained by Kubo et al. (1973).5 A detailed proof of these results is given in the electronic supplementary material. This electronic supplementary material also includes an alternative derivation using stochastic calculus rather than asymptotic methods.

3.1 Evolution equations for the VEGFR model

The results obtained in §3 are valid, in general, as long as the transition rates in the ME fulfil the homogeneity condition stated by equation (3.6). In this section, we apply these results to the particular case of Model 3 described in §2, i.e. we use equations (3.9) and (3.11) to formulate the systems of ODEs for the leading-order contributions to the first and second cumulants (i.e. the first and second moments, respectively). Model 3 is the most general of the three described in §2. Models 1 and 2 can be obtained as particular cases of Model 3 as detailed below.

The conservation laws Embedded Image and SF=nssB (sF refers to the fraction of unbound SH2 domains) are used and the quantities NNR+NS, nRNR/N, nsNs/N and Embedded Image have been defined. L stands for the concentration of ligand, in this particular case VEGF. We assume the so-called excess of ligand regime, namely the concentration of (free) ligand is not affected by binding (Sulzer et al. 1996); N(t=0)=2.3×105 (table 2). The evolution equations corresponding to Model 2 can be obtained from equations (3.12)–(3.25) by setting Embedded Image, Embedded Image, kre=0, kd=0 and ks=0. The equations for Model 1 are equations (3.12) and (3.13) with Embedded Image, kre=0 and ks=0.

By substituting the corresponding values of a(x,r,t) and r from Model 3 (table 1) into equation (3.9) we obtainEmbedded Image(3.12)

Embedded Image(3.13)

Embedded Image(3.14)

Embedded Image(3.15)

Embedded Image(3.16)

Embedded Image(3.17)

Embedded Image(3.18)

Embedded Image(3.19)

Embedded Image(3.20)

Embedded Image(3.21)

Embedded Image(3.22)

Embedded Image(3.23)

Embedded Image(3.24)

Embedded Image(3.25)

Likewise, a system of ODEs can be written for q21(t). This quantity is a symmetric 13×136 matrix and, therefore, has 91 independent components. Hence, the corresponding ODE system has 91 equations. Furthermore, this system of 91 ODEs is coupled to equations (3.12)–(3.25).7 As a result, a full analysis of the fluctuations for the proposed model implies a system of 104 ODEs. In general, for a system with dimension d, the system of ODEs determining the dynamics of the first- and second-order cumulants at leading order has d(d+3)/2 equations. In fact, the most serious shortcoming of the method presented here is that the size of the resulting system of ODEs makes the analysis painstaking, even for modestly complex models like the one given in table 3. With some degree of uncertainty in the parameter values and 104 equations, further simplification seems necessary. Consequently, only the behaviour of the mean value of the full model of table 3 (equations (3.12)–(3.25)) is analysed. However, if we restrict ourselves to the receptor model described in table 1, we obtain a system that we can easily handle.

Using equation (3.11) and table 1, the equations for the fluctuations of u and b corresponding to the receptor model read asEmbedded Image(3.26)Embedded Image(3.27)Embedded Image(3.28)where Embedded Imageij=Embedded Imageji. Q11 and Q22 correspond to the variance of u and b, respectively. These equations are to be solved together with equations (4.2) and (4.3) with Embedded Image, kre=0 and ks=0. Using the conservation law Embedded Image, we obtain the following expression for the variance of x, Q33, in terms of the dependent variables of equations (3.26)–(3.28):Embedded Image

It is important to bear in mind that, while the dynamics of the mean value of the variables corresponding to the receptor model (u, b and x) is unaffected by the dynamics of the intracellular SH2 domains, the dynamics of the fluctuations (q21) of the receptor variables depends on fluctuations and mean values of SH2 variables. Therefore, whereas our restricted analysis may provide useful insights, its conclusions should not be applied to the dynamics of the fluctuations of the whole system.

3.2 Parameter values

Many of the parameter values used in our simulations are available from the literature (see table 2 for a summary). Some of them, however, could not be found in the literature. Such is the case for the parameters related to cross-linking (i.e. dimer formation). We proceed by fixing the value of Embedded Image and then fitting the value of Embedded Image to obtain model results that are consistent with experiments; in particular, we will try to reproduce the detailed experimental data provided by Park et al. (2003) regarding the short-time behaviour of the PDGF/PDGFR system, which is very similar to the VEGF/VEGFR system, in fibroblasts.

There are three aspects of the experimental results of Park et al. (2003) against which we benchmark our model. For ligand concentrations between 10 and 0.01 nM, the number of phosphorylated receptors (which we compare to the number of surface dimers) achieves a maximum within 20 min (the peak of VEGFR phosphorylation is typically reached within 5–10 min of stimulation with VEGF; Mac Gabhann & Popel 2005b). Moreover, there is virtually no response to ligand concentrations below 0.01 nM. Finally, the peak level of receptor activation as a function of ligand dose can be fitted to a Hill curve with Hill exponent n>1, thus exhibiting positive cooperativity.

We have found that for all of these three conditions to be satisfied by Model 3, as shown in figures 3 and 4, for Embedded Image, Embedded Image has to be greater than or equal to 4.6×103 s−1. For smaller values, the last condition is not satisfied, as high-dose inhibition of receptor activation is observed within the aforementioned interval of ligand concentration (results not shown). Thus, values of Embedded Image and Embedded Image will be used in our simulations (table 2). Smaller values of Embedded Image yield a response that is too slow compared with the experimental estimates. According to Mac Gabhann & Popel (2005b), the VEGFR is activated within 5 min of stimulation with VEGF. Values of Embedded Image produce the opposite effect: activation of the VEGFR occurs on an unrealistically short time-scale.

Figure 3

 Simulation results corresponding to Model 3 with Graphic and Graphic, respectively. This plot shows the time course of the proportion of surface dimers for different values of L. Key: solid line corresponds to L=10 nM, dashed line to L=1 nM, dot-dashed line to L=0.1 nM and dotted line to L=0.01 nM. Ligand is introduced at t=0.

Figure 4

(a) Simulation results corresponding to Model 3 with Graphic and Graphic. The squares in this plot show the maximum value achieved by the proportion of surface dimers (figure 3) as a function of ligand concentration. Solid line corresponds to a Hill curve Graphic with n=1.2. (b) Experimental results by Park et al. (2003). Solid line corresponds to a Hill curve Graphic with n=1.5. In both panels, dashed lines correspond to a Hill curve with n=1, which has been included for comparison. Ligand is introduced at t=0.

The value of the parameter Δ, i.e. the radius of the ‘cross-linking surface’ (Alarcón & Page 2006), has been estimated from the average diameter of a typical protein (5 nm) and set to be Δ=2.5 nm.

In order to estimate the number of SH2 domains, we use the value of the SH2 domain affinity to phosphorylated domains within a dimerized receptor, As, provided by Helmreich (2001). If As=100 nM, it seems sensible to assume that this is similar to the concentration of SH2 domains in physiological conditions. A quick calculation shows that, for the value of Vc provided and assuming that the concentration of SH2 domains, Cs, is 100 nM, the number of SH2 domains is of the order of Ns=180 000 (number of domains=NA×Cs×Vc, where NA is Avogadro's number and Cs is the concentration given in moles per litre). In our simulations, NS will be taken as the total number of SH2 monomers.

We consider that each of the kinase domains within a receptor dimer has only one tyrosine residue that is cross-phosphorylated under ligand-induced dimerization, thus providing four high-affinity docking sites for SH2 domains (Cross et al. 2003). Actually, this number is thought to be larger: six or more according to Cross et al. (2003). We have performed simulations with up to eight docking sites (results not shown) and the results do not change substantially. The reason for this can be seen in figure 5. We can see that the proportion of receptor dimers bound to SH2 domains consistently decreases with the number of SH2 domains bound to it. Figure 5 shows that the proportion of receptor dimers with four bound SH2 domains is two orders of magnitude smaller than the proportion of receptor dimers with only one bound SH2 domain. Thus, adding more phosphorylation sites only adds more complexity to the model without providing extra insight.

Figure 5

Simulation results corresponding to Model 3. These plots show the time course of the proportion of surface dimers bound to SH2 domains (x1, x2, x3 and x4) for L=10−8 M.

3.2.1 Non-dimensionalization

To perform our simulations, we have re-written our models in terms of dimensionless variables. We have defined a dimensionless time Embedded Image and, accordingly, all the rate constants in tables 1, 3 and 4 can be expressed in non-dimensional terms by dividing them by koff. Additionally, we have defined A=kon/koff such that the quantity AL is dimensionless.

4. Overexpression of surface VEGFR: effects on anti-angiogenic therapy

Here, we analyse and simulate the evolution equations obtained for Models 1–3. The effects of anti-angiogenic therapy (in the form of an anti-VEGF drug) are simulated for Models 2 and 3. For the reasons pointed out in §1, we restrict our analysis of Models 2 and 3 to the equations for the first moment.

Model 1, being the simplest model discussed in this paper, allows us to analyse the corresponding system for both the first and the (central) second moments.

Most of the simulations presented next concern the effect of anti-angiogenic therapy, in particular, an anti-VEGF drug. In practice, such drugs are specific antibodies against VEGF that bind and inactivate VEGF molecules. We are not going to give a detailed analysis of the VEGF/anti-VEGF drug interactions, rather we implement this therapy in our model by reduction in the concentration of VEGF, L.

4.1 Response of Model 3 to anti-VEGF treatment

We have performed simulations of the effect of anti-VEGF on the dynamics of Model 3. Model 3 together with the set of parameter values shown in table 2 will be taken as our base-line or physiological scenario. Our first step is to show the effect of an anti-VEGF drug on Model 3 (i.e. on a physiological setting). Our results are shown in figure 6 and table 5, where we show two indicators of the effect of the drug: the peak value of active receptors (i.e. surface receptor dimers) and the integral of the number of active SH2-carrying domains (i.e. those bound to surface receptor dimers; Park et al. 2003),Embedded Image(4.1)where Embedded Image denotes the duration of the experiment or measurement.

Figure 6

Simulation results corresponding to Model 3. (a) corresponds to results showing the proportion of surface dimers (upper plot) and the proportion of bound SH2 domains (lower plot) without anti-VEGF treatment for L=10 nM. (b) Results from a simulation in which anti-VEGF therapy was implemented by reducing the level of VEGF from L=10 nM to L=0.01 nM at time t=10−2. (c): idem at t=10−3. Parameter values have been taken from table 2.

View this table:
Table 5

Graphic for Model 3.

We observe that the receptor dimerization peak is very sensitive to the schedule of application of the drug: application of the anti-VEGF drug at t=10−2 (where t=0 corresponds to the time at which ligand is introduced) bears virtually no effect. If the drug is applied earlier at t=10−3, the effect on the activation peak is much more dramatic (figure 6). According to table 5, Embedded Image is less sensitive to the schedule, but there still seems to be some optimal time for drug administration.

4.2 Overexpression of surface receptors by inhibition of endocytosis: Model 2

Our main concern in this work is to provide plausible mechanisms for resistance to (or poor performance of) anti-angiogenic therapy. There is growing experimental evidence supporting scenarios in which the normal functions of ECs are disturbed upon engulfment by a tumour mass (Holash et al. 1999; Ferrara 2002; Jain 2005; Zhang et al. 2005). There is evidence of overexpression of VEGFR and PDFGR in tumour vessels (Ferrara 2002; Zhang et al. 2005) and extra evidence is provided by experiments on retinal ECs stimulated with oestrogen (Suzuma et al. 1999).

In the light of this evidence, we analyse the behaviour of our models upon overexpression of surface receptors. Such overexpression can be the consequence of two different processes. One way of achieving larger numbers of surface receptors is by upregulation of receptor synthesis (Suzuma et al. 1999; Zhang et al. 2005). The incorporation of this effect in our model requires some caveats regarding our model formulation and will be tackled in §5.

Another way of achieving an increase in the number of surface receptors is by inhibition of receptor endocytosis. Polo et al. (2004) present evidence of such a mechanism in cancer cells in which endocytosis of epithelial growth factor receptor was inhibited in a Src-dependent manner upon receptor binding. Within the context of our models, Model 2 can be understood as the non-endocytosis limit of Model 3. Thus, we study Model 2 as a limiting case of endocytosis inhibition in Model 3.

The equations for the first moment of Model 2, as obtained from equations (3.12)–(3.25), as prescribed in §3.1, are given byEmbedded Image(4.2)

Embedded Image(4.3)

Embedded Image(4.4)

Embedded Image(4.5)

Embedded Image(4.6)

Embedded Image(4.7)

Equations (4.2) and (4.3) are decoupled from the rest of the system and, therefore, can be studied in isolation. Using the dimensionless quantities (§3.2.1), they can be rewritten asEmbedded ImageEmbedded Image(4.8)where we have dropped the hats for the dimensionless parameters. It has been shown that in models of receptor binding with cross-linking (i.e. receptor oligomerization), the cellular response to the presence of the ligand is negligible for both very low and very high concentrations of ligand (Lauffenburger & Linderman 1993; Sulzer et al. 1996). In fact, from the steady-state equations corresponding to equations (4.8),Embedded Image(4.9)

Embedded Image(4.10)

We can see that if AL≫1, then u*=O(ϵ) and b*=nRO(ϵ) with ϵ=(AL)−1. If AL≪1, u*=nRO(ϵ) and b*=O(ϵ) with ϵ=AL. In both cases, the number of dimerized receptors is so small that it is impossible for the cell to produce a response. equations (4.9) and (4.10) yield the following solution for the steady state of the receptor model:Embedded Image(4.11)

Embedded Image(4.12)

Embedded Image(4.13)

Figure 7 represents x*(AL) and shows that there is an intermediate range of values of the quantity AL such that the number of dimerized receptors satisfies Embedded Image. It follows that u*=O(ϵ) and b*=O(ϵ). In fact, x* is an even function of log(AL).

Figure 7

x* as a function of the dimensionless quantity AL for different values of the dimensionless quantity Graphic. Parameter values have been taken from table 2. Squares, Graphic; triangles, Graphic; and circles, Graphic.

In spite of the fact that the steady-state response of the system is symmetric with respect to log(AL) and quite insensitive to its value over an interval of values of log(AL) around log(AL)=0, the dynamical behaviour of the model strongly depends on the value of AL. It is straightforward to see that the relaxation time, τ, i.e. the time the system takes to settle down to a close-to-equilibrium state, is essentially determined by AL, the relaxation process being faster for larger values of AL. Two illustrative examples of this are shown in figure 8a,b for log(AL)=−3 and 1, respectively. Numerical solution of the systems of ODEs equations (4.2)–(4.7) has been obtained using Matlab's stiff solver ode23s.

Figure 8

Numerical solution for the first cumulant equations corresponding to Model 2 (table 3). (a) corresponds to log(AL)=−6, panel (b) to log(AL)=1 and (c) to log(AL)=6. Parameter values have been taken from table 2.

From these numerical results and the model equations, we can see that Embedded Image. This implies that, in the pathological situation, this model is aimed to reproduce, the stationary response of the system for L=10−7 M has the same intensity as that for L=10−13 M, but whereas in the former case this response built-up in a time of the order of Embedded Image, in the latter the time required is Embedded Image, which in dimensional terms corresponds to 0.0167 and 16 667 min, respectively. The reader should note that values of L=10−7 M (100 nM) for the VEGF concentration are possibly unrealistic as such high concentrations are unlikely to be found in either physiological or pathological circumstances. We are only using these extreme values here to illustrate our point.

An observation with respect to the physiological situation described by Model 3 is that blocking receptor endocytosis appears to induce an important change in the dynamical behaviour of the system. While Model 3 supports a scenario in which there is a fast, transient (a peak in receptor) activation followed by a decay to a stationary state, Model 2 suggests that, in the pathological setting, the response is slow and sustained.

This observation may have deep therapeutic significance. We have run simulations of anti-VEGF therapy to compare the responses of the physiological (Model 3) with the pathological (Model 2) situations. As in §4.1, we have performed simulations in which at some point during the evolution of the system the concentration of VEGF has been reduced from L=10−8 M (10 nM) to L=10−11 M (0.01 nM). The results are shown in figure 9 and table 6.

Figure 9

Simulation results corresponding to Model 2. (a) corresponds to results showing the proportion of surface dimers (upper plot) and the proportion of bound SH2 domains (lower plot) without anti-VEGF treatment for L=10 nM. (b) Results from a simulation in which anti-VEGF therapy was implemented by reducing the level of VEGF from L=10 to 0.01 nM at time t=10−3. (c): idem at t=104. Parameter values have been taken from table 2.

View this table:
Table 6

Graphic for Model 2.

From these results, it is clear that this therapeutic intervention, in spite of reducing the concentration of VEGF by three orders of magnitude, has had virtually no effect on the peak or steady state of the surface dimers or bound SH2 in the pathological system. Moreover, comparing figure 6 and table 5 to figure 9 and table 6, respectively, we see that the anti-VEGF therapy has done much worse in the pathological than in the physiological case. A far more efficient clearance of active VEGF by the anti-VEGF drug is needed in order to produce better outcomes (see fifth entry in table 6).

Furthermore, as the dynamics of the system has been slowed down by three orders of magnitude (from time-scales of the order of magnitude of minutes to hours or days) by the reduction of ligand concentration, even if a decrease in the angiogenic activity is initially observed, angiogenic activity could resume simply because VEGF clearance by the drug was not effective enough and pass undetected by the clinicians depending on how the follow-up is done.

These results, thus, point towards overexpression of surface VEGFR by inhibition of endocytosis as a possible factor related to resistance to anti-angiogenic therapy. The fact that the resulting stationary response curve is symmetric with respect to log(AL) means that the system can afford a reduction of several orders of magnitude in the concentration of VEGF without producing a significant effect on the long-term behaviour, as the system slowly relaxes to a steady state of very similar ‘stationary’ angiogenic potency. Only if the anti-VEGF is efficient enough to drive the system out of the responsive region, will the treatment produce significant results.

From this discussion, we can infer that the therapeutic outcome of an anti-VEGF drug could be improved by some strategy aimed to narrow the width of the bell-shaped stationary response curve. Figure 7 shows that decreasing the value of the dimerization rate Embedded Image yields such an effect. Simulation results of simultaneously reducing L and Embedded Image are shown in figure 10 and table 6 (see sixth entry). We can see a substantial improvement in anti-VEGF performance by comparing the sixth and the fourth entries of table 6, corresponding to anti-VEGF with and without Embedded Image reduction, respectively.

Figure 10

Simulation results corresponding to Model 2 where the effect of combining anti-VEGF therapy and reduction of Graphic is demonstrated. Both therapies have been applied at t=104. VEGF concentration has been reduced from L=10 to 0.01 nM. The dimerization rate has been reduced from Graphic to Graphic. The rest of the parameter values have been taken from table 2.

5. Overexpression of surface receptors by upregulation of receptor synthesis

So far, we have analysed Model 2 as a limiting case of Model 3 in which endocytosis is downregulated, thus providing a model for pathological overexpression of surface receptors by endocytosis inhibition. The case in which such overexpression is due to upregulation of receptor synthesis is, strictly speaking, beyond the scope of our model formulation as detailed in §2. The reason for this is that a requirement of our stochastic model formulation is that the total number of particles in the model should be constant. When formulating Model 3, we have introduced the additional requirement that the rate of receptor degradation must be precisely balanced by receptor synthesis. This model is, therefore, not appropriate to analyse the effects of upregulation of receptor synthesis.

Here, using the Law of Mass Action, we formulate a generalization of equations (3.12)–(3.25) that explicitly accounts for receptor synthesis,Embedded Image(5.1)Embedded Image(5.2)Embedded Image(5.3)Embedded Image(5.4)

Embedded Image(5.5)Embedded Image(5.6)Embedded Image(5.7)Embedded Image(5.8)Embedded Image(5.9)Embedded Image(5.10)Embedded Image(5.11)Embedded Image(5.12)

Embedded Image(5.13)Embedded Image(5.14)Embedded Image(5.15)where ks is the rate of receptor synthesis and y stands for the total proportion of receptors. The proportions are still measured with respect to N=NR+NS, where NR is now the initial number of receptors, i.e. Embedded Image. This system of ODEs will be referred to as Model 4.

The physiological value of the constant ks, i.e. the base-line case scenario which we will consider as modelling normal EC function, has been chosen, fixing all the other parameter values as in table 2, to obtain values of the total number of receptors of the order of 104, which seems to be a reasonable estimate. Figure 11 shows that ks=9×10−5 s−1 satisfies this requirement for a range of biologically realistic values of the concentration of VEGF.

Figure 11

Simulation results for Model 4 (equations (5.1)–(5.15)) for a physiological situation. This plot shows the proportion of receptors. Dashed line corresponds to L=10−8 M, solid line to L=10−9 M, dot-dashed line to L=10−10 M and dotted line to L=10−11 M. Parameter values have been taken from table 2.

Simulations of Model 4 with parameter values taken from table 2 for different values of the ligand concentration L show an interesting property of the model: it exhibits perfect adaptation. Figure 12 shows that for a wide range of values of L, both the number of surface dimers (not shown) and, consequently, the number of bound SH2 domains relax to the same stationary value regardless of the value of L. The transient, of course, depends on the particular value of L.

Figure 12

Simulation results for Model 4 (equations (5.1)–(5.15)) for a physiological situation. This plot shows the number of SH2 domains bound to receptor dimers on the surface of the cell. Solid line corresponds to L=10−8 M, dashed line to L=10−9 M, dot-dashed line to L=10−10 M and dotted line to L=10−18 M. Parameter values have been taken from table 2.

This property is typical of chemotactic systems which respond to an abrupt change in the concentration of a chemoattractant substance but then adapt to some sort of background state, i.e. the steady state (see Tyson et al. (2003) and references therein). In fact, this property fits into the general picture of the cell responses triggered by VEGF or PDGF, as both substances act as chemoattractants. ECs are known to respond chemotactically to VEGF as part of the angiogenic process (Terranova et al. 1985). PDGF, in turn, is known to act as a chemoattractant for some normal cells (Grotendorst et al. 1982) and some cancer cell lines (Klominek et al. 1998).

The perfect adaptation property exhibited by Model 4 appears to fit well with the picture of the physiological response to VEGF, i.e. a transient activation followed by a relaxation, and consequent (self-)deactivation of the cell response, towards the steady state. The same behaviour was observed in Model 3.

Increasing the rate of receptor synthesis in order to reproduce the pathological situation described by Ferrara (2002) and Zhang et al. (2005) yield several interesting results, which can be observed in figure 13. The perfect adaptation behaviour is robust to an increase in the value of ks, with the difference that the steady-state value of active dimers and, consequently, of bound SH2 domains increases with ks. This, in turn, leads to a dynamical behaviour in which the initial, transient activation tends to disappear (in the sense that the relative height of the peak with respect to the steady-state value is much smaller), in favour of a slow, sustained response much like the one observed in Model 2. The difference with respect to the behaviour of Model 2 is that the slowing down shown by Model 4 is not as dramatic as the one exhibited by Model 2.

Figure 13

Simulation results for Model 4 (equations (5.1)–(5.15)). This plot shows the proportion of SH2 domains bound to surface receptors dimers. Solid line corresponds to L=10−8 M, dashed line to L=10−9 M, dot-dashed line to L=10−10 M and dotted line L=10−18 M. (a) corresponds to ks=4.5×10−4 and (b) to ks=9×10−4. Other parameter values have been taken from table 2.

5.1 Response of Model 4 to anti-VEGF therapy

5.1.1 Physiological case

We first analyse the effect of anti-VEGF on Model 4 in physiological conditions (i.e. with parameter values as given in table 2). The results are qualitatively similar to those obtained for Model 3. We observe a strong dependence of the efficacy of the therapy on the time at which the anti-VEGF is administered, especially regarding the activity peak (figure 14). Such dependence is also observed in Embedded Image, but to a lesser extent (see table 7, entries 1–3).

Figure 14

Simulation results of the effect of anti-VEGF therapy on Model 4 (equations (5.1)–(5.15)) in a physiological situation. The levels of VEGF have been reduced from L=10 to 0.01 nM. (a) corresponds to drug administration at t=10−4 and (b) to administration time t=10−2. Parameter values have been taken from table 2.

View this table:
Table 7

Graphic for Model 4.

5.1.2 Pathological case

Our investigation of the response of Model 4 with upregulated receptor synthesis to anti-VEGF therapy also yields similar results to those obtained using Model 2.

As in the physiological case (base-line rate of receptor synthesis), early administration of the drug leads to an efficient suppression of the activation peak (figure 15). However, owing to the perfect adaptation behaviour exhibited by Model 4, the system recovers to the same steady state, regardless of treatment with the anti-VEGF drug. This has an important implication. Owing to the increased rate of receptor synthesis, the steady-state activation level is higher than in the physiological case (compare figures 12 and 13). In fact, the steady-state activation level corresponding to ks=4.5×10−4 (pathological) is not much smaller than the peak activation level for L=1 nM in the physiological case (ks=9×10−5). This implies that, regardless of the suppression of the activity peak in response to anti-VEGF treatment, the angiogenic response may not be inhibited. The results summarized in table 7 regarding the integrated response Embedded Image lead to the same conclusion: the values obtained for Embedded Image in the pathological case, even in the presence of treatment (entries 5–12), are substantially bigger than those obtained for the physiological case (entries 1–4).

Figure 15

Simulation results of the effect of anti-VEGF therapy on Model 4 (equations (5.1)–(5.15)) in a pathological situation. The value of the rate of receptor synthesis has been increased to ks=4.5×10−4, i.e. fivefold its physiological value. The levels of VEGF have been reduced from L=10 to 0.01 nM. (a) corresponds to an untreated case and (b) to drug administration at t=10−4. Other parameter values have been taken from table 2.

The response of Model 4 to anti-VEGF treatment at later stages in the evolution of the system is shown in figure 16. Owing to the perfect adaptation property of this model, the treatment has no substantial effect other than a short, transient decrease in activity.

Figure 16

Simulation results of the effect of anti-VEGF therapy on Model 4 (equations (5.1)–(5.15)) in a physiological situation. The levels of VEGF have been reduced from L=10 to 0.01 nM. (a) corresponds to drug administration at t=104 and (b) to a close-up of (a). Parameter values have been taken from table 2.

6. Effects of fluctuations for Model 1

Model 1, i.e. the model that describes only ligand/receptor binding and receptor dimerization, is the simplest model proposed in this paper and as such can be thoroughly studied, including the behaviour of the fluctuations, without having to carry out painstaking computations involving hundreds of equations. Here, we carry out such an analysis together with a series of caveats around a commonly sustained notion, namely that the cellular response triggered by a bivalent ligand binding monovalent or bivalent receptor is symmetric with respect to log(AL) (Lauffenburger & Linderman 1993; Sulzer et al. 1996).

Regarding the latter issue, we remark that this is a steady-state result and the dynamical behaviour of the system is quite different depending on whether log(AL)>0 or log(AL)<0. Whereas the steady state of the equations for the first moments is independent of the sign of log(AL), the dynamics of the response to ligand binding strongly depends on the sign of log(AL): the time-scale on which the first moments of the model variables reach the steady state is O((AL)−1). Therefore, the dynamics is much slower for log(AL)<0.

A second point we would like to raise is that the models on which the claims of a symmetric response are based do not include receptor internalization (Sulzer et al. 1996). From figure 4, we can see that the peak in receptor activity (dimerization) when endocytosis is taken in to account is moved towards values of AL>1.

The third factor we think is missing in the usual picture of the cellular response to receptor dimerization is the presence of fluctuations. This is an issue which, in general, is difficult to analyse. Usually, the only way forward is performing simulations of the stochastic dynamics which is computationally intensive. In the particular case of Model 1, which has only two (independent) variables, significant analytical progress can be made using the results obtained using the WKB approximation. In particular, we have analysed the steady state of the system of ODEs consisting of equations (4.2), (4.3) and (3.26)–(3.28). The system of equations corresponding to the steady state has been solved in Matlab using an iterative method (GMRES). The results for Embedded Image are shown in figure 17.

Figure 17

Variance of x for Model 1 (table 1) calculated at the steady state shown on a logarithmic scale. Circles correspond to Graphic and triangles to Graphic. Parameter values have been taken from table 2.

Figure 17 shows that, actually, the steady-state fluctuations around the average value are non-symmetrical. However, given that the total number of particles is NR=50 000, the magnitude of the fluctuations in the present case is very small compared with the steady-state average.

This means that, in the present model, the fluctuations, in spite of not having the symmetry properties of their first moment counterparts, are not expected to have a major effect on the behaviour of the system, in particular, they are not likely to seriously affect the parity of x* with respect to log(AL). However, in other cases where the fluctuations play a more relevant role, their asymmetry could compromise the prediction of the Law of Mass Action models regarding the parity properties of the cellular response.

7. Discussion and conclusions

We have analysed a number of stochastic models describing several parts of the initiation of the signalling cascade triggered upon VEGF binding to a VEGFR molecule. These models include descriptions of VEGF/VEGFR binding, VEGFR dimerization, VEGFR endocytosis and early signalling events (i.e. activation of enzymes carrying SH2 domains). These phenomena correspond to the early steps in the angiogenic process and, therefore, a thorough understanding of all these processes and the interactions between them might help to improve existent therapies targeting angiogenesis. We have also formulated and analysed a deterministic model which allows us to study the effects of upregulation of receptor synthesis that appears to arise in tumour ECs. Our main aim is to produce plausible mechanisms for resistance to anti-angiogenic treatment.

The models analysed in this paper are formulated in terms of a Markov process and, therefore, mathematically described in terms of the corresponding master equation. This equation becomes difficult to solve for models of moderate complexity, like the ones considered here. However, for some systems, including those involving chemical reaction-like kinetics, an approximate solution can be found using the WKB method when the size of the system is large enough. This result, first proved by Kubo et al. (1973), has been extended and adapted to deal with our models. The result is a series of systems of ODEs for the cumulants of order n, qn (with (q1)i=〈xi〉, Embedded Image and so on) that can be analysed using standard methods. The equations for (q1)i=〈xi〉, i.e. for the system's average behaviour, correspond to the Law of Mass Action.

Using this method, we have analysed our base-line or physiological model (Model 3). Unknown parameter values have been estimated by benchmarking the model with experimental results on dose-dependent studies of the PDGFR, which belongs to the same family as the VEGFR and is equivalent to the VEGFR in almost every significant aspect. Model 3, for physiologically relevant VEGF concentrations, exhibits a dynamical behaviour in which an initial transient peak in activation is observed and then a decay to some steady-state value.

Once this physiological scenario has been fixed, we have analysed the effects of overexpression of surface receptors, including its potential effects on response to anti-VEGF treatment. We have first studied the effect of receptor overexpression by inhibition of endocytosis (Model 2).

A second way of achieving overexpression of surface receptors is by upregulation of receptor synthesis. This scenario cannot be studied within the framework of our stochastic models, as our model formulation requires a system with a constant number of particles. Instead, we have formulated a deterministic model based on the Law of Mass Action that allows us to study this situation (Model 4). The physiological scenario, characterized by the value of the rate of receptor synthesis corresponding to this model, yields a dynamical behaviour similar to the one exhibited by Model 3: an initial transient activation followed by a relaxation to a steady state. However, there is an important difference with respect to Model 3, namely Model 4 shows perfect adaptation to the VEGF concentration. Within the framework of Model 4, the pathological case is characterized by an increased rate of receptor synthesis.

The main result of our analysis is that both mechanisms of overexpression of surface receptors lead to a substantially increased resistance to treatment with an anti-VEGF drug. In both cases, the dynamical mechanism appears to be similar: the transient activation exhibited by the physiological case is replaced by a slower and more sustained response. Moreover, in both cases, there is an increased sensitivity to low values of the concentration of VEGF. Model 2 exhibits close-to-full activation for concentrations as low as 10−5 nM (figure 7), whereas physiological activation occurs in the proximity of 1 nM (see figure 4b and Park et al. (2003)).

In the case of Model 4, increasing the rate of receptor synthesis to pathological levels leads to a larger steady-state activation value than the one observed for physiological receptor synthesis, which means that the angiogenic response may not be shutdown after the initial transient is over. In addition, the fact that this system exhibits perfect adaptation with respect to ligand concentration makes this system very sensitive to low VEGF concentrations and also extraordinarily resilient to anti-VEGF treatment. This property may sound unrealistic and could have been introduced in the model by some over-simplistic hypothesis (see below). We have to remark, though, that this increased sensitivity appears in a pathological situation, i.e. in a regime in which the system was never meant to operate. Within the physiological regime, our model produces results that are compatible with experimental results. Furthermore, we recall that perfect adaptation is a property typical of chemotactic systems and that both VEGF and PDGF are chemoattractants for some cell types, chemotactic migration up VEGF gradients being instrumental for successful angiogenesis.

In fact, Model 4 exhibits perfect adaptation as a consequence of one of our model formulation hypotheses, namely internalized receptor dimers are degraded at a much higher rate than internalized receptor monomers. Relaxing this hypothesis leads to a ligand concentration-dependent steady state (result not shown). Further exploration of this issue, i.e. how varying the relative values of the degradation rates can lead to different types of cellular response, will be the subject of future work.

Model 1, being the simplest of the models studied in the present paper, can be more thoroughly analysed. In particular, the resulting system of ODEs for both first and second moments has a reasonable number of equations, which allows us to study the effect of fluctuations on the cellular response. We have shown that the steady-state fluctuations of x are not symmetric with respect to log(AL), although in this case, owing to the size of the fluctuations relative to the average value, this is not expected to affect the predictions made by models formulated in terms of the Law of Mass Action. In other cases, where fluctuations are relevant, this might compromise the predictions of the deterministic models regarding the behaviour of the cell response with respect to the ligand concentration.

A feature that Models 2–4 inherit from Model 1 is the inhibition of the cell response for high ligand concentrations. Although the experimental evidence for such behaviour may not be extensive, the work by Cai et al. (2006) appears to point in that direction. Figure 18b shows their experimental results (data extracted from figure 5c of Cai et al. (2006)), in which the proliferative activity induced by VEGF on retinal microvascular ECs is measured with respect to the activity of unstimulated cells. Figure 18a shows the peak of activated (dimerized) surface receptors as a function of L for Models 3 and 4 (in physiological conditions). We can see that both experiments and theory point to an inhibition of the cellular response for high ligand concentration. However, our models do not predict the bimodal response curve found by Cai et al. (2006).

Figure 18

(a) Simulation results corresponding to Models 3 and 4 (with ks=9×10−5 s−1). The squares (Model 3) and triangles (Model 4) in this plot show the maximum values achieved by the proportion of surface dimers (figure 3) as a function of ligand concentration when the models were simulated until t=1.2, i.e. 20 min in dimensional terms. (b) Experimental results by Cai et al. (2006).

Our main aim in this work is, specifically, to study the properties of the VEGFR system in relation to issues regarding resistance to anti-angiogenic therapy. In spite of this, our models may have wider application. They should certainly be applicable to receptors within the PDGF-like family, to which the VEGFR itself belongs. However, with the necessary caveats on issues such as parameter values, our models should be generic enough to be applicable to the study of other families of RTK, as the ingredients we have included in our models (i.e. receptor dimerization, binding of SH2-carrying enzymes and receptor endocytosis) are generic elements of the function of all the RTKs. In spite of this, we remark that application of our model, for example, to the fibroblast growth factor receptor may require deeper changes to the model. Fibroblast growth factor molecules are actually monomers that first bind heparan sulphate proteoglycans, forming multimers, thus allowing receptor cross-linking. We should extend our model to account for ‘n-mers’ binding to monovalent receptors. Other types of system to which our model could not be applied in a straightforward manner are the ephrins, as these are not diffusible but membrane-bound ligands.

A problem to which our models could be applied is the response to epithelial growth factor in some tumour cells, where endocytosis is known to be inhibited (Polo et al. 2004). Another interesting problem to which our model could be relevant is the angiogenic rescue that, in some types of tumours, takes place when the tumour has degraded the pre-existing vasculature (Holash et al. 1999; Yancopoulos et al. 2000). Such rescue is thought to be mediated by upregulation of angiopoietin-2 and VEGF but the precise mechanisms are not known. VEGFR upregulation might have a role to play in this process. Together with the incorporation of further details on the regulation of receptor synthesis, these problems will be the subject of future research.


T.A. and K.M.P. thank the EPSRC for financial support under grant GR/509067. K.M.P. thanks the Joint Research Councils (EPSRC, BBSRC and MRC) for support under grant GR/R47455. The authors would like to thank three anonymous referees for their valuable comments on an earlier version of this paper, which have made a major contribution to its final form.


  • Electronic supplementary material is available at or via

  • These three types of VEGFR are surface receptors. There is also a soluble form of VEGFR-1.

  • Throughout the paper, we use the same convention: an upper-case letter represents numbers of molecules of a given type, whereas the corresponding lower-case letter represents the proportion of molecules of that particular kind with respect to the total number of molecules. An exception to this rule is L, whose meaning is explained in the text.

  • Teis & Huber distinguish between active and inactive RTKs. We will assume that ‘active’ refers to dimerized receptors, which seems to be pretty clear from the context, and that ‘inactive’ refers to both unbound receptors and non-dimerized ligand/receptor complexes.

  • Within the statistical physics community, this approximation is often referred to as the eikonal approximation.

  • Kubo et al. state the multidimensional result without a proof.

  • The system of equations (3.12)–(3.25) has 14 equations, but the conservation law Embedded Image allows us to reduce the dimensionality of the system by one unit

  • In fact, this method produces a hierarchy of ‘kinetic’ equations where the cumulants of order n depend on all the cumulants of order up to n−1.

    • Received August 14, 2006.
    • Accepted September 24, 2006.


View Abstract