Skip to main content
  • Research article
  • Open access
  • Published:

Computational modeling reveals molecular details of epidermal growth factor binding



The ErbB family of receptors are dysregulated in a number of cancers, and the signaling pathway of this receptor family is a critical target for several anti-cancer drugs. Therefore a detailed understanding of the mechanisms of receptor activation is critical. However, despite a plethora of biochemical studies and recent single particle tracking experiments, the early molecular mechanisms involving epidermal growth factor (EGF) binding and EGF receptor (EGFR) dimerization are not as well understood. Herein, we describe a spatially distributed Monte Carlo based simulation framework to enable the simulation of in vivo receptor diffusion and dimerization.


Our simulation results are in agreement with the data from single particle tracking and biochemical experiments on EGFR. Furthermore, the simulations reveal that the sequence of receptor-receptor and ligand-receptor reaction events depends on the ligand concentration, receptor density and receptor mobility.


Our computer simulations reveal the mechanism of EGF binding on EGFR. Overall, we show that spatial simulation of receptor dynamics can be used to gain a mechanistic understanding of receptor activation which may in turn enable improved cancer treatments in the future.


Amplification of genes for the ErbB family of receptors is associated with poor outcome in women's cancers, including breast, ovarian and endometrial cancer. Under non-pathological conditions, epidermal growth factor (EGF) receptor (EGFR) or ErbB1 is activated by ligand-induced receptor dimerization, resulting in autophosphorylation and phosphorylation of various cellular substrates [1]. However, while it is clear that overexpression is a factor leading to ligand-independent signaling via these receptors, the mechanism by which functional dimerization and activation occurs is unknown. Since EGF binding represents the initial step for activating EGFR, considerable work has been devoted to elucidating the mechanisms of ligand binding and dimerization [17]. However, molecular details of ligand-induced receptor dimerization are not as well understood.

Apart from in vitro biochemical experiments to study mechanisms of EGFR activation [1], recent developments in microscopy have made it possible to visualize protein dynamics in living cells [8]. The current imaging methods either have a high spatial resolution, such as electron microscopy experiments using immunogold labeling [9] and covalent linking to chemical conjugates like ferritin [10], or high temporal resolution, such as fluorescence confocal microscopy [11], single particle tracking [5] and more recently quantum dots based ligands [12]. However, with the currently available imaging technologies, combined high temporal and spatial resolution (of multiple receptors) has not been achieved.

Computational efforts devoted to understanding the extracellular mechanisms leading to EGFR activation are mostly equilibrium studies [7, 1315] or continuum reaction-diffusion models; see references in [3, 16]. Continuum partial differential equation based models have also been used to represent signaling processes in the plasma membrane assuming a continuum distribution of receptors [17].

While such studies have provided useful insights, they are not ideally suited for describing cell surface heterogeneities, such as microdomains and anomalous diffusion of surface receptors [18], which are important to capture the spatiotemporal receptor dynamics, and lack spatial correlations, known to arise from bimolecular reaction events [19], such as dimerization. Monte Carlo (MC) techniques have proven powerful for systems biology modeling [2022]. In the past, spatial MC approaches have provided mechanistic understanding in other biological systems; see for example [20, 2333].

In this study, we have used a spatial MC framework which not only enables a realistic representation of the plasma membrane, but also facilitates integration of different types of biological data produced from biochemical and microscopy studies to gain insight into the mechanistic details of the underlying biological process. We have developed a general kinetic, lattice MC modeling framework to model the ligand (EGF) binding and dimerization of the EGFR. We compare our simulation results with single particle tracking experiments and analyzed the dominant mechanism of ligand binding and dimerization.


Comparison of stochastic and deterministic models

Microscopic events modeled in this work are shown in Figs. 1 and 2. In order to test the MC algorithm and explore possible differences between stochastic and deterministic models, we have performed a number of simulations for various parameters. Results from the hybrid null-event MC algorithm were compared with an ordinary differential equations (ODEs) model with a set of parameters for which the process is reaction limited, i.e., diffusion is fast compared to reaction, in order to test the validity of the MC algorithm. Specifically, we simulated a dimerization reaction in the absence of ligand considering a high receptor number density of 11,000 per μm2 and assuming no dimers initially. Fig. 3(a) compares the concentration trajectories of dimerized EGFR. This comparison confirms that the hybrid null-event MC algorithm captures the time scales of the system resulting in the correct transient concentration profile. Additional validation carried out under diffusion control has again demonstrated the accuracy of our MC method (Mayawala et al., in preparation).

Figure 1
figure 1

Schematic of simulated microscopic events. Each receptor can diffuse to an empty neighboring site, react with a neighboring receptor to form a dimer, and bind ligand. All events are reversible.

Figure 2
figure 2

Reactions events considered in our model as given in [14].

Figure 3
figure 3

Comparison of hybrid null-event MC and ODE models in terms of (a) dimerized EGFR in the absence of ligand at a high receptor density and diffusivity (11,000 per μm2, D = 2 × 10-14 m2sec-1) and assuming no initial dimers, and (b) EGF bound EGFR in the presence of ligand (160 nM) at low receptor density (125 per per μm2) and D = 2 × 10-15 m2sec-1. The reactions on the figure indicate the dominant processes responsible for the concentration trajectories. Error bars indicate 2 standard deviations obtained from 10 independent MC simulations.

Next in Fig. 3(b) we compared the MC and ODE concentration profiles of EGF bound EGFR monomer in the presence of ligand (160 nM), with a receptor number density of 125 receptors per μm2 and a low diffusivity of 2 × 10-15 m2s-1. The low values of receptor density and diffusivity result in a diffusion controlled case. Corresponding to these parameters, the receptor dimerization rate in the spatial MC model was slower compared to that of the ODE model. The diffusion limited dimerization of EGF bound EGFR monomer leads to a higher concentration of unbound receptor in the spatial MC model than in the ODE model. Thus, spatiotemporal MC simulations are required to capture the transient concentration profiles of the signaling species under diffusion limited conditions. Overall, low receptor densities and low diffusivities may render the system diffusion limited. Under such conditions, well-mixed simulations do not provide accurate dynamics. Use of spatial MC bypasses the question whether the system is diffusion or reaction limited. In a forthcoming communication, we will address quantitatively the conditions for which spatial MC simulations are needed.

Partial differential equations (PDEs) have traditionally been used to model diffusion-reaction processes when spatial effects become important. However, accurate representation of receptor-receptor reactions typically requires MC simulation due to the spatial inhomogeneous distribution of receptors stemming from spatial correlations [19, 28]. Aside from spatial correlations, realistic representation of the plasma membrane microdomains and anomalous diffusion make MC simulation indispensable [18]. Due to these limitations, PDE models have not been employed here.

Comparison of hybrid null-event MC simulations with single particle tracking experiments

The dynamics of the ligand binding events were compared with the single particle tracking experiment of Sako et al. [5] at an EGF concentration of 0.16 nM in the 0–60 sec time interval. To compare simulation results with experimental data, EGF was assumed to be associated with Cy3 dye. A dimerized receptor with two EGF molecules was taken to fluoresce twice as intensely as a receptor (single or dimerized) bound to one EGF molecule. The predicted initial increase of low intensity spots (monomers plus dimers having one EGF bound) followed by a slower increase in high intensity spots (dimers with 2 molecules of EGF) is qualitatively consistent with the experimental data (see Fig. 4(a)). The initial increase in the low intensity signal was due to the rapid binding of EGF on predimerized EGFR. Furthermore, the increase in the total density of Cy3-EGF spots (total bound EGF on all receptors), shown in Fig. 4(b), is also consistent with experimental data.

Figure 4
figure 4

(a) Evolution of intensity of dimerized receptors with two ligands (high intensity spots) and of monomer plus dimerized receptors with a single ligand bound (low intensity spots) along with the data of single particle tracking experiments by Sako et al. over time intervals of 20 sec. The simulations were performed for a receptor number density of 5500 per μm2, a diffusivity of D = 2 × 10-14 m2sec-1, and 18% dimers initially. The simulation intensity has been normalized with the experimental data. (b) Comparison of predicted density of Cy3-EGF spots with experimental data of Sako et al. The densities are normalized with the value at 60 sec. Good agreement of simulations with experimental data is found. In both panels, error bars indicate 2 standard deviations obtained from 10 independent MC simulations.

