This Article Abstract Full Text (PDF) supplemental All Versions of this Article: biophysj.105.061564v1 89/4/2806    most recent Alert me when this article is cited Alert me if a correction is posted Services Similar articles in this journal Similar articles in PubMed Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via HighWire Citing Articles via Google Scholar Google Scholar Articles by Skupsky, R. Articles by Nossal, R. J. Search for Related Content PubMed PubMed Citation Articles by Skupsky, R. Articles by Nossal, R. J.

Distinguishing Modes of Eukaryotic Gradient Sensing

R. Skupsky ^* , W. Losert  and R. J. Nossal ^*

^* Laboratory of Integrative and Medical Biophysics, National Institute of Child Health and Human Development, National Institutes of Health, Bethesda, Maryland 20892; and Department of Physics, University of Maryland, College Park, Maryland 20742

Correspondence: Address reprint requests to Ron Skupsky, E-mail: skup{at}helix.nih.gov or Ralph Nossal, E-mail: nossalr{at}mail.nih.gov.

	ABSTRACT

TOP ABSTRACT INTRODUCTION MODEL DEVELOPMENT NUMERICAL TECHNIQUES MODEL CHARACTERIZATION RESULTS SUMMARY AND DISCUSSION SUPPLEMENTARY MATERIAL ACKNOWLEDGEMENTS REFERENCES

 
We develop a mathematical model of phosphoinositide-mediated gradient sensing that can be applied to chemotactic behavior in highly motile eukaryotic cells such as Dictyostelium and neutrophils. We generate four variants of our model by adjusting parameters that control the strengths of coupled positive feedbacks and the importance of molecules that translocate from the cytosol to the membrane. Each variant exhibits a qualitatively different mode of gradient sensing. Simulations of characteristic behaviors suggest that differences between the variants are most evident at transitions between efficient gradient detection and failure. Based on these results, we propose criteria to distinguish between possible modes of gradient sensing in real cells, where many biochemical parameters may be unknown. We also identify constraints on parameters required for efficient gradient detection. Finally, our analysis suggests how a cell might transition between responsiveness and nonresponsiveness, and between different modes of gradient sensing, by adjusting its biochemical parameters.

	INTRODUCTION

TOP ABSTRACT INTRODUCTION MODEL DEVELOPMENT NUMERICAL TECHNIQUES MODEL CHARACTERIZATION RESULTS SUMMARY AND DISCUSSION SUPPLEMENTARY MATERIAL ACKNOWLEDGEMENTS REFERENCES

 
Many eukaryotic cells respond with directional movement to spatial and/or temporal gradients of small molecules that bind to cell surface receptors. This process, called chemotaxis, is important in phenomena as diverse as the immune response of higher animals, wound healing, neuronal patterning, vascular, and embryonic development (1–5), as well as the food gathering and social behavior of some ameboid cells (6). Gradient sensing is the aspect of chemotaxis whereby cells transduce the external distribution of chemotactic ligand into an internal distribution of signaling molecules that effect the morphological and mechanical changes necessary for movement (7). In this article we focus on gradient sensing in highly motile cells such as neutrophils and the aggregating slime-mold, Dictyostelium discoideum. These cells respond with directional movement to differences in ligand concentration of a few percent across their length with cellular velocities of order of 10 µm/min (approximately one cell length per minute), over several orders of magnitude in absolute concentration (8,9).

Phosphoinositides in gradient sensing
Recent experiments in both Dictyostelium and neutrophils have suggested that phosphoinositide (PI) signaling at the plasma membrane mediates gradient sensing in these cells by localizing molecules that relay signals from receptor activation to cytoskeletal rearrangements (10,11). PIs are signaling lipids that are phosphorylated by kinases, and dephosphorylated by phosphatases, at different positions on their inositol headgroup. Depending on their phosphorylation state, PIs can recruit specific molecules from the cytosol to the membrane, including those that affect cellular movement and those that affect their own interconversion (12,13). The feedback regulation and membrane localization of PIs make them well suited to mediate cell surface processes that require highly localized and amplified signals in space and time, such as chemotaxis. In particular, 3'PIs (PIs that are phosphorylated in the 3' position) are thought to play a causal role in nucleating the actin-based protrusions necessary for cell movement, and markers for their production show qualitatively similar dynamics to actin polymerization in chemotaxing cells (14–16). These markers show similar dynamics as well in cells that are round and cannot form protrusions due to treatment with actin depolymerizing agents (17,18), suggesting that aspects of gradient sensing may be decoupled from the morphological and mechanical events involved in chemotaxis.

Recent models relating to PI-mediated gradient sensing
The suggestion that gradient sensing can be decoupled from motility, and that it is mediated by a feedback scheme such as those implicated in PI signaling, has inspired several recent mathematical models. Each accounts for characteristic behaviors in a different way. Levchenko and Iglesias (21) have analyzed a general model that maps onto a scheme of receptor-mediated production of phosphatidylinositol (3,4,5) tris-phosphate (PI(3,4,5)P₃) with feedback through small GTPases (henceforth referred to as the LI model). Narang et al. (22) have analyzed a model abstracted from a scheme of receptor-mediated regulation of phosphatidylinositol (4,5) bis-phosphate (PI(4,5)P₂) levels, modulated by phospholipase C (PLC) activity and feedback through substrate delivery from other membrane compartments (henceforth, NSL model). Postma and Van Haastert (23) have analyzed a general model in which a cytosolic effector molecule enhances receptor-mediated production of a lipid second messenger that, in turn, recruits the effector molecule from the cytosol to the membrane (henceforth, PvH model). These models share important features with our model; other recent models take substantially different approaches (see, for example, Rappel et al. (19) and Haugh and Schneider (20)).

