### 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

${\Gamma}_{}^{d}=\frac{4D}{{a}_{}^{2}},\phantom{\rule{0.5em}{0ex}}\left(1\right)$

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*

${\Gamma}_{i\to j}^{d}=\frac{1}{4}{\Gamma}_{}^{d}{\sigma}_{i}(1-{\sigma}_{j})\phantom{\rule{0.5em}{0ex}}\phantom{\rule{0.5em}{0ex}}\phantom{\rule{0.5em}{0ex}}\phantom{\rule{0.5em}{0ex}}\phantom{\rule{0.5em}{0ex}}j\in {B}_{i}\phantom{\rule{0.5em}{0ex}}\left(2\right)$

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 $\frac{1}{4}{\Gamma}_{}^{d}$, 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

$\text{A}\to \text{C,}\phantom{\rule{0.5em}{0ex}}\phantom{\rule{0.5em}{0ex}}\phantom{\rule{0.5em}{0ex}}{\Gamma}_{i}^{r}=k{\sigma}_{i}\phantom{\rule{0.5em}{0ex}}\left(3\right)$

where ${\Gamma}_{i}^{r}$ 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:

$\text{forthereaction:A+B}\to C,\phantom{\rule{0.5em}{0ex}}{\Gamma}_{i}^{r}=\frac{k}{4}{\sigma}_{i}{\sigma}_{j},\phantom{\rule{0.5em}{0ex}}\left(4\right)$

$\text{andforthereaction:2A}\to \text{C,}{\Gamma}_{i}^{r}=\frac{k}{2}{\sigma}_{i}{\sigma}_{j}.\phantom{\rule{0.5em}{0ex}}\left(5\right)$

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

${p}_{i}^{x}=\frac{{\Gamma}_{i}^{x}}{{\Gamma}_{\mathrm{max}\phantom{\rule{0.5em}{0ex}}}},\phantom{\rule{0.5em}{0ex}}\left(6\right)$

where Γ_{max} is a normalization constant to ensure that the selection probability of each event is always less than or equal to 1 and ${\Gamma}_{i}^{x}$ 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

$\begin{array}{l}{\Gamma}_{\mathrm{max}\phantom{\rule{0.5em}{0ex}}}=4\left(\frac{{\Gamma}^{d}}{4}+\mathrm{max}\phantom{\rule{0.5em}{0ex}}\left\{{\displaystyle \sum _{all\phantom{\rule{0.5em}{0ex}}forward\phantom{\rule{0.5em}{0ex}}reaction\phantom{\rule{0.5em}{0ex}}events}{\Gamma}^{r}}\right\}\right),\\ +\mathrm{max}\phantom{\rule{0.5em}{0ex}}\left\{{\displaystyle \sum _{all\phantom{\rule{0.5em}{0ex}}backward\phantom{\rule{0.5em}{0ex}}reaction\phantom{\rule{0.5em}{0ex}}events}{\Gamma}^{r}}\right\}\end{array}$

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

${\Gamma}_{\mathrm{max}\phantom{\rule{0.5em}{0ex}}}={\Gamma}^{d}+4\left(\frac{{k}_{1f}}{2}+\frac{{k}_{2f}}{4}+\frac{{k}_{3f}}{2}\right)+{\displaystyle \sum _{i=4}^{6}{k}_{if}}+{\displaystyle \sum _{i=1}^{6}{k}_{ib}},\phantom{\rule{0.5em}{0ex}}\left(7\right)$

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 $\sum _{x=1}^{m}{p}_{i}^{x}}>r$ 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

$\Delta t=\frac{1}{{\displaystyle \sum _{i=1}^{No.ofsites}\left({\sigma}_{i}{\displaystyle \sum _{j\in {B}_{i}}^{}{\Gamma}_{i\to j}^{d}\left(1-{\sigma}_{j}\right)}\right)}}=\frac{1}{\frac{1}{4}{\Gamma}^{d}{\displaystyle \sum _{occupiedsites}^{}\left({\displaystyle \sum _{j\in {B}_{i}}^{}\left(1-{\sigma}_{j}\right)}\right)}}.\phantom{\rule{0.5em}{0ex}}\left(8\right)$

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 ~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 diffusivity 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–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].