PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
J Phys Chem B. Author manuscript; available in PMC 2013 March 22.
Published in final edited form as:
PMCID: PMC3321252
NIHMSID: NIHMS363864

The membrane environment can promote or suppress bistability in cell signaling networks

Abstract

Many key biochemical reactions that mediate signal transduction in cells occur at the cell membrane, yet how the 2-dimensional membrane environment influences the collective behavior of signaling networks is poorly understood. We study models of two topologically different signaling pathways that exhibit bistability, examining the effects of reduced protein mobility and increased concentration at the membrane, as well as effects due to differences in spatiotemporal correlations between the membrane environment and 3-dimensional cytoplasm. The two model networks represent the distributive enzymatic modification of a protein at multiple sites and the positive feedback-mediated activation of a protein. In both cases we find that confining proteins to a membrane-like environment can markedly alter the emergent dynamics. For the distributive protein modification network, increased concentration promotes bistability through enhanced protein-protein binding, while lower mobility and membrane-enhanced spatiotemporal correlations suppress bistability. For the positive feedback-mediated activation network, confinement to a membrane environment enhances protein activation, which can induce bistability or stabilize a monostable, active state. Importantly, the influence of the membrane environment on signaling dynamics can be qualitatively different for signaling modules with different network topologies.

Keywords: signal transduction, chemical master equation, two-dimensional

Introduction

Diverse cell types detect environmental stimuli by means of protein receptors associated with the cell membrane. This results in signal transduction, a process which translates receptor engagement to a cellular response via a network of biochemical reactions. Many of the earliest signaling events occur at the cell membrane, rather than in the cytoplasm. Although the influence of the cell membrane on signal transduction has been difficult to study by traditional experimental methods, recent experiments have begun to reveal that the spatiotemporal dynamics of molecules at the membrane play an important role in shaping the qualitative features of signaling.15

A mechanistic understanding of how the membrane environment affects the spatiotemporal dynamics of typical signaling reactions is not available. Seeking such an understanding requires confronting interesting theoretical questions that are pragmatically important. At the membrane, diffusion occurs in an effectively two-dimensional environment, the mobility of proteins is low, and recruited proteins can have a higher effective concentration than when in the cytoplasm. Early studies exploring consequences of the membrane on simple reactions were carried out by Adam and Delbrück, who recognized that cytoplasmic proteins could potentially decrease the expected time to find a membrane-bound target protein by first associating with the membrane and then diffusing in two dimensions.6,7 Although subsequent studies have explored effects of concentration and two-dimensional diffusion on signaling reactions at the membrane,815 most theoretical and computational studies of signal transduction have neglected the presence of the membrane, and a comprehensive framework for understanding the spatiotemporal dynamics of signaling at the membrane remains to be established. In particular, it is not understood how stochastic effects, spatiotemporal correlations, and differences in protein concentration and mobility influence signaling modules with different topologies.

In an effort to gain insight into the influence of the membrane on signal transduction, we use computational methods that explicitly account for stochastic fluctuations to study two reaction networks that exhibit bistability, focusing on the roles of diffusion, dimensionality, and confinement on the presence and characteristics of bistability. Bistability is an important feature of many signaling networks, and can be found in signaling pathways that allow cells to make digital choices, such as those relevant to cell division, apoptosis, and T cell activation against foreign antigen.16 Bistable signaling networks can reside in one of two distinct, stable steady states. Switching between the states can occur by external influences that perturb the signaling network (e.g., receptor stimulation) or by spontaneous, fluctuation-driven transitions.

The two reaction networks considered here are (i) the distributive modification of a substrate protein with two phosphorylation sites, and (ii) a network describing the activation of Ras, which includes a positive feedback loop involving the guanine exchange factor called Son of sevenless (SOS). Both of these networks may be relevant in early T cell signaling events that lead to the decision of whether a T cell will activate against a particular peptide. The first network is a caricature of the Src kinase-mediated phosphorylation of multiple ITAMs in the cytoplasmic tail of the T cell receptor. The second is an important early event in T cell signaling as active Ras proteins can initiate numerous downstream signaling pathways. Signaling networks in diverse cell types also have similar signaling modules.

In both networks, stochastic effects, spatial inhomogeneities that influence protein-protein interactions, and protein mobility may play an important role in shaping collective properties of the dynamics. Traditionally, the effects of diffusion in reaction-diffusion systems have been treated by modifying the kinetic rates in well-mixed kinetic equations to account for the time it takes molecules to diffuse into contact with reaction partners.17 This framework, originally formulated for infinite three-dimensional systems, becomes more difficult in two dimensions, although much effort has been put into understanding simple reaction schemes.18,19 The difficulty in two dimensions, using the language of enzyme-substrate binding, arises from the fact that the rate at which substrate molecules arrive at an enzyme due to diffusion cannot keep up with the rate at which the substrate is bound and transformed, leading to a time-dependent effective reaction rate that decreases with time. Existing theoretical methods are not readily applied to finite systems of arbitrary geometry,20 and recent work has demonstrated that spatiotemporal correlations in particle motion not accounted for by these methods can affect the collective behavior of signaling networks.21 Features of the membrane environment are likely to influence effective kinetic rates and correlations in molecular motions, making the qualitative behavior of the signaling networks difficult to capture using traditional methods. Indeed, our results demonstrate that for the two reaction networks we study, the membrane environment can exhibit qualitatively different behavior than a cytoplasmic environment due to the interplay between protein concentration, diffusion, and dimensionality. Moreover, a particular feature that differs between the membrane and cytoplasmic environments can have qualitatively different effects, depending upon the topology of the signaling network.

Methods

