A bistable mechanism for directional sensing

We present a generic mechanism for directional sensing in eukaryotic cells that is based on bistable dynamics. As the key feature of this modeling approach, the velocity of trigger waves in the bistable sensing system changes its sign across cells that are exposed to an external chemoattractant gradient. This is achieved by combining a two-component activator/inhibitor system with a bistable switch that induces an identical symmetry breaking for arbitrary gradient input signals. A simple kinetic example is designed to illustrate the dynamics of a bistable directional sensing mechanism in numerical simulations.


Introduction
Chemotaxis is the directed movement of a cell towards a chemical source. It is of fundamental importance for many biomedical processes including wound healing, cancer metastasis and morphogenesis of the nervous system [1]- [3]. Directed locomotion of bacteria has been intensively investigated and detailed models of prokaryotic chemotaxis are available [4]. On the other hand, directional motion of eukaryotic cells is more complex. Here, less is known about the chemotactic signaling pathways that link the membrane receptor input to rearrangements of the cytoskeleton and directed actin polymerization [5]. Eukaryotic cells like neutrophils or the social amoeba Dictyostelium discoideum can detect chemoattractant gradients as shallow as a 1% difference in concentration between the front and the back of the cell and exhibit robust directional motion over a large range of different gradient steepnesses [6]- [8].
This response behavior cannot be achieved by simple amplification or any other linear mapping of the external gradient into the cell. Instead, it can be expected that highly nonlinear interactions govern the early stages of chemotaxis, generally referred to as directional sensing. During the first 10 s after gradient exposure, an intracellular symmetry breaking occurs that is reflected in asymmetric spatial distributions of numerous proteins across the cell [9]. These subcellular reorganizations have been experimentally observed by fluorescence microscopy imaging of various green fluorescent protein (GFP)-tagged constructs. Among these redistributing components, the small GTPase Ras, pleckstrin homology (PH) domain proteins, and the p21-activated kinase A (PAKa) have been identified, see [10] and references therein. Although a complete picture of the precise role of these individual players is still missing, it is generally assumed that the asymmetric rearrangement of proteins during directional sensing controls the downstream events that lead to actin driven membrane protrusions and, ultimately, cell locomotion [11].
Numerous efforts have been undertaken to develop mathematical descriptions of directional sensing. Since many molecular details of the chemotactic signaling network remain unknown, most models proposed to date consequently focused on abstract, low-dimensional descriptions. The most common approach is based on two-component activator/inhibitor kinetics. Combining slow diffusing, locally acting excitatory components with a rapidly diffusing, global inhibitor (local excitation/global inhibition (LEGI) models) amplification of the external gradient signal is obtained [9,12,13]. Other variants have been proposed, like the 'first hit' [14] and the 'intermediate depletion' model [15]. While the former is unable to adapt to subsequent changes in the direction of the initial gradient stimulus, the latter suffers from a strong dependence on the average chemoattractant concentration that is not observed in experiments. Intracellular symmetry breaking via a Turing instability in response to external gradient stimuli has been considered by Meinhardt [16]. This approach was later refined and adapted to specifically describe the dynamics of asymmetric localization of phosphoinositides such as PIP 3 at the leading edge of migrating cells [17,18]. As more and more elements of the chemotactic pathway are identified, first attempts were made to capture directional signaling in high-dimensional realistic models [19]. Here, the authors simulate a complex signaling network to describe temporally multiphasic responses in membrane translocation events following gradient stimulation. Current modeling efforts focus on stochastic effects that influence the distribution of occupied receptors across the membrane and become important for directional sensing in shallow gradients [20]- [23]. In [8], we performed a quantitative study of chemotactic motion of the eukaryotic micro-organism D. discoideum. Chemotaxis of Dictyostelium cells was characterized in wellcontrolled linear chemoattractant gradients using microfluidic technology. In this study, we observed chemotactic responses for gradient signals that range over three orders of magnitude in gradient steepness (between 10 −2 and 10 nM µm −1 ). Within this interval, chemotactic motion was robust and showed only a weak dependence of the chemotactic velocity on the gradient steepness. This has prompted us to look for models that respond with an identical symmetry breaking to gradient stimuli of very different steepness.
The models summarized in the previous paragraph do not meet this condition. Gamba et al have designed a model of directional sensing that is based on phase separation of membrane phospholipids into PIP 2 -and PIP 3 -rich domains, mediated by phosphoinositide diffusion and the enzymatic activity of phosphatidylinositol 3-kinase (PI3K) and its counteracting phosphatase PTEN [24,25]. Their diffusion-limited phase separation model accounts for a strong symmetry breaking even at shallow gradients. But it does not show the initial transient uniform response to both isotropic and gradient stimuli that is observed in experiments on phosphoinositide signaling in chemotactic cells [26]. Recently, Levine et al proposed a model that takes the switch-like nature of gradient sensing into account and, at the same time, also reproduces the transient response dynamics correctly [27]. Their 'balanced inactivation model' is composed of a two-component activator/inhibitor system. Under gradient stimulation, this model exhibits a clear symmetry breaking into an activated front and a back, where activation is vanishing. Note, however, that activation levels at the front depend on the external signal and scale linearly with the gradient steepness.
Here, we propose a mechanism where the activated and quiescent states at the front and back of a cell are independent of the external gradient signal as long as a critical gradient value is exceeded. We achieve this property by building our model on the simple pattern formation paradigm of bistability: an underlying LEGI-type activator/inhibitor system is combined with an autocatalytic step that results in an identical symmetry breaking event for arbitrary gradient inputs. Note that our model is restricted to the initial symmetry breaking during directional sensing. Down-stream cytoskeletal responses like the protrusion and bifurcation of pseudopods [28] or the dynamics of complex actin structures in the cell cortex [29] are not described by our model. In the following sections, we first introduce the underlying idea of a bistable model for directional sensing and later present numerical simulations of a specific kinetic example.

