Abstract
Many mathematical models have been proposed for the process of cell polarization. Some of these are ‘functional models’ that capture a class of dynamical behaviour, whereas others are derived from features of signalling molecules. Some mechanistic models are detailed, and therefore complex, whereas others are simplified. Each type contributes to our understanding of cell polarization. However, the huge variety at different levels of detail makes comparisons challenging. Here, we provide examples of both elementary and more detailed models for polarization. We also display how a recent mathematical method, local perturbation analysis, can provide an appropriate tool for such comparisons. This technique simplifies and speeds up the model development process by revealing the effect of model extensions, parameter variations and in silico manipulations such as knockout or overexpression of key molecules. Finally, simulations in both one dimension and two dimensions, and particularly in deforming twodimensional ‘cells’, can highlight behaviour not captured by traditional simulation methods.
1. Introduction
The topic of eukaryotic cell polarization has attracted the attention of many modellers, resulting in a large and growing literature. While each group has its own motivation and paradigms, some questions, on which we focus here, are of universal relevance. These include:

— How can we hope to understand the complex networks involved in cell polarization using quantitative (modelling) approaches?

— Can we be confident that the simplified models capture some properties well?

— What is the relationship between detailed and strippeddown models? Do we have to use ‘functional’ models to capture the phenomena, or can we preserve some biological properties?

— What kinds of tools can we use to study and understand models at various levels of complexity?

— Are there reasonable methods that we can use to compare distinct models at various levels of detail?