Spatial and temporal dynamics of the signaling networks are modeled as discrete-space, continuous-time, stochastic reaction-diffusion processes. Space is discretized into a cubic lattice, with each lattice site having a size commensurate with a few proteins. Each site is assumed to be well mixed, and reactions occur in a stochastic manner between molecules occupying the same lattice site. We assume diffusive motion of proteins on scales larger than a lattice site, and each molecule hops with a rate determined by its diffusion coefficient to a randomly chosen nearest neighbor lattice site. The lattice spacing is 0.01 μm, and reflecting boundary conditions are used. The state of the system is completely specified by the number and type of molecules at each spatial location. Given an initial state of the system, the time evolution of the state of the system is governed by the spatially resolved chemical master equation,22 which provides an equation for the time evolution of the probability distribution over all possible states. Simulations are performed using our previously developed Stochastic Simulation Compiler,23 which exactly samples trajectories for systems of arbitrary geometry using rule-based methods (rule based methods for homogeneous systems are used in Ref. 24).

To explore the effects of confining reactions near the membrane, we simulate reactions in a slab in which one of the dimensions is systematically varied while keeping the other dimensions and the number of molecules fixed. The system dimensions are 1μm × 1μm × l, with l, the “confinement length,” allowed to vary. Confining the system to a single layer of lattice sites gives an effectively two-dimensional system. We additionally vary the diffusion coefficient to mimic the effects of the membrane on diffusion. Typical membrane diffusion coefficients, which reflect membrane-associated protein motion in a cellular environment, range from 0.01 to 0.1 μm2/s,15,17 while typical cytoplasmic diffusion coefficients range from 1 to 10 μm2/s.25

The slab-like system geometry provides a conceptual framework by which to gain insight into how physical features of the membrane environment affect signaling networks. By independently controlling protein mobility and the confinement length, we can make the system either more cytoplasm-like or more membrane-like. Additionally, the confined geometry provides a model for cells such as T cells and B cells, in which the cell nucleus can occupy a large fraction of the cell interior. In such cells, the gap between the nuclear and plasma membranes can be small (< 0.4 μm), leading to slab-like environments in which the spatial extent of one dimension is much smaller than the extent of the other dimensions.

Results

Distributive reaction network

A distributive reaction network is one in which a substrate molecule is catalytically modified multiple times by an enzyme that must unbind from the molecule between each catalytic step. Here we consider the simple example of a substrate protein containing two residues that can each be either phosphorylated or unphosphorylated. A kinase protein catalyzes the phosphorylation of each site, and a phosphatase protein catalyzes the dephosphorylation of each site. We denote the state of the substrate protein by distinguishing between three states: no sites phosphorylated (S0), one of the sites phosphorylated (S1), or both sites phosphorylated (S2). Denoting the kinase by E, the phosphatase by P, and inactive forms of the enzymes by E* and P*, the reaction network considered is:

S0+Ek2k1S0Ek3S1+ES1+Ek5k4S1Ek6S2+ES2+Pk2k1S2Pk3S1+PS1+Pk5k4S1Pk6S0+PEkaE,PkaP

This generic distributive reaction network has been used by Takahashi et al. as a model for a portion of the mitogen-activated protein kinase cascade.21 We choose kinetic parameters to be consistent with their work. We vary the value of ka, which reflects the time required for a kinase to exchange ADP with ATP, over several orders of magnitude as it is unknown for most protein kinases. We find its value to be qualitatively important when considering the effects of the membrane environment, which highlights the importance of experimentally determining its value for protein kinases.

This network model is sufficiently general to provide insight into signaling pathways that contain opposed distributive modification cycles, provided appropriate kinetic parameters are chosen. For example, the phosphorylation of ITAM regions on the ζ-chain of the T cell receptor by the tyrosine kinase Lck involves multiple kinase modifications. It is unknown whether this reaction occurs by a distributive or the alternative processive mechanism, in which an enzyme does not have to unbind from a substrate protein between catalytic events. Additionally, it is unknown how the membrane environment might affect qualitative cellular-level behavior of distributive reaction networks, so principles gleaned from simulations such as those in this paper may prove useful in designing experiments to distinguish whether a phosphorylation cascade occurs via a distributive mechanism.

It was demonstrated by Markevich et al. that a distributive reaction network can exhibit bistable behavior.26 The physical argument is based on the sequestration of enzymes by the substrate molecule, and a necessary condition for the presence of bistability is that the number of substrate molecules must exceed the number of at least one type of enzyme.27 Here we consider the case in which the substrate molecules outnumber each type of enzyme. When bistability is present, in one steady state, fully phosphorylated proteins (S2) outnumber the unphosphorylated proteins (S0), while in the other steady state, the relative numbers are reversed. Consider a system that resides in the steady state with mostly unphosphorylated proteins. Most kinases will be bound and sequestered by S0 proteins. Since there are relatively few S2 proteins, there will be an excess of free phosphatase. Hence, when a kinase creates a singly phosphorylated substrate protein, the new S1 protein is much more likely to bind to a phosphatase than a kinase in a well-mixed system. This provides a driving force to return the protein to the unphosphorylated state, thereby stabilizing the S0 state. An analogous argument holds for the other steady state.

In three dimensions, Takahashi et al. recently showed that a smaller diffusion coefficient can change the characteristics of the distributive reaction network due to enhanced rebinding of substrate protein and enzyme after the first catalytic reaction.21 Rebinding of enzyme and substrate can result in a second catalytic event quickly following the first, leading to an effective processive-like step taking S0 directly to S2. Processive enzyme catalysis cannot exhibit bistable behavior. It was also noted that increasing ka served to weaken features associated with the distributive reaction network. Recent experimental results provide evidence that molecular crowding in mammalian cells, which presumably increases the likelihood of rebinding between proteins, converts the double phosphorylation of MAP kinase proteins from a distributive to a processive mechanism.28

Our focus is on the effects of confining the distributive reaction network to a 2-dimensional membrane-like environment. Features of the membrane include increased concentration and smaller distances between proteins, lower protein mobility, and differences in the properties of random walks in two and three dimensions. Our simulation results parse the consequences of the interplay between these effects, which are difficult to intuit.