The possible sequences of events leading to the formation of EGF bound dimerized EGFR at 60 sec are shown in Fig. 5. Sako et al. [5] suggested sequence 1 as being dominant. However, the experimental study alone cannot unambiguously determine the sequence due to its limited spatial resolution and the fact that only ligand bound receptors can be tracked. Our simulations showed that 95–100% of the receptors follow sequence 1, 0–4.9% sequence 2, and the remaining receptors follow sequence 3. Our results are consistent with the hypothesis of Sako et al. [5]. This comparison serves as a model validation step. Small adjustments (20–30%) in the equilibrium and kinetic parameters tabulated in Table 1, which are well within the margins of error, lead to nearly proportional changes in intensity, i.e., no dramatic differences in the simulation profiles are seen (see appendix for details).

Figure 5
figure 5

Sequence of reactions resulting in dimerized receptors with both receptors bound to ligand for simulations of Fig. 4. All reactions are reversible.

Table 1 Kinetic (reaction events given in Fig. 2) and transport (monomer and dimer diffusion) parameters used in hybrid null-event MC model (factors of 1/2 and 1/4 discussed in the Methods section have to be considered).

Effect of ligand concentration on signaling reaction mechanism in A-431 cells (high receptor density)

Single particle tracking experiments [5] are typically limited to low ligand concentrations. High concentration of ligand would lead to fluorescence of a large number of EGFRs making it impossible to visualize individual particles. However, simulations can be used to elucidate the influence of extracellular EGF concentration on EGFR dimerization. Our simulations indicated that the relative contributions of sequences 1–3 at 60 sec change with ligand concentration (Fig. 6(a)). At low ligand concentration, sequence 1 dominates, whereas at higher ligand concentration, a significant fraction of dimers form via sequence 2. Furthermore, sequence 3 also occurs to appreciable extent at high concentration of EGF. At low ligand concentration, most of the ligand gets bound to dimerized receptors, which have a higher ligand affinity; however, the extent to which free EGFR dimerization can occur is limited. At higher ligand concentration, when a significant fraction of ligand is attached to monomers, the coupling between ligand attached monomer and free or ligand attached monomer gives rise to dimers. The relative contribution of the sequences also changes with time. Specifically, initial ligand binding occurs on predimerized receptors, and hence, the relative contribution of sequence 1 is higher at short times. At longer times, after binding of ligand on monomers, sequences 2 and 3 start contributing. With an increase in ligand concentration, the contributions of sequences 2 and 3 increase at a faster rate. The contribution of sequence 3 is higher at longer times after accumulation of ligand bound monomers. As a final note, the time needed to reach equilibrium substantially decreases as the concentration of ligand increases (not shown), e.g., to a total of a few sec at 160 nM. As a result, high ligand concentrations may challenge single particle tracking experiments also in terms of temporal resolution.

Figure 6
figure 6

Contributions of the different reaction mechanisms at 60 sec for different concentrations of EGF with (a) a receptor number density of 5500 receptors per μm2 and D = 2 × 10-14 m2sec-1, (b) a receptor number density of 125 receptors per μm2 and D = 2 × 10-14 m2sec-1, and (c) a receptor number density of 125 receptors per μm2 and D = 2 × 10-15 m2sec-1.

Support for the suggested mechanisms also comes from biochemical studies. The experimental study of [34] reported that at low doses of EGF, inhibition of high affinity binding by mAb108 can kill almost 50–100% of EGF binding, indicating that most of the early binding takes place by sequence 1 at low EGF concentration. However, this inhibition is overcome at higher concentration (~20–50 times) of EGF, which is indicative of substantial formation of EGF bound dimerized EGFR via sequence 2, consistent with the results of our simulations. A larger scale simulation with variable receptor densities in different regions of the plasma membrane will be developed in the future for quantitative comparison with such biochemical experiments. A recent equilibrium based study [13] has shown that such spatial heterogeneities have strong influence on the amount of EGF binding on EGFR, motivating a more detailed analysis of EGFR on the plasma membrane.

Effect of ligand concentration and receptor mobility on signaling reaction mechanism in cells with normal receptor density

Two important factors influencing ligand binding and dimerization are the receptor density and receptor mobility. The receptor density can significantly influence the mechanism of EGF binding as shown in Fig. 6(b). At lower receptor density (125 receptors per μm2) sequence 1 occurs to a much lower extent as compared to the A-431 cells. For this lower receptor density, at lower EGF concentration sequence 2 is dominant, whereas at higher EGF concentrations, sequence 3 is dominant. Sequence 1 is not important at low receptor density, because of the low amount of EGF free dimers (negligible at the low receptor density considered in this work).

A tenfold decrease in receptor mobility (from 2 × 10-14 m2/s to 2 × 10-15 m2/s) leads to a very small increase in the extent of sequence 3, at the expense of sequences 1 and 2 (compare Figs. 6(b) and 6(c)). This small increase is observed only at low EGF concentration. At higher EGF concentration this increase is even smaller. Sequence 3 occurs to a larger extent at slower diffusion because dimerization is slowing down and so more monomers associate with ligand. At higher EGF concentration, this effect is not as prominent because EGF binding is faster leading to more EGF bound EGFRs, thereby increased dimerization occurs among EGF bound monomer EGFRs even with a higher receptor diffusivity.

Several studies have indicated inhomogeneities in the plasma membrane and excellent reviews have been published on this topic including [18, 3538]. These studies have suggested localization of receptors within small regions, called microdomains, in the plasma membrane. An implication of the containment of receptors in the microdomains is the observation of lower macroscopic diffusivity as has been discussed in [39]. As a result, the microscopic diffusivity can potentially be at least 1–2 orders of magnitude faster than the diffusivity reported in literature. Therefore, we have also studied the effect of a higher diffusivity. In contrast to decreasing diffusivity from 2 × 10-15 m2/s to 2 × 10-14 m2/s mentioned above, larger changes are observed at high ligand concentration (e.g., 1600 nM) and a receptor density of 125 receptors per μm2 for a change in diffusivity from 2 × 10-14 m2/s to 2 × 10-13 m2/s. Specifically, the contribution of sequence 2 increases from ~15% to ~30% at the expense of sequence 3 which decreases from ~85% to ~70%. An increase in receptor diffusivity leads to an increased rate of dimerization between an occupied and a free receptor in comparison to ligand binding on a free receptor. Overall, a faster diffusivity can lead to an overall increase in the dimerization rate but this effect is not dramatic under our simulation conditions.


Our simulation results suggest future single particle tracking experiments or related microscopy experiments. It may be difficult to perform the single particle tracking experiments of [5] at higher ligand concentration in A-431 cells due to the difficulty in visualization of single EGFR and possibly to the short time scales over which transients are over. However, such experiments can potentially be performed in cells with a lower average receptor density. On such cells, the increased contributions of sequences 2 and 3 should be observed to further validate our model. Possible discrepancies between experiments and model could provide new insights to enhance our current understanding of the underlying signaling processes.

The variation in receptor density and receptor mobility can stem from different cell types as well as different spatial features/locations in the plasma membrane (see Methods section for references). Future microscopy experiments should be designed to observe the reaction events and transients of low and high intensity spots, as reported by [5], in different domains of the plasma membrane in the same cell. Such data can then be used to estimate the local density of the receptors which in turn can help in understanding the receptor distribution in the plasma membrane.

This work shows the influence of receptor density and receptor mobility as a biophysical control of signaling processes over the inflexible thermodynamic and biochemical properties. A key suggestion from this work is that it is not adequate to treat the receptor-receptor interactions based only on their kinetic and thermodynamic parameters. Instead, their spatial properties, such as local receptor densities and local lateral mobility, can play significant roles in determining the intracellular signaling. Herein, we have used a simplistic representation of the plasma membrane. This model can be extended to include experimentally reported anomalous hop-like diffusion [18] and other spatial features, such as clathrin coated pits, lipid rafts, caveolae and other microdomains [35, 40]. Such features not only facilitate an enhanced control at the level of plasma membrane but can also be important for the wide diversity of signaling outcomes from limited varieties of ligands and receptors.