In addition to being based on different biochemical schemes, the above models demonstrate qualitative differences in behavior, suggesting that they represent different modes of gradient sensing. For example, in the LI model the steady-state response of the cell always reflects the current stimulus, whereas in the NSL model, once elicited, a cellular response can persist independently of the external stimulus under some conditions. The PvH model requires a high baseline concentration of translocating molecule on the membrane for efficient gradient sensing. Qualitative comparisons of these models, addressing some of these differences, are published in several recent reviews (7,24).

Our model
If the gradient sensing machinery of the cell is modeled as a reaction-diffusion system, we expect to find qualitative differences in systems that include different spatial couplings and/or exhibit different types of bifurcations. Indeed, a general picture of PI signaling suggests that the existence of coupled positive feedbacks and/or cooperative interactions can lead to bifurcations; molecules that translocate to the membrane from a shared pool in the cytosol can affect spatial couplings. To our knowledge, a systematic and quantitative analysis of how these elements lead to qualitative differences in gradient sensing behavior (such as those noted above), has not been done.

To this end we develop a mathematical model of PI-mediated gradient sensing at an intermediate level of detail. Our model consists of a set of reaction-diffusion equations for the spatiotemporal patterns of 3'PIs on the plasma membrane, as well as for the kinase that generates them and the phosphatase that deactivates them. It reproduces much of the translocation dynamics observed experimentally in Dictyostelium.

By appropriately adjusting parameters, we generate four variants of our model. Our model variants differ in whether coupled positive feedbacks lead to multiple steady states, and whether redistribution of translocating molecules plays an essential role in amplifying responses to gradients; they illustrate the qualitatively different modes of gradient sensing that result in each case.

Distinguishing modes of gradient sensing
Each of our model variants demonstrates three characteristic gradient sensing behaviors, which are enumerated in the next section. We find that differences between the variants become evident in simulated dose-response experiments that highlight transitions between efficient and inefficient gradient sensing. These results are used to define criteria that distinguish between the "modes" of gradient sensing illustrated by our model variants. Applying these criteria to analyze the parameter space of our model suggests that boundaries between different types of behavior can be sharp, and regions that display a given behavior can be narrow with respect to variation of some combinations of parameters. Thus, efficient gradient sensing might require homeostatic mechanisms that regulate combinations of parameters to be within specified ranges. Further, because cells in a given population will have a distribution of biochemical parameters, we expect that subpopulations might function in different regions of parameter space, making use of different modes of gradient sensing. Finally, because biochemical parameters can vary during the course of development (e.g., through changes in gene expression), any given cell might transition between efficient and inefficient gradient sensing, and between different modes of gradient sensing, to suit its needs.

	MODEL DEVELOPMENT

TOP ABSTRACT INTRODUCTION MODEL DEVELOPMENT NUMERICAL TECHNIQUES MODEL CHARACTERIZATION RESULTS SUMMARY AND DISCUSSION SUPPLEMENTARY MATERIAL ACKNOWLEDGEMENTS REFERENCES

 
Characteristic behaviors
A mathematical model of eukaryotic gradient sensing should reproduce three characteristic behaviors, which can be seen in the dynamics of markers for 3'PIs on the plasma membrane:

1. Cells adapt to the average concentration of ligand in solution. This allows sensitivity to relative gradients over many orders of magnitude in absolute concentration (9,25). In the particular case where a uniform dose of ligand is applied suddenly, the cellular response is transient and returns to baseline (26).
   
2. When exposed to a shallow and static gradient of ligand, the cell responds with a sharp and persistent internal gradient of signaling molecules (27). This permits immune cells to migrate toward a source of infection and amoeboid cells to move toward a food source over long distances.
   
3. If the gradient of ligand changes direction, the distribution of signaling molecules will follow with some fidelity (18). A neutrophil thereby can capture a moving bacterium, and Dictyostelium cells can form dynamic patterns during development in response to changing gradients of cAMP.
   

Conceptual framework: local activation/global inhibition
The above behaviors suggest that gradient sensing is mediated by a balance of molecules that act locally, and those that act globally. The former generate a strong response to the greater stimulus at the front of a cell in a gradient (a local quantity), whereas the latter mediate adaptation to the average stimulus (a global quantity). This observation is often summarized as "local activation, global inhibition" (31,32). The equations written to describe such systems (including our own) generally share important features with the activator/inhibitor models developed by Meinhardt and Geirer to explain patterning in development (29), which were extended by Meinhardt to explain aspects of chemotactic response (30).

Geometry and space/timescales
We treat the cell as a disk with the cytosol as its interior and the plasma membrane as its perimeter, reflecting the geometry of a rounded cell where actin has been depolymerized. X marks the position along the membrane and is normalized so that the circumference of the cell is 1. Cytosolic molecules translocate to this boundary, along which 3'PIs diffuse (Fig. 1). Although more pattern forming possibilities would be available if, for example, we treated the membrane as two-dimensional, this simplified geometry adequately accounts for the gradient sensing possibilities that we investigate.

View larger version (32K): [in this window] [in a new window]   FIGURE 1  Schematic of the model. The plasma membrane is modeled as the perimeter of a circle, parameterized by the normalized coordinate, X, along which 3'PIs diffuse and to which cytosolic molecules translocate. Translocation dynamics depicted at the top of the circle occur at the leading edge of a cell in a gradient of stimulus, or transiently in response to a uniform step stimulus. Dynamics depicted at the bottom occur at the back of a cell in a gradient, or during adaptation to a uniform stimulus. Dynamics of molecules labeled in bold type are calculated explicitly in our model; those in normal font are set to steady-state values with respect to those modeled, and those in gray are held fixed. Abbreviations are listed in the "Abbreviations used" footnote at the beginning of the article. The schematic does not include all processes modeled.