In our simulations, we use a typical membrane protein concentration of 100 molecules/μm2 for the substrate protein, using initial conditions of NS2 = NS0 = 50 and NE = NP = 25, with all other species and complexes initially unpopulated. We simulate many independent dynamic trajectories in which the position and state of all proteins are tracked in time. Initially, the molecules are placed uniformly at random in space. Each trajectory corresponds to an experiment with an isolated patch of membrane. Starting simulations with NS2 = NS0 and all enzymes unbound gives an unbiased starting point: for a bistable system, half of the trajectories (on average) should go to each of the steady states.

Confining the distributive network promotes bistability through increased concentration

Results in Figure 1 show the effects of decreasing the confinement length. For each simulation, we hold constant the number of proteins, the diffusion coefficient, and the rate of enzyme activation. The diffusion coefficient is characteristic of cytoplasmic diffusion and the value of ka is of the same order of magnitude as that estimated by recent experiments measuring the half-life of ADP release after phosphorylation for a serine-threonine protein kinase.29 At large l, the distribution of the number of S2 molecules at steady state is unimodal, and as the system is confined, the distribution at first broadens. At sufficiently small l, the distribution becomes bimodal, which implies that a population of isolated experiments under these conditions would partition into two subpopulations, one containing cells with most substrate proteins fully phosphorylated and one containing cells with most substrate proteins unphosphorylated. This trend continues for values of l smaller than those considered in Figure 1, as is demonstrated in the Supporting Information. Stochastic fluctuations at short times determine which state a trajectory reaches.

Figure 1
Decreasing the confinement length promotes the emergence of bistability in the distributive reaction network

A dominant feature as the system becomes confined is that the concentration of molecules increases. This favors reactions involving the binding of enzyme and substrate proteins, leading to a higher fraction of bound enzyme-substrate pairs. This enhances the sequestration effect that is necessary for bistability in distributive reaction networks. This results in the emergence of bistability as the system is confined, with more confinement resulting in less frequent fluctuation-driven transitions between states.

Decreasing protein mobility and/or increasing the rate of enzyme activation suppresses bistability

The membrane environment can also suppress bistability because of smaller diffusion coefficients. A transition from bistable to monostable behavior can be seen by decreasing the diffusion coefficient while holding the confinement length (l) and other parameters fixed. At sufficiently small values of D, the system is monostable, with a unimodal distribution of NS2. Increasing D, when l is sufficiently small, eventually leads to bistable behavior, with larger D giving rise to more separated states that are less likely to undergo fluctuation-driven switching. The overall dependence of the bistability on D and l, with ka = 0.7 s1, is summarized in the “bistability diagram” in Figure 2. The effect of increasing ka, and hence reducing the refractory time of enzymes, is shown in Figure 3. Experimental studies measuring the rate of ADP release for the protein kinase PKA estimate a rate 30 to 50 times larger than ka = 0.7 s1, suggesting that enzyme turnover rates may vary over orders of magnitude, depending on the kinase.30 In general, as ka increases, greater confinement and higher mobility are necessary to support bistability.

Figure 2
Bistability diagram for the distributive reaction network
Figure 3
Approximate dividing lines between monostable and bistable regions for various values of ka

Protein mobility affects features of bistability in part because the transport of molecules influences the reaction kinetics.17 As mobility decreases, effective diffusion-influenced kinetic rates are reduced, leading to a smaller fraction of bound enzyme-substrate pairs. This weakens the sequestration effect necessary for bistability. At sufficiently high mobility, effective kinetic rates approach their well-mixed values, yielding the maximum sequestration effect.

Both protein mobility and the rate of enzyme reactivation (or refractory time of the enzyme) can affect features of bistability through their effect on the rebinding between enzyme and substrate molecules. The half-life of an inactive enzyme associated with ka = 0.7 s1 is approximately 1 s. In this time, the characteristic distance traveled by the enzyme in an unbounded 2-dimensional domain is (4Dt)1/2 = 2 μm. This distance is of the same order as the system size simulated, which implies that the location of the enzyme upon reactivation will be weakly correlated with its position immediately after catalysis. Thus, a processive-like step is unlikely to occur as spatial correlation between the enzyme and substrate molecule does not persist on time scales associated with enzyme reactivation. For ka = 0.7 s1, rebinding plays a minimal role compared with the influence of transport on the effective kinetic rates for the values of D considered. The effects of rebinding become more prominent as ka increases and D decreases, since spatial correlation between the enzyme and substrate molecule is more likely to persist for times associated with the reactivation of the enzyme. For a given value of l and D, the effects of diffusion on the transport of molecules are expected to be similar. This suggests that differences in the line of conditions separating bistable and monostable responses as ka is varied (Figure 3) are due primarily to rebinding effects. Re-binding plays a more prominent role with increasing ka, converting the distributive mechanism to an effectively processive one (which does not support bistability). A larger value of ka, at a given value of l, requires higher protein mobility to minimize rebinding and support bistability, which is similar to the results in Ref. 21.

Decreasing protein mobility enhances rebinding for multiple reasons: (i) the enzyme is more likely to become active before hopping occurs; (ii) once an enzyme is active and in contact with S1, a binding event is more likely to occur; and (iii) if hopping occurs before enzyme activation, the enzyme and substrate protein are closer together on average when the enzyme becomes active. To gain insight into the importance of rebinding in the distributive reaction network, we developed a quantitative estimate of the probability of rebinding and catalysis between one E* and one S1 protein, given they just reacted and occupy the same lattice site (see the Appendix for details). Figure 4 shows that the effects of rebinding are much more pronounced with ka = 700 s1, in both two and three dimensions, and that rebinding is more significant in two dimensions.

Figure 4
Rebinding is more pronounced in two dimensions and for smaller enzyme refractory times (larger ka)

At fixed concentration, shape of the signaling volume influences bistability

