Coupled stochastic spatial and non-spatial simulations of ErbB1 signaling pathways demonstrate the importance of spatial organization in signal transduction

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.

Figure 1. Parameter optimization and summary of reaction network.

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. . Red line: results obtained using the ODE model of . 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. are shown in blue. Numbering of reactions is arbitrary.

Figure 2. Illustration of the simulated space of the cell, consisting of two distinct domains: the cell membrane and the cytosol.

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.

Figure 3. Comparison of the CSNSA and ODE solutions for receptor phosphorylation, PLCγ and SHC recruitment following EGF stimulation.

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.

Figure 4. The spatial model predicts that receptor clustering enhances signaling efficiency by creating locally high receptor densities.

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.

Figure 5. The spatial kinetic Monte Carlo algorithm, as implemented in the CSNSA.

This algorithm differs from the original algorithm of Mayawala et al 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.

Figure 6. Schematic of CSNSA.

Coupled Spatial Nonspatial Simulation Algorithm, CSNSA, combines the spatial stochastic algorithm depicted in the right branch, with the spatial kinetic Monte Carlo algorithm 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.

Figure 7. Schematic of the SSA algorithm, as coupled to the hybrid algorithm.

This algorithm is used for all cytosolic interactions. Being a rejection free algorithm, a successful event (reaction) is chosen and executed in each iteration.