Bistable model
We assume that the chemotactic pathway is mediated by a membrane-bound species P. It lives on the cytosolic side of the plasma membrane and is diffusively mobile along the inner membrane surface. The formation of P will depend on extracellular chemoattractant signals that are communicated to the inside of the cell via transmembrane receptors. In the presence of P, downstream events are initiated that influence actin polymerization and rapid reorganization of the cytoskeleton. In this simplified picture, directional sensing is defined as the formation of an asymmetric membrane distribution of P in response to an extracellular gradient signal. The asymmetry in P will result in asymmetric stimulation of cytoskeletal activity that controls the formation of membrane protrusions and directional cell movement.
In general, the dynamics of P on the inner side of the plasma membrane follows a reaction-diffusion equatioṅ where p denotes the concentration of P and D P the coefficient of diffusion along the membrane. Note that throughout this work, physical quantities are given in dimensionless variables. They have been non-dimensionalized by the reference quantities that can be found in appendix A. Quantities that carry a dimension are marked with an asterisk. The kinetics of P formation and degradation depends on the external chemoattractant signal and is summarized in the interaction function f . For the model presented here we require that the form of f leads to bistability in the dynamics of P. This typically implies that the local kineticsṗ = f has three fixed points, p 1 < p 2 < p 3 , where p 1 and p 3 are linearly stable and p 2 is unstable [30]. Here, we assume that the two stable fixed points p 1 and p 3 are independent of the external chemoattractant signal S, whereas the unstable fixed point p 2 depends on S (see appendix B for the more general case, in which both p 2 and p 3 are depending on S). Cellular systems are noisy, so that localized transitions from one fixed point to the other can occur. Due to diffusive coupling, such perturbations may induce similar transitions at neighboring locations, so that transitions from one fixed point to the other can spread through the system in the form of a trigger wave. Their shape and velocity are uniquely defined. Let p = p(ξ ) with ξ = x − ct denote a transition from p 1 to p 3 , that spreads with velocity c through the system, i.e. p → p 1 ( p → p 3 ) for ξ → ∞ (ξ → −∞). The sign of the trigger wave velocity c is determined by the sign of the integral [31], Depending on the sign of c, domains of p 3 (c > 0) or p 1 (c < 0) will grow, so that the system eventually converges to one of the two stable fixed points 4 . As the central assumption of our model, it is required that the trigger wave velocity changes sign across a cell that is exposed to an external gradient in the chemoattractant S (s denoting the concentration of species S). In particular, in the front half of the cell, pointing towards higher chemoattractant concentrations, the wave velocity shall be positive, so that p eventually takes the value of the stable fixed point p 3 . In the back half of the cell, c is negative and p converges to p 1 . The interface between the two membrane fractions of high and low P concentration is formed by a frozen trigger wave of zero velocity.
How can such a dependence of the dynamics of P on the external chemoattractant signal arise? We assume that the receptor input is transmitted to the P-signaling system via a twocomponent activator/inhibitor module that shows similar dynamical behavior as previously proposed LEGI-type models of directional sensing [12]. The system shall have the following properties: (i) both activator and inhibitor are produced locally proportional to the external chemical signal. (ii) The activator A is an immobile, membrane bound species, so that its concentration a on the membrane reflects the external gradient across the cell. (iii) The inhibitor B is a cytosolic, fast diffusing component. Assuming that diffusive spreading of the inhibitor occurs instantaneously, its concentration b is proportional to the midpoint concentration of the external gradient everywhere in the cell. This is illustrated in figure 1(b), where the activator and inhibitor concentrations are shown across the one-dimensional projection of a cell exposed to an external gradient stimulus as indicated schematically in figure 1(a). (iv) The sign of the trigger wave velocity depends on the ratio of the activator and inhibitor concentrations a and b and has to change across a cell that is exposed to an external gradient. In the example described below, we have, c = 0 for a = b, c > 0 for a > b, and c < 0 for a < b. Perturbations in p then induce a convergence to p = p 3 in the front and to p = p 1 in the back of the cell as shown in figure 1(c). (v) Inhibitor levels are slightly higher than the averaged activator concentration in stationary cases. Under this condition, the entire membrane will be eventually driven to p = p 1 for uniform, stationary levels of chemoattractant, since c < 0. This situation is schematically illustrated in figures 1(d)-(f). (vi) Production and degradation of the inhibitor are slower than the corresponding timescales of the activator dynamics, i.e. upon uniform stimulation, the activator reaches stationary levels earlier than the inhibitor, see figures 1(g) and (h). For intermediate times, the activator may thus exceed the concentration of the inhibitor, so that a transient response of p = p 3 can be observed uniformly across the cell as displayed in figure 1(i). For longer times, perturbations will decay to p = p 1 as required for stationary cases.

