Figures
Abstract
Background
The ErbB family of receptors activates intracellular signaling pathways that control cellular proliferation, growth, differentiation and apoptosis. Given these central roles, it is not surprising that overexpression of the ErbB receptors is often associated with carcinogenesis. Therefore, extensive laboratory studies have been devoted to understanding the signaling events associated with ErbB activation.
Methodology/Principal Findings
Systems biology has contributed significantly to our current understanding of ErbB signaling networks. However, although computational models have grown in complexity over the years, little work has been done to consider the spatial-temporal dynamics of receptor interactions and to evaluate how spatial organization of membrane receptors influences signaling transduction. Herein, we explore the impact of spatial organization of the epidermal growth factor receptor (ErbB1/EGFR) on the initiation of downstream signaling. We describe the development of an algorithm that couples a spatial stochastic model of membrane receptors with a nonspatial stochastic model of the reactions and interactions in the cytosol. This novel algorithm provides a computationally efficient method to evaluate the effects of spatial heterogeneity on the coupling of receptors to cytosolic signaling partners.
Conclusions/Significance
Mathematical models of signal transduction rarely consider the contributions of spatial organization due to high computational costs. A hybrid stochastic approach simplifies analyses of the spatio-temporal aspects of cell signaling and, as an example, demonstrates that receptor clustering contributes significantly to the efficiency of signal propagation from ligand-engaged growth factor receptors.
Citation: Costa MN, Radhakrishnan K, Wilson BS, Vlachos DG, Edwards JS (2009) Coupled Stochastic Spatial and Non-Spatial Simulations of ErbB1 Signaling Pathways Demonstrate the Importance of Spatial Organization in Signal Transduction. PLoS ONE 4(7): e6316. https://doi.org/10.1371/journal.pone.0006316
Editor: Henrik Jönsson, Lund University, Sweden
Received: February 10, 2009; Accepted: June 17, 2009; Published: July 23, 2009
Copyright: © 2009 Costa et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Funding: The research was partially supported by the Department of Energy (DE-FG02-05ER25702). MC was supported by an NSF IGERT fellowship. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Competing interests: The authors have declared that no competing interests exist.
Introduction
The ErbB family of receptors, under normal physiological conditions, regulate key cellular processes such as growth, proliferation and differentiation [1], [2], [3]. Overexpression of these receptors deregulates normal cellular function and is a contributing factor to tumorigenesis [4]. There are four members of the ErbB family (ErbB1, ErbB2, ErbB3 and ErbB4) and each family member has its own unique ligand specificity [5], kinase activity [2] and spatial organization on the membrane [1], [6]. In our current study, we have focused solely on the epidermal growth factor receptor (typically abbreviated ErbB1 or EGFR) and the ErbB1 activation of ERK, which is a mitogen activated protein kinase [7]. Ligand binding to ErbB1 stabilizes a conformation of the extracellular domain that allows receptor dimerization [8]. Dimerized receptors are active tyrosine kinases, capable of transautophosphorylation [8]. Phosphorylation of receptor cytoplasmic tails results in recruitment of SH2-containing adaptor and signaling proteins, such as Grb2, Sos, and Shc, that form a signaling scaffold to activate ERK [9].
Due to the importance of the ErbB1-activated ERK pathway, several ordinary differential equation (ODE) models have been developed to gain insight into this pathway [10], [11], [12], [13]. While ODE models have provided insight into the dynamics of this pathway, these models assume that the cell is a homogeneous well-mixed system. In other words, the ODE models neglect spatial localization and organization, such as membrane receptor clustering [3], [14]. Over the past decade, ODE models of the ErbB1-induced ERK pathway have evolved in complexity, becoming both larger and having more experimentally constrained parameters [15]. The first ErbB1/EGFR model was introduced in 1996 and had 35 reactions [16], whereas the most complete models available contain hundreds of reactions [10], [15].
The question remains whether these well-mixed deterministic models are capable of quantitatively describing the temporal dynamics of signaling, since there is significant evidence that cell membrane organization promotes the formation of localized “signaling platforms” [17], [18], [19], [20]. Major advances in our understanding of the membrane have led to a revision of the original Fluid Mosaic model (Singer and Nicholson, 1972), to a more ordered structure with distinct membrane microdomains of lipids and proteins [21], [22], [23] Advanced microscopy techniques have demonstrated that membrane properties, such as transient confinement zones and corrals, may restrict and govern the spatial-temporal dynamics of lipids and membrane proteins [24], [25], [26], [27], [28], [29]. The challenge is to develop computational approaches that can account for membrane spatial heterogeneity and evaluate the impact on signal propagation.
Spatial modeling has been implemented in many scientific disciplines, including physics, material sciences, chemistry, engineering and biological systems. However, the modeling methodologies used vary, with typical approaches including partial differential equations [30], agent-based modeling [31] and spatial Monte Carlo (MC) methods [32], [33], [34]. Spatial MC platforms are particularly powerful numerical simulation tools for studying the dynamics of membrane components [35], [36], [37], [38]. The application of spatial MC methods has been implemented by our group [36] to study ErbB reaction/diffusion and herein to study the effect of spatial heterogeneity on signal propagation. We report the development of a new computational framework that merges a spatial kinetic Monte Carlo (SKMC) algorithm for modeling reaction and diffusion events on the membrane with a stochastic simulator algorithm (SSA) [39] for modeling cytosolic reactions. This new algorithm, the Coupled Spatial and Non-spatial Simulation Algorithm (CSNSA), has enabled us to determine the effects that receptor clustering has on the initiation of signaling.
Results
Establishing Parameters for the Spatial Model
One goal of our study was to evaluate whether simulation results from a spatial stochastic model would differ significantly from a deterministic solution that assume all components are well-mixed. As a starting point, we began with the original ODE model developed by Kholodenko and colleagues [12]. We noted, however, that the ODE model produced results that deviated from the same group's experimental data [12]. We performed a sensitivity analysis to identify the most important enzymatic reaction parameters in the system. Based upon this analysis, we determined that incorporation of receptor degradation mechanisms results in a better fit to the experimental data (Figure 1A) and we fit the new parameters using the PottersWheel MatLab toolbox [40]. Additional reactions added during our model development are illustrated in blue within Figure 1B and the entire set of reaction parameters are summarized in Table 1. Our model modifications are consistent with other models that include negative feedback reactions [10], [11], [13]. In addition, it is noteworthy that the new parameters fit using the ODE model were not explicitly dependent on receptor diffusion. Appendix S1 describes our analytical approach to demonstrate the validity of this fit.
A) Optimization of modeling parameters based upon sensitivity analysis and ODE solution. Green line: Kinetics of Shc phosphorylation in EGF-stimulated hepatocytes (20 nM EGF) as determined by Kholodenko et al. [12]. Red line: results obtained using the ODE model of [12]. Blue line: improved fit of ODE solution to experimental data after incorporation of receptor degradation reactions. B) Summary of reaction network in the ODE and CSNSA models. Note that, in the spatial CSNSA model, stars mark membrane reactions handled by the spatial stochastic Monte Carlo algorithm. All remaining reactions are governed by the Gillespie algorithm. Additional reactions that were added to the original ODE model from Kholodenko et al. [12] are shown in blue. Numbering of reactions is arbitrary.
Validating the CSNSA hybrid approach
The novelty of the CSNSA approach lies in its computationally efficient framework that considers receptor diffusion and reaction in the 2-dimensional confines of the plasma membrane, while cytosolic reactions occur stochastically under well-mixed conditions. The simulated space is illustrated in Figure 2, with a full description of the CSNSA algorithm in the Methods section below. As an initial test, results were compared with the ODE solution (as described in Figure 1) and the experiment results in Kholodenko et al [12]. The simulation space was populated with an initial random distribution of receptor at a density of 141 receptors per µm2, each diffusing at 1×10−14 m2s−1 [41]. In both ODE and CSNSA models, reactions were initiated by addition of EGF ligand (20 nM). Results show that, when receptors are randomly distributed, the two approaches give similar results for the rate and extent of ErbB1 phosphorylation and for the recruitment of PLCγ (Figure 3). The CSNSA model predicts a slightly lower peak value and less sustained recruitment of Shc (Figure 3) when compared to the ODE solution. These results emphasize that the CSNSA hybrid stochastic model is comparable to deterministic solutions in the absence of local concentration gradients or membrane inhomogeneities.
The CSNSA model incorporates a Monte Carlo approach to simulate receptor diffusions and interactions on the cell membrane and couples to a spatial stochastic algorithm (Gillespie) for all cytosol interactions.
Simulated kinetics of ErbB1 phosphorylation (A), PLCγ recruitment (B) and Shc phosphorylation after EGF (20 nM) using the ODE model (dashed lines) or the CSNSA model (solid black line). Results (A,B) from both simulation methods compare well with experimental data (solid circles) reported by Kholodenko et al. [12]
Predicting the Impact of Receptor Density vs. Clustering
We next used the CSNSA to determine the effects of receptor spatial distribution and density on downstream signaling. We defined three different conditions, as shown in the schematic of Figure 4. In the first condition (magenta), the simulation space contained a modest density of dispersed receptors (106 receptors per µm2). In the second condition (dark blue), the simulation space contained a high density of well dispersed receptors (705 receptors per µm2). The final simulation condition (cyan) began with a dense cluster of receptors, which was initially confined to a central region of 705 receptors per µm2 and then permitted to diffuse over time to encompass the entire simulation space for a final density of 106 receptors per µm2. For each regime we examined how initial receptor density and clustering conditions influenced coupling to four of ErbB1's adaptor proteins. The temporal profiles of the cytosolic species Grb2, Sos, and pShc and membrane-bound PLCγ are shown in Figure 4B–E.
A) Schematic illustration of three simulation cases: dispersed (left), high-receptor density (middle), and highly clustered (right). See legend for key to colored lines in each plot. Results predict the kinetics of Grb2 activation (B), PLCγ phosphorylation (C), Shc phosphorylation (D) and Sos activation (E). Active Grb2 is equivalent to: RGrb2+RGrb2Sos+RpShcGrb2+RpShcGrb2Sos+Grb2Sos+pShcGrb2+pShcGrb2Sos; Total phosphorylated PLCγ = RpPLCγ+pPLCγ+pPLCγI; total phosphorylated Shc = RpShc+RpShcGrb2+RpShcGrb2Sos+pShc+pShcGrb2+pShcGrb2Sos; total Sos RGrb2Sos+RpShcGrb2Sos+Grb2Sos+pShcGrb2Sos.
All temporal profiles of the CSNSA were compared with their ODE solutions (shown in purple and red). The most notable differences came from the clustered regime (cyan), which had the same receptor concentration of 106 receptors per µm2 as the non-clustered regime (magenta) but was initially confined to a smaller region. The clustered regime showed a marked increase in the amplitude of signal propagation in comparison to the ODE solution. These data demonstrate that spatial models are needed to accurately predict the consequence of membrane heterogeneity on signal propagation and set the stage for more refined considerations of signaling platforms.
Discussion
In this work, we describe a new, efficient computation framework for evaluating the contributions of spatial organization to important cellular processes. Although applied here to study ErbB1 signal initiation at the plasma membrane, the algorithm should be readily adaptable to other receptor systems, organelle sites and biochemical cascades. We show that, when considering well-mixed systems, solutions obtained using the CSNSA hybrid model and the more traditional ODE solutions are comparable. However, given the growing evidence for membrane compartmentalization at both the plasma membrane and internal organelles [6], [42], [43], we propose that the spatial stochastic model will more accurately predict the outcomes of events that take place between membrane proteins and lipids and their cytosolic binding partners.
As an example, we used CSNSA to demonstrate that receptor clustering creates a more efficient signaling environment. The existence of receptor clusters is well established [23], [44], [45], but the significance of this membrane organization has been approached in only a few recent publications [31], [46]. Our previous work concluded that ligand-independent ErbB1 dimerization is likely to be dependent on two factors: density and the probability of receptor “fluxing” from a closed (dimerization-incompetent) to an open (dimerization-competent) conformation [31], [47]. Because clustering creates locally high receptor concentrations, it enhances the probability for collision between receptors that are transiently in the conformationally “open” state [31]. Here, we show that ErbB1 clustering also enhances the signaling output of receptors, based upon the more efficient recruitment of PLCγ1, Grb2, Sos and Shc.
The importance of spatial effects is emerging as an important topic in systems biology, with technologies such as single particle tracking and electron microscopy demonstrating unique spatial domains [25], [26], [48], [49], [50], [51], [52]. In this work, we applied a novel algorithm to show a direct link between spatial heterogeneity and downstream signaling. We propose that future studies of receptor signaling should seek to gain a fundamental understanding of the spatial interactions and spatial organization of the receptors and to apply these concepts to predictions of signaling output. ErbB receptor clustered domains have been observed in many cancers using different microscopy techniques [6], [44]. Understanding this bigger picture of spatial-temporal protein interactions will drive forth knowledge of cell signaling events and offer the potential to lead towards better drug treatment options.
Methods
Coupled Spatial, Non-spatial Simulation Algorithm (CSNSA)
The Coupled Spatial Non-spatial Simulation Algorithm, CSNSA, is a hybrid model that considers the diffusive behavior and organization of receptors and other membrane components within a 2-D framework, bordered by a well-mixed cytosol. A spatial kinetic Monte Carlo algorithm was employed to capture the spatial-temporal dynamics of receptors on the cell membrane [36] (Figure 5); we used a null-event algorithm that allows ease of implementation and variation of the underlying model. For computational simplicity, the cytosol is treated as a well-mixed solution and modeled with the stochastic simulation algorithm of Gillespie [39]. This assumption is reasonable in the cytosol, given that the diffusivity of proteins in the cytosol (1×10−10 m2s−1) [53] is four orders of magnitude larger than that in the plasma membrane (1×10−14 m2s−1) [41].
This algorithm differs from the original algorithm of Mayawala et al [46] in the time update, which occurs recursively until a successful event is selected. Time is not updated when a null event occurs. A detailed description is provided in the text.
The two algorithms are coupled using the CSNSA, which employs a novel algorithm that selects and executes reactions that allow the molecular species to evolve in space and time. The coupling method takes into account the stochastic nature of biological systems. The first step of the CSNSA is to select a spatial domain (cell membrane or cytosol) and thus the corresponding algorithm for the next event. The selection is made by computing the probabilities of a membrane (SKMC) event or a cytosolic (SSA) event, which are calculated as:where
is defined as,
The total transition rate for the SKMC, , is the sum of all transition rates for all SKMC events, or more specifically the transition rate for diffusion (
) and the sum of the reaction events (
) for all
reaction types,
, where
is the total transition rate for each reaction type defined over all lattice sites
,
.
is defined as the sum of the transition diffusion rate
over all lattice sites
,
. Thus,
is defined as:
The SSA only accounts for stochastic variations in species populations and does not consider the spatial organization in the cytosol, and therefore does not contain a diffusion term. The is defined as the sum of
over all reaction types,
.
The combined MC method operates like a single MC method by considering the superposition of all processes. Time is updated in a “combined” manner from with an average time step as,
. Given that the two algorithms are different (null-event vs. rejection free), the CSNSA is a hybrid method. In order to properly match time scales, upon selection of a spatial event, the SKMC model is continuously executed until a successful event is selected, as shown in Figure 6, based on probability theory described in [33]. The complete algorithm, which is shown in Figure 7, was implemented in Fortran 90. Since the algorithm is stochastic, 10 simulations with different seeds for the random number generator were used. The CSNSA was benchmarked by comparison to an ODE model in a reaction-limited system, where the diffusion coefficient in the CSNSA was made fast compared to the reaction rates (Figure 4). The typical CPU time for 50 receptors/lattice is ∼15 min, for 125 receptors/lattice is ∼2880 min, and for 500 receptors/lattice is ∼14400 min on an Intel® Xeon™ CPU 3.2 GHz processor with 8.00 GB of Ram.
Coupled Spatial Nonspatial Simulation Algorithm, CSNSA, combines the spatial stochastic algorithm [39] depicted in the right branch, with the spatial kinetic Monte Carlo algorithm [56] in the left branch. Upon selection of a branch, a successful event has been executed, species populations are updated, transition rates and probabilities are recomputed, and time advances. The CSNSA is described in greater detail within the text.
This algorithm is used for all cytosolic interactions. Being a rejection free algorithm, a successful event (reaction) is chosen and executed in each iteration.
Spatial Kinetic Monte Carlo (SKMC)
Once an algorithm is selected and executed, transition probabilities are computed again at each time step. Computing involves computing the
values for the SKMC over the entire lattice. This computation is the most CPU intensive step in the simulation algorithm. We, therefore, used an optimized computation method. In order to maximize efficiency, a local region that is affected by the previous reaction event is defined [36], and the
for each lattice site is computed for this region both before and after the event has been executed. This eliminates scanning the entire lattice before and after an event is implemented, and the new
is calculated by:
where,
is the total transition probability computed initially or at a previous successful MC event,
is the sum of transition probabilities of all sites affected by an executed event based on the old configuration, and
is the sum of transition probabilities of all sites affected by an executed event based on the new configuration.
The SKMC algorithm is a modified null-event lattice MC method; for further details see Mayawala et al. [36]. All reactions that are on the lattice or reacting with a species on the lattice are handled by the SKMC (see Figure 2, * denotes membrane reactions and ` denotes interfacial reactions). These reactions include ligand association and dissociation, receptor dimerization and decomposition, receptor phosphorylation and dephosphorylation, and phosphorylated receptor associating with and disassociating from cytosolic species. When an interfacial reaction occurs, a molecule of the cytosolic species is subtracted from the cytosolic population and the membrane species is converted to a new species at the same location on the lattice.
The spatial domain is a two-dimensional lattice with periodic boundary conditions. The initial condition of the lattice is dependent on user specifications and can either be randomly populated or clustered in pre-defined domains. The algorithm is implemented by selecting an occupied lattice site, choosing a successful (reaction or diffusion) or unsuccessful (null) event based on the probabilities, and if a successful event was chosen, executing the event.
An event is selected by computing the probability distribution for all events, defined as: , for lattice site i and event x. Table 2 shows the events executed by this algorithm and the equations for computing ΓX for each event. Γmax is defined as
where the multiple of four accounts for events occurring in each of the four directions on the square lattice.
The spatial algorithm is coupled with the Stochastic Simulation Algorithm (SSA); therefore, unlike the original SKMC algorithm [36], the CSNSA is recursive in that it continuously selects an event until a successful event is chosen and executed as shown in Figure 6; therefore time is not updated if an unsuccessful event is selected.
Stochastic Simulation Algorithm (SSA)
The non-spatial SSA developed by Gillespie [39] was used to model protein association reactions in the cytosol. The algorithm begins with initializing species populations and time; then propensities for all reactions are computed, and an event is randomly selected and the time is updated. This is a rejection free method; therefore, a reaction event is chosen and time is updated by an increment whose average is .
Interfacial Reactions
Interfacial reactions occur when a cytosolic species binds to or detaches from a receptor on the square lattice. In the former case, a molecule from the cytosolic species is subtracted from the cytosol population and a new product is produced at the site that was previously occupied by the reacting receptor. In the latter case, the converse procedure occurs. An example is shown in Table 1 (Interfacial Reaction #1), in which a cytosolic species, Shc, binds to a receptor, R, occupying site k producing product R-Shc at site k.
The rate constants for cytosolic reactions are calculated by first computing the cytosolic volume (Vcyt = 1/3 rL2 µm3), where r is the radius of the cell, and L is the lattice dimension. Next we compute the number of molecules per µm3, Nsp. By multiplying the product of Vcyt and Nsp with the rate constant (given in terms of molecules−1 s−1for bimolecular reactions or s−1 for unimolecular reactions), we obtain a transition rate with units of molecules s−1.
Sensitivity Analysis
To elucidate a mechanism that agrees with the experimental results [12] and explains the biological nature of our system, we modified the reaction scheme developed by Kholodenko et al. [12]. A sensitivity analysis was performed on the reaction mechanism, using the decoupled direct method and the backward differentiation formula method, as implemented in the NASA Glenn chemical kinetics and sensitivity analysis code LSENS [54], [55]. In addition to the species concentrations, these methods automatically follow the temporal evolution of the first-order sensitivity coefficients dC/dηj. The vector C contains the concentrations of all biochemical species and η j is a parameter of interest, such as an initial concentration or a rate constant. The parameters of the new system were refined, and fits were performed for the new reactions shown in blue in Figure 1 and for the Michaelis-Menten reactions using PottersWheel. The parameters to refine were determined to be sensitive using the LSENS package.
Acknowledgments
We thank Drs. Kapil Mayawala, Abhijit Chatterjee, Genie Hsieh and Tomas Mazel for helpful discussions on this manuscript. The image in Figure 2 was provided by Ulises Martinez.
Author Contributions
Conceived and designed the experiments: DGV JE. Performed the experiments: MC. Analyzed the data: MC KR DGV JE. Contributed reagents/materials/analysis tools: DGV BW. Wrote the paper: MC DGV JE.
References
- 1. Lajoie P, Partridge EA, Guay G, Goetz JG, Pawling J, et al. (2007) Plasma membrane domain organization regulates EGFR signaling in tumor cells. J Cell Biol 179: 341–356.
- 2. Linggi B, Carpenter G (2006) ErbB receptors: new insights on mechanisms and biology. Trends Cell Biol 16: 649–656.
- 3. Yarden Y, Sliwkowski MX (2001) Untangling the ErbB signalling network. Nat Rev Mol Cell Biol 2: 127–137.
- 4. Britten CD (2004) Targeting ErbB receptor signaling: a pan-ErbB approach to cancer. Mol Cancer Ther 3: 1335–1342.
- 5. Hynes NE, Lane HA (2005) ERBB receptors and cancer: the complexity of targeted inhibitors. Nat Rev Cancer 5: 341–354.
- 6. Yang S, Raymond-Stintz MA, Ying W, Zhang J, Lidke DS, et al. (2007) Mapping ErbB receptors on breast cancer cell membranes during signal transduction. J Cell Sci 120: 2763–2773.
- 7. Santos SD, Verveer PJ, Bastiaens PI (2007) Growth factor-induced MAPK network topology shapes Erk response determining PC-12 cell fate. Nat Cell Biol 9: 324–330.
- 8. Blinov ML, Yang J, Faeder JR, Hlavacek WS (2006) Depicting signaling cascades. Nat Biotechnol 24: 137–138; author reply 138.
- 9. Blinov ML, Faeder JR, Goldstein B, Hlavacek WS (2006) A network model of early events in epidermal growth factor receptor signaling that accounts for combinatorial complexity. Biosystems 83: 136–151.
- 10. Sasagawa S, Ozaki Y, Fujita K, Kuroda S (2005) Prediction and validation of the distinct dynamics of transient and sustained ERK activation. Nature Cell Biology 7: 365–373.
- 11. Schoeberl B, Eichler-Jonsson C, Gilles E, Muller G (2002) Computational modeling of the dynamics of the MAP kinase cascade activated by surface and internalized EGF receptors. Nature Biotechnology 20: 370–375.
- 12. Kholodenko B, Demin O, Moehren G, Hoek J (1999) Quantification of Short Term Signaling by the Epidermal Growth Factor Receptor. The Journal of Biological Chemistry 274: 30169–30181.
- 13. Hendriks BS, Opresko LK, Wiley HS, Lauffenburger D (2003) 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 278: 23343–23351.
- 14. Orton R, Sturm O, Vyshemirsky V, Calder M, Gilbert D, et al. (2005) Computational modelling of the receptor-tyrosine-kinase-activated MAPK pathway. Biochem J 392: 249–261.
- 15. Kiyatkin A, Aksamitiene E, Markevich N, Borisov N, Hoek J, et al. (2006) Scaffolding Protein Grb2-associated Binder 1 Sustains Epidermal Growth Factor-induced Mitogenic and Survival Signaling by Multiple Positive Feedback Loops. The Journal of Biological Chemistry 281: 19925–19938.
- 16. Huang CY, Ferrell JE Jr (1996) Ultrasensitivity in the mitogen-activated protein kinase cascade. Proc Natl Acad Sci U S A 93: 10078–10083.
- 17. Jiang X, Sorkin A (2002) Coordinated traffic of Grb2 and Ras during epidermal growth factor receptor endocytosis visualized in living cells. Mol Biol Cell 13: 1522–1535.
- 18. Bluthgen N, Bruggeman FJ, Legewie S, Herzel H, Westerhoff HV, et al. (2006) Effects of sequestration on signal transduction cascades. FEBS J 273: 895–906.
- 19. Borisov NM, Markevich NI, Hoek JB, Kholodenko BN (2005) Signaling through receptors and scaffolds: independent interactions reduce combinatorial complexity. Biophys J 89: 951–966.
- 20. Kholodenko BN, Sauro HM (2008) Spatio-temporal dynamics of protein modification cascades. SEB Exp Biol Ser 61: 141–159.
- 21. Gallegos AM, Storey SM, Kier AB, Schroeder F, Ball JM (2006) Structure and cholesterol dynamics of caveolae/raft and nonraft plasma membrane domains. Biochemistry 45: 12100–12116.
- 22. Schroeder F, Gallegos AM, Atshaves BP, Storey SM, McIntosh AL, et al. (2001) Recent advances in membrane microdomains: rafts, caveolae, and intracellular cholesterol trafficking. Exp Biol Med (Maywood) 226: 873–890.
- 23. Lillemeier BF, Pfeiffer JR, Surviladze Z, Wilson BS, Davis MM (2006) Plasma membrane-associated proteins are clustered into islands attached to the cytoskeleton. Proc Natl Acad Sci U S A 103: 18992–18997.
- 24. Orr G, Hu D, Ozcelik S, Opresko LK, Wiley HS, et al. (2005) Cholesterol dictates the freedom of EGF receptors and HER2 in the plane of the membrane. Biophys J 89: 1362–1373.
- 25. Andrews NL, Lidke KA, Pfeiffer JR, Burns AR, Wilson BS, et al. (2008) Actin restricts FcepsilonRI diffusion and facilitates antigen-induced receptor immobilization. Nat Cell Biol 10: 955–963.
- 26. Fujiwara T, Ritchie K, Murakoshi H, Jacobson K, Kusumi A (2002) Phospholipids undergo hop diffusion in compartmentalized cell membrane. J Cell Biol 157: 1071–1081.
- 27. Koyama-Honda I, Ritchie K, Fujiwara T, Iino R, Murakoshi H, et al. (2005) Fluorescence imaging for monitoring the colocalization of two single molecules in living cells. Biophys J 88: 2126–2136.
- 28. Ritchie K, Kusumi A (2003) Single-particle tracking image microscopy. Methods Enzymol 360: 618–634.
- 29. Murase K, Fujiwara T, Umemura Y, Suzuki K, Iino R, et al. (2004) Ultrafine membrane compartments for molecular diffusion as revealed by single molecule techniques. Biophys J 86: 4075–4093.
- 30. Mac Gabhann F, Popel AS (2005) Differential binding of VEGF isoforms to VEGF receptor 2 in the presence of neuropilin-1: a computational model. Am J Physiol Heart Circ Physiol 288: H2851–2860.
- 31. Hsieh MY, Yang S, Raymond-Stinz MA, Steinberg S, Vlachos DG, et al. (2008) Stochastic simulations of ErbB homo and heterodimerisation: potential impacts of receptor conformational state and spatial segregation. IET Syst Biol 2: 256–272.
- 32. Chatterjee A, Vlachos D (2005) Temporal acceleration of spatially distributed kinetic Monte Carlo simulations. Journal of Computational Physics 211: 596–615.
- 33. Chatterjee A, Vlachos D (2007) An overview of spatial microscopic and accelerated kinetic Monte Carlo methods. J Computer-Aided Mater Des 14: 253–308.
- 34. Woolf PJ, Linderman JJ (2004) An algebra of dimerization and its implications for G-protein coupled receptor signaling. J Theor Biol 229: 157–168.
- 35. Brinkerhoff CJ, Woolf PJ, Linderman JJ (2004) Monte Carlo simulations of receptor dynamics: insights into cell signaling. J Mol Histol 35: 667–677.
- 36. Mayawala K, Vlachos DG, Edwards JS (2005) Computational modeling reveals molecular details of epidermal growth factor binding. BMC Cell Biol 6: 41.
- 37. Woolf PJ, Linderman JJ (2003) Self organization of membrane proteins via dimerization. Biophys Chem 104: 217–227.
- 38. Woolf PJ, Linderman JJ (2003) Untangling ligand induced activation and desensitization of G-protein-coupled receptors. Biophys J 84: 3–13.
- 39. Gillespie D (1977) Exact Stochastic Simulation of Coupled Chemical Reactions. The Journal of Physical Chemistry 18: 2340–2361.
- 40. Maiwald T, Timmer J (2008) Dynamical modeling and multi-experiment fitting with PottersWheel. Bioinformatics 24: 2037–2043.
- 41. Kusumi A, Nakada C, Ritchie K, Murase K, Suzuki K, et al. (2005) 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. Annu Rev Biophys Biomol Struct 34: 351–378.
- 42. Smith MG, Simon VR, O'Sullivan H, Pon LA (1995) Organelle-cytoskeletal interactions: actin mutations inhibit meiosis-dependent mitochondrial rearrangement in the budding yeast Saccharomyces cerevisiae. Mol Biol Cell 6: 1381–1396.
- 43. Yang HC, Simon V, Swayne TC, Pon L (2001) Visualization of mitochondrial movement in yeast. Methods Cell Biol 65: 333–351.
- 44. Nagy P, Vereb G, Sebestyen Z, Horvath G, Lockett SJ, et al. (2002) Lipid rafts and the local density of ErbB proteins influence the biological role of homo- and heteroassociations of ErbB2. J Cell Sci 115: 4251–4262.
- 45. Wilson BS, Pfeiffer JR, Raymond-Stintz MA, Lidke D, Andrews N, et al. (2007) Exploring membrane domains using native membrane sheets and transmission electron microscopy. Methods Mol Biol 398: 245–261.
- 46. Mayawala K, Vlachos DG, Edwards JS (2006) Spatial modeling of dimerization reaction dynamics in the plasma membrane: Monte Carlo vs. continuum differential equations. Biophys Chem 121: 194–208.
- 47. Ozcan F, Klein P, Lemmon MA, Lax I, Schlessinger J (2006) On the nature of low- and high-affinity EGF receptors on living cells. Proc Natl Acad Sci U S A 103: 5735–5740.
- 48. Hansen M, Prior IA, Hughes PE, Oertli B, Chou FL, et al. (2003) C-terminal sequences in R-Ras are involved in integrin regulation and in plasma membrane microdomain distribution. Biochem Biophys Res Commun 311: 829–838.
- 49. Ritchie K, Iino R, Fujiwara T, Murase K, Kusumi A (2003) The fence and picket structure of the plasma membrane of live cells as revealed by single molecule techniques (Review). Mol Membr Biol 20: 13–18.
- 50. Suzuki K, Ritchie K, Kajikawa E, Fujiwara T, Kusumi A (2005) Rapid hop diffusion of a G-protein-coupled receptor in the plasma membrane as revealed by single-molecule techniques. Biophys J 88: 3659–3680.
- 51. Ritchie K, Shan XY, Kondo J, Iwasawa K, Fujiwara T, et al. (2005) Detection of non-Brownian diffusion in the cell membrane in single molecule tracking. Biophys J 88: 2266–2277.
- 52. Birtwistle MR, Hatakeyama M, Yumoto N, Ogunnaike BA, Hoek JB, et al. (2007) Ligand-dependent responses of the ErbB signaling network: experimental and modeling analyses. Mol Syst Biol 3: 144.
- 53. Morimatsu M, Takagi H, Ota KG, Iwamoto R, Yanagida T, et al. (2007) Multiple-state reactions between the epidermal growth factor receptor and Grb2 as observed by using single-molecule analysis. Proc Natl Acad Sci U S A 104: 18013–18018.
- 54. Radhakrishnan K (1999) LSENS: Multipurpose Kinetics and Sensitivity Analysis Code. AIAA J 41: 848–855.
- 55. Radhakrishnan K (1991) Combustion Kinetics and Sensitivity Analysis Computations,” in Numerical Approaches to Combustion Modeling, ES.Oran and JP Boris. AIAA J 83–128.
- 56. Mayawala K, Vlachos DG, Edwards JS (2005) Heterogeneities in EGF receptor density at the cell surface can lead to concave up scatchard plot of EGF binding. FEBS Lett 579: 3043–3047.