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).

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).

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.

• 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).
• 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.
• 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.

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 (more ...)

A rounded Dictyostelium cell might have a radius (r_c) of 4 μm, and the fastest cellular responses occur on timescales of seconds (17,18). Diffusion coefficients for cytosolic proteins, membrane bound proteins, and membrane lipids, might be of order 20 μm²/s, 0.03 μm²/s, and 0.5 μm²/s, respectively (33–35). On cellular length scales, these estimates result in diffusion times of order 100 ms, 100 s, and 1 s, respectively ( where d is the dimension of the space and D is the relevant diffusion coefficient). Due to these differences in timescales, we simplify our model by treating cytosolic molecules as being uniformly distributed and membrane bound proteins as fixed. Lipid diffusion, however, occurs on the same timescale as cellular response and is calculated. Thus, in our model, cytosolic molecules act globally, coupling reactions at all points on the membrane, proteins act locally, and the spatial characteristics of lipids are context dependent. For a contrasting interpretation, see the model of Rappel et al. (19) where delays due to cytosolic diffusion play an essential role.

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.

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 (more ...)

PI(3,4)P₂ dynamics generally follows PI(3,4,5)P₃ dynamics with a slight lag (47). However, there is evidence in Dictyostelium cells suggesting that disruption of a phosphatase that acts on both (PTEN, discussed below) affects their dynamics differently (48). Thus, we model PI(3,4,5)P₃ and PI(3,4)P₂ separately, and use the scaled variables, P₃ and P₂, to represent their concentrations. Their sum, which we denote P_n, is a primary output of our model.

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.

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.

TABLE 1 Model variables and uniform steady-state values in the unstimulated cell

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).

TABLE 2 Model parameters

3′PI dynamics: a primary model output The following equations describe the dynamics of the 3′PIs, PI(3,4,5)P₃ (P₃), and PI(3,4)P₂ (P₂), on the membrane:

(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).

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 γ₀.

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 (more ...)

The uniform steady-state values of our model variables, which are used to initialize our simulations for each variant, may be found in Table 1; parameter values may be found in Table 2; considerations involved in setting them are discussed in the Supplementary Material. Characteristics of our model variants are summarized in Table 3 and discussed below together with relationships to other models in the literature.

TABLE 3 Characteristics of model variants

Case 1 Coupled positive feedbacks do not result in multiple steady states for any stimulus, and redistribution of signaling molecules does not significantly contribute to response amplification. We expect this variant to share features with the LI model.

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.

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)

(9)

(10)

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.

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 (more ...)

To analyze steady-state responses to nonuniform γ, we fix averages and cytosolic concentrations in our equations at the values determined by the circled points on the solid curves in Fig. 4 a. This corresponds to analyzing the response of the amplification module at a single point on the membrane where concentrations everywhere else are held fixed at uniform steady-state solutions for fixed γ (variations at a single point on the membrane will not affect averages or cytosolic concentration). Solutions are represented by the dashed curves in Fig. 4 a. Their slope, near their intersections with the solid curves, is an important determinant of steady-state cellular response to a gradient in (which generates a gradient in γ). The dashed curves double back on themselves more strongly for Cases 3 and 4 than they do for Cases 1 and 2. This suggests the possibility of multiple steady states in response to nonuniform γ, resulting from the simultaneous activation of coupled positive feedbacks. For Cases 2 and 4, the solid curves differ from the dashed curves more significantly at their intersections—though the uniform response of the amplification module is relatively weak for these Cases, the response at a point can be quite sharp. This difference reflects the greater importance of translocating molecules in amplifying responses to nonuniform γ.

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.

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.

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.

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 (more ...)

In Fig. 5 b, we summarize the results for peak responses to uniform stimuli of magnitude measured by the parameter, S, with a dose-response curve for each variant. Responses are weaker for Cases 2 and 4 than for Cases 1 and 3, due to depletion of cytosolic molecules, consistent with the weaker uniform response of their amplification modules (Fig. 4 a, solid curves). However, whereas all of the curves show quantitative differences, they are qualitatively similar — continuous and monotonically increasing.

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. (more ...)

For Cases 3 and 4, a discontinuity in the dose-response curve indicates the existence of multiple steady states, highlighting the effects of coupled positive feedbacks, and distinguishing these cases qualitatively from Cases 1 and 2 (the discontinuity for Case 4 is at zero gradient). For Case 3, a threshold in stimulus must be overcome to induce a highly polarized response to small gradients. This is achieved if the gradient is applied together with a sufficient uniform stimulus (black curve). Otherwise, the response at small gradients remains weakly polarized (gray curve). For Case 4, the uniform state is unstable to arbitrarily small gradients and no such threshold needs to be crossed; other steady-state solutions are not plotted for Case 4.

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 (more ...)