Kinetic example and numerical simulations
Let us now turn to a simple example that fulfills the general requirements outlined in the previous section. We will first describe the specific model composed of an activator/inhibitor system (LEGI module) that couples to an autocatalytic membrane bound species P (bistable module). We then present a numerical analysis of the one-dimensional model and an example of a two-dimensional simulation.

Model
The external chemoattractant signal S acts on a two-component LEGItype activator/inhibitor system. According to condition (i), we require that both the activator A and the inhibitor B are produced proportional to the external chemoattractant signal S. Furthermore, A and B are part of other metabolic pathways, so that they are continuously formed and degraded independently of the external signal S, The S-independent production and degradation reactions ensure that both species exhibit a residual background concentration in the absence of S. As required by conditions (ii) and (iii), the activator A is an immobile, membrane bound component, while the inhibitor B spreads in the cytosol by diffusion. If cytosolic diffusion of B occurs quickly, we can approximate the production of B to depend on the averaged concentration of S,s = s dx, at each location in the cell. We thus obtain the following kinetics for A and B: Here, we have assumed that the concentrations of species A and B for the endogenous production of A and B are constant (large reservoirs). The numerical values of the rate constants of activator and inhibitor kinetics do not matter as long as conditions (v) and (vi) are fulfilled. The values that were chosen for our numerical simulations can be found in the caption of figure 3. We assume that the action of the activator/inhibitor system on downstream events in the signaling pathway is determined by the ratio of activator and inhibitor concentrations [12]. Such a dependence can occur via a third species Q that is formed and decomposed depending on the presence of the activator and inhibitor, Again, we assume a constant substrate concentration Q (large reservoir), so that the kinetics of Q is If the dynamics of Q is fast compared to the timescales of the LEGI system, the concentration of Q will be determined at all times by the ratio of activator and inhibitor concentrations, For simplicity, we take k 4 = k −4 , so that q = b/a. The concentration of Q is the readout quantity of the LEGI module.