We have developed a computational framework for studying cell surface receptor dynamics that can bridge biochemical data on one hand with various microscopy experiments on the other, which currently lack simultaneous high spatial and temporal resolution. This work provides an important step forward in the era of in vivo imaging based modeling approaches [8, 41]. For example, in this work comparison of MC simulations with single particle tracking experiments reveals how the sequence of receptor-receptor and ligand-receptor reaction events depends on the ligand concentration, receptor density and receptor mobility. Our computer simulations reveal the underlying mechanism on the plasma membrane leading to dimerized and ligand (EGF) bound receptors.

Considering the interest in targeted antibodies for cancer therapeutics [42], a detailed understanding of the biochemical mechanisms involved in signal sensing at the plasma membrane is desired. Future advancements in medicine will ultimately also include the mathematical analysis and modeling of ErbB receptor diffusion, dimerization and clustering, which will increase our understanding of tumorigenesis and lead to medical advances, such as individualized therapy for heterogeneous cancers.


A hybrid null-event algorithm

A coarse, molecular level based computational framework that leaves atomistic details out, e.g., conformations and vibrations of proteins into potential energy surface minima, but still provides the sequence of molecular events at the receptor length scale was employed. Herein, we have utilized a kinetic, lattice MC method for simulating the EGFR dynamics. Microscopic events modeled include receptor dimerization and decomposition, ligand-receptor association and dissociation, and Brownian diffusion of receptors (see Figs. 1 and 2). Formation of high-mers that happens at longer times is not considered in this work.

The existence of multiple timescales in the system and low surface density of receptors make lattice MC simulations computationally prohibitive. We have devised a hybrid between the continuous time MC method [43] and the null-event MC method (see [44] for an overview of spatial MC algorithms) to increase the speed of the null-event algorithm but maintain its flexibility. In our hybrid null-event algorithm, only lattice sites filled with receptors were randomly selected, resulting in two to four orders of magnitude speedup, depending on receptor density, relative to a traditional null-event MC algorithm where all sites are randomly chosen. Additionally, operations that are often involved in a continuous time MC method, such as summation of transition probabilities and searches over the entire lattice, are avoided, resulting in further acceleration of simulations.

The spatial domain was represented using a two-dimensional square lattice that was initially randomly populated with a given density of receptors. Periodic boundary conditions were employed [45]. After initialization, an occupied site was randomly picked and one of the microscopic events is possibly selected to occur based on probabilities described below.

Transition probability of diffusion

The probability of diffusion per unit time in all four directions on a square lattice was calculated using random walk theory from the diffusivity

Γ d = 4 D a 2 , ( 1 ) MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciGacaGaaeqabaqabeGadaaakeaacqqHtoWrdaqhaaqaaaqaaiabdsgaKbaacqGH9aqpdaWcaaqaaiabisda0iabdseaebqaaiabdggaHnaaDaaabaaabaGaeGOmaidaaaaacqGGSaalcaWLjaGaaCzcamaabmaabaGaeGymaedacaGLOaGaayzkaaaaaa@39A2@

where a is the microscopic lattice pixel dimension (taken here to be 2 nm) and D is the corresponding diffusivity of a receptor or a dimer. The transition probability of diffusion per unit time in moving from site i to site j is

Γ i j d = 1 4 Γ d σ i ( 1 σ j ) j B i ( 2 ) MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciGacaGaaeqabaqabeGadaaakeaacqqHtoWrdaqhaaqaaiabdMgaPjabgkziUkabdQgaQbqaaiabdsgaKbaacqGH9aqpdaWcaaqaaiabigdaXaqaaiabisda0aaacqqHtoWrdaqhaaqaaaqaaiabdsgaKbaacqaHdpWCdaWgaaqaaiabdMgaPbqabaGaeiikaGIaeGymaeJaeyOeI0Iaeq4Wdm3aaSbaaeaacqWGQbGAaeqaaiabcMcaPiaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaemOAaOMaeyicI4SaemOqai0aaSbaaSqaaiabdMgaPbqabaGccaWLjaGaaCzcamaabmaabaGaeGOmaidacaGLOaGaayzkaaaaaa@5509@

where B i denotes the set of sites to which diffusion from site i can occur. In our model, this set includes all 4 first-nearest neighboring sites. σ i is the occupancy (discrete) function that is 1, if site i is filled, or 0, if site i is empty (a single index indicating the site is herein used to simplify notation). According to Eq. (2), the transition probability of diffusion per unit time along any direction on the square lattice can be 0 or 1 4 Γ d MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciGacaGaaeqabaqabeGadaaakeaadaWcaaqaaiabigdaXaqaaiabisda0aaacqqHtoWrdaqhaaqaaaqaaiabdsgaKbaaaaa@317F@ , depending on the occupancy of the first-nearest neighboring site in the corresponding direction.

Transition probability of reactions

The transition probability of a reaction was obtained in terms of the macroscopic reaction rate constants, k. For a first-order (e.g., the decomposition of an EGFR dimer) or pseudo-first order (e.g., the EGF binding onto a receptor because EGF is assumed at a constant concentration) reaction (see Fig. 2), one has

A C, Γ i r = k σ i ( 3 ) MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciGacaGaaeqabaqabeGadaaakeaacqqGbbqqcqGHsgIRcqqGdbWqcqqGSaalcaaMc8UaaGPaVlaaykW7cqqHtoWrdaqhaaqaaiabdMgaPbqaaiabdkhaYbaacqGH9aqpcqWGRbWAiiaacqWFdpWCdaWgaaqaaiabdMgaPbqabaGaaCzcaiaaxMaadaqadaqaaiabiodaZaGaayjkaiaawMcaaaaa@43ED@

where Γ i r MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciGacaGaaeqabaqabeGadaaakeaacqqHtoWrdaqhaaqaaiabdMgaPbqaaiabdkhaYbaaaaa@3100@ is the transition probability of reaction at site i. For a bimolecular reaction on a square lattice (e.g., receptor-receptor dimerization), the transition probability per unit time at selected site i was modeled as:

for the reaction: A + B C , Γ i r = k 4 σ i σ j , ( 4 ) MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciGacaGaaeqabaqabeGadaaakeaacqqGMbGzcqqGVbWBcqqGYbGCcqqGGaaicqqG0baDcqqGObaAcqqGLbqzcqqGGaaicqqGYbGCcqqGLbqzcqqGHbqycqqGJbWycqqG0baDcqqGPbqAcqqGVbWBcqqGUbGBcqqG6aGocqqGGaaicqqGbbqqcqqGGaaicqqGRaWkcqqGGaaicqqGcbGqcqGHsgIRieaacqWFdbWqcqGGSaalcqqGGaaicqqHtoWrdaqhaaqaaiab=LgaPbqaaiab=jhaYbaacqGH9aqpdaWcaaqaaiab=TgaRbqaaiabisda0aaaiiaacqGFdpWCdaWgaaqaaiab=LgaPbqabaGae43Wdm3aaSbaaeaacqWFQbGAaeqaaiabcYcaSiaaxMaacaWLjaWaaeWaaeaacqaI0aanaiaawIcacaGLPaaaaaa@5F0B@

and for the reaction: 2A C,  Γ i r = k 2 σ i σ j . ( 5 ) MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciGacaGaaeqabaqabeGadaaakeaacqqGHbqycqqGUbGBcqqGKbazcqqGGaaicqqGMbGzcqqGVbWBcqqGYbGCcqqGGaaicqqG0baDcqqGObaAcqqGLbqzcqqGGaaicqqGYbGCcqqGLbqzcqqGHbqycqqGJbWycqqG0baDcqqGPbqAcqqGVbWBcqqGUbGBcqqG6aGocqqGGaaicqqGYaGmcqqGbbqqcqGHsgIRcqqGdbWqcqqGSaalcqqGGaaicqqHtoWrdaqhaaqaaGqaaiab=LgaPbqaaiab=jhaYbaacqGH9aqpdaWcaaqaaiab=TgaRbqaaiabikdaYaaaiiaacqGFdpWCdaWgaaqaaiab=LgaPbqabaGae43Wdm3aaSbaaeaacqWFQbGAaeqaaiabc6caUiaaxMaacaWLjaWaaeWaaeaacqaI1aqnaiaawIcacaGLPaaaaaa@6145@