To emphasize that the emergence of bistability upon confinement (change in geometry) is not simply a concentration effect, we examined the steady state distribution of NS2 in two slab-shaped systems that have equal volume and contain equal numbers of molecules. However, they differ in shape, with dimensions of 2μm × 2μm × 0.01μm and 1μm × 1μm × 0.04μm. Although both systems have the same concentration of proteins, differences in shape affect the distribution of distances between proteins, with the average distance between neighboring molecules smaller in the second case. This holds since, in general, the expected distance from a molecule to its nearest neighbor is well approximated by setting ρV(r) ~ 1, where V (r) is the volume of the region contained within a sphere of radius r, ρ is the concentration of molecules, and ρV(r) gives the expected number of molecules within the region. For the two systems considered in this section, the average distance between neighboring molecules is larger than the slab thickness, so we can approximate the volume within radius r of a point at the center of the slab as V(r) ~ πr2λ, where λ denotes the slab thickness. Using ρV(r) ~ 1 and solving for r, it follows that the expected distance to a nearest neighbor scales as λ−1/2. Thus, at constant concentration, as the slab thickness increases, the average distance between neighboring molecules decreases. This result holds for confined regions in which the slab thickness is smaller than the average distance between neighboring molecules. At sufficiently large values of λ, the average distance between nearest neighbors approaches a constant value.

For the two systems considered here, the average distance between neighboring molecules in the less confined system is roughly half the distance between molecules in the more confined system. A more complete analysis, presented in the Supporting Information, shows that the average distance from a molecule to its kth nearest neighbor is smaller in the less confined system. The smaller average distance between molecules reduces the effect of diffusion on the effective kinetic rates and reduces the likelihood of rebinding, since molecules must diffuse a smaller distance before they become effectively well mixed. This is reflected in the steady state distributions shown in Figure 5: while bistability is present in both systems, the stable states are more separated when λ = 0.04 μm, with less frequent fluctuation-induced switching between steady states. Notice that, in the case in which concentration is fixed, a difference in shape (less confinement) promotes bistability. This result shows that differences in the time (or distance) required for a protein to encounter other proteins after its last signaling act can result in qualitative differences in the input-output characteristics of signaling modules on cell membranes. In the limit of very high mobility, the systems are well mixed and shape does not affect the dynamics.

Figure 5
The shape of the confining region, at fixed volume and concentration, affects features of the bistability

Positive feedback regulation of Ras activation

Ras is a small G protein (GTPase) that plays an important role in many signaling pathways (e.g., T cell receptor and growth factor receptor signaling).31 Ras is inactive when bound to GDP and active when bound to GTP. In its active state it can enzymatically modify numerous substrate proteins involved in signal transduction, particularly those involved in cell proliferation.32

Ras activity is regulated by the GTPase-activating protein RasGAP and the guanine nucleotide exchange factors RasGRP and SOS. RasGAP, when bound to Ras, enhances the rate of conversion of GTP to GDP, thus promoting the inactivation of Ras. RasGRP and SOS both catalyze the exchange of GDP with GTP, thus activating Ras. SOS has been shown to have both an allosteric and a catalytic pocket that bind to Ras. The activation of Ras leads to more SOS molecules with RasGTP bound to the allosteric pocket, and an enhancement in the rate of catalysis.33 Thus, there is a SOS-mediated positive feedback loop in the Ras activation network, with RasGTP enhancing the production of additional RasGTP. This positive feedback loop can lead to bistability, hysteresis, and digital signaling.34

Ras is a membrane-anchored protein, and RasGRP and SOS are recruited to the membrane by diacylglycerol (DAG) and adaptor proteins, respectively. The degree to which they are recruited to the membrane is regulated by upstream signaling pathways. The importance of the membrane environment in Ras activation by SOS is illustrated by recent experiments in which the rate of Ras activation is greatly enhanced when membrane-bound Ras is used instead of Ras in solution.35 We consider the effects of membrane confinement and diffusion on the following simple representation of the Ras activation network (adapted from Das et al.34). We denote RasGDP by D and RasGTP by T, and we assume the allosteric pocket of SOS must be occupied for catalytic activity.

SOS+Dk2k1SOS·DSOS+Tk4k3SOS·TSOS·D+Dk6k5SOS·D·Dk7SOS·D+TSOS·T+Dk9k8SOS·T·Dk10SOS·T+TGAP+Tk12k11GAP·Tk13GAP+DDAG+GRPk15k14DAG·GRPDAG·GRP+Dk17k16DAG·GRP·Dk18DAG·GRP+T

This network exhibits bistability in appropriate parameter regimes. When the network is bistable, the steady states can be characterized by the number of active Ras molecules. One state is “active” and contains a relatively large number of RasGTP, while the other state is “inactive” and contains relatively few RasGTP. To test for the presence of bistability, trajectories are started from two different sets of initial conditions. In each case, all molecules begin unbound and are uniformly distributed in space. In one case all Ras are in the active form and in the other case all Ras are in the inactive form. If the system is bistable, trajectories starting from the all-RasGDP state will likely fall into the basin of attraction for the inactive state, while trajectories starting from the all-RasGTP state will likely fall into the basin for the active state. If the system does not exhibit bistability, then all trajectories will evolve in time toward a common state. In bistable systems, it is common for one steady state to be more stable than the other. In our stochastic simulations, this is reflected by stochastic switching from the less stable state to the more stable state on the time scale of the simulations. As such, in simulations with steady states of disparate stability, the less stable steady state is only transiently populated by trajectories starting in its basin of attraction.

Membrane confinement of the Ras network and/or decreasing protein mobility promotes Ras activation

As before, we simulate the Ras network in a slab geometry which is confined in one spatial dimension, varying the diffusion coefficient as well as the number of DAG and RasGRP molecules. The simulations are performed with 150 Ras molecules, 80 SOS molecules, 8 RasGAPs, and a variable number of DAG and RasGRP molecules (NDAG = NGRP = 0; NDAG = 6, NGRP = 20; and NDAG = 24, NGRP = 80). The three most confined systems all exhibit bistable behavior (Figure 6), but with l = 0.01μm, trajectories from the inactive state spontaneously switch to the active state, and with l = 0.03μm, trajectories from the active state spontaneously switch to the inactive state. This gives an indication of the relative stability of each steady state. At l = 0.04μm, the system is at the limit of bistability, and most trajectories from the active regime quickly populate the inactive state.