Bistable module.
Following the general outline in section 2, the LEGI module transmits the receptor input signal S to the bistable signaling system of a membrane bound species P.
In general, temporal changes in the concentration of P, d p dt = f ( p), will be determined by production and degradation terms, For our specific example, we assume an autocatalytic production of P that takes place in an enzymatic reaction following Michaelis-Menten-type kinetics. Degradation of P is taken to be linear, so that We take cooperative binding effects into account by introducing a Hill coefficient α (in the present example, we set α = 2). Under the assumption of a quasi-steady-state, the concentration of the enzyme-substrate complex X is constant. With the total enzyme concentration e 0 conserved, the well-known Michaelis-Menten rate law can be derived for the production of P. We thus find the following expression for the interaction function f : Here, v = k 2 e 0 denotes the maximal reaction velocity and K = (k −1 + k 2 )/k 1 the Michaelis constant. The shape of the interaction function f can be seen in figure 2. It has two stable fixed points p 1 and p 3 , and an unstable fixed point p 2 , for which p 1 < p 2 < p 3 . The fixed points are determined by How is the dynamics of the species P linked to the LEGI module and thus to the external signal S? As a key property, our model has to fulfill condition (iv). The latter requires that the trigger wave velocity, i.e. the integral (2), changes sign across a cell that is exposed to an external gradient in S. This can be achieved if we assume that the unstable fixed point p 2 depends on the concentration of Q, i.e. on the ratio of activator and inhibitor concentrations, while the two stable fixed points are independent of Q. To obtain such a dependence, we require K + 1 = v/k 3 , so that equation (12) reduces to In this case, the two stable fixed points become p 1 = 0 and p 3 = 1, and the unstable fixed point is p 2 = K . The dependence of p 2 on the ratio of the activator and inhibitor concentrations is now implemented by choosing K =K q. Since K = (k −1 + k 2 )/k 1 , the q-dependence of K implies a q-dependence of at least one of the rate constants, e.g. k −1 =k −1 q, withk −1 k 2 . Note that from the condition K + 1 = v/k 3 , it also follows that v/k 3 = k 2 e 0 /k 3 has to depend on q. We achieve this by assuming that the total enzyme concentration is e 0 =ẽ 0 (K q + 1) with k 2ẽ0 /k 3 = 1.
The value ofK has to be chosen such that the trigger wave velocity changes sign across a cell that is exposed to a chemoattractant gradient. Under a gradient stimulus, the activator and inhibitor concentrations change from a > b to a < b across the cell, see figure 1(b). We therefore require that the integral (2), and thus the trigger wave velocity, is zero for a = b. This ensures an opposite sign of the trigger wave velocity in the parts of the cell where a > b and a < b, respectively. In figure 2, f ( p) is shown for a > b (black) b = a (dark gray) and a < b (light gray) to illustrate this. Under the condition that A = 0 for a = b, equation (14) leads to the relation √ K arctan 1 √ K = 1 − 1 2(K +1) , which can be solved numerically to yieldK = 0.4357.
In the preceding discussion, several restrictions were imposed on the choice of the model parameters to ensure that p 3 is independent of q and to enforce a zero trigger wave velocity for q = 1. In appendix B, we have analyzed the robustness of our model under conditions where these restrictions, in particular the condition K + 1 = v/k 3 , are relaxed.