Here the reacting species (A and B or A and A) occupy adjacent sites i and j, and the units of k are (molecules/site)-1sec-1. The factor of four in k/4 in Eq. (4) accounts for the fact that all four neighboring sites of site i were randomly chosen to search for the existence of species B. Similarly, the factor of 2 in Eq. 5 was due to the degeneracy of species participating in the homodimerization.

Event selection and time advancement

After an occupied site, say i, was selected, the transition probabilities per unit time of all possible events were computed. The probability for a certain event 'x' at site i was calculated as

p i x = Γ i x Γ max , ( 6 ) MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciGacaGaaeqabaqabeGadaaakeaacqWGWbaCdaqhaaqaaiabdMgaPbqaaiabdIha4baacqGH9aqpdaWcaaqaaiabfo5ahnaaDaaabaGaemyAaKgabaGaemiEaGhaaaqaaiabfo5ahnaaBaaabaGagiyBa0MaeiyyaeMaeiiEaGhabeaaaaGaeiilaWIaaCzcaiaaxMaadaqadaqaaiabiAda2aGaayjkaiaawMcaaaaa@40D7@

where Γmax is a normalization constant to ensure that the selection probability of each event is always less than or equal to 1 and Γ i x MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciGacaGaaeqabaqabeGadaaakeaacqqHtoWrdaqhaaqaaiabdMgaPbqaaiabdIha4baaaaa@310C@ is the transition probability of event x (reaction or diffusion) at site i, as defined above. The maximum values of transition probabilities per unit time of events described by Eqs. (2)-(4) are Γd/4, k, k/4, and k/2, respectively. We defined the normalization constant as

Γ max = 4 ( Γ d 4 + max { a l l f o r w a r d r e a c t i o n e v e n t s Γ r } ) , + max { a l l b a c k w a r d r e a c t i o n e v e n t s Γ r } MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciGacaGaaeqabaqabeGadaaakqaabeqaaiabfo5ahnaaBaaabaGagiyBa0MaeiyyaeMaeiiEaGhabeaacqGH9aqpcqaI0aandaqadaqaamaalaaabaGaeu4KdC0aaWbaaeqabaGaemizaqgaaaqaaiabisda0aaacqGHRaWkcyGGTbqBcqGGHbqycqGG4baEdaGadaqaamaaqafabaGaeu4KdC0aaWbaaeqabaGaemOCaihaaaqaaiabdggaHjabdYgaSjabdYgaSjabbccaGiabdAgaMjabd+gaVjabdkhaYjabdEha3jabdggaHjabdkhaYjabdsgaKjabbccaGiabdkhaYjabdwgaLjabdggaHjabdogaJjabdsha0jabdMgaPjabd+gaVjabd6gaUjabbccaGiabdwgaLjabdAha2jabdwgaLjabd6gaUjabdsha0jabdohaZbqabiabggHiLdaacaGL7bGaayzFaaaacaGLOaGaayzkaaGamaiGaaaqciilaWcabaGaey4kaSIagiyBa0MaeiyyaeMaeiiEaG3aaiWaaeaadaaeqbqaaiabfo5ahnaaCaaabeqaaiabdkhaYbaaaeaacqWGHbqycqWGSbaBcqWGSbaBcqqGGaaicqWGIbGycqWGHbqycqWGJbWycqWGRbWAcqWG3bWDcqWGHbqycqWGYbGCcqWGKbazcqqGGaaicqWGYbGCcqWGLbqzcqWGHbqycqWGJbWycqWG0baDcqWGPbqAcqWGVbWBcqWGUbGBcqqGGaaicqWGLbqzcqWG2bGDcqWGLbqzcqWGUbGBcqWG0baDcqWGZbWCaeqacqGHris5aaGaay5Eaiaaw2haaaaaaa@9C41@

where multiplication by a factor of 4 accounts for the microscopic process in each of the four directions on the square lattice. For the reaction events given in Fig. 2

Γ max = Γ d + 4 ( k 1 f 2 + k 2 f 4 + k 3 f 2 ) + i = 4 6 k i f + i = 1 6 k i b , ( 7 ) MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciGacaGaaeqabaqabeGadaaakeaacqqHtoWrdaWgaaqaaiGbc2gaTjabcggaHjabcIha4bqabaGaeyypa0Jaeu4KdC0aaWbaaeqabaGaemizaqgaaiabgUcaRiabisda0maabmaabaWaaSaaaeaacqWGRbWAdaWgaaqaaiabigdaXiabdAgaMbqabaaabaGaeGOmaidaaiabgUcaRmaalaaabaGaem4AaS2aaSbaaeaacqaIYaGmcqWGMbGzaeqaaaqaaiabisda0aaacqGHRaWkdaWcaaqaaiabdUgaRnaaBaaabaGaeG4mamJaemOzaygabeaaaeaacqaIYaGmaaaacaGLOaGaayzkaaGaey4kaSYaaabCaeaacqWGRbWAdaWgaaqaaiabdMgaPjabdAgaMbqabaaabaGaemyAaKMaeyypa0JaeGinaqdabaGaeGOnaydacqGHris5aiabgUcaRmaaqahabaGaem4AaS2aaSbaaeaacqWGPbqAcqWGIbGyaeqaaaqaaiabdMgaPjabg2da9iabigdaXaqaaiabiAda2aGaeyyeIuoacqGGSaalcaWLjaGaaCzcamaabmaabaGaeG4naCdacaGLOaGaayzkaaaaaa@6577@

where the subscripts f and b refer to the forward and backward reaction events, respectively, and i denotes the corresponding reaction according to Fig. 2. A random number 'r' was finally chosen from a uniform distribution between 0 and 1. The events were randomly ranked. The smallest value of m satisfying x = 1 m p i x > r MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciGacaGaaeqabaqabeGadaaakeaadaaeWbqaaiabdchaWnaaDaaabaGaemyAaKgabaGaemiEaGhaaaqaaiabdIha4jabg2da9iabigdaXaqaaiabd2gaTbGaeyyeIuoacqGH+aGpcqWGYbGCaaa@3A7E@ criterion was chosen. If no value of m satisfied the above criterion, then no event was selected to occur (a null-event) and a new occupied site was again randomly picked. Otherwise, the mthevent was selected, the populations on the lattice were updated accordingly to reflect the stoichiometry of the reaction or diffusion process, and the real time was advanced based on the most frequently selected event as suggested in [44]. In our simulations, real time was advanced based on the diffusion of receptors according to which the average time step after each successful diffusion event was calculated as