Figure 6
Confining the Ras activation network affects properties of the bistability at fixed protein mobility

These results can be understood in terms of concentration effects: as the system is confined, the concentration of proteins increases, promoting binding between proteins. This is the same principle that resulted in enhanced binding between substrate proteins and enzymes in the distributive reaction network. Here, more concentrated conditions favor the binding of Ras molecules to the allosteric and catalytic pockets of SOS, which enhances the production of RasGTP. At sufficiently high levels, RasGTP can initiate and sustain positive feedback, which is necessary for bistability.

Along with the results for D = 1μm2/s (Figure 6), Figure 7 shows that reducing the mobility of molecules decreases the stability of the inactive state relative to the active state. At high mobility, the system is bistable, with the active state more stable than the inactive state, while at low mobility, the system is monostable and resides in the active state. Lower mobility allows more cells to achieve the highly active state in the bistable system. Two effects promote stability of the active state as the mobility of proteins decreases. The first is that RasGAPs, which have the highest binding rate (k11) to Ras, experience effects of lowering diffusion at higher mobility than other binding reactions; i.e., the fastest reaction becomes diffusion-influenced first. This serves to decrease the effective binding rate between RasGAP and RasGTP, thus increasing the survival time of an active Ras protein. This promotes positive feedback through SOS and can enhance the further production of RasGTP. This result highlights that different reactions are influenced differently upon reduction in protein mobility, and predicting the effects of diffusion on the network requires knowledge of rates of the biochemical reactions. Additionally, as the rate of diffusion decreases, rebinding between recently bound proteins becomes more likely. For example, RasGTP bound to the allosteric site of SOS becomes more likely to rebind, increasing its effective lifetime.

Figure 7
Decreasing the diffusion coefficient stabilizes the active Ras state relative to the inactive Ras state

The average time that RasGTP is bound to the SOS allosteric pocket, per binding event, is τ=k41. Two effects increase the effective lifetime: multiple bindings before hopping and multiple diffusive excursions and returns before the molecules become mixed. We treat the probability of returning as in the distributive reaction network and assume that if RasGTP and SOS diffuse beyond a characteristic distance apart, they are effectively mixed. This distance is set roughly by the expected distance to the nearest RasGDP molecule, which competes for the allosteric pocket of SOS. The number of bindings before hopping is geometrically distributed, as is the number of returns before diffusion beyond a cutoff distance. The expected number of each type of event is

E[numberofbindingsbeforehopping]=pon1ponE[numberofreturnsbeforemixing]=preturn1preturn

Note that here pon is the probability that RasGTP binds SOS before either molecule hops away. Putting these results together, the effective lifetime is

τeffective=(1+pon1pon)τ+pon(preturn1preturn)(1+pon1pon)τ=(1+pon1pon)(1+ponpreturn1preturn)τ

With D = 0.01μm2/s, the effective time (including rebinding) that RasGTP is bound to the allosteric pocket of SOS, compared to the case without rebinding, is 1.9 times longer in two dimensions and 1.3 times longer in three dimensions. Additionally, as mobility decreases, it is more likely for a recently produced RasGTP to stay in the proximity of the SOS molecule for long enough to enhance its binding to the allosteric pocket when RasGDP unbinds. Both of these effects serve to enhance the effects of positive feedback.

An additional feature of interest is that within the monostable regime, the expected time to reach the active steady state from the all-RasGDP initial condition varies non-monotonically with the diffusion coefficient, as is demonstrated by the response times reported in Figure 7. As the diffusion coefficient decreases from 0.1μm2/s, the average time to reach the active state first decreases (faster response time); however, for sufficiently small D, the average time increases (slower response time). Above the crossover diffusivity, decreasing mobility decreases the response time by reducing the effectiveness of RasGAPs and enhancing positive feedback through rebinding effects. Below the crossover diffusivity, the response time is dominated by the time it takes for molecules to diffuse into contact. Hence, reducing D in this regime increases the overall response time. The non-monotonic response time is similar to a result from Ref. 21, in which lowering mobility first decreased response times due to enhanced rebinding but eventually increased response times due to longer molecular encounter times.

The bistability diagram for the Ras activation network is shown in Figure 8 and summarizes the influence of confinement, the diffusion coefficient, and the number of DAG and RasGRP on bistability (previous studies have explored the influence of RasGRP on Ras activation in well-mixed systems34,36). There are three dominant effects at play in the results: (i) confinement increases the concentration of molecules, which enhances protein-protein binding and the effects of positive feedback; (ii) decreasing D enhances production of RasGTP through reduced effective RasGAP activity and enhanced positive feedback effects due to rebinding; and (iii) increasing NGRP enhances the activation of Ras, thus biasing the system toward an active state.

Figure 8
Bistability diagrams for the Ras-SOS reaction network

Discussion

Our results reveal how certain general features of the membrane environment can qualitatively change the dynamics and collective behavior of signaling networks. The specific features we studied, compared with a 3-dimensional cytoplasmic environment, are increased protein concentration, decreased protein mobility, and different spatiotemporal correlations between molecules. The character of the influence of these features can depend upon the topology of the signaling network.

For the distributive reaction network, we find that effects of membrane confinement can either enhance or suppress bistability. Increased concentration at the membrane enhances binding between enzymes and substrate molecules, promoting the sequestration of enzymes and supporting bistability. Protein mobility and spatiotemporal correlations are closely intertwined, with correlations arising from the spatial proximity between molecules at times soon after they first interact and mobility influencing the time scale over which spatial correlations persist. The nature of the spatiotemporal correlations is different in the membrane environment because of lower protein mobilities and differences in the character of 2-dimensional and 3-dimensional random walks. The importance of spatiotemporal correlations is highlighted by the strong dependence of bistability on rebinding between enzymes and substrate proteins in the distributive reaction network. Additionally, protein mobility and the distribution of distances between proteins influence the time required for a molecule to diffuse to a reaction partner. This influences the effective kinetic parameters governing the system, with slower diffusion leading to smaller effective binding rates and relatively fast reactions being influenced before slower reactions upon lowering protein mobility. For example, in the Ras activation network, RasGAP binding experiences the effects of diffusion at the highest mobility, biasing the system toward active Ras as mobility decreases.