Equations and parameters.
We have performed numerical simulations of equation (1) for the dynamics of p, where f is given by equation (11). The activator and inhibitor dynamics are based on equations (5). The parameter values can be found in the caption of figure 3. Simulations are performed on a one-dimensional domain of length L = 10 that is discretized with a grid length of x = 0.02. To mimic a two-dimensional cut through a circular cell, periodic boundary conditions are imposed. A linear gradient stimulus that spans froms − s max tos + s max then maps onto the cell perimeter in the form of a sine function, s =s + s max sin( 2π L x). The time step is t = 0.001 and for the diffusion constant of P we choose D P = 0.1, a typical value for membrane diffusion. The system of equations is integrated numerically using a forward Euler scheme and a nearest-neighbor representation of the Laplacian operator. We assume that the cellular system is noisy, so that, on a microscopic level, the quantity p undergoes random fluctuations. From time to time, fluctuations add up and become sufficiently large to induce a local transition between the two stable fixed points of the system. In our simulations, noise is taken indirectly into account by initiating random localized transitions between the two fixed points. The positions of perturbations are randomly chosen from a uniform distribution across the cell. The perturbations are x = 0.08 in size and their amplitudes are uniformly distributed between 0 and 1. Also the waiting times between such perturbations is random and drawn from a uniform distribution that extends from zero to some maximal waiting time. Depending on the sign of the trigger wave velocity, perturbations may grow or decay.

Response to shallow gradients.
Numerical simulations were performed to mimic stimulation with directional and uniform cues. In figure 3(b), a p-profile is shown that results from exposure to a linear gradient in s across the cell. The gradient stimulus is shown in figure 3(a) and corresponds to a sine profile along the perimeter of a circular cell. Figure 3   One of the outstanding characteristics of directional sensing in eukaryotic cells is an extraordinary sensitivity to shallow gradients. Experiments have shown that concentration differences as low as 1% between the front and the back of a cell can be identified [6]- [8]. The bistable mechanism we propose here can account for a reliable performance of directional sensing over a large range of gradient steepnesses. An example for sensing of shallow gradients is shown in figures 3(a) and (b). A gradient in S that is ranging from s = 9.9 to 10.1 across the cell is sufficient to induce a clear symmetry breaking: in the front half of the cell, p converges to the stable fixed point of p 3 = 1, in the back it takes the value of the second stable fixed point p 1 = 0. Although the limit for detection of shallow gradients will be ultimately determined by the noise level of the system, the constant values of the fixed points p 1,3 ensure that our model shows an almost identical response over a wide range of gradient stimuli. Only two features of the response dynamics depend on the gradient steepness and are described in the following.

Threshold for directional responses.
The fraction of the membrane with p = 0 is slightly larger than the p = 1 part. This asymmetry increases for decreasing gradients. It is caused by the inhibitor level that always exceeds the averaged activator concentration as required by condition (v) in section 2. This is schematically shown in figure 4(a). With decreasing gradient, the membrane fraction for which a > b becomes smaller until it disappears for Here,s denotes the average concentration of S over the cell and s max the maximal deviation from the average, i.e. the gradient stimulus spans froms − s max tos + s max across the cell.