Δ t = 1 i = 1 N o . o f s i t e s ( σ i j B i Γ i j d ( 1 - σ j ) ) = 1 1 4 Γ d o c c u p i e d s i t e s ( j B i ( 1 - σ j ) ) . ( 8 ) MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciGacaGaaeqabaqabeGadaaakeaaiiaacqWFuoarieaacqGF0baDcqGF9aqpdaWcaaqaaiab+fdaXaqaamaaqahabaWaaeWaaeaacqWFdpWCdaWgaaqaaiab+LgaPbqabaWaaabCaeaacqWFtoWrdaqhaaqaaiab+LgaPjabgkziUkab+PgaQbqaaiab+rgaKbaadaqadaqaaiab+fdaXiab+1caTiab=n8aZnaaBaaaleaacqGFQbGAaeqaaaGccaGLOaGaayzkaaaabaGae4NAaOMae8hcI4Sae4Nqai0aaSbaaSqaaiab+LgaPbqabaaakeaaaiabggHiLdaacaGLOaGaayzkaaaabaGae4xAaKMae4xpa0Jae4xmaedabaGae4Nta4Kae43Ba8Mae4Nla4Iae4hiaaIae43Ba8Mae4NzayMae4hiaaIae43CamNae4xAaKMae4hDaqNae4xzauMae43CamhacqGHris5aaaacqGF9aqpdaWcaaqaaiab+fdaXaqaamaalaaabaGae4xmaedabaGae4hnaqdaaiab=n5ahnaaCaaabeqaaiab+rgaKbaadaaeWbqaamaabmaabaWaaabCaeaadaqadaqaaiab+fdaXiab+1caTiab=n8aZnaaBaaaleaacqGFQbGAaeqaaaGccaGLOaGaayzkaaaabaGae4NAaOMae8hcI4Sae4Nqai0aaSbaaSqaaiab+LgaPbqabaaakeaaaiabggHiLdaacaGLOaGaayzkaaaabaGae43Ba8Mae43yamMae43yamMae4xDauNae4hCaaNae4xAaKMae4xzauMae4hzaqMae4hiaaIae43CamNae4xAaKMae4hDaqNae4xzauMae43CamhabaaacqGHris5aaaacqGFGaaicqGFUaGlcaWLjaGaaCzcamaabmaabaGae4hoaGdacaGLOaGaayzkaaaaaa@9004@

The summation in the denominator was updated each time a successful event happens by subtracting the previous occupancy values affected by the selected event and adding the new ones. In this way, the summation was carried out only once (at the beginning of a simulation) over all occupied sites. Finally, a new site is picked randomly until the desired real time is reached.

The hybrid null-event MC algorithm explained above was implemented (for details on the null-event algorithm see [44]) using Fortran 90. For each data point, 10 simulations with different seeds of the random number generator were used to collect statistics.

Simulation size and model parameters

In this work, the cell surface was represented using a 2-dimensional square lattice with each pixel being 2 nm × 2 nm in size. The number density of receptors ranges from ~102 receptors per μm2 on normal cells [46] to ~103 receptors per μm2 on human epithelioid carcinoma cells (A-431 cells) which overexpress EGFR [2]. However, the local density of receptors can be much higher in-vivo because of the localization of receptors in certain regions of the plasma membrane, such as in lipid rafts [35, 47, 48]. We simulated a low density (to represent normal cells) of 31 receptors on 500 nm × 500 nm mesh and a high density (to represent A-431 cells) of 55 receptors on 100 nm × 100 nm mesh, which are equivalent to receptor number densities of 125 and 5500 per μm2, respectively. The diffusivity of monomer EGFR has been reported to be around 2 × 10-14 m2s-1 [49, 50]. Lower macroscopic diffusivities for EGFR have also been observed [49], which may be due to containment within cytoskeletal elements [51] or lipid rafts [52]. Based on these suggestions, the effect of a slower diffusivity is also analyzed by considering a diffusivity of 2 × 10-15 m2s-1. Simulations have been performed in the 0–60 sec time interval to capture the initial transients.

The model parameters are summarized in Table 1. We assumed two types of EGF binding on the cell surface: low affinity binding (on monomer EGFR) and high affinity binding (on dimerized EGFR). Experimental information suggests the existence of a high affinity (for receptor-ligand association) receptor population, most of which, if not all, is present in the form of dimers; see review by [1]. A recent equilibrium study has also shown that this interpretation is consistent with the experimentally reported concave-up shape of the Scatchard plot [13].

Experimental studies have provided evidence of predimerized receptors on A-431 cells to different extents [9, 5356]. Consequently, a fraction (~82%) of receptors was initially placed at random locations as monomers and the remaining as dimers on simulated A-431 cells for comparison with single particle tracking data. Corresponding to the dimerization equilibrium constant for this data, we found that at the lower receptor number density of 125 per μm2, there is negligible number of dimers in the absence of ligand. The receptor dimerization constants vary with ligand occupancy. Several experimental studies have shown that dimerization between unbounded receptors occurs with lower affinity than that between one bounded and one unbounded receptor. Finally, dimerization between two ligand bounded receptors occurs with the highest affinity [6, 57].


A sensitivity analysis was performed in which each kinetic parameter (k i , i = 1f, 4f, 5f, 6f, 1b, 4b, 5b, 6b) was increased by 20%, and the change in the high intensity spots was observed at three different times (20, 40 and 60 sec) from the mean of 10 independent MC simulations. The normalized sensitivity coefficient, reported in Fig. 7, is defined as

Figure 7
figure 7

Normalized sensitivity coefficients at 3 different times (20, 40 and 60 sec) calculated by introducing a 20% increase in the kinetic parameter indicated on the x-axis.

( I I o ) / I o ( k i k i o ) / k i o = I I o 0.2 I o , ( A1 ) MathType@MTEF@5@5@+=feaafiart1ev1aqatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciGacaGaaeqabaqabeGadaaakeaadaWcaaqaamaalyaabaGaeiikaGIaemysaKKaeyOeI0IaemysaK0aaSbaaSqaaiabd+gaVbqabaGccqGGPaqkaeaacqWGjbqsdaWgaaWcbaGaem4Ba8gabeaaaaaakeaadaWcgaqaaiabcIcaOiabdUgaRnaaBaaaleaacqWGPbqAaeqaaOGaeyOeI0Iaem4AaS2aaSbaaSqaaiabdMgaPjabd+gaVbqabaGccqGGPaqkaeaacqWGRbWAdaWgaaWcbaGaemyAaKMaem4Ba8gabeaaaaaaaOGaeyypa0ZaaSaaaeaacqWGjbqscqGHsislcqWGjbqsdaWgaaWcbaGaem4Ba8gabeaaaOqaaiabicdaWiabc6caUiabikdaYiabdMeajnaaBaaaleaacqWGVbWBaeqaaaaakiabcYcaSiaaxMaacaWLjaWaaeWaaeaacqqGbbqqcqqGXaqmaiaawIcacaGLPaaaaaa@5543@

where I is the % of high intensity spots (y axis in Fig. 4a) upon perturbing the kinetic parameter, ki, and Io is the nominal value corresponding to the original set of kinetic parameters, kio. Only 8 of the 12 kinetic parameters were independently perturbed because of the two equilibrium constrains reported in [14], i.e.,

K 2 = K 1 K 5 K 4 ,  and  ( A2 ) MathType@MTEF@5@5@+=feaafiart1ev1aqatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciGacaGaaeqabaqabeGadaaakeaacqWGlbWsdaWgaaWcbaGaeGOmaidabeaakiabg2da9maalaaabaGaem4saS0aaSbaaSqaaiabigdaXaqabaGccqWGlbWsdaWgaaWcbaGaeGynaudabeaaaOqaaiabdUealnaaBaaaleaacqaI0aanaeqaaaaakiabcYcaSiabbccaGiabbggaHjabb6gaUjabbsgaKjabbccaGiaaxMaacaWLjaWaaeWaaeaacqqGbbqqcqqGYaGmaiaawIcacaGLPaaaaaa@4213@

K 3 = K 1 K 5 K 6 ( K 4 ) 2 . ( A3 ) MathType@MTEF@5@5@+=feaafiart1ev1aqatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciGacaGaaeqabaqabeGadaaakeaacqWGlbWsdaWgaaWcbaGaeG4mamdabeaakiabg2da9maalaaabaGaem4saS0aaSbaaSqaaiabigdaXaqabaGccqWGlbWsdaWgaaWcbaGaeGynaudabeaakiabdUealnaaBaaaleaacqaI2aGnaeqaaaGcbaWaaeWaaeaacqWGlbWsdaWgaaWcbaGaeGinaqdabeaaaOGaayjkaiaawMcaamaaCaaaleqabaGaeGOmaidaaaaakiabc6caUiaaxMaacaWLjaWaaeWaaeaacqqGbbqqcqqGZaWmaiaawIcacaGLPaaaaaa@4193@