The last result also highlights that the effect of changing the same parameter due to membrane confinement can be different, depending upon the topology of the signaling network. Reduced protein mobility due to membrane confinement abrogates bistability for both reaction networks studied, but leads to qualitatively different outcomes. In the distributive reaction network, loss of bistability due to lower protein mobility results in a new monostable state with lower levels of active molecules compared to the more active bistable state. In sharp contrast, in the Ras activation network, loss of bistability due to lower mobility results in a monostable state with levels of active molecules similar to the previously more active bistable state. In one case, lower mobility leads to fewer active molecules, while in the other case, lower mobility promotes activation.

In our analysis of the two signaling networks, we used a lattice-based model in which proteins, considered to be well-mixed within a lattice site, diffuse by hopping between neighboring lattice sites. Other methods include continuous-space, particle-based methods,13,21 which could potentially give different rebinding properties on molecular length scales. The diffusion coefficients we use to model protein motion are typically reported for motion over length scales large compared to protein size, and treating motion on molecular length scales as diffusive may not be appropriate for membrane-bound proteins. Additionally, it is unlikely that the lattice and continuum methods would give qualitatively different behavior in the distributive reaction network, as the values of ka governing the refractory time of the enzymes are typically much smaller than the rate at which proteins hop to nearest neighbor lattice sites. Hence, it is unlikely that the well-mixed assumption of the lattice site leads to an overestimate of rebinding at short times, as might be the case if ka were sufficiently large. Additionally, in the Ras activation network, the two methods are likely to produce qualitatively similar results, since in both cases the time that Ras is effectively bound to the allosteric site of SOS is enhanced in two dimensions.

Studying the dynamics of signaling networks within a slab-like geometry gave physical insight into how the membrane environment can affect signaling networks. It will be of considerable interest to study systems in which molecules confined to a membrane interact with molecules that diffuse within the cytoplasm. When introducing the distributive reaction network, we discussed the phosphorylation of ITAM regions on the cytoplasmic tail of T cell receptors by the kinase Lck. In this biologically motivated network, the substrate molecule and kinase are confined the membrane, but the phosphatase can diffuse within the cytoplasm. In Figure S5 in the Supporting Information, we demonstrate that such a system can exhibit bistability, although it is expected that details of the signaling output will depend on details of the system geometry.

A deep understanding of the role of the membrane in cell signaling will result from synergy between theory, computation, and experiments. Many of our results can be tested by studying model lipid bilayers with well-defined concentrations of tethered proteins and comparing their behavior with experiments in which the proteins interact in solution (such as in a cuvette) and in live cells. Various other effects of the membrane environment, such as those due to cytoskeletal dynamics and receptor orientation,37 remain to be carefully studied. It will be especially important to use theory to disentangle specific effects from general principles in studying these influences of the membrane. Theoretical studies would benefit from techniques to account for diffusion-influenced kinetic parameters in arbitrary geometries, more sophisticated techniques to account for rebinding and spatiotemporal correlations between molecules, and methods to account for the spatially inhomogeneous and noisy environment in which signaling takes place.

Supplementary Material

1_si_001

Acknowledgments

This research was supported by a NIH Director’s Pioneer award (AKC) and NIH grant 1P01AI091580-01 (SMA, AKC, AW, JPR, JTG). We are grateful to Prof. John Kuriyan for important input and discussions.

Appendix: Probability of rebinding and catalysis

In this section, we present a quantitative estimate of the probability of rebinding and catalysis between one E* and one S1 protein, given they just reacted and occupy the same lattice site. These results also apply to the phosphatase reactions and include dependence on the dimensionality of the system, the rate of enzyme activation, and the concentration of molecules.

The basic idea behind the calculations is that once the initial reaction takes place, the enzyme and substrate molecule must remain “sufficiently close” to avoid becoming well-mixed and sequestered by other molecules in the system. We implement this condition by assuming that the enzyme and S1 cannot diffuse beyond a characteristic distance apart, with the distance set by the expected distance to the nearest particle that can bind (and hence sequester) one of the molecules. Let ρ denote the concentration of such “traps” in this section. For simplicity, we assume the lattice is infinitely extended in all directions.

Active enzyme

We first calculate the probability that an active enzyme and S1 molecule that occupy the same lattice site react before diffusing beyond a characteristic distance apart. In a subsequent section, we extend the result to the case in which the enzyme is initially inactive.

Given that the enzyme and substrate molecule occupy the same lattice site, the probability of binding before either molecule hops to a neighbor is

pon=γonγon+γhop=(k4/a3)(k4/a3)+2n(D/a2),

with γ a total rate, k4 the bulk on-rate between E and S1, a the lattice spacing, D the diffusion co-efficient, and n the number of nearest neighbor lattice sites. The probability of binding depends on the dimensionality, since n = 2d, where d denotes the dimensionality. The probability of hopping before binding is phop = 1 − pon. Once the enzyme and substrate protein are bound, the probability of unbinding before reacting is

poff=γoffγoff+γrxn=k5k5+k6,

The probability of reacting before unbinding is prxn = 1 − poff. For the distributive reaction network, poff ≈ 0.10.

Now assume that a hopping event occurs before a reaction takes place and that the enzyme and substrate protein occupy nearest neighbor lattice sites. Let preturn be the probability that, as they undergo random walks, they co-occupy the same lattice site before diffusing beyond the prescribed cutoff distance. We call this a “re-encounter” or “return” event and construct an approximation for preturn in the next section. Treating preturn as known, we can write a formally exact expression for the probability that E and S1, starting from the same lattice site, react before diffusing beyond the cutoff distance. We denote this probability by Pr{react}.