Timescale of directional responses.
The time it takes to establish a directional response increases with decreasing gradient. The reason for this is a dependence of the trigger wave velocity on the ratio of activator and inhibitor concentrations, q = b/a. This dependence is displayed in figure 4(b). The shallower the gradient, the smaller the difference between a and b, so that the trigger wave velocity is confined to small values around zero. This results in a slower growth of perturbations and thus in longer transient times to establish a directional response. For values of q = 0.5, . . . , 2, a linear relation between q and the trigger wave velocity is found. In this regime, which spans up to gradients of about 50%, the model can be operated safely.
Situations outside this range are not considered here 5 .
We have characterized the dependence of the response time on the gradient steepness. To this end, we define the deviation of the actual system state from a complete directional response as where N is the number of grid points along the cell membrane andp i = 0 or 1 for a < b or a > b, respectively 6 . Following exposure to a gradient, δ will decay as p approaches the response statep. In figure 5(a), the decay of δ is shown for a selection of different gradients. We determine the response time by fitting with an exponential of the form C e −t/τ + D with the fit parameters C, D and τ . The latter is taken as the characteristic timescale of the response. In figure 5(b), the timescale τ is shown for gradients ranging from 1 to 10%. At shallow gradients, increased response times are observed. For gradients above 4%, τ stabilizes at values below 3 and is only weakly dependent on the gradient steepness.

Influence of noise.
The presence of noise is an essential prerequisite for our bistable model of directional sensing. Without noise, the p-signaling system cannot initiate transitions between the two stable fixed points p 1 and p 3 in response to changes in the external signal S. As explained above, noise is indirectly taken into account by imposing localized perturbations that are randomly distributed in space and time (see the paragraph on equations and parameters at the beginning of this section). How does the response dynamics of the model depend on these random perturbations?
The essential parameter that influences the response dynamics is the waiting time between the application of individual perturbations. The waiting times are random and drawn from a uniform distribution between zero and some maximal waiting time. Two limiting cases can be distinguished. For infinitely long waiting times, no response to external stimuli will occur since no transitions from the quiescent state p = 0 to the fixed point at p = 1 are initiated. In the limit of waiting times that are much shorter than the characteristic timescale of the trigger wave velocity, the system is dominated by perturbations, so that no coherent response profile can form. We have performed a series of numerical simulations, in which we systematically investigated the influence of the average waiting time noise t between perturbations on the response dynamics of the model. Again, we recorded the temporal behavior of the deviation δ of p from the target response statep; for a definition of δ see equation (16).
In figure 6(a), the time evolution of δ is shown for a number of different average waiting times noise t. For small noise t, the high frequency of perturbations prevent a relaxation to the complete response state and δ takes a finite value different from zero. As noise t increases, δ is converging to values close to zero indicating a clear directional response of the system. However, for large noise t, perturbations become rare and strongly increased decay times of δ are observed. To quantify the characteristic decay times, an exponential decay function C e −t/τ + D is fitted to the δ(t)-profiles. The resulting values of the characteristic response timescale τ are displayed in figure 6(b). It shows a steep increase for short average waiting times and levels off to a linear increase for growing noise t.

Two-dimensional example.
We have performed two-dimensional simulations of our model. Contrary to the previous computations, where the cytosolic concentration of the inhibitor B was taken proportional to the average stimulus across the cells, we now take cytosolic diffusion of B explicitly into account, with D B the corresponding cytosolic diffusion constant. Simulations were performed in polar coordinates on a circular domain of radius r = 4. Angle and radius were discretized with r = 0.2 and θ = 2π/80, respectively, the time step was t = 10 −6 until t = 4 and t = 10 −4 for larger times. In figure 7, an example of a simulation is displayed that corresponds to a two-dimensional cut through a circular cell. The cell is exposed to an extracellular gradient in S that points in the upward direction and changes by 10% across the cell. On the left-hand side, figure 7(a), the activator and inhibitor concentrations are shown on the membrane and in the cytosol, respectively. Note that due to the finite diffusion constant of B, a slight gradient in the cytosolic inhibitor concentration can be observed. In figure 7(b), the distribution of P along the membrane is displayed. It can be seen that two domains of p = 1 and p = 0 are formed that are connected by narrow interfaces in the regions where the activator and inhibitor concentrations are equal, a = b. These interfaces correspond to stationary trigger waves.