The equilibrium relations determine the changes in the kinetic parameters of the dependent reactions 2 and 3 (see Fig. 2) upon perturbing those of the linearly independent reactions. A change in an equilibrium constant can be associated with a change in the forward, backward, or both rate constants. For simplicity a change in the rate constant of a forward (backward) linearly independent reaction is taken to cause a change in the forward (backward) rate constant of the linearly dependent reactions. Specifically, one has

k 2 f = ( k 2 f ) o f k 1 f f k 5 f f k 4 f , ( A4 ) MathType@MTEF@5@5@+=feaafiart1ev1aqatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciGacaGaaeqabaqabeGadaaakeaacqWGRbWAdaWgaaWcbaGaeGOmaiJaemOzaygabeaakiabg2da9maabmaabaGaem4AaS2aaSbaaSqaaiabikdaYiabdAgaMbqabaaakiaawIcacaGLPaaadaWgaaWcbaGaem4Ba8gabeaakmaalaaabaGaemOzay2aaSbaaSqaaiabdUgaRnaaBaaameaacqaIXaqmcqWGMbGzaeqaaaWcbeaakiabdAgaMnaaBaaaleaacqWGRbWAdaWgaaadbaGaeGynauJaemOzaygabeaaaSqabaaakeaacqWGMbGzdaWgaaWcbaGaem4AaS2aaSbaaWqaaiabisda0iabdAgaMbqabaaaleqaaaaakiabcYcaSiaaxMaacaWLjaWaaeWaaeaacqqGbbqqcqqG0aanaiaawIcacaGLPaaaaaa@4E8B@

k 3 f = ( k 3 f ) o f k 1 f f k 5 f f k 6 f ( f k 4 f ) 2 , ( A5 ) MathType@MTEF@5@5@+=feaafiart1ev1aqatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciGacaGaaeqabaqabeGadaaakeaacqWGRbWAdaWgaaWcbaGaeG4mamJaemOzaygabeaakiabg2da9maabmaabaGaem4AaS2aaSbaaSqaaiabiodaZiabdAgaMbqabaaakiaawIcacaGLPaaadaWgaaWcbaGaem4Ba8gabeaakmaalaaabaGaemOzay2aaSbaaSqaaiabdUgaRnaaBaaameaacqaIXaqmcqWGMbGzaeqaaaWcbeaakiabdAgaMnaaBaaaleaacqWGRbWAdaWgaaadbaGaeGynauJaemOzaygabeaaaSqabaGccqWGMbGzdaWgaaWcbaGaem4AaS2aaSbaaWqaaiabiAda2iabdAgaMbqabaaaleqaaaGcbaWaaeWaaeaacqWGMbGzdaWgaaWcbaGaem4AaS2aaSbaaWqaaiabisda0iabdAgaMbqabaaaleqaaaGccaGLOaGaayzkaaWaaWbaaSqabeaacqaIYaGmaaaaaOGaeiilaWIaaCzcaiaaxMaadaqadaqaaiabbgeabjabbwda1aGaayjkaiaawMcaaaaa@56B4@

k 2 b = ( k 2 b ) o f k 1 b f k 5 b f k 4 b ,  and  ( A6 ) MathType@MTEF@5@5@+=feaafiart1ev1aqatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciGacaGaaeqabaqabeGadaaakeaacqWGRbWAdaWgaaWcbaGaeGOmaiJaemOyaigabeaakiabg2da9maabmaabaGaem4AaS2aaSbaaSqaaiabikdaYiabdkgaIbqabaaakiaawIcacaGLPaaadaWgaaWcbaGaem4Ba8gabeaakmaalaaabaGaemOzay2aaSbaaSqaaiabdUgaRnaaBaaameaacqaIXaqmcqWGIbGyaeqaaaWcbeaakiabdAgaMnaaBaaaleaacqWGRbWAdaWgaaadbaGaeGynauJaemOyaigabeaaaSqabaaakeaacqWGMbGzdaWgaaWcbaGaem4AaS2aaSbaaWqaaiabisda0iabdkgaIbqabaaaleqaaaaakiabcYcaSiabbccaGiabbggaHjabb6gaUjabbsgaKjabbccaGiaaxMaacaWLjaWaaeWaaeaacqqGbbqqcqqG2aGnaiaawIcacaGLPaaaaaa@53F0@

k 3 b = ( k 3 b ) o f k 1 b f k 5 b f k 6 b ( f k 4 b ) 2 , ( A7 ) MathType@MTEF@5@5@+=feaafiart1ev1aqatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciGacaGaaeqabaqabeGadaaakeaacqWGRbWAdaWgaaWcbaGaeG4mamJaemOyaigabeaakiabg2da9maabmaabaGaem4AaS2aaSbaaSqaaiabiodaZiabdkgaIbqabaaakiaawIcacaGLPaaadaWgaaWcbaGaem4Ba8gabeaakmaalaaabaGaemOzay2aaSbaaSqaaiabdUgaRnaaBaaameaacqaIXaqmcqWGIbGyaeqaaaWcbeaakiabdAgaMnaaBaaaleaacqWGRbWAdaWgaaadbaGaeGynauJaemOyaigabeaaaSqabaGccqWGMbGzdaWgaaWcbaGaem4AaS2aaSbaaWqaaiabiAda2iabdkgaIbqabaaaleqaaaGcbaWaaeWaaeaacqWGMbGzdaWgaaWcbaGaem4AaS2aaSbaaWqaaiabisda0iabdkgaIbqabaaaleqaaaGccaGLOaGaayzkaaWaaWbaaSqabeaacqaIYaGmaaaaaOGaeiilaWIaaCzcaiaaxMaadaqadaqaaiabbgeabjabbEda3aGaayjkaiaawMcaaaaa@5688@