The first time course, labeled T = 150 s, represents cellular response under conditions of slow gradient rotation. The shape of the P_n distribution remains relatively steady, and its direction follows that of the gradient with a slight lag (the grayscale pattern keeps its shape and is translated diagonally with a period of 1). Time courses are qualitatively similar for all cases under these conditions (Case 4 is depicted).

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 similar to Case 1 (not shown). For Case 4, the polarized distribution remains stable when gradient rotation becomes too fast to follow, and an oscillatory motion results. The transition from gradient-following to oscillatory behavior occurs abruptly, and is indicated by the sharp drop in at shorter T. Similar to Case 3, efficient gradient following continues at shorter T if a stronger gradient is used (not shown).

These simulations indicate that differences between our model variants are most apparent in the characteristics of transitions in response that occur when gradient rotation becomes too fast to follow. The observed differences do not imply that one variant functions more efficiently than the others; rather, each might have a different utility to any given cell types. In particular, these simulations demonstrate that under conditions where a bistable response is possible and the polarized response could potentially “get stuck” in an initial direction (as in Cases 3 and 4), the direction of polarization can still turn to follow a slowly rotating gradient. Although the simulations analyzed thus far are designed to probe the roles of coupled positive feedbacks and translocation in gradient sensing, we note that other types of dynamic simulations and experiments could be useful to further analyze our model, or to investigate other relevant biochemical/biophysical mechanisms (see, for example, Rappel et al. (19) and Subramanian and Narang (71)).

TABLE 4 Criteria for distinguishing modes of gradient sensing

For each combination of parameters, we generated dose-response curves to uniform stimuli and to gradients applied together with, and after equilibration to, a uniform stimulus. The criteria in Table 4 were applied. The results are summarized in Fig. 8 a, where the numbers label regions on this surface in parameter space according to the mode of gradient sensing that best characterizes the simulated responses. An asterisk (*) in a region marks the combination of parameters used to define the corresponding model variant. Regions that gave ambiguous results with respect to translocation were labeled 3′ or 4′ if they met criteria for Modes 3 or 4 with respect to coupled positive feedbacks, and 1′/2′ otherwise (our criteria relating to coupled positive feedbacks alone do not distinguish between Cases 1 and 2). If the relative polarization in response to a gradient defined by G = 0.05 was <0.4 (about half the average value calculated for our model variants; see Fig. 6 b), the region was labeled as “weakly polarized”. We also simulated responses to a rotating gradient with a period of 150 s. Regions were labeled as “poorly following” if either the steady normalized polarization or angular velocity was below 0.6 (see Fig. 7). In these regions of inefficient gradient sensing, we did not apply our criteria.

FIGURE 8 Partitioning of parameter space. Parameters affecting the roles of coupled positive feedbacks and translocation are varied. Criteria in Table 4 are applied to simulation results, as described in the text, and regions of parameter space are labeled by (more ...)

In Fig. 8 a, we see that decreasing κ_K (which increases the effectiveness of Loop II) leads to transitions to modes with successively higher gain in response. For example, decreasing κ_K at 1/κ_c = 0.5 results in the sequence of transitions: ‘Weakly polarized’ → Mode 2 → Mode 4 → Mode 4′ → Mode 3′ → ‘Poorly following’. Increasing 1/κ_c leads to transition to modes where redistribution of translocating molecules makes coupled positive feedbacks more effective in responses to gradients. For example, decreasing 1/κ_c at κ_K = 3 results in the sequence of transitions: ‘Weakly polarized’ → Mode 1′/2′ → Mode 2 → Mode 4. In general, transitions between modes that differ in the roles played by coupled positive feedbacks are sharp (indicated by a bold line), whereas transitions between modes that differ in the importance of translocation pass through regions of parameter space whose qualitative behavior is ambiguous (labeled 1′/2′, 3′, or 4′). Further, the regions where each variant functions efficiently might be narrow in some directions, suggesting constraints on parameters.

The qualitative results of varying other model parameters can be similarly understood, based on the degree to which they adjust the strength of positive feedbacks and the importance of translocating molecules in our model. As an example, we have varied the parameter ζ_T in Eq. 4. This adjusts the importance of constitutive removal of PTEN from the membrane. Larger ζ_T means a less significant inverse translocation of PTEN from the membrane in response to external stimuli, making Loop II less effective. Beginning with the set of parameters that defines Case 1, ζ_T and 1/κ_c were varied to generate Fig. 8 b, in the same way that Fig. 8 a was generated. To maintain a similar baseline value of T_m when ζ_T was varied, λ_T was adjusted such that = constant, where is the value of K_m in the unstimulated cell for Case 1.

Comparing the structures in Fig. 8, a and b, we see that similar transitions result when κ_K is decreased (panel a) or when ζ_T is decreased (panel b). Thus, we might expect that their simultaneous variations can compensate for each other. This is demonstrated in Fig. 8 c, where we see that decreasing κ_K while appropriately increasing ζ_T often result in the same mode of gradient sensing (other parameters were set as in Case 4). Similar results were found for variations of other parameters, such as κ₃, which adjusts the value of K_m for which binding of PI3K to the PITP becomes saturated (results not shown).

We have developed a model of eukaryotic gradient sensing at an intermediate level of detail. The model is based on a scheme of 3′PI signaling thought to be important in chemotaxis of highly motile cells such as Dictyostelium and neutrophils. By appropriately adjusting parameters, we defined and analyzed four variants of our model, whose qualitative features are summarized in Table 3. Our variants differ in whether coupled positive feedbacks lead to multiple steady states, and whether redistribution of molecules that translocate from the cytosol to the membrane leads to significant differences in responses to uniform stimuli and to gradients. We expect these possibilities to be available to any gradient sensing cell.

Each of our model variants accounts for the characteristic behaviors of gradient sensing enumerated in our Introduction. Experimentally, these are generally elicited by strong stimuli, such as a saturating uniform dose of chemoattractant, or a steep gradient (e.g., generated by a pipette that leaks chemoattractant close to the cell). Our simulations suggest that each variant behaves differently under conditions where the gradient sensing response begins to fail. For example: when a uniform stimulus becomes too small to elicit an observable response; when a gradient is too small to produce a highly polarized response; or when the pipette in a rotating gradient experiment begins to move too quickly for the cell to follow. These transitions, from efficient gradient detection to failure, are best observed in dose-response experiments (Figs. 5–7).

Based on differences in simulated dose-response experiments to uniform stimuli and to gradients of different strengths (Figs. 5 and 6), we have defined criteria that determine whether a cell with unknown biochemical parameters makes use of a mode of gradient sensing illustrated by one of our model variants (Table 4). Using these criteria, we have characterized several surfaces in the parameter space of our model by determining which mode of gradient sensing best describes the behavior at each point (Fig. 8). Because our criteria depend on qualitative differences in response, rather than particular details of our model, we expect them to be applicable as well to other gradient sensing cells that display characteristic behaviors, and to analyzing models based on other biochemical schemes.

Our analysis of cellular response to rotating gradients (Fig. 7) suggests that the differences that define our variants have consequences for cellular behaviors in natural settings, where chemotactic stimuli vary in both space and time. For example, Cases 1 and 2 always follow a rotating gradient perfectly, their polarization gradually weakening for fast rotations. The polarized response for Case 3 can be turned on and off by strong stimuli that change quickly. Case 4 remains persistently polarized in an initial direction when the gradient changes quickly. Each type of behavior might have a different utility for the cells that use them. These differences would influence, for example, the patterns formed during Dictyostelium aggregation, and the population dynamics of neutrophils migrating to a source of infection.

Applying the criteria in Table 4 to real cells could guide the identification of important positive feedbacks and translocating molecules that shape essential features of gradient sensing responses; our model could be used to suggest the qualitative results of perturbing these elements. For such applications, cell to cell variability will have important consequences. The values of biochemical parameters in each cell, even in a clonal population, will generally be different (different concentrations of each protein will be expressed, etc.). Responses to a given stimulus will then vary, and transitions in the response of a given cell will occur at different stimuli. Thus, multiple experiments must be done on single cells to apply the criteria in Table 4, and analysis of many cells will be necessary to characterize the distribution of behaviors in a population. The need for quantitative measurements on single cells to observe transitions in behavior has been emphasized in connection to other cellular systems (see, for example, Cluzel et al. (72) and Ferrell and Machleder (73)); such measurements are only recently being done in connection to gradient sensing (for example, Postma et al. and others (15,17,35,74).

In developing our model, we found that many parameter variations had similar effects and could compensate for each other. Further, certain combinations of parameters were more important than others in determining qualitative behaviors. For example, varying any of the several parameters that control the strength of the positive feedbacks in our model might lead to similar bifurcations. This is illustrated in Fig. 8, where the boundaries that separate Cases 1 and 2 (which have unique steady states) from Cases 3 and 4 (which have multiple steady states) are highlighted in bold. Combinations of parameters near these boundaries often resulted in more efficient gradient sensing responses. For example, Modes 1 and 2 showed enhanced responses to small gradients; Modes 3 and 4 followed rotating gradients more efficiently. As with our analysis of constraints in Fig. 4, these observations suggest that some combinations of parameters must be set quite carefully for efficient gradient sensing, and that homeostatic mechanisms exist to maintain their balance in cells.

The partitioning illustrated in Fig. 8 suggests that we think of the different modes of gradient sensing, represented by our model variants, as existing within volumes of parameter space and functioning efficiently within subvolumes. These subvolumes might have sharp boundaries and may be narrow with respect to variations of some combinations of parameters. This presents the intriguing possibility that different cells in a population, which differ in the values of their biochemical parameters, might employ different modes of gradient sensing and function at varying levels of efficiency. Each mode may have a different utility to a cell of a given type. Thus, during the course of development, a cell might transition between efficient gradient detection and nonresponsiveness, or between different modes of gradient detection, by modulating its biochemical parameters to suit its needs.

An online supplement to this article can be found by visiting BJ Online at http://www.biophysj.org.

We thank Erin Rericha for a critical reading of the manuscript and helpful discussion. This work was done in partial fulfillment of the requirements for R. Skupsky's doctorate from the University of Maryland. This study utilized the high-performance computational capabilities of the Biowulf PC/Linux cluster at the National Institutes of Health (http://biowulf.nih.gov).