Discussion
We have proposed a mechanism for directional sensing based on bistable dynamics of a membrane bound autocatalytic component P. Our model is constructed of two modules: the external signal S acts on a two-variable, LEGI-type activator/inhibitor system, which controls the autocatalytic kinetics of a bistable species P. In section 3, we presented a kinetic example that fulfills the requirements for bistable directional sensing that were explained in section 2.
Let us summarize the predictions of our model. They do not depend on the specific choice of the kinetic equations (5) and (11) (15) is not fulfilled. The reason for this is that the two stable fixed points p 1 and p 3 are independent of the external stimulus S. (4) There is a threshold to trigger a response in p, i.e. uniform stimuli have to exceed a certain minimal increase in s and directional stimuli have to exceed a minimum gradient steepness, see condition (15). This is due to an inhibitor concentration that exceeds the averaged activator levels as illustrated in figure 4(a) for the gradient case. The properties (1) and (2) have been experimentally established and form the basis of previous LEGI-type directional sensing models [26]. To assess the validity of our present model, experiments are required that focus on the predictions (3)-(6). In particular, systematic mapping of the intracellular response as a function of the external gradient signal is needed. The onedimensional simulations can then be compared to two-dimensional confocal cuts through the cell body. Systematic experiments rely on techniques to produce well-controlled gradient signals with high spatiotemporal resolution. In the recent literature, several experimental efforts are reported to perform quantitative experiments on chemotaxis and on the intracellular rearrangement of proteins related to directional sensing [32,33]. In particular, advances in microfluidic technology have opened up new perspectives to construct devices for quantitative migration experiments [8,34,35] and for well-controlled single cell stimulation, see [36] and references therein. In a recent publication, we have presented a technique for controlled directional stimulation of chemotactic cells based on the photochemical release of signaling substances in a microfluidic chamber [37]. We suggest that this approach may provide the basis to experimentally test the modeling predictions made here.
In earlier modeling work, it has been speculated that autocatalytic steps are involved in membrane phospholipid kinetics [17]. In the more recent literature, the role of phospholipid signaling is controversially discussed in the light of current experimental findings [38,39]. It has been suggested that PI3K is not essential for directional movement under strong stimuli but plays an important role for directional sensing in shallow gradients [40]. However, no ultimate consensus has been reached about the role of phosphoinositide signaling in the chemotactic pathway. In the present work, we therefore do not propose any explicitly biochemical assignment to the model components A, B and P. The specific reaction scheme that leads to equation (11) is a simplistic example, intended to illustrate the dynamics of a bistable directional sensing mechanism. It is not supported by direct biochemical evidence. Other kinetic schemes can be found that fulfill the assumptions of section 2 equally well.
We have proposed a generic mechanism of directional sensing that is based on bistable dynamics. In particular, our model combines a LEGI-type activator/inhibitor module with an additional bistable reaction. The activator/inhibitor part of the model does not involve any autocatalytic terms or other nonlinear interactions. Thus it will not generate any selforganized structures on its own. Only in combination with the additional bistable component are the characteristic properties of a pattern forming system introduced. This allows us to combine the advantages of a LEGI model-adaptation after uniform stimuli, independence of the midpoint concentration-with the switching dynamics of a bistable system. To show the desired dynamical behavior, only a number of general conditions have to be fulfilled (listed in section 2). As a result, we obtain a model that shows a clear and almost identical symmetry breaking over a large range of gradient input signals. This is in contrast to earlier models, where a dependence of the directional response on the gradient steepness [27] or the midpoint concentration is observed [15,16]. In such models, nonlinear behavior is typically introduced directly at the level of the activator/inhibitor system. In comparison, our combined model allows for more flexibility in creating the desired response dynamics.
Similar to previous models of directional sensing, our model is restricted to the initial stages of gradient sensing. Note that this is only a limited part of the entire chemotactic process, which involves complex cytoskeletal mechanics and shape deformations leading to the formation of pseudopods and actual cell movement. An integral model that includes a description of cell polarization and motility remains the ultimate future goal of modeling efforts in this field. Parameter