Coupling to outside stimuli
In both Dictyostelium and neutrophils, ligand binding activates receptors, which activate the heterotrimeric G-proteins (HTGs) to which they are coupled, in a pattern that closely reflects that of ligand in solution (18,34,36). Both receptor and HTG activation drive recruitment and activation of type I PI3' kinases (PI3Ks, which are enzymes that phosphorylate PIs in the 3' position), though many details are unknown (37,38). For these reasons and because, at least in Dictyostelium, receptor and HTG desensitization do not seem to drive adaptation on timescales considered (39), we let a single variable, , represent ligand-mediated receptor and HTG activation together, which drive PI3K recruitment and activation. In our model, this defines the external stimulus at each point on the membrane.

Biochemical observations, model variables, and network topology
3'PIs
The 3'PIs thought to be relevant in gradient sensing are PI(3,4,5)P₃ and PI(3,4)P₂, both of which comprise on the order of 0.02% of total plasma membrane lipid in resting cells and specifically act to recruit a similar set of cytosolic molecules to the membrane (40–42). Production of PI(3,4,5)P₃ by PI3K acting on PI(4,5)P₂, and of PI(3,4)P₂ by dephosphorylation of PI(3,4,5)P₃, are thought to be the relevant production pathways (40,43). In addition, there is evidence in neutrophil-like cell lines for a positive feedback from 3'PIs to delivery of PI(4,5)P₂ to PI3K, involving small GTPases of the Arf and Rho family (44–46). This is highlighted by Loop I in Fig. 2.

View larger version (21K): [in this window] [in a new window]   FIGURE 2  Biochemical scheme. Numbered feedback loops highlight a suggested network topology. Solid arrows denote chemical conversions, whereas dashed arrows denote stimulation/activation. Loops I and II, representing positive feedbacks, constitute an amplification module; represents the outside stimulus that drives this module by stimulating PI3K translocation to the membrane. Negative feedback Loop III, representing PI3K phosphorylation and inactivation, accounts for response adaptation. Subscripts m and c denote membrane-bound and cytosolic, respectively, and * denotes phosphorylated species, which are catalytically inactive in our model. Abbreviations and fonts are used as in Fig. 1.

Membrane bound PI3K and PTEN
The enzyme that generates 3'PIs, a PI3K, and an enzyme that dephosphorylates them, the PTEN phosphatase, are also essential components of our model. As mentioned, PI3K localization and activation are thought to be coupled to outside stimuli. Further, PI3K localization in Dictyostelium seems to parallel 3'PI localization upon cellular stimulation (37). If we consider 3'PIs to be the primary signal that localizes other molecules in gradient sensing, this suggests a positive feedback from 3'PIs to the enzymes that produce them, represented by Loop II in Fig. 2. PTEN translocation to the membrane occurs in a pattern inverse to that of PI3K in Dictyostelium (37,49). This amplifies the effects of Loop II. In our model, the variables K_m and T_m represent scaled concentrations on the membrane of PI3K and PTEN, respectively.

Cytosolic/inactive PI3K and PTEN
Because cellular response eventually adapts to the average stimulus, as do PI3K and PTEN activities (37,49), some form of integral feedback regulation must exist (50). We represent this by negative feedback Loop III in Fig. 2. Both PI3K and PTEN activity are known to be controlled by phosphorylation in some cell types (51,52). We use the scaled variables K_c* and T_c* to represent the fractional concentrations of PI3K and PTEN, respectively, which are cytoplasmic and phosphorylated. These are catalytically inactive in our model.

Proposed biochemical mechanisms
Loop I: positive feedback through substrate delivery
A possible mechanism for the feedback in Loop I, many elements of which have been studied in neutrophils or neutrophil-like cells, is depicted in Fig. 2. 3'PIs recruit GTP exchange factors (GEFs) to the membrane, where they catalyze the exchange of GDP for GTP in specific small GTPases (g-P) of the Arf and Rho family (53,54). These GTPases are then activated and stabilized on the membrane (55,56), and play roles (together with their regulators) in remodeling the membrane and actin network (57–59). Some have been shown to stimulate PI(4)P5' kinases (PIPKs) to make additional PI(4,5)P₂ (45,46).

Experimental observations do not indicate an accumulation of free PI(4,5)P₂ upon cellular stimulation (60,61), suggesting that the PI(4,5)P₂ generated by the feedback in Loop I is used immediately. This observation could be explained if the generated PI(4,5)P₂ was bound to a transfer protein (PITP) and passed directly to PI3K for conversion to PI(3,4,5)P₃ (there is evidence in neutrophils for PITP involvement in PI(3,4,5)P₃ production (62)). In our model, we assume this mechanism and do not include the dynamics of free PI(4,5)P₂. This simplification is consistent with the lack of clear evidence for spatial gradients of free PI(4,5)P₂ in gradient sensing cells (60,63).

Loop II: positive feedback through regulation of enzymatic activity
For the feedback in Loop II, we propose that 3'PIs recruit an as-yet unidentified molecule to the membrane, which stabilizes membrane-bound PI3K. Membrane-bound PI3K then produces more PI(3,4,5)P₃. To account for a PTEN dynamic inverse to that of PI3K, we propose that PI3K in its capacity as a protein kinase (64,65), or another molecule whose dynamics parallels PI3K dynamics, phosphorylates PTEN (to our knowledge an interaction between PI3K and PTEN has not yet been directly investigated experimentally). In our model, phosphorylated PTEN is cytosolic and inactive.