One way to construct Pr{react} is to consider all possible sequences of events leading to a reaction. It is clear that the final two events must be the sequence in which the two molecules bind and then react (this occurs with probability pon prxn). Before this event, any sequence of events leading to E and S1 occupying the same lattice site without the molecules diffusing too far apart is allowed. The two molecules, given they start at the same lattice site, return to the same configuration after binding and then unbinding (this occurs with probability pon poff) or after hopping and then undergoing a diffusive excursion that returns to a common lattice site (this occurs with probability phop preturn). Accounting for all possible sequences ending with catalysis gives the probability of reaction

Pr{react}=(ponprxn)l=0k=0l(lk)(phoppreturn)k(ponpoff)lk=(ponprxn)l=0(phoppreturn+ponpoff)l=ponprxn1phoppreturnponpoff

The following section is devoted to developing an approximation for preturn.

Probability of return

Begin with E and S1 occupying nearest neighbor lattice sites. We wish to estimate preturn, the probability that they eventually occupy the same lattice site without diffusing beyond a fixed distance apart. We begin by making the assumption that space is continuous, with molecules distributed uniformly at random. Using properties of spatial Poisson processes, it is straightforward to show that the expected distance to a nearest “trap” molecule (distributed with concentration ρ) is

E2[r]=12ρE3[r]=Γ(1/3)3((4/3)πρ)1/3(0.55)ρ1/3,

where the subscripts denote two and three dimensions (the dimensions of the concentrations are not explicitly distinguished, but are assumed to be consistent with the dimensionality).

As an approximation to the lattice random walk problem, consider a particle diffusing in continuous space between two concentric, absorbing spherical boundaries. The particle diffuses with diffusion coefficient Dtot = 2D and represents the relative distance between the two molecules. We wish to compute the probability that the particle reaches the inner sphere before it reaches the outer sphere, which corresponds to the enzyme and S1 molecule reaching a common lattice site before diffusing beyond a cutoff distance. Let the radius of the inner sphere be half the lattice spacing, a/2, and let the particle start at position r = a to mimic the nearest neighbor condition on the lattice. We use the expected distance to the closest trap as the radius of the outer sphere. This calculation falls into the category of first passage problems, a thorough review of which can be found in Ref. 38.

Noting the spherical symmetry in the problem, let Φ(r) denote the probability that the diffusing particle, starting at radial position r, reaches the inner boundary before reaching the outer boundary. With this notation, preturn = Φ(a). The function Φ is known to satisfy the Laplace equation, [nabla]2Φ = 0, subject to the boundary conditions Φ(R) = 1 and Φ(R+) = 0 (here, R = a/2 and R+ = Ed[r]).38

It is straightforward to solve these equations in two and three dimensions.38 The return probability can be written