— What can we learn from the much more computationally intensive simulations of twodimensional deforming cells that may not be evident in onedimensional simulations?
All types of models have their uses and their limitations. Ultimately, we would like to be in a position to faithfully account for all the components in the complex signalling networks in the cell, and how they work together in space and time. Practically speaking, this is still a challenge for the future, as much remains to be determined experimentally before such detailed models can be constructed.
Models for cell polarity are traditionally sets of reaction–diffusion (RD) equations. These partial differential equation (PDE) models are very sensitive to changes in kinetic terms, and no universal method for classifying or categorizing their behaviour exists. This means that fully analytic comparison of distinct models is a difficult process. Simulations provide a simpler (but less satisfying) option: we can simulate various proposed models using similar protocols to see how they compare. (An example of such onedimensional protocols was given by Jilkine & EdelsteinKeshet [1].) One issue is that each model has a variety of possible behaviours in various parameter regimes, and comparing across all models and regimes becomes combinatorially difficult. For such systems, specifically those with slow and fast diffusing components, a recent tool of local perturbation analysis (LPA) proves very helpful.
Which models are more useful, those that aim to be allencompassing and detailed, or those that are simple and mathematically tractable? Are both types relevant to biology? Here, we argue that it is possible to preserve some biological realism in strippeddown abstract models and still refer to specific biological components. An example from our work is the activity cycling of GTPases, reviewed briefly below. Chronologically, the detailed versions of such models led to simplifications that could be analysed mathematically. In turn, the simple models were ideal for understanding the underlying mechanism at play in the tangle of interactions of the larger networks. How can we be sure that simple and detailed versions of the models preserve some common essence? Aside from simulations, it is desirable to consider the dynamic regimes that models can display as certain key parameters are varied. This type of ‘signature’ provides a summary of model dynamics that no one set of simulations adequately encompasses. We illustrate a type of bifurcation analysis that provides such a signature.
2. Models: from detailed to simple
Modelling biological systems requires a balance between detail (which implies complexity) and mathematical tractability (which motivates simplicity). These two extremes are at opposite poles of our ability to understand the implications of a model analytically. An eminent applied mathematician, Lee Segel, taught his students to seek to identify ‘which of many variables are important, and which can be disregarded in simplifying the problem’ (in the recent words of his son, Joel Segel).
Understanding the regulatory machinery responsible for cell polarity is a prime example where this approach has been useful. While this machinery is complex and multifaceted, in many cases conserved elements and motifs emerge, as shown by the examples in figure 1. Mathematical modelling is well suited for identifying these motifs and deciphering their functions. One example is the central regulatory role of small GTPases such as Cdc42, Rac and Rho. Subsets of these regulators are ubiquitously involved in control of actin nucleation and growth, and myosin contraction in most eukaryotic cells. While their specific interactions can vary across cell types, their core function is conserved and central to eukaryotic cell motility and (in some cases) cell polarization.
(a) Lessons we learned from small GTPase models
Models for networks of small GTPases interacting with one another and with other components have been described [5–8]. At present, there is virtually no kinetic or biochemical detail that allows each binding and regulation step to be modelled mechanistically. Moreover, attempting to do so can lead to highly complex models for small parts of the system, which are not ideal. Most such models are at least partly phenomenological. We started with the model for three interacting GTPases, Rac, Cdc42 and Rho, with mutually inhibitory interactions between the latter two and input to the central regulator Cdc42 [6, fig. 3c]. The mutual inhibition resulted in multiple steady states (bistability) with regions of high Cdc42 (front) and low Cdc42 (rear), which were at the time associated with reciprocal (low and high) Rho. However, RD models of these three interacting proteins result in waves of activity that sweep across and take over the whole cell, so no polarization can be sustained. Our first lesson from this disappointment was that not all reasonable verbal models in fact have the expected behaviour when implemented in a spatial setting. While the model displays bistable kinetics in its timedependent version, it fails to account for the existence of multiple regions in a single spatial domain.
A second lesson we learned is that some types of biological details are important. For example, GTPases have both active and inactive forms. Transitions between these forms are on the timescale of seconds, much faster than protein synthesis, and hence the total amount is roughly constant on the timescale of cell polarity. Active GTPases are membranebound, and so diffuse more slowly than their cytosolic inactive counterparts. When this biological detail is incorporated [5], the model becomes a sixvariable RD system with three active and three inactive forms and conservation of each of the proteins overall. While the inactive form does not participate in signalling directly, its availability affects activation: when it is depleted below some level, no further activation can take place. Incorporation of these conservative (in)activation kinetics leads to sustained polarization. Bistability initiates a wave of activity, and depletion of the inactive substrate stops this wave before it traverses the cell. This type of behaviour has been denoted wavepinning (WP).
Another lesson gained from experience is that some modelling detail is important. For example, for simple exchange between active and inactive forms of the same protein, the basic structure of the wellmixed model (ignoring space at present) must be of the following type: 2.1where a, b are rates of activation and inactivation that could depend on u and on other variables. In the specific case of small GTPases, a = a(u,…) is the rate of guanine nucleotide exchange factor (GEF)mediated activation and b = b(u,…) is that of GTPase activating protein (GAP)mediated inactivation. Both might be functions of the level of active GTPase and/or other active components that are known or hypothesized to affect (in)activation. (As v is inactive and does not participate in signalling, neither a nor b should depend on it.) In general, either a or b must be nonlinear to account for switchlike cellular responses. Moreover, it is customary to assume some saturation in the dependence of these rates on components that affect them.
One more observation that emerges from our experience is that some modelling detail is less important at this stage of exploration. Indeed, it has been shown that from a modelling perspective, numerous biochemical regulatory circuits (with positive or negative feedbacks) are functionally the same in the sense of having similar dynamic behaviour [9]. Until details about signalling systems are determined experimentally, it is therefore impossible to distinguish between (say) certain types of positive feedback to GEFs versus a specific reciprocal negative feedback to GAPs in small GTPase models. Furthermore, the switchlike response could be due to zeroorder ultrasensitivity or to cooperativity (Hill functions). In future experiments, it may be possible to test for such details using knockout experiments whereby either the GEF or the GAP is gradually silenced. At present, the choice of where the feedback acts and how it acts is, to a great extent, arbitrary in the construction of preliminary models, as they can be formally shown to be dynamically equivalent [9].
(b) Simplifying the small GTPase models
Mathematical analysis of the phenomenon of WP that we observed in the sixvariable PDE model for Cdc42, Rac and Rho was challenging. In such cases, one is often constrained to exhaustive simulations (but see appendix in [6]). Consequently, it is hard to fully understand how the system works, and we were motivated to simplify and abstract this model. To do so, we formulated a prototype consisting of a single GTPase in active (u) and inactive (v) forms (resulting in a twovariable RD system). There, nonlinear positive feedback from u to its own activation replaced the previous double negative feedback. We preserved the details that had proven to be important: large differences in the rates of diffusion, constant total amount of the protein and nonlinear saturating feedback. We found that the resulting simple singleGTPase system is capable of polarizing via an apparently similar WP behaviour. The strippeddown system was then sufficiently simple to be analysed mathematically [2,10]. Interestingly, the same model was also adopted to describe polarization of PAR proteins in the early embryo of the nematode Caenorhabditis elegans [11,12]. We refer to this model henceforth as the WP model.
The significance of the WP model is as follows: first, it is a simple prototype system to describe robust polarization that has a threshold behaviour (shown below). Several recent cell motility simulations such as (e.g. [13,14]) have adopted this model to enable polarization in prototype motile cells. By construction, the WP model has aspects that specifically describe known proteins (GTPases), so it is not simply a ‘functional’ model in the sense of having appropriate dynamic behaviour. As noted in several works (e.g. [11,12]), it generalizes to signalling proteins that cycle on and off the membrane (for example, PAR proteins). Finally, its simplicity allows for extensive mathematical analysis [10], generalization to a stochastic version [15] and extension to include the effect of Factin [3], described below.
3. Tools for analysis of polarization and other reaction–diffusion models
The vast majority of cell polarity models are based on systems of RD equations (but see Houk et al. [16], where membrane tension plays an important global inhibitory role). A unifying feature of these models is a dichotomy between fast and slow diffusing components. A typical twovariable model can be described in general by the PDEs 3.1We will take u to be slow diffusing (D_{u} < D_{v}) and, wherever relevant, it will represent the active form of a polarity regulator. Classical analysis of such systems, dating back to Turing [17], is based on linear stability analysis (LSA) to check for instability of a homogeneous steady state (HSS) to small amplitude noise. This reduces the PDE problem (equations (3.1)) to a linear algebraic problem of determining whether eigenvalues of the linearized system can have positive real parts. When this is the case, the HSS can be destabilized and some type of pattern (spatially nonuniform solution) may result. LSA gives clues about pattern initiation and (in some cases) of the types of pattern to expect, but it cannot predict the full evolution of the system because nonlinear effects generally come into play once the perturbation grows. While LSA is an important tool in the applied mathematician toolbox, it has limitations. First, it can only check for instabilities that are provoked by arbitrarily small noise. Second, it does not easily scale up to larger nonlinear networks, where finding the HSS and eigenvalues can be extremely challenging.
For this reason, complementary tools have proved useful in our work on polarization. We review one such method, LPA, in the next section.
(a) Brief survey of local perturbation analysis
The technique described below was invented jointly by AFM Marée and V. Grieneisen. It was initially described in the appendix of Veronica Grieneisen's PhD [18], Utrecht University. It was then used extensively to parametrize models in Marée et al. [5,8] and more recently in the context of Holmes et al. [3,4]. Informal justification of this method will be presented here with detailed mathematical justification pending.
Before describing LPA, we note for future reference the typical form of kinetic equations describing the wellmixed (spatially uniform) system: 3.2We will refer to this wellmixed system and to the full RD system (equations (3.1)).
The basic idea of the LPA is to take advantage of the diffusion discrepancy and consider the limit D_{u} → 0, D_{v} → ∞ in the RD system (equations (3.1)). We probe stability of the HSS with respect to a ‘local’ perturbation in the form of a narrow peak of u at some location in the domain. (The magnitude of the peak need not be small, but its ‘total mass’ should be negligible.) In this limit, the localized peak of u (denoted u_{L}) does not influence the surrounding background levels of u (denoted u_{G}). As v diffuses ‘infinitely fast’, it is uniform in space and can be described by a single global variable (v_{G}). Further, since the mass of the perturbation is negligible, u_{L} does not influence the evolution of v_{G}. So, u_{G} and v_{G} interact according to the wellmixed chemical reaction kinetics and the growth or decay of u_{L} is influenced by v_{G}.
The dynamics of this peak and the global variables describing the broader domain can then be described by a collection of evolution ordinary differential equations (ODEs) for u_{L}, u_{G} and v_{G}. 3.3This system of ODEs describes the growth or decay of the stimulus peak, providing stability information for the underlying spatial system.
We now apply this method to the polarization models described above. A common feature of polarization models is interconversion between two forms with the total average concentration w conserved. We show two such examples below. In this case, g(u,v) = −f(u,v). We use conservation, to eliminate one variable (v_{G} in this case), leaving the reduced LPA system 3.4The original problem consisting of PDEs has been reduced to a set of ODEs describing the growth of a narrow pulse. Dynamics of these ODEs under a wide range of parameter variation can be explored readily using bifurcation software.
Figure 2 depicts the results of such an analysis. For the given parameter set, the wellmixed WP model has a unique stable steadystate value of u for all values of the parameter k_{0}, as shown in figure 2a. As k_{0} increases, the steadystate level of the active form u increases monotonically. We refer to this curve as the ‘global branch’ in what follows. As shown in figure 2b, the LPA diagram includes both the global branch and an additional (grey) curve hereafter denoted the ‘local branch’. The interpretation of this LPA diagram will be explained in the next section, where a comparison with a second model is made.
LPA has several benefits. It is fairly easy to formulate the LPA ODEs by identifying the slow and fast variables. It is then easy to explore the parameter regimes of growth/decay and/or threshold behaviour of the perturbation using bifurcation software. Below, we further illustrate several advantages of LPA. (i) It provides a standardized tool for comparison of distinct RD models for polarization. (Models with similar mechanisms have a similar LPA ‘signature’. Models that look superficially similar but have distinct mechanisms have different LPA signatures.) (ii) LPA helps to point to what changes when a model is extended or modified in some way. For example, if new feedback in introduced or another component is added to a model, then LPA can show how the behaviour changes over a whole range of parameter variation. (iii) LPA is relatively easy to scale up to larger systems of interacting components. This makes LPA much more ‘userfriendly’ than LSA.
4. Model comparisons
As previously mentioned, the literature on models for cell polarization is large, with a variety of models of many types [19–21]. How should we compare such models and their functionality? In principle, we can do so by simulating the response of models to similar protocols [1], but this can at best probe a limited set of parameters in each case. Here, we show that LPA can help to inform such comparisons.
(a) Example 1: comparison of wavepinning and Otsuji models
Here, we use LPA to compare the polarization model of Otsuji et al. [19] (abbreviated OT) to the WP model of Mori et al. [2]. Both OT and WP share the form of equations (3.1) with g(u,v) = −f(u,v), implying conservation of the total amount of protein. In both models, u and v represent active and inactive forms of GTPases, respectively. The kinetics for these models are given by, for the WP model [2] 4.1
and for the OT model [19] 4.2
In the kinetics of equation (4.1), there is a rate of activation of u from v (basal value k_{0}) and an inactivation rate δ. (This follows the form of equations (2.1).) A Hill function depicts positive feedback from u to its own activation. The total average concentration, w, is determined by initial conditions, not by other parameters. While OT kinetics (equation (4.2)) can be rewritten similarly in the form of equations (2.1), both activation and inactivation rates depend explicitly on the total average amount of protein w.
Figure 3 compares the LPA ‘signatures’ of WP and OT. Here, the steadystate value of u_{L} is plotted against the bifurcation parameter w. In each case, there is a monotonically increasing ‘global branch’ (thick dark curve), together with an additional ‘local branch’ that stems from spatial features captured by the LPA equations (thinner grey curve). The OT model in figure 3a has a strictly decreasing local branch. For small values of w, only the global branch is stable. The intersection point is a transcritical (branch point) bifurcation at which the two branches exchange stability. Thereafter, for larger w values, we have a Turing regime, whereby a perturbation of arbitrarily small size can lead to patterning. To the left of that bifurcation, LPA suggests that only perturbations of sufficiently large amplitude might lead to destabilization. (In practice, because LPA is an asymptotic approximation, and as the local branch climbs very steeply upwards towards the left, it is hard to cause patterning in that regime.)
Figure 3b for the WP model is significantly different. First, we see four rather than two major regions (I–IV, left to right, see legend of figure 3). In Regions I and IV, a stable HSS can never be perturbed. In Region II, the HSS is stable (solid curve) but can be destabilized by a sufficiently large local perturbation (past the threshold of the local branch, dot–dashed). For values of w between the two intersections (transcritical bifurcations), the HSS is unstable to any perturbation (Turing regime, Region III). Between the rightmost intersection and the ‘knee’ to its right, the system can only be prodded away from its elevated HSS by a pulse that locally lowers the value of u by a sufficient amount (Region IIb). It turns out that both Regions II and IIb correspond to wavepinning behaviour, a fact that can only be determined by fullscale simulations. Examples of full PDE simulations in one spatial dimension are shown in figure 3c for OT and figure 3d for WP. In general, OT tends to form steep peaks, whereas WP has plateaus with broad shoulders.
In short, WP has distinct regimes of (i) complete stability, (ii) Turing instability and (iii) two distinct regimes of threshold response. OT shares the Turing regime, and a threshold response that in practice proves hard to trigger. Overall, aside from providing insight into parameter regimes for a given model, LPA produces a graphical summary of the range of distinct dynamic properties, providing a ‘signature’ that can be compared across RD polarization models.
(b) Example 2: the effect of model extension and additional feedback
The WP model has been studied earlier in Mori et al. [2,10] and by LPA and related methods in Walther et al. [15]. Here, we show that LPA can be used to understand the effects of model extension or modification. In Holmes et al. [3], we considered an extension of the WP model that included interaction with Factin (see figure 1b). This was motivated by the following observations: (i) waves of Factin and nucleation promoting factor activity have been observed in cells [22–24] and (ii) waves of small GTPases and edge protrusion (presumably powered by actin polymerization) have also been observed [25]. Accordingly, a modified WP model was linked to nucleation of Factin. Feedback from Factin was assumed to accelerate inactivation of the small GTPase (figure 1b), that is, to act as negative feedback channelled through GAPs. The resulting WPF model is among the simplest descriptions of Factin waves driven by nucleation promoting factors.
To account for the additional influence of Factin (F) on the system, the following minimal extension of the WP model is used 4.3where the GTPase kinetics includes an Factindependent inactivation term 4.4The model was explored by simulations guided by LPA analysis [3]. Figure 4 compares the bifurcation diagrams for the WP model alone (where F is treated as a constant parameter, i.e. ε = 0) and the fully dynamic A, I, F model. We see the appearance of several new Hopf bifurcations that were absent in WP, arising as a result of the dynamic feedback from Factin. These regimes lead to a variety of wavelike behaviours [3], including static polarization, oscillatory polarization, travelling wave trains, reflecting waves and other exotic patterns as shown in figure 5. Note that such wavelike patterns are absent in the WP model alone.
5. Comparing simple and detailed models for polarization
Ultimately, we are interested in learning about real signalling networks and their dynamic behaviours. To do so, we must return to the more complex biological scenarios, where GTPases interact with each other and with other signalling modules. A benefit of this approach is that the model components are linked to known molecules which could, in principle, be manipulated in experimental tests.
Examples of this type of modelling include Dawes & EdelsteinKeshet [7], Marée et al. [8], Kholodenko [9] and Marée et al. [26]. In Kholodenko [9] and Marée et al. [26], models for active and inactive Cdc42, Rac and Rho were simulated spatially in one dimension and two dimensions, in both static and deforming ‘cells’. In Dawes & EdelsteinKeshet [7], Marée et al. [8], GTPases were linked to phosphoinositides PIP, PIP_{2} and PIP_{3} via positive feedback to/from Rac and Cdc42. (The role of PIPs has been controversial in the biological literature, and we sought to explore how they affect the small GTPase layer.) These models were based on mutual inhibition between Cdc42 and Rho supplemented by positive feedback from Cdc42 to Rac and from Rac to Rho.
In recent experimental work described by Lin et al. [27], HeLa cells cultured in microfluidic chambers were highly constrained to polarize along the direction of narrow channels. Biochemical constructs in these cells were engineered to convert a controlled external gradient of a small molecule (rapamycin) into a robust internal gradient of Rac, bypassing all the innate receptorbased signalling systems. Cells in several states were tested for polarization responses, and results were quantified and compared to a sequence of models for internal signalling. We noted immediately that the model based on mutual antagonism between Cdc42 and Rho could not account for the Racmediated polarization response: in fact, that model polarizes in a direction opposite to the gradient of Rac activation. Rather than simply add supplemental interactions to an existing set of proposed interactions, we instead sought the simplest connectivity between these GTPases capable of recapitulating this data. Figure 1c shows the schematic diagram of this ‘minimal circuit’. See Holmes et al. [4] for details of the models and analysis.
The level of complexity of the detailed network in figure 1c is significantly greater than that of the WP model in figure 1a: the model consists of nine PDEs (for C, R, ρ active/inactive forms, and P_{1}, P_{2}, P_{3} representing the three forms of PIP lipids.) However, carrying out LPA analysis is relatively straightforward, because the resulting 15 LPAODEs (PIP_{1,2,3} are considered slow variables) are simpler to analyse with bifurcation software than most systems of even two to three nonlinear PDEs.
A typical LPA bifurcation plot with respect to the basal activation rate for Rac (I_{r}) leads to the diagram shown in figure 6a. A direct comparison can be made with figure 2b which provides the LPA diagram for WP with respect to a similar basal activation rate parameter k_{0}. We observe a number of similarities: both models share a monotonically increasing global HSS branch (dark curves). Both share similar regimes: a wavepinning polarization regime where the HSS is sensitive to a sufficiently large stimulus peak (II), a Turing noisesensitive regime (III) and a regime where no perturbation induces patterning and only the HSS is stable (IV). It is interesting to note that, even though the models are significantly different, they share certain key features. (This suggests that lurking behind a facade of details is the core mechanism that was identified in the far simpler WP model.) Furthermore, these similarities are easily detected using the LPA.
A second feature of figure 6 illustrates the usefulness of the LPA diagram in understanding the effect of specific experimental manipulation. Figure 6b shows the effect of PI3K knockout, which reduces the feedback between PIP_{3} and Rac (marked with f_{1} in figure 1c). The series of increasingly thick curves represent increasing levels of feedback (from none at f_{1} = 0 to essential at f_{1} = 1). As feedback increases, the grey loop (‘local branch’) shifts to the left. This means that at a given value of the activation rate k_{0}, the cell is more sensitive to stimulation: a pulse of smaller amplitude is able to breach the threshold and lead to polarization. In this sense, LPA also illustrates how feedbacks of various sorts affect cell behaviour, with minimal computational effort.
6. Cell shape and model comparisons
The ultimate test of polarization models is how they perform in a deforming cell, and how they are linked to cell motility. In real cells, assembly of the cytoskeletal polymer actin prods the cell membrane forward at the leading edge, while contraction of the cytoskeleton by myosin molecular motors pulls up the rear. Anchorage by adhesion molecules provides requisite traction. Cell polarity is an additional layer of chemistry that keeps cell migration oriented.
The simplest possible model for a deforming cell requires only a recipe that couples intracellular chemistry to membrane velocity. To focus on the interaction between cell shape and cell polarity, it makes sense to bypass a detailed description of molecular force generation by coupling the chemical kinetics of cell polarity directly to membrane displacement. This phenomenological approach allows for a direct comparison of different polarity mechanisms, without additional complicating factors.
(a) Level set methods
Investigating the influence of cell geometry on cell polarity requires a scheme for describing arbitrary cell shapes. One possibility would be a Lagrangian description of points along a parametrized cell boundary. Level set methods provide an Eulerian alternative. As a simple example, in two dimensions, consider a circular cell that grows bigger, without changing shape. Level set methods describe this scenario in terms of a cone that points downwards as it descends along a vertical axis. Intersection of this moving cone with a fixed horizontal plane then represents the edge of the growing cell. At any point in the plane, the corresponding cone elevation ψ gives the shortest distance to the cell boundary, with negative elevations inside. This defines a ‘distance map’ for which the zero level set (ψ = 0) provides an implicit description of the cell boundary. In general, once intracellular physics sets the velocity V for membrane displacement, advection of ψ then moves the boundary appropriately (∂ψ/∂t = − V · ∇ψ), using standard (though nontrivial) numerical methods [28,29]. Level set methods make it easy to compute the outward normal and curvature κ of a cell boundary ( and ).
(b) The moving boundary node method
Given a level set description for the edge of a cell, it remains to describe cell polarity within the contained region. The moving boundary node method [30] provides a means of solving the requisite differential equations. In a level set representation, the cell boundary usually falls between the nodes of a computational grid, with the boundary location determined by interpolation. However, nodes near the boundary can be moved onto the boundary without difficulty; taken together, ψ and give exactly the required displacements. This yields dynamic, boundaryfitted coordinates and a set of distorted volume elements that exactly fit the shape of a cell. Computational nodes retain simple Cartesian connectivity, which makes it easy to compute fluxes between adjacent volume elements. This leads to a finitevolume discretization of reaction–advection–diffusion equations inside moving, deforming cells. The method is accurate to second order and conserves mass.
(c) Comparison of wavepinning and Otsuji polarity
For both the Turing regime of OT (figure 3a) and region II of WP (figure 3b), perturbation of an unstable, spatially homogeneous steady state triggers polarization of a static circular cell. Disrupting the homogeneity of active protein u with a 10% spatial gradient across the cell provides sufficient stimulus. Compared with WP (figure 7a), polarization for OT (figure 7f) is more localized and has a peak of greater magnitude. These polarized states, for static cells, provide initial conditions for cells with moving boundaries.
Perhaps the simplest possible assumption is a linear relationship between the concentration of active protein, u, and velocity, v, along outward normal, , at each point on the perimeter of a cell 6.1where V_{0} = 0.5 µm s^{−1} sets a realistic velocity for typical crawling cells, while u_{SS} is the steadystate concentration for a polarized, static, circular cell, and u* is an outward growth threshold, chosen here as the value of u_{SS} at which ∇u_{SS} attains a maximum. Points on the cell boundary where u − u* is positive move outward, along , whereas the cell contracts at points where u − u* is negative. For the same parameter values considered in one dimension (figure 3), WP yields convex, translating cells (figure 7b–d), whereas OT yields constricting, concave cells (figure 7g–i).
Constriction of cells under OT can be understood as dominance of reactions over diffusion. This can be expressed as a ratio of timescales 6.2where τ_{R} is the characteristic time for OT reactions while R = 5 µm is the chosen cell radius, so τ_{D} is the characteristic time for intracellular diffusion. Recall that a_{1} is the reaction rate for OT (equation (4.2)), whereas k_{0}, δ and γ are WP reaction rates (equation (4.1)). Substituting the maximum of k_{0}, δ or γ for a_{1} in τ_{R} gives the characteristic time for WP reactions instead. For either OT (figure 7g–i) or WP (figure 8d,e), a ratio of τ_{D}/τ_{R} = 2500 yields constriction. In this case, relatively fast reactions cause localized depletion of u that diffusion fails to replenish so u drops below u* to give inward velocity at the point of concavity. By contrast, a ratio τ_{D}/τ_{R} = 250 yields convex translating cells for both WP (figure 7a–d) and OT (figure 8d,e).
Notably, constriction under WP (figure 8d,e) highlights some limitations of working in one dimension. In two dimensions, the cell develops a narrow bottleneck (figure 8e, inset) between two compartments. The width of the bottleneck and the average radius of each compartment give a constricted cell three inherent length scales, compared to the single length scale for polarity simulations in one dimension (figure 3). This is significant because the tradeoff between reaction and diffusion timescales (equation (6.2)) depends on the characteristic length of the region in question. These geometrical effects somehow combine to manifest as emergence of two separate regions in which u is elevated (figure 8e), starting from just one region (figure 7a). Over time, the concentration profile for constriction under WP develops a local maximum, towards the rear of the cell (figure 8c, grey curves). This does not happen in one dimension (figure 3d), for a comparable initial condition.
7. Discussion
In this paper, the main emphasis has been on GTPasebased models for polarity initiation in motile eukaryotic cells. The models surveyed here consist of sets of RD equations. We have shown several methods for comparing across distinct models (from simple abstract ones to more detailed and complex ones). One method that is of particular utility is the LPA, a useful tool for systems with slow and fast rates of diffusion. We have shown that seemingly related models for polarization in fact have very different LPA ‘signatures’ (as shown in figure 3 for two types of small GTPase models), whereas models at different levels of detail can be related by a common or similar LPA signature (e.g. figures 2b and 6a.)
LPA can complement other methods of analysis but is not a replacement for fullscale simulations, such as those of figures 3c,d and 5. Those are still needed to explore the patterns that form after the shorttime initiation phase. A first step is to determine the behaviour of such models in onedimensional spatial domains where we can ask whether robust polarization results from various stimuli.
However, as shown in §6, even fullscale onedimensional simulations are not always predictive of the evolution of polarization in a deforming twodimensional ‘model cell’. While the latter defines a challenging moving boundary problem, such simulations are an ultimate test of models for polarization and cell motility.
Funding statement
The authors acknowledge support from NSERC discovery grant to L.E.K. and NSERC postgraduate fellowship to M.D. L.E.K. and W.R.H. have also been partly supported by National Institutes of Health under grant no. R01 GM086882 to Anders Carlsson.
Acknowledgements
L.E.K. organized review, wrote §§1–5; W.R.H. developed LPA, edited paper and provided figures; M.Z. created, simulated, analysed twodimensional deforming cell simulations and wrote §6 and M.D. provided figures 2,3a,b and 5.
Footnotes
One contribution of 17 to a Discussion Meeting Issue ‘Cellular polarity: from mechanisms to disease’.
 © 2013 The Author(s) Published by the Royal Society. All rights reserved.