Loop III: negative feedback for adaptation
To account for response adaptation, we propose that PI3K on the membrane is phosphorylated by an as-yet unidentified kinase, which is constitutively active on the membrane. Phosphorylated PI3K is cytosolic and inactive in our model. PI3K phosphorylation acts as a mechanism of global inhibition because it depletes the cytosolic pool of active/unphosphorylated PI3K, which is a shared pool for recruitment to the entire membrane.

Assessing the biochemical scheme
Several difficulties, often encountered when modeling cellular signal transduction, need to be confronted in developing our biochemical scheme. First, the regulatory mechanisms that are necessary to develop our model have not all been studied in, and may not all exist in, a single cell type. In particular, most of the biochemical interactions suggested for Loop I, as well as regulation of PI3K and PTEN by phosphorylation, have been studied in neutrophils, but have not been sufficiently addressed in Dictyostelium. On the other hand, the translocation dynamics that our model reproduces have been studied primarily in Dictyostelium. Further, several of the relevant biochemical mechanisms, in particular those involved in Loops II and III, are in fact unknown, and the mechanisms that we have proposed are uncertain. In light of these difficulties, we will be primarily concerned in our analysis with the overall roles of coupled positive feedbacks and translocation, and with the differences in gradient sensing mechanism to which they give rise. These types of results should depend on qualitative features of our model, such as network topology, rather than biochemical details. In addition, our model can be used to investigate the consistency of the proposed mechanisms with current data, and to suggest experiments that probe their validity. An example of such an application is outlined below.

Consider the proposed role of PI3K in our model. Recent observations of Dictyostelium cells in which 3'PI dynamics have been perturbed suggest that the qualitative features of PI3K and PTEN translocation, in response to a saturating uniform stimulus, remain unchanged under these conditions (37,48,66). We have simulated the response of our model under conditions of PI3K inhibition and still found a substantial PI3K translocation to the membrane in response to large uniform stimuli, in accordance with these experimental observations (though the magnitude is slightly less and the time course slightly faster; results not shown). More significant differences in translocation dynamics, between conditions where PI3K is inhibited and where it is not, are found in responses to gradients and/or weaker stimuli. Responses to these types of stimuli should be investigated experimentally to further probe a possible involvement of PI3K activity in regulating translocation dynamics.

To reproduce the PTEN translocation observed in Dictyostelium cells where PI3K is inhibited, however, our model must be extended to include PI3K-independent regulation of PTEN (or we must assume that another kinase besides PI3K phosphorylates and inactivates PTEN). Binding to PI(4,5)P₂, which might be depleted by phospholipase C (PLC) activity upon cellular stimulation (47,67), is a potential alternative mechanism for regulating PTEN (66,68). Such an extension of our model is the subject of current work.

From the above discussion, it is clear that forms of regulation not included in our biochemical scheme must be relevant for gradient sensing. Thus, our model should be viewed in a more general sense as a module for PI3K-mediated gradient sensing. Under normal conditions, it reproduces the enumerated characteristic responses and the discussed translocation dynamics. However, it must be integrated with models of other biochemical processes to account for cellular response under a wider range of conditions.

Steady-state assumptions and intermediate level of detail
In our model equations, which are presented in the next subsection, we only explicitly consider the variables summarized in Table 1. The spatiotemporal dynamics of the other regulatory molecules in our proposed biochemical scheme are generally less well characterized. Further, to our knowledge, there are no observed delays in cellular response specifically associated with their activation. For these reasons, we introduce the following simplifying procedure concerning these regulatory molecules: kinetic equations are written for their dynamics; time derivatives are set to zero; the resulting steady-state equations are used to express their concentrations in terms of the variables of our model (see Supplementary Material for more details). This procedure has the virtues of preserving steady-state solutions as well as many effects of translocation. Many uncertain biochemical details no longer appear explicitly in our equations, which are often consistent with alternatives that give rise to similar qualitative behaviors. On the other hand, our model parameters are related to real cellular biochemical parameters. Thus, we can investigate how qualitative behaviors depend on these quantities, and whether our biochemical mechanisms are consistent with current data (e.g., our discussion of PI3K inhibition, above). In this way, our model is developed at an intermediate level of detail.

Generally, the parameter is used to represent scaled forward rate constants, to represent backward rate constants, and to represent saturation concentrations and/or concentrations at which a term becomes effective. The parameter is used to represent constitutive processes acting in parallel with the regulated processes of our model. In a lowest order approximation, these constitutive terms account for some of the many processes not included in our model, and for the observation that inhibition of any single process in our model does not generally lead to complete inactivation of the signaling pathway. The biochemical quantities represented by our model parameters are summarized in Table 2. Note that many of our model parameters are combinations of concentrations and rate constants (see Supplementary Material).

(1)

(1a)

(1b)

(2)

Because the feedbacks in our model generally depend on the sum of the 3'PIs (P_n), which is a primary output of our model, we have made the substitution Eq. 1b wherever possible.

In Eq. 1, the term accounts for P₃ production due to PI3K acting on PITP-bound PI(4,5)P₂. This term couples Loops I and II, and saturates at large K_m because the PITP is depleted. The factor is proportional to the concentration of PITP-bound PI(4,5)P₂. Its form (Eq. 1a) is obtained in the Supplementary Material by writing kinetic equations for the molecules in Loop I (see Fig. 2), many of whose spatiotemporal dynamics are not well characterized, and setting their concentrations to steady state with respect to the variables of our model. The denominator of includes a local term (), which accounts for saturation of this feedback due to depletion of membrane bound molecules, and a global term (), which accounts for saturation of this feedback due to depletion of cytosolic molecules that translocate to the membrane (in particular, the GEF in Loop I). If the latter term dominates, then under conditions where redistribution of translocating molecule keeps fixed, will vary approximately linearly with P_n, and hence contain a term linear in P₃. The degradation term in Eq. 1 is linear in P₃, as well. Then, under these conditions, we might expect a sharp change in P₃ production as K_m and T_m vary, and the balance between production and degradation shifts. The terms, in Eq. 1 account for constitutive production of P₃, independent of the various molecules in Loop I (see Table 2). The loss terms in Eq. 1 account for dephosphorylation of P₃ at the 3' position by PTEN (T_m) and for conversion to P₂ by a phosphatase, such as SHIP (43), whose dynamics are not included in our model ().