where subscript 'o' denotes the nominal value, and f k i MathType@MTEF@5@5@+=feaafiart1ev1aqatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciGacaGaaeqabaqabeGadaaakeaacqWGMbGzdaWgaaWcbaGaem4AaS2aaSbaaWqaaiabdMgaPbqabaaaleqaaaaa@3122@ is 1, if the kinetic parameter (k i ) is not perturbed, and 1.2, if the parameter is increased by 20%.


  1. Jorissen RN, Walker F, Pouliot N, Garrett TPJ, Ward CW, Burgess AW: Epidermal growth factor receptor: mechanisms of activation and signalling. Experimental Cell Research. 2003, 284: 31-53. 10.1016/S0014-4827(02)00098-8.

    Article  CAS  PubMed  Google Scholar 

  2. Wiley HS: Anomalous binding of epidermal growth factor to A431 cells is due to the effect of high receptor densities and a saturable endocytic system. J Cell Biol. 1988, 107: 801-810. 10.1083/jcb.107.2.801.

    Article  CAS  PubMed  Google Scholar 

  3. Wiley HS, Shvartsman SY, Lauffenburger DA: Computational modeling of the EGF-receptor system: a paradigm for systems biology. Trends in Cell Biology. 2003, 13: 43-50. 10.1016/S0962-8924(02)00009-0.

    Article  CAS  PubMed  Google Scholar 

  4. Schlessinger J: Allosteric regulation of the epidermal growth factor receptor kinase. J Cell Biol. 1986, 103: 2067-2072. 10.1083/jcb.103.6.2067.

    Article  CAS  PubMed  Google Scholar 

  5. Sako Y, Minoghchi S, Yanagida T: Single-molecule imaging of EGFR signalling on the surface of living cells. Nature Cell Biology. 2000, 2: 168-172. 10.1038/35004044.

    Article  CAS  PubMed  Google Scholar 

  6. Lemmon MA, Bu Z, Ladbury JE, Zhou M, Pinchasi D, Lax I, Engelman DM, Schlessinger J: Two EGF molecules contribute additively to stabilization of the EGFR dimer. The EMBO Journal. 1997, 16: 281-294. 10.1093/emboj/16.2.281.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  7. Klein P, Mattoon D, Lemmon MA, Schlessinger J: A structure-based model for ligand binding and dimerization of EGF receptors. PNAS. 2004, 101: 929-934. 10.1073/pnas.0307285101.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  8. Weijer CJ: Visualizing signals moving in cells. Science. 2003, 300: 96-100. 10.1126/science.1082830.

    Article  CAS  PubMed  Google Scholar 

  9. Vanbelzen N, Rijken PJ, Hage WJ, Delaat SW, Verkleij AJ, Boonstra J: Direct visualization and quantitative-analysis of epidermal growth factor-induced receptor clustering. Journal Of Cellular Physiology. 1988, 134: 413-420. 10.1002/jcp.1041340312.

    Article  CAS  Google Scholar 

  10. Haigler HT, McKanna JA, Cohen S: Direct visualization of the binding and internalization of a ferritin conjugate of epidermal growth factor in human carcinoma cells A-431. J Cell Biol. 1979, 81: 382-395. 10.1083/jcb.81.2.382.

    Article  CAS  PubMed  Google Scholar 

  11. Verveer PJ, Wouters FS, Reynolds AR, Bastiaens PIH: Quantitative Imaging of Lateral ErbB1 Receptor Signal Propagation in the Plasma Membrane. Science. 2000, 290: 1567-1570. 10.1126/science.290.5496.1567.

    Article  CAS  PubMed  Google Scholar 

  12. Lidke DS, Nagy P, Heintzmann R, Arndt-Jovin DJ, Post JN, Grecco HE, Jares-Erijman EA, Jovin TM: Quantum dot ligands provide new insights into erbB/HER receptor–mediated signal transduction. Nature Biotechnology. 2004, 22: 198 -1203. 10.1038/nbt929.

    Article  CAS  PubMed  Google Scholar 

  13. Mayawala K, Vlachos DG, Edwards JS: Heterogeneities in EGF receptor density at the cell surface can lead to concave up scatchard plot of EGF binding. FEBS Letters. 2005, 579: 3043-3047. 10.1016/j.febslet.2005.04.059.

    Article  CAS  PubMed  Google Scholar 

  14. Wofsy C, Goldstein B, Lund K, Wiley HS: Implications of epidermal growth factor (EGF) induced EGF receptor aggregation. Biophys J. 1992, 63: 98-110.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  15. Chamberlin SG, Davies DE: A unified model of c-erbB receptor homo- and heterodimerisation. Biochimica et Biophysica Acta (BBA) - Protein Structure and Molecular Enzymology. 1998, 1384: 223-232. 10.1016/S0167-4838(97)00203-3.

    Article  CAS  Google Scholar 

  16. Lauffenburger DA, Linderman JJ: Receptors models for binding, trafficking, and signaling. 1993, New York, Oxford University Press

    Google Scholar 

  17. Haugh JM: A unified model for signal transduction reactions in cellular membranes. Biophysical Journal. 2002, 82: 591-604.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  18. Kusumi A, Nakada C, Ritchie K, Murase K, Suzuki K, Murakoshi H, Kasai RS, Kondo J, Fujiwara T: Paradigm Shift of the Plasma Membrane Concept from the Two-Dimensional Continuum Fluid to the Partitioned Fluid: High-Speed Single-Molecule Tracking of Membrane Molecules. Annual Review of Biophysics and Biomolecular Structure. 2005, 34: 351-378. 10.1146/annurev.biophys.34.040204.144637.

    Article  CAS  PubMed  Google Scholar 

  19. Ziff RM, Gulari E, Barshad Y: Kinetic phase transitions in an irreversible surface-reaction model. Physical Review Letters. 1986, 56: 2553-2556. 10.1103/PhysRevLett.56.2553.

    Article  CAS  PubMed  Google Scholar 

  20. Slepchenko BM, Schaff JC, Carson JH, Loew LM: Computational Cell Biology: Spatiotemporal Simulation of Cellular Events. Annual Review of Biophysics and Biomolecular Structure. 2002, 31: 423-441. 10.1146/annurev.biophys.31.101101.140930.

    Article  CAS  PubMed  Google Scholar 

  21. Takahashi K, Arjunan SNV, Tomita M: Space in systems biology of signaling pathways – towards intracellular molecular crowding in silico. FEBS Letters. 2005, 579: 1783-1788. 10.1016/j.febslet.2005.01.072.

    Article  CAS  PubMed  Google Scholar 

  22. Lemerle C, Ventura BD, Serrano L: Space as the final frontier in stochastic simulations of biological systems. FEBS Letters. 2005, 579: 1789-1794. 10.1016/j.febslet.2005.02.009.

    Article  CAS  PubMed  Google Scholar 

  23. Saxton MJ: Single-particle tracking: effects of corrals. Biophys J. 1995, 69: 389-398.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  24. Andrews SS, Bray D: Stochastic simulation of chemical reactions with spatial resolution and single molecule detail. Physical Biology. 2004, 1: 137-151. 10.1088/1478-3967/1/3/001.

    Article  CAS  PubMed  Google Scholar 

  25. Goldman J, Andrews S, Bray D: Size and composition of membrane protein clusters predicted by Monte Carlo analysis. European Biophysics Journal. 2004, 33: 506-512. 10.1007/s00249-004-0391-6.

    Article  CAS  PubMed  Google Scholar 

  26. Levin MD, Shimizu TS, Bray D: Binding and Diffusion of CheR Molecules Within a Cluster of Membrane Receptors. Biophys J. 2002, 82: 1809-1817.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  27. Mahama PA, Linderman JJ: A Monte Carlo study of the dynamics of G-protein activation. Biophys J. 1994, 67: 1345-1357.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  28. Shea LD, Omann GM, Linderman JJ: Calculation of diffusion-limited kinetics for the reactions in collision coupling and receptor cross-linking. Biophys J. 1997, 73: 2949-2959.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  29. Shimizu TS, Aksenov SV, Bray D: A Spatially Extended Stochastic Model of the Bacterial Chemotaxis Signalling Pathway. Journal of Molecular Biology. 2003, 329: 291-309. 10.1016/S0022-2836(03)00437-6.

    Article  CAS  PubMed  Google Scholar 

  30. Stiles JR, Jr TMB, Salpeter EE, Salpeter MM: Monte Carlo simulation of neurotransmitter release using MCell, a general simulator of cellular physiological processes. Computational Neuroscience. Edited by: Bower JM. 1998, New York, Plenum, 279-284.

    Chapter  Google Scholar 

  31. Woolf PJ, Linderman JJ: Untangling Ligand Induced Activation and Desensitization of G-Protein-Coupled Receptors. Biophys J. 2003, 84: 3-13.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  32. Woolf PJ, Linderman JJ: An algebra of dimerization and its implications for G-protein coupled receptor signaling. Journal of Theoretical Biology. 2004, 229: 157-168. 10.1016/j.jtbi.2004.03.012.

    Article  CAS  PubMed  Google Scholar 

  33. Zhong H, Wade SM, Woolf PJ, Linderman JJ, Traynor JR, Neubig RR: A Spatial Focusing Model for G Protein Signals. Regulator Of G Protein Signaling (Rgs) Protein-Mediated Kinetic Scaffolding. J Biol Chem. 2003, 278: 7278-7284. 10.1074/jbc.M208819200.

    Article  CAS  PubMed  Google Scholar 

  34. Bellot F, Moolenaar W, Kris R, Mirakhur B, Verlaan I, Ullrich A, Schlessinger J, Felder S: High-affinity epidermal growth factor binding is specifically reduced by a monoclonal antibody, and appears necessary for early responses. J Cell Biol. 1990, 110: 491-502. 10.1083/jcb.110.2.491.

    Article  CAS  PubMed  Google Scholar 

  35. Laude AJ, Prior IA: Plasma membrane microdomains: organization, function and trafficking (Review). Molecular Membrane Biology. 2004, 21: 193-205. 10.1080/09687680410001700517.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  36. Parton RG, Hancock JF: Lipid rafts and plasma membrane microorganization: insights from Ras. Trends in Cell Biology. 2004, 14: 141-147. 10.1016/j.tcb.2004.02.001.

    Article  CAS  PubMed  Google Scholar 

  37. Kusumi A, Sako Y: Cell surface organization by the membrane skeleton. Current Opinion in Cell Biology. 1996, 8: 566-574. 10.1016/S0955-0674(96)80036-6.

    Article  CAS  PubMed  Google Scholar 

  38. Vereb G, Szollosi J, Matko J, Nagy P, Farkas T, Vigh L, Matyus L, Waldmann TA, Damjanovich S: Dynamic, yet structured: The cell membrane three decades after the Singer-Nicolson model. PNAS. 2003, 100: 8053-8058. 10.1073/pnas.1332550100.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  39. Ritchie K, Shan XY, Kondo J, Iwasawa K, Fujiwara T, Kusumi A: Detection of Non-Brownian Diffusion in the Cell Membrane in Single Molecule Tracking. Biophys J. 2005, 88: 2266-2277. 10.1529/biophysj.104.054106.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  40. Roy CL, Wrana JL: Clathrin- and non-clathrin-mediated endocytic regulation of cell signalling. Nature Reviews Molecular Cell Biology. 2005, 6: 112-126. 10.1038/nrm1571.

    Article  PubMed  Google Scholar 

  41. Phair RD, Misteli T: Kinetic Modeling Approaches to In Vivo Imaging. Nature Reviews Molecular Cell Biology. 2001, 2: 898-907. 10.1038/35103000.

    Article  CAS  PubMed  Google Scholar 

  42. Roskoski JR: The ErbB/HER receptor protein-tyrosine kinases and cancer. Biochemical and Biophysical Research Communications. 2004, 319: 1-11. 10.1016/j.bbrc.2004.04.150.

    Article  CAS  PubMed  Google Scholar 

  43. Gillespie DT: Exact stochastic simulation of coupled chemical reactions. Journal of Chemical Physics. 1977, 81: 2340-2361. 10.1021/j100540a008.

    Article  CAS  Google Scholar 

  44. Reese JS, Raimondeau S, Vlachos DG: Monte Carlo algorithms for complex surface reaction mechanisms: Efficiency and accuracy. Journal of Computational Physics. 2001, 173: 302-321. 10.1006/jcph.2001.6877.

    Article  CAS  Google Scholar 

  45. Jakobsson E: Computer simulation studies of biological membranes: progress, promise and pitfalls. Trends in Biochemical Sciences. 1997, 22: 339-344. 10.1016/S0968-0004(97)01096-7.

    Article  CAS  PubMed  Google Scholar 

  46. Benveniste M, Livneh E, Schlessinger J, Kam Z: Overexpression of epidermal growth factor receptor in NIH-3T3- transfected cells slows its lateral diffusion and rate of endocytosis. J Cell Biol. 1988, 106: 1903-1909. 10.1083/jcb.106.6.1903.

    Article  CAS  PubMed  Google Scholar 

  47. Simons K, Toomre D: Lipid rafts and signal transduction. Nature Reviews Molecular Cell Biology. 2000, 1: 31-39. 10.1038/35036052.

    Article  CAS  PubMed  Google Scholar 

  48. Pike LJ: Lipid rafts: bringing order to chaos. J Lipid Res. 2003, 44: 655-667. 10.1194/jlr.R200021-JLR200.

    Article  CAS  PubMed  Google Scholar 

  49. Kusumi A, Sako Y, Yamamoto M: Confined lateral diffusion of membrane receptors as studied by single particle tracking (nanovid microscopy). Effects of calcium-induced differentiation in cultured epithelial cells. Biophys J. 1993, 65: 2021-2040.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  50. Rees AR, Gregoriou M, Johnson P, Garland PB: High affinity epidermal growth factor receptors on the surface of A431 cells have restricted lateral diffusion. The EMBO Journal. 1984, 3: 1843-1847.

    PubMed Central  CAS  PubMed  Google Scholar 

  51. Wiegant FA, Blok FJ, Defize LH, Linnemans WA, Verkley AJ, Boonstra J: Epidermal growth factor receptors associated to cytoskeletal elements of epidermoid carcinoma (A431) cells. J Cell Biol. 1986, 103: 87-94. 10.1083/jcb.103.1.87.

    Article  CAS  PubMed  Google Scholar 

  52. Pralle A, Keller P, Florin EL, Simons K, Horber JKH: Sphingolipid-Cholesterol Rafts Diffuse as Small Entities in the Plasma Membrane of Mammalian Cells. J Cell Biol. 2000, 148: 997-1008. 10.1083/jcb.148.5.997.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  53. Martin-Fernandez M, Clarke DT, Tobin MJ, Jones SV, Jones GR: Preformed Oligomeric Epidermal Growth Factor Receptors Undergo an Ectodomain Structure Change during Signaling. Biophys J. 2002, 82: 2415-2427.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  54. Yu X, Sharma KD, Takahashi T, Iwamoto R, Mekada E: Ligand-independent Dimer Formation of Epidermal Growth Factor Receptor (EGFR) Is a Step Separable from Ligand-induced EGFR Signaling. Mol Biol Cell. 2002, 13: 2547-2557. 10.1091/mbc.01-08-0411.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  55. Moriki T, Maruyama H, Maruyama IN: Activation of preformed EGF receptor dimers by ligand-induced rotation of the transmembrane domain1. Journal of Molecular Biology. 2001, 311: 1011-1026. 10.1006/jmbi.2001.4923.

    Article  CAS  PubMed  Google Scholar 

  56. Gadella TWJ, Jovin TM: Oligomerization of epidermal growth factor receptors on A431 cells studied by time-resolved fluorescence imaging microscopy. A stereochemical model for tyrosine kinase receptor activation. J Cell Biol. 1995, 129: 1543-1558. 10.1083/jcb.129.6.1543.

    Article  CAS  PubMed  Google Scholar 

  57. Sherrill JM, Kyte J: Activation of epidermal growth factor receptor by epidermal growth factor. Biochemistry. 1996, 35: 5705-5718. 10.1021/bi9602268.

    Article  CAS  PubMed  Google Scholar 

  58. Schlessinger J: The epidermal growth factor receptor as a multifunctional allosteric protein. Biochemistry. 1988, 27: 3119-3123. 10.1021/bi00409a002.

    Article  CAS  PubMed  Google Scholar 

  59. Kawamoto T, Sato JD, Le A, Polikoff J, Sato GH, Mendelsohn J: Growth stimulation of A431 cells by epidermal growth factor: identification of high-affinity receptors for epidermal growth factor by an anti-receptor monoclonal antibody. PNAS. 1983, 80: 1337-1341.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  60. Lin CR, Chen WS, Lazar CS, Carpenter CD, Gill GN, Evans RM, Rosenfeld MG: Protein kinase C phosphorylation at Thr 654 of the unoccupied EGF receptor and EGF binding regulate functional receptor loss by independent mechanisms. Cell. 1986, 44: 839-848. 10.1016/0092-8674(86)90006-1.

    Article  CAS  PubMed  Google Scholar 

  61. Hendriks BS, Opresko LK, Wiley HS, Lauffenburger D: Quantitative Analysis of HER2-mediated Effects on HER2 and Epidermal Growth Factor Receptor Endocytosis: Distribution Of Homo- And Heterodimers Depends On Relative Her2 Levels. J Biol Chem. 2003, 278: 23343-23351. 10.1074/jbc.M300477200.

    Article  CAS  PubMed  Google Scholar 

Download references


This work was supported by grants from the US Department of Energy (DE-FG02-05ER25702) and the National Science Foundation (CTS-0312117). KM thanks Abhijit Chatterjee for useful discussions.

Author information

Authors and Affiliations


Corresponding authors

Correspondence to Dionisios G Vlachos or Jeremy S Edwards.

Additional information

Authors' contributions

KM carried out the simulations and drafted the manuscript. DGV and JSE edited the manuscript. All authors participated in the analysis of the data. All authors read and approved the final manuscript.

Authors’ original submitted files for images

Rights and permissions

Open Access This article is published under license to BioMed Central Ltd. This is an Open Access article is distributed under the terms of the Creative Commons Attribution License ( ), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Reprints and permissions

About this article

Cite this article

Mayawala, K., Vlachos, D.G. & Edwards, J.S. Computational modeling reveals molecular details of epidermal growth factor binding. BMC Cell Biol 6, 41 (2005).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: