Deterministic Versus Stochastic Cell Polarisation Through Wave-Pinning

Cell polarization is an important part of the response of eukaryotic cells to stimuli, and forms a primary step in cell motility, differentiation, and many cellular functions. Among the important biochemical players implicated in the onset of intracellular asymmetries that constitute the early phases of polarization are the Rho GTPases, such as Cdc42, Rac, and Rho, which present high active concentration levels in a spatially localized manner. Rho GTPases exhibit positive feedback-driven interconversion between distinct active and inactive forms, the former residing on the cell membrane, and the latter predominantly in the cytosol. A deterministic model of the dynamics of a single Rho GTPase described earlier by Mori et al. exhibits sustained polarization by a wave-pinning mechanism. It remained, however, unclear how such polarization behaves at typically low cellular concentrations, as stochasticity could significantly affect the dynamics. We therefore study the low copy number dynamics of this model, using a stochastic kinetics framework based on the Gillespie algorithm, and propose statistical and analytic techniques which help us analyse the equilibrium behaviour of our stochastic system. We use local perturbation analysis to predict parameter regimes for initiation of polarity and wave-pinning in our deterministic system, and compare these predictions with deterministic and stochastic spatial simulations. Comparing the behaviour of the stochastic with the deterministic system, we determine the threshold number of molecules required for robust polarization in a given effective reaction volume. We show that when the molecule number is lowered wave-pinning behaviour is lost due to an increasingly large transition zone as well as increasing fluctuations in the pinning position, due to which a broadness can be reached that is unsustainable, causing the collapse of the wave, while the variations in the high and low equilibrium levels are much less affected.


Introduction
Many eukaryotic cell types undergo directed movement in a variety of scenarios. Such motility is important in embryogenesis (Charest and Firtel 2007), wound healing, immune surveillance (Ridley et al. 2003), and cancer metastasis (Ridley et al. 2003). As a first step in this process, cells polarise, forming a distinct front and rear distinguished by biochemical profiles of signalling molecules that regulate lamellipodial extension (Ridley 2006). An important part of that internal polarizing biochemistry is based on the activity and distribution of Rho GTPases. These switch-like signalling proteins exhibit a distinct active (GTP, membrane-bound) form and an inactive (GDP) form that is largely cytosolic. Only the active, GTP-form is able to interact with downstream effectors to exert its biological function. Interchange between these two forms is mediated by GTPase-activating proteins (GAPs), which augment inactivation, and guanine nucleotide exchange factors (GEFs), which facilitate activation. It has been established that the active form increases its own rate of activation via various self-recruitment mechanisms (Raftopoulou and Hall 2004;Li et al. 2003). While the active form binds the plasma membrane, the inactive form can be both in the membrane or released to the cytoplasm, a process which is positively regulated by binding to guanine nucleotide dissociation inhibitors (GDIs).
When a cell is stimulated, some Rho GTPase activity (notably, Cdc42 and Rac1) is focused at the leading edge (Ridley et al. 2003), inducing localized actin polymerization that generates protrusive forces propelling the cell (Raftopoulou and Hall 2004). Here, we are concerned about the onset of polarity and its maintenance, thus focusing only on the polarization of the Rho pattern, and not on the downstream remodelling of the cytoskeleton (or possible feedbacks that this might generate).
Based on general interactions between Cdc42, Rac, and Rho, and taking into account known parameters for the kinetics and diffusion of the active and inactive forms, we have shown earlier that sustained polarization within a cell is possible, even when the well-mixed system has only one equilibrium, and in the spatial setting this equilibrium is stable against both homogeneous and small non-homogeneous perturbations (Marée et al. 2006;Jilkine et al. 2007). Later, we determined the mathematical essence of the mechanism by studying a reduced deterministic model of cell polarization, coining it "wave-pinning" (Mori et al. 2008). It remains, however, unclear to what extent stochasticity at low molecule numbers can influence the potential of the mechanism to initiate and sustain polarity within the cell. We therefore compare and contrast the deterministic and the stochastic version of the core model for wave-pinning. A simple 1D geometry in which this generic Rho GTPase can be studied is shown in Fig. 1, where the organelles and nucleus are omitted, L is a cell diameter, and the chemical system is modelled by a two-component reaction with distinct rates of diffusion D a D b across L, since proteins diffuse much more slowly in the lipid membrane than in the cytosol. Here, A is the active and B the inactive small GTPase (with concentrations a(x, t) and b(x, t)). The height H and width W of the compartment are assumed to be reasonably small, so gradients are described in the x direction for x ∈ [0, L], t ≥ 0. The system of reaction-diffusion equations of the deterministic GTPase model in Mori et al. (2008) is Here, f (a, b) is the rate of GTPase interconversion. Equations (1a), (1b) are taken with no-flux boundary conditions: a x (0) = a x (L) = 0, b x (0) = b x (L) = 0 and the inter-conversion rate is modelled to include auto-activation of A (through a positive feedback of A onto its own production): where k 0 and δ denote the basal rates of activation and inactivation of A, respectively, γ is the rate of maximal feedback strength, and K is the concentration of A leading to a half-maximal feedback level. We first briefly describe the deterministic aspects of this model, and build on the previous analysis by introducing a local perturbation analysis that leads to insights on how the initiation of polarisation depends on the parameters and on the total amount of molecules. We then explore how the polarisation mechanism reacts when only a limited number of molecules is available and stochasticity starts to impact on the polarisation state of the cell. To do so, we describe and analyse an analogous stochastic version (low copy number regime) of the same model. We confirm our stochastic implementation by showing that simulating large molecule and lattice numbers approaches the thermodynamic limit. To analyse the equilibrium behaviour of the stochastic system, we introduce statistical tools which provide us with intriguing insights regarding the dynamics in the low copy number regime, namely it being dominated by spatial fluctuations of the transition zone rather than temporal fluctuations in the activity level, and loss of polarity due to the region of high activity, through stochastic fluctuations, reaching a broadness that is unsustainable, causing the sudden collapse of the whole wave. Bifurcation analysis of a simplified model of a pinned wave provides us with a straightforward rationale for the behaviour of the stochastic system close to the point where the wave is lost due to stochastic fluctuations.

Wave-Pinning
Given appropriate conditions, within model (1a), (1b) a stimulus-pulse of GTPase located within an otherwise homogeneous domain, for example at one end of the cell, leads to the formation of a travelling wave of activation that slows down and stalls, delimiting a spatial region of activation, i.e. creating a robustly polarised cell (Mori et al. 2008). In our simulations, this initial pulse is captured by a first-order reaction converting b to a within a small domain of the cell. This is done through the term k s b(x, t), which is both added to Eq. (1a) and subtracted from Eq. (1b). The wave-pinning regime depends on the relative rates of diffusion and the total amount T = L 0 a(x, t) + b(x, t) dx in the system. The mechanism of wave-pinning can be attributed to the following: relatively rapid diffusion of b (D b D a ) leads to a moreor-less constant level of b over the cell, while the existence of three roots of f (a, b) for the fixed well-mixed equilibrium b level allows for a sufficiently large local perturbation in a to locally reach a distinct activity level (a process that we have coined "Δ-perturbability", see below). Mass conservation ensures that, while this peak of increased a levels expands its domain over the cell with its front propagating like a wave, the more-or-less homogeneous level of b drops. This global decrease of b slows down and eventually limits the spatial propagation of the wave, pinning it at an equilibrium position (Mori et al. 2008). Even though wave-pinning requires the existence of three roots of f (a, b) for fixed b level, it is important to realize that it is not a consequence of bi-stability and subsequent front propagation between two stable states. (Note that in reaction-diffusion systems, the terminology bi-stability is used to denote cases in which the corresponding well-mixed system has two distinct stable steady states.) The well-mixed ODE system has only one equilibrium, and in the PDE this equilibrium is stable against both homogeneous and small nonhomogeneous perturbations. Nevertheless, in the spatial setting a sufficiently large local perturbation can trigger the travelling wave, which subsequently stalls, giving rise to sustained polarity.

Wave-Pinning Versus Propagation Failure
Because we compartmentalize space in this study to perform stochastic simulations, it is relevant to introduce yet another mechanism, coined propagation failure (Britton 1985;Keener 1987). As a possible source of confusion, propagation failure has previously also been referred to as "pinning of waves" (Fáth 1998), thus evoking the need to emphasize its clear distinction from "wave-pinning" as defined in Mori et al. (2008).
Propagation failure describes a specific phenomenon that can be observed in bistable systems in which travelling waves fail to propagate when space is discrete. This may occur when both the wave velocity is low and the discretisation of the space is coarse (relative to the diffusion coefficient) (Keener 1987;Fáth 1998). Under such conditions, propagation failure can manifest itself if, at the location of the wave front, the diffusive flux from one sub-domain into the next becomes insufficient to bring the levels of that sub-domain above the threshold required for the amplification and subsequent propagation of the wave. In contrast, the phenomenon of wave-pinning does not require a discretised space. Instead, when the triggered wave spreads over the domain, the velocity of the wave decreases, because of the drop in the available inactive form that is used up by being converted into the active form. Nevertheless, we here find that both phenomena become coupled to one another when space is discretised. Due to the slowing down of the wave during the wave-pinning process, inevitably the velocity of the wave eventually becomes sufficiently low that propagation failure will occur within coarse grids. Consequently, when we discretise space in this study, which we do in both numerical PDE simulations and in Gillespie simulations, propagation failure occurs for large sub-domain sizes as well as low diffusion rates.
Given that the sub-division into compartments is a computational method, but does not represent a biological property of the cell, we will ensure below that propagation failure does not play a role in the dynamics presented in this paper nor influences the biological insights we derive here. This brings us to the next issue, which is how to distinguish propagation failure from wave-pinning, given that in both cases the wave stalls.

Analysis of Polarity Initiation
The full bifurcation analysis of any system of partial differential equations (PDEs) is a challenging undertaking. While Mori et al. (2011) focused on the requirements of the travelling wave to stop, we will here discuss an analysis regarding the potential to initiate polarity and a travelling wave, in which we probe the homogeneous state of the cell with a local perturbation. In short, we ask what happens if a local perturbation is introduced to a resting cell (being at a uniform steady state), by observing whether such a perturbation will diverge to a distinct local equilibrium (eventually causing polarisation through wave-pinning), or alternatively dampen out, returning to the rest state corresponding to the global state of the cell. This analysis provides a straightforward test whether a (sufficiently large) perturbation can "invade" the initially uniform steady state solution. We refer to this reduced model as the "local perturbation analysis" (LPA) model, or system, as it allows us to study invasion criteria for a local perturbation of any given amplitude. (Note that such Δ-perturbability does not directly imply sustained polarisation through wave-pinning, see below.) To address the onset of polarisation without having to deal with the full complexity of the PDEs, we break down the spatial system into two effective compartments, one corresponding to the levels of the active and inactive form at the site of the local perturbation, and another corresponding to the global values over the rest of the cell. The simplified representation of the deterministic system is given in Eqs. (4a), (4b), used to predict the total amount of small GTPase T for which to expect initiation of polarisation and wave-pinning. This analysis allows us to compare results from the deterministic and the stochastic version of the model later on. For the local perturbation analysis, we make the following assumptions and approximations to the PDE model: • We ask whether the value of A at a site of the localized pulse a L (0) will diverge from the uniform global concentration of active GTPase a G (t). Since this active form has a very low rate of diffusion, we consider the limit D a ≈ 0 and treat a L (t) as a purely local variable, that can vary independently from a G (t). This is equivalent to assuming that any perturbation in A will be spatially confined to the site of the perturbation and will initially evolve independently of the rest of the domain. • Since the inactive GTPase B has a relatively fast rate of diffusion, we take the limit at which D b ≈ ∞ and consider b(t) to be a purely global variable (b L (0) = b G (t) ≡ b(t)). Restated, any local perturbation in B caused by the local perturbation in A will be instantly adjusted to the global, homogeneous concentration profile. This leads to the following LPA model: Furthermore, given that we consider only a narrow initial pulse of activation that hardly affects the overall cell levels, it is reasonable to approximate Eliminating b(t) by conservation leads to a system of two ODEs: The approximations required for the polarity-invasion analysis are depicted in Fig. 2. Note that since the activating pulse is confined to a sufficiently small subsection of the domain its variation will not affect the total amount of inactive form Fig. 2 Local perturbation analysis. Schematic representation of the deterministic local pulse analysis (LPA). In our reduction, we introduce a local perturbation Δa, and consider the evolution of the local variable a L , the inactive protein (represented by the dotted lines, evolves to the straight black line, such that b L = b G = b), and the active protein in the rest of the cell, a G and, therefore, the dynamics of B solely depends on the global level of A in the extended domain. We can use such approximations to address the following questions: under what circumstances would a localized pulse of activation grow in magnitude compared to the surrounding levels? How large should the amplitude of the stimulus be to trigger a new state (e.g. depicting an initial polarisation)? And, if a new, bounded state exists, what values do we expect it to have? These answers depend on the parameters and on the total number of molecules (i.e. amount of small GTPase) in the system. Note that we do not restrict attention to small amplitude perturbations (as is common in linear stability analysis, e.g. of reaction-diffusion equations).
The LPA system (4a), (4b) can be studied in several ways, e.g. via the a G a L phase plane and ODE bifurcation analysis. Performing such analysis reveals five qualitatively distinct regions of polarity behaviour, depicted in Fig. 3. The behaviour in each region can be captured with a sequence of corresponding a G a L -phase plane diagrams. On the bifurcation diagram, we plot the steady state value a L versus the bifurcation parameter T , representing the total amount of the small G-protein. The five distinct regimes of behaviour (labelled I-V) are separated by two saddle-node and two transcritical bifurcations. The shape of the curve traversing the diagram from lower left to upper right is shared with the bifurcation curve that would be obtained when the PDE system is taken to be well mixed (D a = D b = ∞). It represents a steady state where a G = a L (unpolarised cell, "rest state"). The two actually coincide, meaning that the equilibrium level in the local patch is neither lower nor higher than the uniform global background level of a G . Note that for any equilibrium found in the well-mixed case, there should exist a corresponding equilibrium a G = a L in the LPA model. The stability of the equilibrium, however, can change (from stable to unstable). For example, while the equilibrium in the well-mixed case is always stable, the portion of that curve in Region III of the LPA model is unstable. We indicate the line a G = a L in the a G a L planes, using a dashed grey line, representing the absence of any local perturbation.
We now explain how to interpret the diagram and its implications for polarity behaviour. (I) In region I the total amount of molecules is low (T < 19.09) and a regime is found in which only one steady state value a L = a G < 0.2 exists. That state of low activation is unresponsive to stimuli, and no pulse can "invade". The timespace plot of a(x, t) stays at, or rapidly returns to, a low uniform level (solid curve) showing steady-state local activated form, a L , versus the total amount of material in the domain (T ). Five distinct regimes of behaviour are found, as explained in the text. Top row: colour plots of the solutions of the PDEs, starting close to the uniform steady state in the given region. In these space-time plots, time axis is horizontal and space axis is vertical. Activity level is depicted using a colour gradient, with red indicating the highest and blue the lowest activity levels. Middle row: phase plane behaviour in the a G a L planes, showing the number and stability of steady states of the reduced system in each of the regions. In regions I and V, only a uniform level of global activity is stable, and no pulse or stimulus can grow. In the intermediate regions, only a sufficiently large pulse (II), or a sufficiently low dip (IV) can grow. In region III, the homogeneous state is unstable to any perturbation, and a variety of patterns can form, depending on initial conditions. Bottom graph: Wave-pinning position as a function of T . In regions i and iii no sustained polarised state is possible, while in region ii wave-pinning occurs, with the position of pinning, L 0 * monotonically increasing with T (Color figure online) no matter how large the amplitude of an applied stimulating pulse. That is, a cell will not polarise when stimulated, it remains unpolarised, in a rest state. (II) This is the region corresponding to the deterministic regime of cell polarisation (wave-pinning). Here, for an intermediate level of substance 19.09 < T < 23.0, there are three coexisting equilibria of Eqs. (4a), (4b), the outer two of which are stable. A value of a G corresponding to the lower branch can be "invaded" by a local pulse, provided its amplitude is large enough to surpass a threshold depicted by dashed curve in region II. The lower branch corresponds to the a G = a L equilibrium, i.e. to a cell that is in a homogeneous rest state. Thus, the rest state is stable against small perturbations, but sufficiently large perturbations can polarise the cell. Note that as the total amount increases, the required amplitude to trigger polarisation decreases sharply, so that close to but just short of T = 23.0 a pulse of very small size can lead to polarisation. We see that the full PDE solutions (top panels in Fig. 3) show the invasion of such a pulse in this regime, which becomes established at a finite amplitude over some fraction of the domain. (III) Other patterned states (e.g. with one or more patches of active GTPase) occur in this region. For 23.0 < T < 25.99, the global steady state a G = a L is un-stable to any perturbation. Here, small amplitude noise or a pulse of small magnitude will disrupt the global state leading to other patterned states. This kind of behaviour is typical of a Turing instability. Indeed, the full PDE solution (with random noise initial conditions such that the total amount falls in this range) produces patterns with multiple peaks. (IV) For even higher values, 25.99 < T < 35.58, the total amount of GTPase is so high that the global level of activation is at an elevated steady state level (highest solid branch of the diagram). Here, an invading "pulse" has to locally deactivate a region in order to "invade" (i.e. the pulse is a dip below the uniform global level). The amplitude of that "dip" must cross the threshold (dashed portion of curve) to trigger the polarisation, as otherwise it decays back to the uniform activation level. As shown in the solution of the PDEs, a dip of sufficiently large amplitude leads to a stable local patch of depressed activity in an otherwise high global level of activity. (V) Finally, above T > 35.58 the potential of polarisation is lost again. That is, no pulse or dip can invade, and the uniform global state is one of high activity everywhere in the domain. The LPA does not address the question at which position the wave will be pinned, but rather if a wave can be triggered and how high it will become. The next question therefore is at which position along the cell length the wave stalls. We indicate the wave position by L 0 , and the equilibrium value of L 0 at which the wave stalls L 0 * . In Mori et al. (2008), the wave-pinning position has been derived mathematically for the limiting case of an infinite difference in diffusion rate between the active and inactive form (i.e. using a sharp front approximation). In the bottom panel of Fig. 3, we show the steady state value L 0 * as a function of T . Regarding the wave-pinning itself, three regions of qualitatively different behaviour can be discriminated. In regions i and iii, no stable polarity can be found, because any wave would completely retract or expand over the whole domain, respectively. In contrast, in region ii we find the possibility of a stable co-existence of a high and a low state is found. Note that the pinning position L 0 * depends on the value of T . Importantly, the figure shows that the interval of sustained polarity is smaller than the interval of regions II-IV. It illustrates that even when a wave can be triggered, it does not always follow that it can also be sustained. The membrane is divided into N equally wide well-mixed sub-domains. Diffusion within the membrane is modelled as first-order reactions while we assume that B diffuses fast enough to warrant modelling it as occupying a well-mixed cytoplasmic pool. Dimensions H and W associated with the domain of length L are necessary for conversion of concentrations to numbers of molecules but do not play any other role: all stochastic simulations are one-dimensional Note that both the results of the LPA and of the wave-pinning position act as an approximation to the actual PDE behaviour, where actual rates of diffusion are finite and initial conditions can affect whether initially a single peak or several peaks emerge. However, it correctly captures the basic boundaries that determine its potential to polarise, the minimum perturbation amplitude required to do so, the expected values to be reached at the local perturbation, and the position at which the wave pins. We will show how this brings valuable insights when interpreting the role of stochasticity in cell polarisation by wave-pinning.

Stochastic Version
Next, we ask how the same polarisation mechanism would behave in the low copy number regime. We ask under what conditions a stochastic equivalent of the deterministic model still presents wave-pinning, i.e. after triggering the formation of sustained regions within the cell with respectively low and high levels of the active form, and if our approach predicts biologically relevant conditions under which wave-pinning may be unsustainable in live, noisy cells.
For stochastic simulations of our system, we resort to the stochastic formulation of chemical kinetics (McQuarrie 1967), and use the stochastic simulation algorithm (SSA) due to Gillespie (1976)  In the deterministic case, we choose experimentally-supported (Postma and van Haastert 2001;Postma et al. 2004) relative diffusion values D b = 100D a that effectively render B homogeneous. In our stochastic simulations, we approximate B as a homogeneous cytoplasmic pool (D b → ∞) to reduce the computational cost, while maintaining a full computation for the heterogeneous A distribution. Diffusion of A in the membrane is treated as a series of first order reactions (spatial SSA )), each with propensity where A i denotes the current number of molecules of A in membrane lattice point i. The above propensity equals zero for diffusion to the left when i = 1, and diffusion to the right when i = N , respectively (no-flux boundary conditions). Similarly, we choose the remaining propensities according to Gillespie (1976): background inactivation out of membrane lattice point i: δA i , pulse activation into membrane lattice point i: k s B N .
We note that B, the current number of molecules of the inactive species in the well-mixed cytoplasmic pool, is rescaled with 1/N because each membrane lattice point only senses this fraction of the available total number of molecules of B. For each propensity p, the probability that the corresponding event occurs within the next dt units of time equals p · dt + o(dt), where o(·) denotes terms that converge to zero quicker than its argument (little o notation, lim dt↓0 o(dt) dt = 0). In the above propensity expressions, most kinetic constants are equal to those used in the deterministic system, Eqs. (1a), (1b), since they are independent of the units of a(x, t) and b(x, t), or A i and B, correspondingly. However, the Michaelis-Menten constant K of Eq. (2) has the same units as a(x, t) and needs to be rescaled for our stochastic simulations: where N A is Avogadro's constant, V a denotes the volume associated with each membrane lattice point, V a = L N W H 2 , and the factor 10 −21 is required to re-scale units (V a has units of µm 3 , and K units of µM). Homogeneous initial concentrations of A and B, a(x, 0) = a 0 and b(x, 0) = b 0 , are rescaled to numbers of molecules equivalently: Note that in our stochastic simulations we need to associate volumes with each lattice point (both for the membrane and for the cytoplasm), since it reformulates a concentration-based model, Eqs. (1a), (1b), as a molecule-based stochastic model. The natural choice for conversion between the two is through proportionality to the dimensions (volume) of the system. Even though we specify a volume for each lattice point in this conversion, the stochastic simulations are effectively one-dimensional, as are our deterministic simulations. That is, in both the deterministic and the stochastic system, we focus on radial polarisation along the diameter of a cell.

Results
For our simulations and comparison between the stochastic and deterministic case, we use parameter values based upon Mori et al. (2008), and summarized in Table 1: is linearly stable (a 0 = 0.2683 µM). We fix the slab height H = 0.2 µm and length L = 10 µm (see Fig. 1) and simulate our stochastic system for varying numbers of molecules by varying the width W of our cell slab. Varying the width W (Fig. 4) allows us to change the number of molecules without changing the initial concentrations of A and B.

Propagation Failure Does Not Affect the Deterministic Simulations
We first determined that the discretisation of space utilized does not cause propagation failure around the parameter values used for our analysis on stochasticity. To do so, we make use of the fact that the initial perturbation that triggers a wave does not necessarily have to be small in width. A wide (but not too wide) perturbation of sufficient amplitude can also trigger a wave. Perturbations that are wider than the final pinning position, however, trigger waves with a negative velocity, i.e. waves that retract until they come to halt at the pinning position. If indeed propagation failure plays a role, both the extending and the retracting wave are expected to halt before their velocities would have become zero in the continuous case. Thus, propagation failure would cause both waves to a halt at distinct positions. We therefore devised the following numerical experiment: We use a total amount of T = 22.68, our default value, from Mori et al. (2008), and within the wave-pinning regime, in which a Δ-perturbation is required to trigger a wave that gives rise to sustained polarity; see the black arrow in Fig. 3. We then start simulations, for varying values of the active form diffusion coefficient D a and number of lattice sites N , with two different wave-shaped initial conditions (see Fig. 5): one narrower and higher, one broader and lower, both indicated by dashed lines in the figure, where the former lies to the left and the latter to the right of the equilibrium wave-pinning position. Figure 5 shows the equilibrium wave profiles for the narrowly initiated waves (thick red lines) and broadly initiated waves (thin blue lines). Indeed, sufficiently large box sizes (low N ) and sufficiently low diffusion coefficients show a discrepancy in the final position of the pinned wave between the retracting and the extending waves, illustrating how too slow diffusion combined with too coarse subcompartmentalization leads to propagation failure. For the default values used in this study (N = 50; D a = 10 −1 µm 2 /s), however, no propagation failure can be observed. A 10-fold coarser grid or 100-fold slower diffusion would be needed for propagation failure to occur.

Equilibrium Behaviour of the Stochastic System
Individual SSA runs (Fig. 6, top) resemble the deterministic behaviour for sufficiently large numbers of molecules (≈6800), and the averaged concentration profile clearly shows wave-pinning after a stimulus is applied (Fig. 6, bottom). With fewer molecules (≈700), individual SSA runs become more erratic and the averaged profile The corresponding mean behaviour in the bottom image confirms this. The bottom inset depicts loss of wave-pinning from the mean behaviour when too few molecules are present in the system is homogeneous across the domain, resulting in loss of wave-pinning (Fig. 6, bottom, inset). As we increase the numbers of molecules by increasing W , we observe convergence of the stochastic system to the deterministic one, indicated by a decreasing observed variance and smoothing of the black solid curve. To study the distribution of pinning positions, we fit, well after wave-pinning, the concentration profile of the active species to the symmetric Richards model (Richards 1959), which very well captures sigmoidal profiles and contains a parameter, c, that defines the position of the inflection point, which we consider to be the pinning position of the wave: where M and m are the concentrations at the high and low plateau, respectively, h determines the slope at the inflection point, and c is the position of the inflection point. The insets in Fig. 6, top, show the distribution in pinning position for the indicated number of molecules, based upon concentration profiles from 100 simulation runs collected between 150 s and 200 s and recorded every 0.1 s (i.e. based upon 50, 100 concentration profiles). We fit it to a normal distribution, indicated by the black line. At high molecule numbers, the distribution closely follows a Gaussian distribution, with the mean corresponding to the predicted value from the PDE analysis (indicated by the dashed line). At low molecule numbers, however, we find that the distribution becomes leptokurtic, illustrating that the wave is more confined to its pinning position than is to be expected from the level of variation observed. Also, at small molecule numbers, the pinning position becomes slightly biased to smaller L 0 * (i.e. to the left), due to the fact that stochasticity also creates an effectively lower level of auto-activation; see the discussion below. At a high molecule number, fluctuations in the pinning position become small and closely spaced around the position within the deterministic model (indicated by the dashed line). At lower molecule numbers, the profile broadens.
To quantify convergence of the stochastic case, A i (t), to the deterministic case, a i (t), at time t and lattice point i, we use the normalized Euclidean distance We compute this distance for J = 100 SSA runs, d(t) (j ) (j = 1, . . . , J ), and report mean observed distances d(t) = (1/J ) J j =1 d(t) (j ) . The observed mean Euclidean distance in Fig. 7 shows the expected convergence of the stochastic case to the deterministic case for increasing copy number. In Fig. 7, disappearance of the wave profile becomes apparent by the increasing value of d(200) for fewer molecules. These results show the stochastic limit to the deterministic behaviour, confirming the qualitative intuition of a stochastic lower limit to the wave pinning mechanism. What is not clear a priori, however, is the dynamics by which wave-pinning is lost at low molecule numbers and how this is affected by other parameter changes.

Stochastic Simulations Are Not Affected by Propagation Failure
Before we discuss in detail the loss of wave-pinning at low molecule numbers, we revisit the phenomenon of propagation failure. We repeated the procedure described Fig. 7 Convergence of stochastic to deterministic system. Mean Euclidean distance over J = 100 SSA runs, d(200) at time t=200 s, after the stimulus. Simulation parameters as in Fig. 6, except for W . (W is varied to obtain different total copy numbers.) This plot shows convergence of the stochastic system to macroscopic predictions as we increase the number of molecules above to determine potential propagation failure within the stochastic description, but now for varying molecule numbers rather than diffusion rates (Fig. 8). Indeed, as before, when the grid becomes rougher, secondary effects of propagation failure manifest themselves and affect the pinning positions through this independent mechanism. Again however, as we showed before in Fig. 5, those effects are not noticeable in the simulation regimes we focused on, and only to a very small extent under extreme conditions of a very coarse grid.
We then studied how stochasticity influenced this relationship. We found that stochasticity reduces rather than increases the parameter regime for which propagation failure can be observed (compare the column corresponding to N = 5 in Fig. 8 with the upper left graph in Fig. 5). This can be understood by realizing that stochasticity can help overcome the threshold to propagate the wave. This means that for coarse grids and low diffusion rates the stochasticity at low molecule numbers can even increase the precision of the pinned wave. Besides deviations at low box numbers, we also observe deviations from the stationary solution of the deterministic model (as indicated in each frame with a black line) for a combination of low molecule and high box numbers (see lower right graph in Fig. 8). This is an artefact attributed to an effective change in the auto-activation function when the number of molecules in a box becomes very small, which will be discussed further below.

Comparison of Deterministic and Stochastic Predictions
We test our predictions from the LPA system, Eqs. (4a), (4b), by classifying individual SSA runs (using the same parameter values as in Fig. 6) as either homogeneous in A (i.e. uniform in A well after the stimulus, e.g. at t = 200 s) or inhomogeneous in A (i.e. where a local pulse invaded the global concentration profile of A, creating at least one high plateau or peak in A).
Comparing the predictions of the deterministic version of the model with the SSA simulations we find interesting differences. With parameter values as in Fig. 6, specifically with D a = 0.1 µm 2 /s (central inverted cup in Fig. 9), we observe a somewhat . This discrepancy may be explained by the relatively fast diffusion of A, compared with the limiting rates used in the theoretical treatment, which destabilizes local perturbations in A (through high curvature in the concentration profile of A), making it harder for perturbations to stabilize and invade the remaining profile of A. Indeed, as we decrease D a by orders of magnitude (Fig. 9, broadening inverted cups), increasingly wider ranges of T allow for stabilization of perturbations in A and their invasion of the global homogeneous A level. Our observations for decreasing D a , and certainly for D a = 0 µm 2 /s (Fig. 9, dashed line), indicate increasing agreement between our stochastic simulations and the predicted wave-pinning criteria.

Loss of Wave-Pinning in Small Number Regimes
Loss of wave-pinning for fewer molecules may either result from increased randomness along the concentration axis (vertical noise, Fig. 6), or greater random fluctuations in the pinning position of a travelling wave (horizontal noise, Fig. 6): With a Fig. 9 Dependence of wave-pinning stability on diffusion and total amount. Fraction of SSA runs (out of 100 runs) observed at t = 200 s to be stably inhomogeneous (wave-pinned, or with multiple peaks in A). System parameters as in Fig. 6 with varying D a . Each inverted-cup-shaped curve corresponds to a different D a value, decreasing from D a =0.1 µm 2 /s to D a = 10 −5 µm 2 /s from the innermost to the outermost cup, using order of magnitude changes in the parameter value (as indicated by the arrows). Dashed line: D a = 0 µm 2 /s. This plot shows for what parameter values (changing total amounts T and diffusion constant D a ) the stochastic system shows stabilization of perturbations in the uniform concentration profile of A (high fraction of inhomogeneous runs). Stabilization of perturbations is difficult for high diffusivity of A (D a = 0.1 µm 2 /s innermost, narrowest cup) since spikes in the concentration profile of A are more readily smoothed out before they grow to sufficiently high levels. For decreasing values of D a , our observations match progressively better with the predictions made using the LPA in the main text. That is because the LPA assumes zero diffusivity for A and, therefore, ignores any effects due to the diffusion of the active form low copy number, individual runs may show wave-pinning but the steepest point of the concentration profile (pinning position) may fluctuate along the horizontal axis due to inherent stochasticity.
To distinguish the effects caused by these distinct types of noise, we quantify the horizontal noise at lattice point i, using a sample of J = 100 SSA runs, with a function similar to auto-correlation, here denoted "spatial auto-correlation": Our estimator of spatial auto-correlation then becomes: This function centres observations A i (m) about the spatial mean and gauges the correlation between observations k time steps apart. We also normalize our estimator ( RS k (i) ∈ (−1, 1)) by dividing by the observed variance about the spatial mean. A high-value RS k (i) denotes a lattice point i where A i is stably far away from the spatial mean, while a low value suggests that A i repeatedly comes near the spatial mean.
We report estimates of auto-correlation (spatial, and temporal further on) as sample means RS k (i) = (1/J ) J j =1 RS k (i) (j ) of the corresponding estimates for J SSA runs, RS k (i) (j ) . We further compute these estimates in time periods when we expect our system to be stable, long after the application of a stimulus and wave-pinning (stimulus applied between 50 s and 70 s and observations between 500 s and 1500 s used for computations).
In a regime where the mean behaviour presents wave-pinning, we observe a dip in spatial auto-correlation in a section of the domain that includes the pinning position (Fig. 10, top). The width of this dip (shaded area in Fig. 10, top) increases with decreasing numbers of molecules in the system (Fig. 10, bottom), and spans the entire domain (10 µm) when wave-pinning is lost.
For quantifying vertical noise, we use an auto-correlation function denoted "temporal auto-correlation" for clarity: where E[A i ] is the expected value and Var[A i ] 2 is the variance of concentration A i . We use the common sample mean and sample variance as estimators of these statistics: Using these, our estimator of temporal auto-correlation becomes  figure) over 100 SSA runs, as a function of total number of molecules. When the dip width approaches the length of the domain of the cell (10 µm), the pinning position stops conferring information that can be used by the cell. We note a sudden and drastic increase in width at approximately 2,000 molecules. Simulation parameters as in Fig. 6, except for W which is varied Temporal auto-correlation measures the randomness of the time evolution of A, which is governed by a continuous-time Markov process. (If we know the distribution of A(t) at present, the future distribution of A(s), s > t, only depends on A(t) and is independent of observations of A before t.) Given the Markovian character of A(t), we expect the temporal auto-correlation to be greatly dependent on the lag k: a small lag (k = 1 s) results in relatively high temporal auto-correlation (black area, Fig. 11 top right and bottom right), while a bigger k = 20 s yields small temporal auto-correlation (light grey area, Fig. 11 top right and bottom right). We also observe that the Markovian character of A(t) is comparable for small and large numbers of molecules since the magnitude of temporal auto-correlation does not change significantly when altering the number of molecules (Fig. 11 right panels). This means that in both small and large copy regimes, temporal auto-correlation decreases comparably fast. While the vertical noise does not seem to decrease when increasing the number of molecules, wave-pinning still shows up as a marked peak in temporal auto-correlation in the transition zone of the wave (Fig. 11, bottom right). Since overall vertical noise is relatively constant, this peak seems to be caused by the wave fluctuating about its pinning position: concentration A i (t) in the transition zone fluctuates between high In the bottom row and for a sufficiently high number of molecules, spatial auto-correlation reveals a clear pattern in the equilibrium state with the dip in autocorrelation robust across different lags k. In the same copy-number regime, temporal auto-correlation shows a pattern for small lags k which, however, vanishes for increasing values of k. Simulation parameters as in Fig. 6, except for W and low values (high and low plateau of wave) causing A i (t) to be far away from its average repeatedly (high temporal auto-correlation).
We find spatial auto-correlation to behave markedly different from temporal autocorrelation for increasing numbers of molecules (Fig. 11, left panels). While spatial auto-correlation is comparable in magnitude to temporal auto-correlation in a small copy number regime (Fig. 11 top panels), we observe much greater spatial autocorrelation than temporal auto-correlation for various lag k values in a large copy number regime (Fig. 11 right panels). The stably high spatial auto-correlation, even for large lags k, far away from the pinning position suggests that the high and low plateau of the wave are persistent in wave-pinning regimes (Fig. 11 bottom right). The random fluctuation of the pinning position is highlighted by the decreasing spatial auto-correlation for increasing lags k (dark to light grey area in Fig. 11 bottom left) in this part of the domain.
To elucidate the source of the horizontal noise and the impact it has on the sustainability of a polarised cell, we make use of another spatially simplified representation of the deterministic system. Without loss of generality, once a wave has been formed, the concentration profile of a can be split into a plateau of length L 0 and approximate concentration a L on the left, and a plateau of length L − L 0 and approximate level a R on the right (Fig. 12 left) where Here, a 1 > a 2 are two roots of f (a, b wp ) = 0, where b wp is the uniform concentration of b for this wave-shaped profile, while it is assumed that D b = ∞ and D a = 0. Then the total amount T of protein in the domain can be approximated as T wp = Lb wp + L 0 a 1 + (L − L 0 )a 2 . In the limit D a →0 (no direct communication between plateaus), a L (t) and a R (t) evolve independently. Eliminating b(t) using mass conservation leads to: which also implies that the position of the wavefront is System (12) is a second deterministic reduction that leads to a way of comparing the stochastic and deterministic model versions. In the deterministic case, we expect that the wave-pinned configuration, Fig. 12, left panel, is stable over time when L 0 = L 0 * (pinning position). As Eq. (13) indicates, noise in the concentration variables will propagate to the width of the high plateau of the pinned wave. We expect that this propagation of noise is qualitatively the same in the full spatial system and conjecture that noise in the pinning position (horizontal noise) is the result of noise in the concentration levels (vertical noise). Given that the total amount is fixed at T wp , a L and a R adjust to varying widths of L 0 : increasing L 0 will typically decrease a L and increase a R accordingly, and vice versa. For system (12), we plot the steady-state concentration of a L as a function of L 0 (Fig. 12, right panel) and observe two critical values for L 0 , L (1) 0 < L (2) 0 , at which saddle-node bifurcations occur. The bifurcation plot in Fig. 12 shows that if we initialize our simplified system in a wave-pinned configuration (high plateau on left, low plateau on right, as in Fig. 12 We expect that the saddle-node bifurcation at L (2) 0 and hysteresis explain the sharp increase in dip width (Fig. 10, bottom) observed in the full spatial system: decreasing the number of molecules in the system increases vertical noise which propagates (2) 0 = 6.75 µm. The corresponding bifurcation plot of a R is symmetric to this one of a L . In our stochastic simulations, for equivalent parameter values and number of molecules = 6,820, we did not observe pinning positions greater than 6 µm (top inset, Fig. 6), as is predicted by this analysis into horizontal noise. As horizontal noise increases, the likelihood that the pinning position randomly overshoots the critical value increases and we are more likely to observe wave collapse (dip width approaching 10 µm, Fig. 10, bottom). This prediction is further supported by our observation that under conditions equivalent to those of Fig. 12, and with 6,820 molecules, we do not observe any pinning positions greater than 6 µm (Fig. 6, top inset). Figure 12 is closely linked to Fig. 3. While the latter shows the impact of a Δperturbation of infinitesimally small width, corresponding to an L 0 = 0 µm, as a function of T , the former shows the impact of perturbations of varying width, for a fixed value of T . The link between the analysis on the loss of wave-pinning with system (12) and the analysis on polarity initiation (LPA, system (4a), (4b)) is very direct. By allowing L 0 ↓ 0, system 12 becomes equivalent to system (4a), (4b): a L in Eq. (12) describes the active level on a vanishingly small domain and, therefore, ceases to affect the inactive species b (i.e. a L -active left-becomes the local perturbation a L -active local-in Eqs. (4a), (4b)), while a R starts occupying the entire domain of length L and, therefore, only the presence of a R affects b (i.e. a R -active right-becomes the global active form a G -active global-in Eqs. (4a), (4b)). Figure 13 brings both pieces of information together, showing a two-parameter bifurcation plot, with L 0 along the x-axis, and T along the y-axis. The right panel of Fig. 12 corresponds to a horizontal cross-section at T = 22.68 (Fig. 13, zone VI), while the bifurcation diagram of Fig. 3 corresponds to a vertical cross-section through this bifurcation diagram at L 0 = 0, along which line we have indicated the two fold bifurcations (FB) and two transcritical bifurcations (TB) that can be seen in Fig. 3, indeed reproduced at exactly the same parameter values.
Increasing L 0 away from zero, the critical total amounts of both fold bifurcations in Fig. 3 change while the critical values of the transcritical bifurcations are unaffected (solid grey lines in Fig. 13). It implies that at levels of T that are less favourable for sustained polarity, narrow Δ-perturbations are still able to trigger a  12)) as a function of L 0 for a fixed value of T within that zone wave, while broad ones can not do that any more. Moreover, it illustrates that when the well-mixed equilibrium is unstable against spatially inhomogeneous bifurcations (zones III-V), the width of the perturbation becomes irrelevant. We also observe four cusp bifurcation points (CP), describing where two fold bifurcation lines merge. These points imply the possible coexistence of waves with different amplitudes. Finally, there are two bifurcations where a fold bifurcation and a transcritical bifurcation collide (FT). They are linked to the symmetry of the two-parameter bifurcation plot about L 0 = 5 µm, which is due to the lack of inherent bias for either a leftoriented or a right-oriented polarisation in system (12): when L 0 < 5 µm, the pinned wave has its high plateau on the left, while for L 0 > 5 µm the high plateau is on the right. Due to this symmetry, we for example observe a bifurcation plot equivalent to Fig. 3 when plotting a R and setting L 0 = 10 µm (data not shown), confirming that initiation of a wave from the left is equivalent to initiation from the right. Together, Fig. 13 reveals that there are seven qualitatively different zones of levels of T , each presenting diverse requirements on wave initiation and maintenance. It allows us to predict the potential to trigger (or sustain, when L 0 = L 0 * ) a wave through a perturbation of any possible width and height, as well as the expected height that such a wave will reach while it travels, stalls, or stochastically fluctuates. Specifically, it sets the boundaries for horizontal fluctuations to trigger a collapse of the polarised cell state.

Validity of Stochastic Model at Large Compartment Numbers
Sub-dividing the domain into N compartments is a computational method which should not influence the biological insights derived here. Above we discussed that the coarseness of the lattice could introduce artefacts due to propagation failure, but that the simulations in this study are sufficiently fine-grained. Figure 8, however, also presented deviations from the expected wave profile when N was large. When testing the behaviour of the stochastic model as N → ∞, we realized that this is due to the fact that the effective auto-activation function starts to change as the number of molecules in a box becomes very small. Given that at high N only small variations between the mean values of neighbouring boxes are expected (as the flux between boxes goes to infinity), we tested if the change in auto-activation could be due to the expected stochastic variations within each box, assuming a Poisson distribution for the number of active molecules within each box. We therefore calculated the predicted mean activation rate as a function of the concentration a for different box sizes, by taking the kernel of the Poisson distribution with the auto-activation function itself. The insets of Fig. 14 show the resulting auto-activation rates for different box sizes, to which subsequently the auto-activation function γ a 2 K 2 +a 2 was fitted, with both K and γ as the fitting parameters. The figure shows that when the number of active molecules becomes small the functional response effectively shifts to the right (i.e. the fitted value of K becomes larger), while the plateau hardly changes (i.e. γ remains more or less 1.0). This phenomenon is due to the plateau in the sigmoidal activation function: while the high spectrum of the Poisson distribution cannot further increase the activation rate, as it is capped, the low spectrum decreases it. Figure 14 repeats the analysis of Fig. 12 for the estimated effective values of K. It shows that consequently at lower molecule numbers the wave becomes lower and is lost more easily due to horizontal fluctuations until, at very low molecule numbers per box, the wave cannot be sustained any longer. Also, the lower critical value of L 0 , L (2) 0 , shifts to lower values (i.e. to the left). This shift to lower L 0 and drop in height of the wave closely corresponds to the observed changes in the stochastic simulations as can be seen in Fig. 8, suggesting that the side-effects observed when increasing N are due to this modification of the auto-activation function only. This reduction in the number of molecules per box when increasing N can be overcome in our modelling set-up by increasing W . Thus, N can be made arbitrarily large, as long as, through modifications of W , the number of molecules per box is maintained above around 14 (below which value the effective K becomes too large to warrant sustainable wave-pinning).

Discussion
In this paper, we compared deterministic and stochastic aspects of a model for cell polarisation of the wave-pinning class (Mori et al. 2008). This work was motivated by recent interest in the influence of stochastic noise in biological systems. In considering stochastic noise, we account for possible effects due to low-copy numbers of signalling proteins, as has been done for instance by Isaacson et al. (2011). The different bifurcation diagrams are made by varying the parameter values γ and K of the Hill-type auto-activation γ a 2 /(K 2 + a 2 ). The effective parameter values γ and K for small box sizes are determined as follows: For the mean number of molecules of active A being equal to 0, 1, 2, . . . molecules per box (indicated by the black dots in the insets), a Poisson distribution for that number of molecules is assumed in order to calculate an expected rate of auto-activation. These observed values are then fitted to the auto-activation term itself using a least-squares fit. The fitted auto-activation functions are shown as grey lines in the insets, together with the residual sums of squares, R 2 . We observe that while the fitted γ hardly differs from the originally used value (γ ≈ 1), the fitted K increases as box volume decreases. The bifurcation diagrams show that this change in effective auto-activation lowers the higher plateau and narrows the range of permissible L 0 values. Note that the difference between the lower panels is virtually undetectable. Indeed, both are almost equivalent to the bifurcation diagram for the deterministic system, which was shown in Fig. 12. It illustrates that a volume of V a = 0.1 µm 3 (bottom left) forms the upper bound for observing this stochastic effect. This corresponds to 137 molecules per box, given that the total concentration is 2.268 µM. The lower bound for wave-pinning to occur is reached when the effective K becomes 1.07, at a box size of V a = 0.009 µm 3 , corresponding to 14 molecules per box Recent studies indicate that noise can have either constructive or detrimental effects in biological systems. For example, noting beneficial effects, Paulsson et al. (2000) observed that stochastic focusing increased sensitivity of cascades, Rao et al. (2002) found that noise-induced population heterogeneity improves fitness, Howard and Rutenberg (2003) argued that biologically relevant oscillations in a twocomponent dynamical system are more robust in the stochastic case than the deterministic one, and Gamba et al. (2005) showed that stochasticity could play a role in chemotactic responses to shallow gradients. On the other hand, detrimental effects were noted by, for example, McAdams and Arkin (1997) who showed that gene expression in a noisy regime resulted in bursts, rather than constant levels of gene expression. For our stochastic model of cell polarisation, we observe that at critically low molecule numbers stochastic noise has an impact on the behaviour of the system in a detrimental manner, eventually destroying polarisation.
As shown in Fig. 6, we verify correctness of our stochastic implementation by comparing our stochastic simulation results in large copy number regimes with the deterministic system (approaching the macroscopic limit). The compartmentalized spatial Gillespie model used here, however, does present deviations when the lattice number N becomes very small or very large. Alternatively, an off-lattice Brownian dynamics model could have been developed (Andrews and Bray 2004;van Zon 2005), to independently confirm the results based upon space discretisation presented here, given that both approaches present their own limitations Chapman 2007, 2009). However, to re-formulate the auto-activation in terms of massaction kinetics only, which is required to use off-lattice alternatives, would require a replacement of the current parsimonious non-linear term by, for example specific enzyme kinetics. Such a description would involve the introduction of new variables and new assumptions as well as relatively ad hoc choices regarding their reactions and behaviour, obfuscating the comparison between the PDE model and the stochastic model.
Here, we studied wave-pinning in a one-dimensional slab of cell material representing radial positional information available to a single cell. Our results suggest that for small copy numbers, radial information within a cell regarding its front and back that is available through the wave-pinning process decreases in quality because of fluctuations of the wave around its pinning position. If we extrapolate our results to a spherical 15 µm-diameter cell and assuming micro-molar small G-protein concentration (corresponding to 10 6 molecules), stable polarisation should typically be observed. However, in vivo effective reaction compartment size is often restricted due to macro-molecular crowding and resulting volume exclusion (Schnell and Turner 2004;Grima 2010): the volume encompassed by a cell is generally occupied by a range of macro-molecules which do not participate in any of the relevant chemical reactions. As these macro-molecules span the cell in a mesh-like fashion, individual effective reaction compartments may emerge which have a potentially small volume. Our results indicate that sufficiently small effective reaction compartments (those that hold 10 3 molecules, bottom Fig. 10), will produce inaccurate positional information (fluctuation of the pinning position). Hence, a cell subject to great amounts of macro-molecular crowding may integrate inaccurate positional information from its individual effective reaction compartments and, therefore, lose a global sense of directionality.
Instead of a gradual loss of wave-pinning, we observe a threshold number of molecules (between 2,000 and 3,000 molecules) below which the wave is suddenly lost. Analysis of a limiting deterministic case shows that sudden loss of wave-pinning is due to saddle-node bifurcations with hysteresis. We conjecture that this also explains sudden loss of polarisation in our spatial stochastic system. Altschuler et al. (2008) studied a related positive feedback model of Cdc42 with a homogeneous cytoplasmic pool of the inactive form. Their model includes selfrecruitment of the active form, but the conditions for polarisation and wave-pinning as defined through the LPA method described in this paper are not fulfilled, hence no stable polarisation can be observed. Nevertheless, by defining polarisation as a transient situation where 10 % of the domain holds more than 50 % of the molecules, they were able to show that within their model decreasing the numbers of molecules, which increases the random fluctuations, could trigger polarisation. This is opposite to what has been shown in our study, in which increasing fluctuations shuts off the polarity. For a fixed positive feedback strength, they found that 1,000 molecules yields the maximum probability (fraction of simulation runs) for polarisation. For 3,000 molecules or more, they observed that fewer than 50 % of the runs would polarise. In contrast, our model predicts that polarisation fails below 2,000 molecules (Fig. 10, bottom panel). Moreover, in our model, the behaviour of the stochastic system more closely resembles that of the deterministic system when the molecule number increases (bottom panel of Fig. 6), i.e. unlike Altschuler et al. (2008) with more molecules the polarity becomes increasingly more robust. Khain et al. (2011) recently discussed stochastic travelling waves in a spatial onedimensional model of spruce budworm populations, in which the high plateau of the wave corresponded to parts of the environment with a great number of budworms (outbreak state) and the low plateau denoted few budworms (refuge state). Note that in their model wave-pinning does not occur. Nevertheless, they compared a deterministic version (thermodynamic limit) of their model with a stochastic version of it and observed differences in wave propagation velocity. They explained that these differences are caused by random jumps, possible within the stochastic system, from the high plateau to the low plateau and vice versa, similar to our explanation for the fluctuations observed in the pinning position. In their case, however, the stochasticity affects the velocity of the travelling wave, while in our study it affects the pinning position.
In future efforts, it would be interesting to study the stochastic model for cell polarisation in higher dimensions, where effects of geometry are non-trivial (e.g. see Strychalski et al. 2010), as well as in off-lattice Brownian dynamics models.
Finally, studies such as this one can be extended to tackle the intriguing question of how cells communicate polarisation within a wider tissue context, which is a subject of ongoing work. We envision that stochasticity within the coupling of cell polarities between cells could play an important role.