Equation 2 describes P₂ dynamics that follows P₃ dynamics with a slight lag. P₂ is generated from P₃ (1st term), as well as from other sources (₂). The loss terms account for dephosphorylation of P₂ at the 3' position by PTEN (T_m), and for dephosphorylation by other phosphatases (represented by ).

In Eqs. 1 and 2, the relative values of _m and _c determine the importance of molecules in Loop I that translocate from the cytosol to the membrane; the ratio ₃/₃ is important in determining the strength of Loop I; the ratios D/₃ and D/₂ control the effects of diffusion on spatial responses. D is the coefficient of lipid diffusion, in units where the circumference of the cell is 1, assumed to be the same for P₃ and P₂.

Membrane bound PI3K: coupling to outside stimuli
To describe the dynamics of membrane bound PI3K (K_m), we write the following equation:

(3)

(3a)

(3b)

The term in Eq. 3 (defined in Eq. 3a) represents recruitment of cytosolic PI3K (K_c) to the membrane in response to receptor activation by outside stimuli (). is scaled by unregulated recruitment of PI3K to the membrane; the factor + 1 thus accounts for the sum of receptor-mediated and unregulated recruitment. Equation 3b expresses conservation of total PI3K and is used to eliminate K_c from our equations (total PI3K is scaled to 1); A_cell is the area of our two-dimensional model cell, equal to 1/4 in units where the circumference of the cell is 1. The loss term in Eq. 3 represents PI3K phosphorylation and removal from the membrane by a constitutively active kinase on the membrane, whose concentration is assumed to be fixed and has been absorbed into the scaled parameter, _K. The factor in the denominator accounts for the feedback in Loop II by decreasing the rate of removal of PI3K from the membrane with increasing P_n. Such a factor can be derived by assuming that 3'PIs recruit molecules to the membrane that bind stochiometrically to PI3K, preventing PI3K from returning to the cytosol (by stabilizing PI3K on the membrane and inhibiting the kinase that phosphorylates it). The concentration of these molecules, and their complex with PI3K, are set to steady state with respect to our model variables. The fraction of PI3K on the membrane that is free to return to the cytosol is then proportional to (derived in the Supplementary Material). This factor becomes important only when Thus, the magnitude of determines the effectiveness of Loop II in signal amplification.

PTEN dynamics: amplifying the effects of PI3K
The following equations describe the dynamics of PTEN on the membrane (T_m), and of the fraction of total PTEN concentration that is cytosolic and phosphorylated (T_c*):

(4)

(4a)

(5)

The first term in Eq. 4 accounts for constitutive recruitment of cytosolic PTEN (T_c) to the membrane with rate constant, _T. Equation 4a expresses conservation of total PTEN and is used to eliminate T_c from our equations (total PTEN is scaled to 1). The loss term represents PTEN removal from the membrane by PI3K-mediated phosphorylation (K_m), as well as constitutive removal (). In Eq. 5, the first term accounts for phosphorylation of PTEN over the entire membrane and subsequent return to the cytosol. The loss term represents dephosphorylation of PTEN in the cytosol, which is constitutive with rate constant _T^*. These equations reproduce a T_m dynamics inverse to K_m dynamics, enhancing the effects of Loop II. Other regulatory mechanisms would have a similar effect if their kinetics paralleled that of PI3K translocation.

Cytosolic/inactive PI3K: adaptation
We write the following equation for the fractional concentration of PI3K that is cytosolic and phosphorylated/inactive (K_c*):

(6)

The first term in Eq. 6 represents phosphorylation of PI3K over the entire membrane. We have chosen parameters such that the reaction that dephosphorylates PI3K in the cytosol (2nd term) is saturated (K_c* >> _K^*). This ensures that, at steady state, we have the following relation (which results from averaging Eq. 3 over the entire membrane, combining with Eq. 6, and replacing the integral with an equivalent spatial average):

(7)

Thus, adaptation to the average stimulus occurs in our model because PI3K is phosphorylated over the entire membrane, depleting the shared cytosolic pool of unphosphorylated PI3K (K_c), such that always return to ₀ (discussed further below).

Modular structure and driving parameter
Equations 1–5 describe positive feedback (Loops I and II), and can be considered to constitute an amplification module; Eq. 6 represents negative feedback (Loop III), and constitutes an adaptation module. The parameter (defined Eq. 3a) is interpreted as a driving parameter. It serves both to couple these modules to each other and to couple the entire system to outside stimuli. The LI model is based on a similar modular structure, as is a more recent variation (69).

This signaling network might function in gradient sensing as follows: receptor activation increases the local value of , and hence the local value of (Eqs. 3 and 3a), recruiting PI3K to the membrane and driving the amplification module (Eqs. 1–5). PI3K on the membrane is then phosphorylated and returns to the cytosol (Eq. 6). The pool of PI3K that is free to return to the membrane (K_c), and hence the value of , drops. If is uniform, returns to ₀ everywhere (see Eq. 7) and the cellular response returns to baseline. Thus, we interpret ₀ as setting the baseline state of the cell. If there is a gradient in , will remain elevated at the front of the cell, where is above its average value, If ₀ is set appropriately, the amplification module will still be driven in this region, but not at the back of the cell where is below and has dropped below ₀.

Model variants and parameters
A highly nonlinear response of our amplification module, and the possibility of multiple steady states, is expected if both Loops I and II become strongly activated upon cellular stimulation. If depletion of translocating molecules saturates either of these loops we expect qualitative differences in responses to uniform stimuli and to gradients, where these molecules can redistribute between the front and back of the cell. To investigate these possibilities, we generate four variants of our model by adjusting the parameter _K (Eq. 3), which controls the effectiveness of Loop II in signal amplification, and the parameter _c (Eqs. 1 and 1a), which controls the degree to which cytosolic depletion saturates Loop I. The resulting differences in the amplification modules of our variants are schematized in Fig. 3. Other model parameters have been set empirically to reproduce characteristic responses seen in Dictyostelium, as most have not been directly measured.

View larger version (24K): [in this window] [in a new window]   FIGURE 3  Model variants. Depicted elements of the model's amplification module are adjusted to define our model variants; thickness of arrow indicates strength of feature. When Loops I and II are sufficiently activated upon cellular stimulation, coupled positive feedbacks can result in multiple steady states (as in Cases 3 and 4). If translocation (Trans.) is important for amplification, then redistribution of molecules between the front and back of the cell can enhance responses to gradients (as in Cases 2 and 4).

Case 2
Multiple steady states are absent. Upon global stimulation, the amplification module depletes cytosolic molecules, and the response saturates; in response to gradients, however, these molecules are redistributed from the back to the front of the cell, stabilizing and enhancing the polarized response. This variant's amplification module shares features with the PvH model.

