Computational modeling reveals molecular details of epidermal growth factor binding.

BACKGROUND
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.


RESULTS
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.


CONCLUSION
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.


Background
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 ligandinduced 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 [1][2][3][4][5][6][7]. 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,[13][14][15] 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].
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 nullevent MC algorithm were compared with an ordinary differential equations (ODEs) model with a set of parame-ters 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 µm 2 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 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.

Schematic of simulated microscopic events
Reactions events considered in our model as given in [14]

Figure 2
Reactions events considered in our model as given in [14]. 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 µm 2 and a low diffusivity of 2 × 10 -15 m 2 s -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.

Next in
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.
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 adjust-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 µm 2 , D = 2 × 10 -14 m 2 sec -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 µm 2 ) and D = 2 × 10 -15 m 2 sec -1 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 µm 2 , D = 2 × 10 -14 m 2 sec -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 µm 2 ) and D = 2 × 10 -15 m 2 sec -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. ments (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).

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. 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 µm 2 ) 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 m 2 /s to 2 × 10 -15 m 2 /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,[35][36][37][38]. 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 m 2 /s to 2 × 10 -14 m 2 /s mentioned above, larger changes are observed at high ligand concentration (e.g., 1600 nM) and a receptor density of 125 receptors per µm 2 for a change in diffusivity from 2 × 10 -14 m 2 /s to 2 × 10 -13 m 2 /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.

Discussion
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.

Conclusion
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 nullevent 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  Fig. 2 Transport parameters D monomer 2 × 10 -14 -2 × 10 -15 m 2 sec -1 [49,50] over the entire lattice, are avoided, resulting in further acceleration of simulations.

) 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).
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 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 where B i denotes the set of sites to which diffusion from site i can occur. In our model, this set includes all 4 firstnearest 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 , 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 where is the transition probability of reaction at site i.
For a bimolecular reaction on a square lattice (e.g., recep- tor-receptor dimerization), the transition probability per unit time at selected site i was modeled as: 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) -1 sec -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 where Γ max is a normalization constant to ensure that the selection probability of each event is always less than or equal to 1 and 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 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 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 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 m th event 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 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 2dimensional square lattice with each pixel being 2 nm × 2 nm in size. The number density of receptors ranges from 10 2 receptors per µm 2 on normal cells [46] to ~10 3 receptors per µm 2 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 µm 2 , respectively. The diffusivity of monomer EGFR has been reported to be around 2 × 10 -14 m 2 s -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 diffu-for the reaction: and for the reaction: 2A sivity of 2 × 10 -15 m 2 s -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,[53][54][55][56]. 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 µm 2 , 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].

Appendix
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 where I is the % of high intensity spots (y axis in Fig. 4a) upon perturbing the kinetic parameter, k i , and I o is the nominal value corresponding to the original set of kinetic parameters, k io . Only 8 of the 12 kinetic parameters were independently perturbed because of the two equilibrium constrains reported in [14], i.e., 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. where subscript 'o' denotes the nominal value, and is 1, if the kinetic parameter (k i ) is not perturbed, and 1.2, if the parameter is increased by 20%.

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.