preturn={ln(R+/a)ln(R+/R)d=21a1R+1R1R+d=3

Inactive enzyme

Now consider the full case in which the inactive enzyme E*and S1 co-occupy the same lattice site. The rate of enzyme activation, ka, is an important parameter, since it controls how long the enzyme is expected to remain inactive. The longer it is inactive, the less likely it is for rebinding to occur. The probability that the enzyme becomes active in the time interval from t1 to t2 is

Pr{activationin(t1,t2)}=t1t2dtkaekat=ekat1(1eka(t2t1))

We will use this result later.

We begin by estimating τreturn, the average time for a diffusive excursion in which the two molecules diffuse apart but return together before diffusing beyond the cutoff distance. To do this, we return to the problem of a particle diffusing between two concentric spheres and compute the expected time for the particle to reach the inner sphere, conditioned on not hitting the outer sphere. This time, starting from position r, is denoted by t(r). Following Ref. 38, this quantity satisfies

Dtot2[Φ(r)t(r)]=Φ(r)

subject to the boundary conditions Φ(r)t(r) = 0 for r [set membership] {R, R+}. We have the solutions for Φ from above. In three dimensions, the solution for τreturn = t (a) is

τreturn=(aR)(2R+Ra)6Dtot

In two dimensions, the solution is more complicated but well-behaved:

τreturn=R+2ln(a/R)+R2[1+ln(R+/R)]ln(R+/a)a2[1+ln(R+/a)]ln(R+/R)4Dtotln(R+/R)ln(R+/a)

These equations are strictly increasing as a function of increasing R+ (decreasing ρ).

For each excursion taken by the molecules, we have the probability that they return and the average time for them to return, given they do not diffuse beyond a cutoff distance. Using these results, along with the distribution of times that the enzyme activates, we can estimate the probability that the proteins react before they diffuse apart.

Let Δt = τhop + τreturn, where τhop = (2nD/a2)−1 is the average time for the molecules to hop apart once they co-occupy a site. For most values of ka and D, the average time for the enzyme to become active is much longer than Δt, which implies that the two molecules are likely to diffuse apart many times before they can bind and react. Enzymes that activate earlier than average are more likely to contribute significantly to rebinding, since the molecules must undergo fewer excursions and returns. Once the enzyme is active and occupies the same lattice site as S1, we can use the probability of binding and catalysis from the previous sections.

In the remainder of the calculation, we break time into increments of size Δt. For each increment, we have the probability that the enzyme activates as well as the probability that the enzyme and substrate molecule have remained sufficiently close throughout the process. Let α denote the probability that E and S1 react, given that they start at the same lattice site (computed previously). Then the probability that E* and S1 react can be approximated as

Pr{S1S2}=(1ekaτhop)α+k=0[eka(kΔT)ekaτhop(1ekaΔt)][(preturn)k+1]α

The first term accounts for cases in which the enzyme activates before any hops have occurred. Terms within the sum account for the probability of k + 1 successful returns without diffusing beyond the cutoff distance (second term in square brackets), as well as the probability that the enzyme activated within the time interval from τhop + kΔt to τhop + (k + 1)Δt (first term in square brackets). Summing the geometric series and substituting the expression for α, we obtain the desired result

Pr{S1S2}=[(1ekaτhop)+ekaτhop(1ekaΔt)(preturn1preturnekaΔt)]×(ponprxn1phoppreturnponpoff)

This probability is plotted as a function of ρ, D, ka, and d in Figure 4.

Footnotes

Supporting Information Available

The Supporting Information includes a full description of the kinetic parameters used in the simulations, a discussion of the identification of bistability in the distributive reaction network, the bistability diagrams used to construct Figure 3 in the text, numerical results augmenting the discussion of Figure 5, and additional simulation results for systems with greater confinement and a system with coupled membrane and cytoplasmic regions. This information is available free of charge via the Internet at http://pubs.acs.org.

References

1. Groves JT, Kuriyan J. Nat Struct Mol Biol. 2010;17:659–665. [PMC free article] [PubMed]
2. Fooksman DR, Vardhana S, Vasiliver-Shamis G, Liese J, Blair D, Waite J, Sacristán C, Victora G, Zanin-Zhorov A, Dustin ML. Annu Rev Immunol. 2010;28:79–105. [PMC free article] [PubMed]
3. Kholodenko BN, Hancock JF, Kolch W. Nat Rev Mol Cell Biol. 2010;11:414–426. [PMC free article] [PubMed]
4. Singleton KL, Roybal KT, Sun Y, Fu G, Gascoigne NRJ, van Oers NSC, Wülfing C. Sci Signal. 2009;2:ra15. [PMC free article] [PubMed]
5. Lee K, et al. Science. 2003;302:1218–1222. [PubMed]
6. Adam G, Delbrück M. In: Structural Chemistry and Molecular Biology. Rich A, Davidson N, editors. Freeman: San Francisco; 1968.
7. McCloskey MA, Poo MM. J Cell Biol. 1986;102:88–96. [PMC free article] [PubMed]
8. Haugh JM, Lauffenburger DA. Biophys J. 1997;72:2014–2031. [PubMed]
9. Shea LD, Omann GM, Linderman JJ. Biophys J. 1997;73:2949–2959. [PubMed]
10. Kholodenko BN, Hoek JB, Westerhoff HV. Trends Cell Biol. 2000;10:173–178. [PubMed]
11. Haugh JM. Biophys J. 2002;82:591–604. [PubMed]
12. Elf J, Ehrenberg M. Syst Biol. 2004;1:230–236. [PubMed]
13. Monine MI, Haugh JM. Biophys J. 2008;95:2172–2182. [PubMed]
14. Alonso S, Bär M. Phys Biol. 2010;7:046012. [PubMed]
15. Dushek O, van der Merwe PA, Shahrezaei V. Biophys J. 2011;100:1189–1197. [PubMed]
16. Ferrell JE. Curr Opin Cell Biol. 2002;14:140–148. [PubMed]
17. Lauffenburger DA, Linderman JJ. Receptors: Models for Binding, Trafficking, and Signaling. Oxford University Press; USA: 1996.
18. Torney DC, McConnell HM. Proc R Soc London, Series A. 1983;387:147–170.
19. Szabo A. J Phys Chem. 1989;93:6929–6939.
20. Bénichou O, Chevalier C, Klafter J, Meyer B, Voituriez R. Nat Chem. 2010;2:472–477. [PubMed]
21. Takahashi K, Tanase-Nicola S, ten Wolde PR. Proc Natl Acad Sci USA. 2010;107:2473–2478. [PubMed]
22. Hattne J, Fange D, Elf J. Bioinformatics. 2005;21:2923–2924. [PubMed]
23. Lis M, Artyomov MN, Devadas S, Chakraborty AK. Bioinformatics. 2009;25:2289–2291. [PMC free article] [PubMed]
24. Faeder JR, Blinov ML, Hlavacek WS. Methods Mol Biol. 2009;500:113–167. [PubMed]
25. Luby-Phelps K. Int Rev Cytol. 1999;192:189–221. [PubMed]
26. Markevich NI, Hoek JB, Kholodenko BN. J Cell Biol. 2004;164:353–359. [PMC free article] [PubMed]
27. Ortega F, Garcés JL, Mas F, Kholodenko BN, Cascante M. FEBS J. 2006;273:3915–3926. [PubMed]
28. Aoki K, Yamada M, Kunida K, Yasuda S, Matsuda M. Proc Natl Acad Sci USA. 2011;108:12675–12680. [PubMed]
29. Keshwani MM, Harris TK. J Biol Chem. 2008;283:11972–11980. [PubMed]
30. Zhou J, Adams JA. Biochemistry. 1997;36:15733–15738. [PubMed]
31. Genot E, Cantrell DA. Curr Opin Immunol. 2000;12:289–294. [PubMed]
32. Campbell SL, Khosravi-Far R, Rossman KL, Clark GJ, Der CJ. Oncogene. 1998;17:1395–1413. [PubMed]
33. Margarit S, Sondermann H, Hall BE, Nagar B, Hoelz A, Pirruccello M, Bar-Sagi D, Kuriyan J. Cell. 2003;112:685–695. [PubMed]
34. Das J, Ho M, Zikherman J, Govern C, Yang M, Weiss A, Chakraborty AK, Roose JP. Cell. 2009;136:337–351. [PMC free article] [PubMed]
35. Gureasko J, Galush WJ, Boykevisch S, Sondermann H, Bar-Sagi D, Groves JT, Kuriyan J. Nat Struct Mol Biol. 2008;15:452–461. [PMC free article] [PubMed]
36. Riese MJ, Grewal J, Das J, Zou T, Patil V, Chakraborty AK, Koretzky GA. J Biol Chem. 2011;286:5254–5265. [PubMed]
37. Wu Y, Vendome J, Shapiro L, Ben-Shaul A, Honig B. Nature. 2011;475:510–513. [PMC free article] [PubMed]
38. Redner S. A Guide to First-Passage Processes. Cambridge University Press; 2001. This material is available free of charge via the Internet at http://pubs.acs.org/