Appendix B. Robustness
In the main text, several restrictions are imposed on the choice of the model parameters. Here, we will investigate the robustness of our model by relaxing the most severe of these restrictions.

B.1. Relaxing the condition K
In the main text, we required K + 1 = v/k 3 . This condition ensured that the fixed point p 3 became p 3 = 1, independent of the external gradient. We will now drop this restriction. The Michaelis constant K still depends on the external gradient, K = (k −1 q + k 2 )/k 1 , but no relation with the maximal reaction velocity v is assumed. The parameter v may take a constant value, independently of K . Solving equation (12) gives the general expressions for the two fixed points where we have introduced the notation e 1 = v/k 3 . As K depends on the external gradient, p 3 will now also change as a function of the gradient signal. Optimally, parameters are chosen such that the trigger wave velocity is zero in the middle of the cell. This leads to a clearly defined, equally sized front and back part of the cell under gradient stimuli. The sign of the trigger wave velocity is determined by the integral in equation (2) With f ( p) as given in equation (11), we find as the general form of equation (14) in the case of p 3 = 1. We now substitute the solution (B.2) for p 3 into equation (B.3) and introduce the notation γ = K /e 2 1 , where K ≈K q with K = (k −1 + k 2 )/k 1 ,k −1 k 2 . From equation (B.3) we then obtain as the condition for zero wave velocity. This equation can be solved numerically to yield γ ≈ 0.21. In the middle of the cell, we have q = 1 and K =K . Equation (B.4) thus restricts the choice ofK and e 1 , so thatK /e 2 1 ≈ 0.21 is fulfilled.

B.2. Dependence of p 3 on the external gradient
Relaxing the condition K + 1 = v/k 3 introduces a dependence of p 3 on the external gradient signal s =s[1 + α sin( 2π L x)], where x denotes the spacial coordinate running along the cell perimeter. How strong is this dependence? From equation (B.2) we can see that the variation of p 3 with the external gradient is proportional to 7 .

B.3. Fluctuations in the choice of parameters: robustness
After relaxing the condition K + 1 = v k 3 , equation (B.4) led to a constraint in the choice of the ratioK /e 2 1 . This constraint arose from the requirement that the trigger wave velocity should be zero in the middle of the cell. However, for a directional response to occur we do not need to fulfill this condition precisely. As long as the trigger wave velocity becomes zero somewhere inside the cell, we will observe a biphasic membrane response with one part of the membrane showing p = p 1 and the other part showing p = p 3 , although now the two parts will not be of 7 We maintain the assumptions thatk −1 k 2 and k i + k i + for i = a, b. equal size any more. Thus, we may allow for fluctuations in the parametersK and e 1 as long as equation (B.4) is satisfied at some position x * along the cell membrane, with 0 < x * < L and L the length of the cell perimeter.
Let us first assume that e 1 is fixed andK may fluctuate, so thatK = γ e 2 1 , i.e. the trigger wave velocity will not be zero in the middle of the cell. Away from the cell center, we have q = 1. Here, a zero trigger wave velocity may result ifK q = γ e 2 1 . Note however that the ratio of activator and inhibitor concentrations q = b/a is determined by the external gradient s =s[1 + α sin( 2π L x)]. The tolerable range of fluctuations inK will thus depend on the steepness α of the external gradient. In particular, we obtain 8 .
(1 − α)γ e 2 1 <k and (1 + α)γ e 2 1 >k From this estimate we conclude that we will observe a biphasic response as long as fluctuations inK remain on the order of the gradient steepness. Similar relations are found for fluctuations in e 1 , i.e. v and/or k 3 , while keepingK constant.