Case 3
Coupled positive feedback in Loops I and II is sufficient to produce multiple steady states. The amplification module efficiently amplifies responses to uniform stimuli as well as to gradients, but some redistribution of translocating molecules between the front and back of the cell is necessary to stabilize polarized responses against diffusion. In a shallow gradient of stimulus, the cell can be either in a slightly or highly polarized state; switching between these states requires overcoming a threshold in stimulus. Thus, the steady-state response of the cell depends on the history of the applied stimulus, as well as on its current value. We expect this variant to share features with the NSL model (as well as Meinhardt's formulation, (30)), where strong polarization requires overcoming a stimulus threshold and, under some conditions, cellular responses can persist after the stimulus is removed.

Case 4
As in Case 2, responses to uniform stimuli are weak due to depletion of cytosolic molecules. However, redistribution of translocating molecules in response to a gradient enhances amplification due to coupled positive feedbacks, and results in multiple steady states. The uniform state of the cell is unstable. A slight gradient will induce a highly polarized state, which is stable and can persist when the stimulus is removed. This variant shares features with the NSL and Meinhardt models, as well.

	NUMERICAL TECHNIQUES

TOP ABSTRACT INTRODUCTION MODEL DEVELOPMENT NUMERICAL TECHNIQUES MODEL CHARACTERIZATION RESULTS SUMMARY AND DISCUSSION SUPPLEMENTARY MATERIAL ACKNOWLEDGEMENTS REFERENCES

 
We converted our system of PDEs to a system of ODEs by discretizing our spatial coordinate with equally spaced points and using a finite difference approximation for diffusion terms. We solved the system using an Euler method, treating diffusive terms in our equations implicitly and reaction parts explicitly (70). The accuracy of our results was checked by varying the time step (h) and number of spatial points (N) by a factor of 2. Deviations between simulations with different h and N were always <0.2% (generally significantly less) for all the patterns of stimulus that were tested.

To characterize polarized distributions of signaling molecules such as P_n, we define the relative polarization as

(8)

where and We define the direction of polarization of the distribution as

(9)

n (n = integer) is added or subtracted appropriately so that changes continuously during the course of a simulation. We then define the angular velocity of rotation of the distribution by

(10)

	MODEL CHARACTERIZATION

TOP ABSTRACT INTRODUCTION MODEL DEVELOPMENT NUMERICAL TECHNIQUES MODEL CHARACTERIZATION RESULTS SUMMARY AND DISCUSSION SUPPLEMENTARY MATERIAL ACKNOWLEDGEMENTS REFERENCES

 
Amplification modules
The different features of each variant's amplification module (Eqs. 1–5) are depicted in Fig. 3, where the strength or importance of an element is indicated by the thickness of the corresponding arrow. To analyze these differences, we calculate steady-state solutions of Eqs. 1–5 at fixed values of the driving parameter, (see Eq. 3a), for each variant. This is equivalent to analyzing the steady-state response of our model without adaptation; the results obtained are suggestive of the differences in response that will be found later when the full system of equations is considered. To further simplify our analysis, we neglect diffusion and make the assumptions that follow concerning the integrals in our equations. As mentioned, these integrals are equivalent to spatial averages.

If the stimulus is uniform then will vary uniformly, as will the response of the amplification module; then averages, corresponding to the integrals in our equations, will be equal to their values at any point on the membrane. These solutions are represented by the solid curves in Fig. 4 a. Notice that for Cases 2 and 4, due to depletion of cytosolic molecules at higher , the steady-state response of the amplification module is less sharp than it is for Cases 1 and 3. The curve for Case 3 shows multiple steady states for a small range of , due to the simultaneous activation of Loops I and II. These solid curves suggest differences in the responses of each variant to uniform stimuli, where increases transiently and then returns to baseline.

View larger version (24K): [in this window] [in a new window]   FIGURE 4  (a) Amplification modules. Equations 1–5, representing our model's amplification module, are solved for steady state at fixed values of the driving parameter, , neglecting diffusion. Solid curves represent uniform responses. Dashed curves represent responses at a single point on the membrane, where other concentrations are held fixed at uniform steady-state values determined by the circled solutions on the solid curves. is the normalized driving parameter, where i is the value of 0 that sets the baseline for Case 1 (see Eq. 7). (b and c) Restrictions on parameters. Steady-state solutions of Eqs. 1–6, calculated with the stimulus given by: where S = 2 and G = 0.05. is the relative polarization of the 3'PI distribution (defined Eq. 8), normalized by the value calculated when all parameters are set at the fixed values used in subsequent simulations (given in Supplementary Material). To investigate restrictions on baseline, is varied (b). To investigate restrictions on diffusion length, is varied; is the fixed value of used in subsequent simulations (which is the same for all cases) (c).

Constraints on baseline
From the curves in Fig. 4 a it is clear that, in order for each variant to reproduce the characteristic behaviors of gradient sensing, there must be different constraints on the baseline state of the cell, as determined by ₀ (see Eq. 7). If the cell is to adapt to all uniform stimuli, ₀ should not be set in the bistable region of the solid curves in Fig. 4 a. Additionally, if there is to be a sharp cellular response to small gradients in , ₀ should be set to a value where response to nonuniform , which is approximated without diffusion by the dashed curves in Fig. 4 a, is sharp.

To investigate constraints on ₀, we simulate cellular response to a static, spatially linear gradient in , using the full set of equations (1–6); ₀ is varied by appropriately adjusting _K^* and _K (see Eq. 7). In Fig. 4 b, we plot the normalized polarization that results (see caption), as a function of the normalized baseline parameter, where _i is the fixed value of ₀ for variant i, which is used in subsequent simulations. Note that the values of _i have been chosen such that for each curve the polarized response at is nearly optimal.

For Cases 1 and 3, the peak response in Fig. 4 b is constrained to a narrow range of . This indicates that for strong responses to relatively weak gradients, ₀ must be set very near an amplification threshold for the uniform response of the amplification module (compare the positions of the peaks here to the solid curves in Fig. 4 a), but robust to perturbations in many other model parameters. The restriction of ₀ to a narrow range of values suggests that in the development of cells that sense gradients by these mechanisms there must exist homeostatic mechanisms that maintain the baseline state near such a threshold. Levchenko and Iglesias, as well, have found responses to gradients to be sensitive to variations in parameters that set the baseline state of the cell.

Cases 2 and 4, on the other hand, rely on redistribution of signaling molecules between the front and back of the cell to amplify responses to gradients. These cases merely require that ₀ be high enough for a large fraction of these molecules to be on the membrane in the unstimulated cell (this fraction increases with which appears in Eqs. 1 and 1a; see also Supplementary Material). Note that for Case 3, the depicted jump in response represents a discontinuity, corresponding to a bifurcation, though only the more polarized steady state is depicted. Postma and Van Haastert have also demonstrated enhanced responses to gradients with a high baseline concentration of translocating molecule on the membrane.

Constraints on diffusion length
Under conditions where Fig. 4 a suggests a possible bistable response to some patterns of , diffusion can destabilize nonuniform steady-state solutions, and cellular response might depend sharply on the relative value of the coefficient of lipid diffusion (D) in our equations. More specifically, steady-state profiles depend on the ratio where is a degradation rate constant and ² can be thought of as the squared length that a lipid may diffuse before being degraded, per unit concentration of degrading enzyme. Postma and Van Haastert have also analyzed the importance of diffusion length in their model, though it does not include multiple steady states.

In Fig. 4 c we vary D, thereby varying is the value of ² used in subsequent simulations, which was the same for all variants (corresponding to on the order of 5% of the cell's circumference). Again, we plot the normalized polarization of the 3'PI distribution in response to a static gradient, applied such that the more polarized steady state results in cases where multiple steady states are possible. Notice that is in the regime where, for each variant, diffusion has a similarly significant effect on the polarized response of the cell. For Case 3, however, the highly polarized state becomes unstable and is abruptly lost as increases (the discontinuity in the curve represents a bifurcation). Although multiple steady states exist for Case 4 as well, the highly polarized state remains stable, as reflected by the continuous curve.

	RESULTS

TOP ABSTRACT INTRODUCTION MODEL DEVELOPMENT NUMERICAL TECHNIQUES MODEL CHARACTERIZATION RESULTS SUMMARY AND DISCUSSION SUPPLEMENTARY MATERIAL ACKNOWLEDGEMENTS REFERENCES

 
How well does each model variant account for the characteristic behaviors of gradient sensing, and how might the variants be distinguished? To investigate this we simulate cellular response to uniform stimuli, to static gradients, and to rotating gradients. The full system of equations (Eqs. 1–6) is used, together with the fixed parameter values that define each variant (given in Table 2 and Supplementary Material). Differences between the model variants are seen most clearly in dose-response curves for these simulations, which highlight transitions in cellular response. We use these results to define criteria that distinguish between the modes of gradient sensing illustrated by our model variants.

Uniform stimulus
Fig. 5 a illustrates a transient response to a uniform step stimulus applied at t = 0. Time courses were qualitatively similar for all cases (Case 4 is depicted). is increased suddenly, producing an increase in , which then returns to baseline. The concentrations of 3'PIs and PI3K on the membrane (represented by the scaled variables P_n and K_m) increase transiently, whereas the concentration of membrane bound PTEN (T_m) decreases. The peak response, quantified by (P_n)_max, increases as the size of the step in increases, and occurs on a similar timescale to adaptation of . The timescale for adaptation increases as the step size decreases, the effect being more pronounced for Cases 1 and 3 than for Cases 2 and 4 (data not shown).

View larger version (16K): [in this window] [in a new window]   FIGURE 5  Response to a uniform step stimulus applied at (a) Sample time course for Case 4, S = 2 (time courses for other cases were qualitatively similar). All quantities marked by a tilde are normalized by their values before t = 0. (b) The normalized peak 3'PI concentration, (Pn)max, recorded for each variant as a function of the step size, S, in a dose-response curve.

Static gradients
Fig. 6 a illustrates the simulated steady-state profile of a cell in a spatially linear gradient of stimulus (Case 4 is depicted, though all of the cases show qualitatively similar profiles). The steady-state relative polarization of the 3'PI distribution for each variant (see Eq. 8), in response to gradients whose strength is measured by the parameter G (see caption), is given in Fig. 6 b. Gradients are applied together with, or after adaptation to, a uniform stimulus. All of the variants show a strong relative polarization in response to relatively weak gradients.

View larger version (18K): [in this window] [in a new window]   FIGURE 6  Steady-state response to a static, spatially linear gradient. The stimulus is defined by: where is a measure of the relative strength of the gradient, is the coordinate across the cell in the direction of the gradient, and rc is the radius of the cell. Sf = 2; Si is taken as 0 or 2, depending whether the gradient is applied together with, or after equilibration to, a uniform stimulus (this is only labeled for Case 3, as responses of other cases are independent of Si). (a) Sample profile in a 5% relative gradient (G = 0.05). Case 4 is depicted (profiles for all cases were similar). (b) The relative polarization of the 3'PI distribution, {Pn} (defined Eq. 8), recorded for each variant as a function of the relative gradient, G, in a dose-response curve.

Rotating gradients
Investigation of responses to rotating gradients (qualitatively investigated experimentally by Parent and Devreotes (32)) may provide insight into the functioning of each variant in natural settings, where stimuli vary in both space and time. To simulate such responses we first applied a static linear gradient together with a uniform step stimulus, and allowed the polarized distribution of signaling molecules to equilibrate. We then began to rotate the gradient at t = 0 with different periods of rotation, T. The space/time plots in Fig. 7 a record P_n as a grayscale value in sample time courses. The initial gradient is in the direction marked by X = 0 and 1 (the normalized spatial variable, X, is periodic); time is measured relative to T (in revolutions); P_n values are normalized by the peak P_n value before gradient rotation.

View larger version (72K): [in this window] [in a new window]   FIGURE 7  Responses to rotating gradients. Unstimulated cells are first polarized in a static gradient defined by where S = 2 and G = 0.075. After cellular equilibration occurs, the gradient begins to rotate (at t = 0) and the stimulus is given by where T is the period of rotation. (a) Sample time courses of the Pn distribution are plotted as grayscale values, normalized by the peak Pn value before gradient rotation. The spatial variable X, is periodic, and t/T measures time in gradient rotations. T = 150 s (Case 4) illustrates efficient following of a slowly rotating gradient, as indicated by the steady diagonal translation of the grayscale pattern, with a period of 1 (responses for other cases were similar). (b) Quantities that characterize the response are recorded as a function of T. and are long-time averages of and which represent normalized polarizations and angular velocities of the 3'PI distribution, respectively (see text and Eqs. 8–10). L is the long-time value of and measures the lag of the polarization angle behind that of the gradient for Cases 1–3; o measures half the range of during the time after oscillations become steady for Case 4; measures the direction of polarization of the Pn distribution (see Eq. 9).

Sample time courses are also given for shorter T, where the model variants demonstrate different kinds of failure in gradient following. For Cases 1 and 2, the polarized distribution becomes gradually washed out; the weakly polarized steady-state distribution that results eventually follows the direction of gradient rotation. Case 3 suddenly becomes depolarized when gradient rotation becomes too fast to follow. Case 4 remains polarized near its initial direction, turning toward the direction of the gradient whenever it is close to the direction of polarization; an oscillatory steady-state behavior results.

For each time course, we calculated the polarization of the 3'PI distribution (Eq. 8) normalized by its value at and the angular velocity of its direction (Eqs. 9 and 10) normalized by the angular velocity of gradient rotation (). indicates that the polarized distribution remains stable during the gradient rotation; indicates that the direction of polarization follows the gradient perfectly (with a slight lag). For each T, we characterize the cellular response by recording the steady-state quantities, and which are long-time averages of and respectively. If no significant oscillations occur, the long-time value of the lag in the direction of polarization, behind that of the gradient, is recorded (_L, measured in revolutions, for Cases 1–3). When steady oscillations in and do occur, the long-time value of the amplitude of oscillations in polarization direction about the average motion is recorded (denoted _o, measured in revolutions, for Case 4). Resulting dose-response curves to gradients rotating with different periods are given in Fig. 7 b.

For Cases 1 and 2, at shorter T, the P_n distribution becomes increasingly depolarized during an initial transient, after which following becomes perfect (). For Case 3, the highly polarized distribution becomes destabilized if it is not sufficiently aligned with the direction of the gradient, as occurs when the rotation becomes too fast to follow. The sharp drop in at shorter T is a discontinuity, and represents a bifurcation where the highly polarized distribution can no longer rotate stably at the frequency of gradient rotation; if the simulation is initialized in the weakly polarized state (that is, the gradient is applied without a sufficient uniform stimulus; see previous subsection), there is no discontinuity in the dose-response (not shown). Applying a weaker gradient resulted in destabilization of the highly polarized steady state at longer T, whereas a stronger gradient results in a highly polarized state that is stable at shorter T. For a sufficiently strong gradient, Case 3 no longer demonstrates a bistable response to static gradients (see Fig. 6 b), and responses to rotating gradients become