Order and Stochastic Dynamics in Drosophila Planar Cell Polarity

Cells in the wing blade of Drosophila melanogaster exhibit an in-plane polarization causing distal orientation of hairs. Establishment of the Planar Cell Polarity (PCP) involves intercellular interactions as well as a global orienting signal. Many of the genetic and molecular components underlying this process have been experimentally identified and a recently advanced system-level model has suggested that the observed mutant phenotypes can be understood in terms of intercellular interactions involving asymmetric localization of membrane bound proteins. Among key open questions in understanding the emergence of ordered polarization is the effect of stochasticity and the role of the global orienting signal. These issues relate closely to our understanding of ferromagnetism in physical systems. Here we pursue this analogy to understand the emergence of PCP order. To this end we develop a semi-phenomenological representation of the underlying molecular processes and define a “phase diagram” of the model which provides a global view of the dependence of the phenotype on parameters. We show that the dynamics of PCP has two regimes: rapid growth in the amplitude of local polarization followed by a slower process of alignment which progresses from small to large scales. We discuss the response of the tissue to various types of orienting signals and show that global PCP order can be achieved with a weak orienting signal provided that it acts during the early phase of the process. Finally we define and discuss some of the experimental predictions of the model.


Introduction
Epithelia in diverse tissues, in addition to their apico-basal polarization, acquire a polarization within the two-dimensional layer of cells -a phenomenon called planar cell polarity (PCP) [1][2][3][4][5]. In the developing wing of Drosophila, PCP determines the growth direction of small hairs that extend radially from cell boundaries. In a wild-type wing, where cells are approximately hexagonal and form a regular honeycomb lattice, all of these hairs point to the distal direction.
A series of recent experiments show that several key proteins [6], including the transmembrane proteins Frizzled (Fz) and Van-Gogh (Vang) and the cytosolic proteins Dishevelled (Dsh) and Prickled (Pk), localize asymmetrically on cell boundaries [7][8][9][10][11][12] defining a direction in the plane within each cell and forming a characteristic zig-zag pattern of protein localization on the lattice (Fig. 1A).
Other experiments show that local PCP orientation depends on inter-cellular signaling. First, mutant clones in which fz or Vang activity is suppressed or amplified, cause characteristic and reproducible inversion of polarity in large patches of cells that are proximal or distal to the clone [13]. These observations are summarized in Figs. 1 C,D. Second, in fat mutant clones [14,15] hairs do not all point correctly in the distal direction, yet, their orientation is strongly correlated between nearby cells and varies gradually across the tissue creating a characteristic swirling pattern.
Thus the experimental evidence suggests that an interaction between neighboring cells tends to locally align their polarity [1,3,14]. This local polarity need not point distally unless, in addition, there is a global orienting signal that picks out the distal direction throughout the wing (most likely originating with the Dpp morphogen gradient which defines the Anterior-Posterior axis of the wing in the larval stage of development [16]). Yet, aside from a clear involvement of protocadherin fat [17,18] the molecular details of this pathway remains for now unknown. The swirling patterns in fat mutants [14] and recent evidence [15,19], suggest that the orienting field is related to the presence of a ''gradient'' in the fat, four-jointed, and dachs pathway.
These observations evoke an analogy between PCP and the behavior of ferromagnets, extensively studied in physics and well understood in terms of statistical mechanics of relatively simple models [20]. In these models each atomic site is assigned a magnetic dipole -spin -which can assume a different orientation (analogous to the direction of polarization in an epithelial cell). The salient properties of ferromagnets arise from the opposing influence of an interaction between neighboring spins, which tends to co-align their orientation, and the influence of thermal fluctuations, which tend to randomize the spin direction. Ferromagnets typically exhibit two phases of behavior: a high temperature phase, where spins are disordered and a low temperature ferromagnetic phase, where the interactions dominate over thermal fluctuations -leading to a spontaneous polarization in an arbitrary direction. In this state even a small external magnetic field has a big effect on magnetic polarization as the spontaneous polarization aligns itself with the external field, yet the dynamics leading to global alignment can be quite slow.
An essential lesson from statistical mechanics is that the ordered and disordered states exist in a broad class of models and can be discussed in a general context, focusing on a classification of the different regimes as a function of a few parameters. We follow this lesson by focusing the study on the competition between the intercellular interaction and the disordering influence of the fluctuations introduced by the noisy molecular interactions. As in statistical mechanics we define a phase diagram which identifies different regimes of behavior in the space of the most relevant parameters. We then address the role of the global directional signal in the dynamics of global alignment.
A molecular model for PCP formation was recently proposed in Ref. [21], and was shown to reproduce a number of experimental findings. This model involves 38 parameters that were adjusted to successfully reproduce a set of wild-type and mutant phenotypes.
Here we pursue an alternative approach and instead of moving on to more and more complex models develop a model with a smaller number of degrees of freedom and a smaller number of parameters. Instead of fixing a particular set of parameters by fitting the data we explore the generic behavior of the model as a function of parameters defining quantitative features characteristic of the different phases. In formulating the model we identify several essential ingredients, required to obtain the characteristic zig-zag pattern and the non-autonomy of fz and Vang mutant clones. We expect our simplified model to capture important properties of PCP, although it does not incorporate all the molecular details.
After discussing the essential ingredients of the model, we obtain a phase diagram describing its steady state properties. We then consider the dynamics of local polarization strength and orientation in the absence and in the presence of a global orienting signal. We show that global alignment can be achieved with a weak global orienting signal provided it is present throughout the tissue at the earliest stage of PCP dynamics. Finally we discuss the experimental predictions coming out of the model and the tools required to test these predictions.

Model ingredients
Three essential ingredients are included in the model, to account for the characteristic zig-zag patterns of protein localization and for the non-autonomy of fz and Vang mutant clones.
Two membrane proteins form complexes across the inter-cellular interface. As in Ref. [21] we assume that two membrane-bound proteins, a and b -standing for Fz and Vangform complexes across inter-cellular interfaces. This is the source of intercellular interaction in the model.
Complex formation across cell interfaces accounts in a simple way for the non-autonomous effect of clones in which either a or b are mutated. However to account for the observed localization of Fz and Vang proteins on the opposite sides of the cell interface there must be a mechanism which prevents a, b (or Fz and Vang) from mingling with each other on the same side of the interface. Thus the next two assumptions introduce molecular interactions acting inside each cell, leading to spontaneous segregation of the complexes and driving the protein distribution towards a nonuniform state.
Complex formation on a single inter-cellular interface is bistable. We assume that complexes of one polarization (a=b) Figure 1. Summary of experimental observations. (A) Protein localization pattern in wild-type wing: Fz (green) localizes on the distal membrane, together with Dsh, while Vang (red) localizes on the proximal membrane, together with Pk. (B) Key PCP proteins localize apically in the adherens junction area, within a strip of about 1m from the top [7,11,12,39]. (C,D) Mutant fz (C) and Vang (D) clones influence the polarity of wild-type cells bordering the clone such that it points towards the clone (fz, C) or away from it (Vang, D). This effect is propagated to a large patch of wild-type cells that are distal to the clone (fz) or proximal to it (Vang) [40]. Over-expression of fz causes an effect similar to that of Vang mutant clones, and over-expression of Vang causes an effect similar to fz mutants. doi:10.1371/journal.pcbi.1000628.g001

Author Summary
Epithelial tissues are often polarized in a preferred direction which determines, for example, the direction of hair growth on mammalian skin, the orientation of scales in fish, the alignment of ommatidia in the fly eye and of sensory hair cells in the vertebrate cochlea. This in-plane polarization, known as planar cell polarity, is one of the morphogenetic fields that play a role in tissue patterning during development. Here we focus on planar cell polarity in the fly wing, where protein localization and inter-cellular ligand-receptor interactions combine with an unknown orienting signal to establish planar cell polarity of the wing epithelium. We demonstrate an analogy between this process and models of ferromagnetism in physical systems that have been studied extensively using the tools of statistical mechanics. The analogy helps in understanding how local interactions between cells can lead to global polarization order and elucidate the role of global orienting signals and the dependence of the dynamics of the process on parameters. We demonstrate that in the absence of an external orienting signal swirling patterns should emerge due to random noise. We propose ways to test this prediction and ways to quantify the magnitude and spatial variation of the unknown external orienting signal.
inhibit formation of complexes of the opposite polarization (b=a), and that this inhibition leads to bistability at the level of a single interface between two cells ( Fig. 2A). As a simple example, consider the following dynamics of complex binding and unbinding on a single, planar interface (Fig. 2B), where u 1 and u 2 represent concentrations of interfacial complexes with two possible polarizations (respectively a 1 b 2 and a 2 b 1 ) and a 1,2 , b 1,2 are concentrations of free (unbound) proteins on the two sides of the interface. The positive and negative feedback on complex formation is represented through the dependence of the rate coefficients K on u 1,2 (see Methods). E.g., enhancement of K(u 2 ; u 1 ) in Eq. (1) with increasing u 1 or suppression with increasing u 2 . If this dependence is sufficiently non-linear, the dynamics lead to two stable steady states: one with u 1 wu 2 , the other with u 1 vu 2 , as illustrated in Fig. 2C. Feedback effects could be equally well modeled by an opposite modulation of K' and in reality quite likely involve modulation of both K and K'. As an example, consider the case where there is only negative feedback through the dependence of K or K' on u 2 . If the free a and b diffuse sufficiently rapidly, a 1,2~at {u 1,2 and b 1,2~bt {u 2,1 , where a t and b t are the total available concentrations of a and b proteins. It is then easy to see that for bistability K' or 1=K must be a convex increasing function of u 2 .
Inhibition acts non-locally within each cell. While bistability of the complex formation would suffice to explain localization of Fz and Vang on the opposite sides of each interface between cells, in order to explain segregation of Fz and Vang to the opposite sides of each cell we assume that the mutual inhibition of u 1 and u 2 complexes acts non-locally within a cell. Hence instead of making K in Eqs. (1-2) be a local function of u 1 and u 2 we assume that K is a function of c Ã , the concentration of a messenger molecule which is itself a non-local function of the u 1 and u 2 distribution over the surface of a given cell. The messenger molecule thus mediates an interaction between u 1 and u 2 complexes, i.e., c Ã diffuses within each cell creating an effective repulsion between u 1 and u 2 complexes on adjacent interfaces. The non-local repulsion will for a broad range of parameters result in a dipole-like distribution of a and b (and hence of the u 1 and u 2 complexes) over the surface of each cell.
A plausible and quite general mechanism for generating such a non-local inhibitory signal involves the modification of a diffusible protein, as illustrated in Figs. 2 D,E where we denote the unmodified and modified protein by c and c Ã , respectively. The rate of modification c?c Ã at a given point on the membrane depends on the local density of u 2 complexes and information is transmitted within the cell by diffusion of the modified protein.
Many variations on this general theme are possible and are discussed in detail in the supporting analysis (Text S1, Part I). Below we follow the scheme shown in Fig. 2E, where the membrane-bound protein a serves the role of the messenger protein c. Modification of a corresponds to the binding of a cytoplasmic protein, and this process is inhibited by u 2 complexes. The fraction x of unmodified proteins then obeys the equation where r is the location on the membrane. The parameters k 2 and a, related to the rate constants for modification of a are discussed in the supporting analysis (Text S1, part I). Note that increase of u 2 increases x and that the influence of u 2 is non-local, with a characteristic range set by k {1 . Finally, we assume that only modified a proteins can form complexes, hence the rate coefficient K is proportional to 1{x (Eq. 7). Other details of the non-local inhibition mechanism are described in Methods. Interestingly, we find that to maintain a non-local signal in a steady state an energy flux is necessary (Text S1, part I).

Stochastic dynamics
There are several reasons why the dynamic equations are not deterministic. Even in the steady state, interfacial complexes not only bind and unbind due to thermal fluctuations, but like nearly everything else inside the cell are being constantly recycled and reassembled. Stochastic fluctuations arise from the molecular noise of reactions and the variability in the state of the cell defining the ''intrinsic'' and ''extrinsic'' noise [22]. It will suffice however to describe stochasticity of complex binding and unbinding as if it were a Poisson process. Equation (1) is thus replaced by a stochastic equation, [and a similar modification applies to Eq. (2)] where the noise j can be approximated as white Gaussian noise if the number of molecules per cell is not too small. Assuming that the dominant contribution comes from the finite number of molecules participating in the binding/unbinding dynamics, the variance of j is inversely proportional to N 0 (see Methods), where N 0 is defined as the number of a molecules per interface: N 0~at A where a t is the total concentration of a molecules (bound and unbound) and A is the area of an interface (about 5m|1m -see Fig. 1B). Since the variance of j decreases with increase of N 0 , 1=N 0 plays a role similar to temperature in a ferromagnet. If there are *10 3 Fz molecules per cell [23], N 0 is of order *10 3 resulting in the root-mean-square fluctuations of the order of 3% (i.e. 1= ffiffiffiffiffiffi N 0 p ) of the mean. Other sources of intrinsic noise, in addition to the stochasticity of binding and unbinding events, may increase the noise variance beyond the above estimate. These additional noise sources include, for example, stochasticity in the signaling pathway that generates the non-local inhibition within each cell, or fluctuations in a t and b t . Such sources of intrinsic noise, acting upstream of u 1 and u 2 , are propagated to the PCP signaling dynamics through the dynamics of complex formation, and can thus be described qualitatively by the noise term in Eq. (4), with an effective value of N 0 that is possibly smaller than predicted from the number of a and b molecules alone.

Phase diagram
What are the consequences of the model defined above when cells are arranged on a hexagonal lattice? Let us first consider the steady state in the deterministic limit. Fig. 3A shows a typical phase diagram on a two-dimensional plane dissecting our five dimensional parameter space (see Methods): the y axis is the range of the non-local interaction in units of the cell lattice spacing, and the x axis the coefficient a which controls inhibition (see Methods).
In the region labeled U there is a unique steady state in which there is no polarization of the protein distribution. In contrast, in region S the stable steady state has the symmetry shown in Fig. 3B: Both a and b distributions carry a vector dipole moment that points towards the center of a side, and due to the lattice symmetry there are six equivalent states of this type. A uniform steady state exists as well, but it is unstable. Region V differs from S in the direction of the dipole, which points towards a vertex instead of pointing towards and edge (Fig. 3C).
The transition from the uniform state, U, to the edge state, S in the phase diagram is continuous: the dipole moment tends to zero when approaching the phase boundary from the S side. A similar transition from a U state to a V state can exist as well, and is present on another two dimensional ''slice'' through the parameter space of our model. This transition is also continuous.
We next consider the effects of stochasticity, which were ignored in the discussion above by setting N 0 ??. When N 0 is finite (similar to a non-vanishing temperature in a spin model), we ask whether the steady state maintains long-range order: i.e. whether a particular orientation is singled out throughout the lattice and the dipole moment has a non-zero average. In the language of the analogy with magnetic systems this would be a ferromagnetic state. The latter disappears as the temperature increases above a certain critical value, giving way to a paramagnetic state where dipole moments point in random directions and the average polarization vanishes (an intermediate state with quasi-long range order may exist as well, in similarity to 2-dimensional clock models [24][25][26][27]). Hence, we expect an ordered state to be stable only when N 0 is sufficiently large, and this is indeed observed in our simulations (Fig. 4). Yet with a realistic number of molecules per cell, in the order of several thousands, the vertex and side states in our model are typically ferromagnetic.
It may thus appear that when N 0 takes realistic values the system is in an ordered state and stochasticity is altogether unimportant. However, as we discuss next, the steady state is not necessarily reached within the time scales of wing development, and stochasticity plays an important role in the dynamics of ordering.

Dynamics of ordering in the absence of a global orienting signal
Let us consider the dynamics of PCP formation, first in the absence of a global orienting signal. Fig. 5 shows results from a stochastic simulation, starting from a state where a and b are uniformly distributed in all cells.
We can identify two stages of the process. The first stage corresponds to a gradual build up of a dipolar polarization on the cellular level. The dipole initially points in a random direction, but as its amplitude increases with time ( Fig. 5B) local polarization begins to re-orient. At the end of this stage, when amplitude saturates, there is no global choice of PCP direction, but the orientation of nearby cells is strongly correlated: as an example, Fig. 5A shows the configuration of dipoles shortly after saturation. The second stage, which follows amplitude saturation, exhibits slow coarsening dynamics [27]: polarity direction is approximately aligned within discrete domains, the size of which gradually expands by movement of their boundaries. Note also the existence of vortex-like defects [28] (Fig. 5A and Fig. S1). Coarsening ultimately leads to a spatially uniform steady state, but this process occurs over a long time scale compared to that of amplitude growth.
A quantitative theory of the early dynamics is obtained from the linear instability of the uniform steady state (described in detail in Text S1, part II). The variance of the local dipole amplitude increases exponentially in time with a characteristic time scale l, where for simplicity numeric prefactors of order unity are omitted (see Text S1, part II). In this equation j 0 is the amplitude of noise in the unstable uniform steady state, and both l and L are found from the instability analysis (Text S1, part II). This prediction is shown in Fig. 5B (dashed line) for comparison with the simulation. Two additional insights come from the analysis of early dynamics (Text S1, part II). First, PCP is initially isotropic, despite the discrete 6-fold symmetry of the hexagonal cell lattice. Consequently, the dipole moment initially has equal probability to point in any direction in the interval ½0,2p). Second, the spatial correlation established during the early dynamics typically has a longer range in the direction parallel to the dipole, compared to the perpendicular direction. These two properties of the dynamics lead to a characteristic swirling pattern before non-linearities set in. The range of correlation at this stage depends on the location in the phase diagram and increases logarithmically as a function of N 0 .

Effect of global orienting signals
We next consider how various types of symmetry-breaking orienting signals influence PCP dynamics.
Boundary orienting signal. For example, a row of cells that do not express b (or, alternatively, a) can serve as a boundary orienting signal. Can such a signal orient a tissue as large as the wing? In a deterministic model without any stochasticity, the boundary is the only cause for symmetry breaking, and will necessarily set polarity orientation throughout the tissue. In the presence of noise a local choice of polarity is established in the bulk of the wing, and competes with the boundary signal.
During amplitude growth a moving front separates two regions of the tissue: between the boundary and the front all cells point in the orientation set by the boundary, whereas beyond the front cells point in all possible directions. The position x of the front increases sub-linearly with time, x!t 1=2 , with a prefactor that depends on the position in the phase diagram and increases logarithmically with increase of N 0 (Fig. 6 and Text S1, part II). After amplitude saturation, when an independent choice of polarity is established in the bulk, front propagation is arrested (more precisely the front continues to diffuse, but this occurs over a much longer time scale than the initial propagation). Our simulations of this process for different parameters (see Fig. 6) suggest that a boundary induced polarization would not reliably spread across hundreds of cells on a plausible time scale.
Bulk orienting signal. Any perturbation that breaks the symmetry of forming a2b versus b2a pairs can potentially act as a bulk signal. Symmetry breaking can occur, for example, through a graded expression of a or b proteins in the tissue, or alternatively, another protein with a graded distribution might sequester or hyper-activate either a or b. Such graded distributions may be expected to arise from morphogen gradients. Yet, the asymmetry on the level of a single cell, due to such an effect is expected to be weak because the concentration gradient of the protein is small on the scale of a single cell. On the other hand, a bulk magnetic field of any magnitude will eventually orient an ordered ferromagnet. In the PCP context, with the developmental time scale of *10 hours corresponding to the PCP amplitude growth stage, an important question to ask is whether a weak bulk field can orient the whole wing within this limited time frame.
To address this question we focus on a particular type of a bulk orienting signal that can be easily quantified. A graded expression of a (or b) within the wing acts as a signal that orients polarity in parallel to the gradient direction. The effect of such a field can be analyzed analytically and is described in the supporting analysis (Text S1, part II).
In our model, a gradient in a expression corresponding to a *10% change across the wing (assuming that the wing is *250 cells across and a 0:05% change between adjacent cells) yields full orientation in the distal direction before amplitude saturates in state B of the phase diagram (Fig. 3A) with N 0~5 000. A signal ten times larger, which corresponds to a two-fold change in  Fig. 5) is initiated with all dipoles pointing in the downwards direction. Stochastic dynamics are then followed to assess whether long range order in the cell array is maintained, and this is done for several different values of N 0 , the average number of molecules per interface. While long range order is maintained for N 0~5 000 and 500 (left and center panels), long range order is destroyed by the stochastic fluctuations for N 0~5 0 (right panels), as quantified by the center panels which track the dynamics of the polarization averaged over all cells (Center panels, amplitude: ½vp x w 2 zvp y w 2 1=2 ; Bottom panels, orientation). The top panels show a snapshot of a subset of cells at the end of the simulation. doi:10.1371/journal.pcbi.1000628.g004 concentration across the wing, is sufficient to achieve full orientation in state A, which is far from the phase transition. These results suggest that a weak orienting signal (e.g. a 0.1% per cell variation in protein level) can effectively orient the wing within the time scale of about 10 hours.
It is possible that the orienting signal is not derived directly from a protein gradient: An early polarization may exist in each cell before the asymmetric localization of the key PCP proteins develops. For example, an early polarization exists in the distribution of Widerborst [29,30]. In addition, recent evidence [19,31] suggests that proteins in the dachs-fat-daschous pathway are asymmetrically distributed as well. Such an early polarization, of a protein other than Fz and Vang, may be a rough readout of a morphogen gradient, and may couple to the dynamics of PCP proteins to establish an orienting signal. In this case the cell-cell interaction in PCP may serve to smooth such a signal, creating a readout that is more spatially-uniform and accurate than the input present in each cell alone. This is demonstrated in Fig. 7, where a noisy orienting signal (yellow arrows) is compared to polarity response (white arrows).

Discussion
Under-expression of fz and Vang. The effect of changing particular parameters of the model may depend on the position within the phase diagram. An example can be seen in the lowerleft part of the diagram in Fig. 3, where increasing the parameter k {1 could either switch from a vertex state to a side state, or vice versa. However, we find quite generally that decreasing a t or b t causes a transition to the uniform state. Reducing protein numbers in the cell corresponds in our model to a simultaneous decrease in a t or b t , and in N 0 which tends to destabilize the ordered state. Hence within our model a decrease in Fz and Vang concentrations increases deviations of hair polarity from the correct distal direction and eventually destroys the ordered state altogether. While our model is in agreement with the broad effect of fz mutation or under-expression, it also generates new and quantitative predictions, as discussed below.
PCP dynamics. The observed asymmetry in distribution of Fz and Vang builds up gradually over a time scale of about ten hours, between 18 and 32 hours after puparium formation [7][8][9][10][11][12]. The simplest interpretation of this observation is that PCP formation takes place during the first stage of the dynamics, before amplitude saturation. The characteristic time scale l should then be of order *10 hours. An alternative scenario is that a local dipole moment builds up in each cell on a much shorter time scale, and that PCP dynamics occurs mostly within the second stage of domain growth. In this latter scenario we expect to observe domains in which polarity points to directions other than the distal one. Since existence of such domains in not reported experimentally, the evidence appears to support the first scenario. Experimental observations were made mostly from static images in which proteins on the two sides of each interface could not be resolved. It will thus be extremely useful to quantify the dynamics of PCP amplitude and orientation, in order to distinguish unambiguously between the two scenarios. Such quantification would make it possible to test the detailed predictions on dynamics.
Swirling patterns in the absence of an orienting signal. The model predicts that swirling patterns should emerge in the absence of an orienting signal. These patterns are consistent with those observed in large fat mutant clones, at least qualitatively. The prediction thus supports the hypotheses that fat mutants lack a coupling with the orienting signal [14]. However, fat mutants differ from wild type tissues in an important way, namely, that their cell arrangement is less ordered than in wild type tissues [15]. There are thus two possible mechanisms leading to disorder in fat mutants: one arising from the role of stochasticity in the absence of an orienting signal, and the other arising from lattice disorder. These two mechanisms are not necessarily mutually exclusive.
Regardless of the mechanism at work in fat mutants, our model predicts that even in an ordered lattice without excess defects, swirling patterns will appear in the absence of an orienting signal, followed by slow coarsening dynamics. We envision three potential ways to test this hypotheses. First, if fat is necessary for the coupling with an orienting signal, but also plays a separate role in lattice repacking, disabling fat activity at a sufficiently late stage of the dynamics, after lattice repacking [32], may inhibit coupling with the orienting signal without influencing lattice order. Second, it may be possible to find other mutations in which lattice order is not disrupted, but the coupling with the orienting signal is absent. Third, it may be possible to negate the effect of the endogenous distally orienting signal by inducing an orienting signal in the proximal direction.
An artificial bulk orienting signal. We predict that graded expression of a or b will act as an orienting signal. This prediction is consistent with experiments in which a gradient in expression of fz was induced using heat-shock promoters, causing inversion of hair-growth direction [33]. A similar effect is expected with a graded expression of Vang.
Inducing a gradient of Fz or Vang protein concentration may provide a way to cancel the endogenous signal in order to test the predictions discussed above. Further, inducing such a gradient in a tunable, quantifiable manner could be a realistic experimental objective, e.g., using light-switchable promoter systems that allow precise spatio-temporal control of gene expression [34]. In addition to testing the prediction that a and b gradients can act as an orienting signal, measuring the magnitude of gradients that induce a significant perturbation in the PCP pattern can provide a way to quantify the magnitude of the endogenous orienting signal and its spatial variation within the wing. Another prediction that could potentially be tested along these lines is that applying an orienting signal only at a late stage of the ordering dynamics will have only weak influence on polarity (see Fig. S2). Finally, inducing an orienting signal in fat mutants could help distinguish between the role of lattice disorder and the role of uncoupling from the orienting signal.
While gradients in fz or Vang expression could be used as a tool to perturb PCP in a controlled manner, experimental evidence suggests that the endogenous orienting signal is not due to a gradient in fz expression [23,33,35,36]: First, a graded expression has not been observed experimentally (although our model suggests that very weak gradients may be sufficient to select the distal orientation). Second, uniform expression of a fz transgene with a heat-shock promoter is sufficient to rescue a null fz genotype.
Bulk vs. boundary signals. Our results demonstrate that a weak concentration gradient within the tissue can produce a reliable response, although the concentration change on the scale of individual cells is very small. The reliable response is achieved by the collective dynamics of the network, which effectively integrates the orienting signal over a region of the tissue larger than the size of an individual cell. Hence our results suggest that in the PCP pathway inter-cellular interactions within the network of cells serve to increase the fidelity of response to a morphogenetic field.
In contrast to the precise readout of a weak bulk signal, a boundary signal cannot effectively propagate in our model over a large number of cells. This result is expected to hold in any model that shares a fundamental aspect with our model, namely, that the uniform state is unstable and gives rise spontaneously to a patterned state driven by noise even in the absence of a global  ordering field. In contrast, in an excitable system where the uniform state is stable, it may be possible to achieve patterning by a propagating front -as observed, for example, in the morphogenetic furrow during drosophila eye development [37].
Phenomenological models in Biology. Modeling in Biology tends to emphasize molecular detail. Yet in biological networks that involve more than a few components the typical situation is that many details are unknown, and it is imperative to devise an approach that can be insightful and predictive even in the absence of complete knowledge. Our strategy was based on building a semi-phenomenological model which attempts to identify the key microscopic aspects (e.g. formation of transcellular heterodimer complexes), build a simple model which parameterizes the many unknowns and systematically identify different regimes of behavior as a function of parameters (e.g. via a phase diagram). We then focus on identifying the observable effects that can help to discriminate between different regimes of the model. For example, the dynamics of intracellular polarization and the ''coarsening'' dynamics that extends local correlations into a global order, are identified as informative quantitative phenotypes deserving careful experimental study. The study is obviously incomplete, as it does not explicitly identify all relevant genes and molecules, but it provides a useful framework allowing to classify phenotypes and accordingly group observed genetic perturbations, and eventually refine the model at an increased level of molecular precision.

Methods
Model. We considered two mechanisms for establishing a non-local inhibitor field in each cell, and results were similar in the two. The examples used in the manuscript use one realization of the model which is summarized below, while a full description of both mechanisms is provided in the supporting analysis (Text S1, part I).
Fast diffusion. We assume that protein diffusion is fast. The relevant time-scale in this context is the typical time for diffusion of a membrane protein from one side of a cell to the opposite side, which (assuming a diffusion coefficient D^0:5m 2 =sec) is of the order of 10 minutes. In comparison, the asymmetric pattern of protein localization arises on a time scale of several hours, so that separation of time scales appears to be a reasonable approximation. Assuming fast diffusion, free a and b concentrations are uniform in each cell. Similarly, diffusion of bound a/b complexes equilibrates their concentration on any given interface. Complexes, however, cannot diffuse from one interface to another without unbinding first, and reforming with new constituents -processes that we assume are slow. Hence complex concentrations can vary between interfaces belonging to the same cell and the dynamics do not necessarily lead to a uniform steady state. For a regular hexagonal array of cells one then needs to keep track of six variables per cell representing the total numbers of interfacial complexes u a i (with i labeling the sides of cell a).
where u 1,i are concentrations of complexes on the six sides of cell 1, having an a within that cell and u' 2,i are concentration of complexes in cell 2 with a b in that cell, and where u 1,0~u ' 2,0~u1 (Fig. S3). The rate K9 is a constant that we set to unity by rescaling time, and K is given by where the term su 1 describes self-excitation, and vxw is the average on the interface of the non-local field x, which represents the fraction of unmodified a proteins. This field obeys Eq. (3), in which r s is a one-dimensional coordinate ranging from 0 to 6S, S~L= ffiffi ffi 3 p is the length of a cell side, and u 2 (r s ) is a step-wise uniform function equal to u 2,i on each of the sides of cell 1.
All the parameters in this equation are dimensionless: we rescale all concentrations by the total concentration of a proteins, so that a t~1 . Lengths are rescaled by setting L~1. The independent parameters in the model are thus K 0 , b t , s, a, k, and N 0 . The parameters of states A and B, as well as the value of l in these states are summarized in Table 1.
The noise term in Eq. (4) is Gaussian with covariance To derive this relation, recall that all concentrations were rescaled so that a t~1 . The total concentration a t~1 of a proteins corresponds to having N 0 molecules per interface, by definition of N 0 . Hence the number of u complexes on the interface is given by N 0 u 1 (u 1 is thus the fraction of a proteins that participate in a complex). Assuming Poisson statistics, the variance in the number of reactions per unit time is given by N 0 (Ka 1 b 2 zK'u 1 ) which, after division by N 2 0 yields Eq. (8). Stochastic simulations. We used a forward explicit Euler method for simulating the stochastic equations on a lattice of cells. In each step a set of 12 linear equations are solved in each cell to obtain the field x on the membrane and its average vxw on each of the six sides. Eq. (4) is then used to update the two complex concentrations on each interface. The time step was 10 {3 . A typical simulation such as the one in Fig. 5 requires *10 hours on an Intel Core 2 processor.
Global field. When analyzing the effect of graded a expression we use an equivalent constant field, projected onto the dipole carrying modes, as described in Text S1, part II (Eqs. S54-S55), rather than include explicitly a graded expression of a, in order to avoid boundary effects. However, we also ran simulations on large lattices with direct gradients of a to verify the applicability of Eq. (S55). In Fig. 7 a local field was associated with each cell. The dynamic equations at each interface involved a local field taken as the average of the fields on the two cells separated by the interface. Dipole moment. We define the magnitude of the PCP dipole in each cell as P~(M a =a t {M b =b t )=2L where M a,b is the dipole moment of a/b protein distribution, M a~P i a i (r i {r cm ) and the sum is over all sides of the cell, a i is the a concentration on side i, r i is side i's center, and r cm is the cell's center. A similar equation holds for M b . The only contribution to M comes from the complexed proteins because the free proteins are uniformly distributed in the cell. where K~10 5 , a 1~b2~1 {u 1 , K'~au 2 =(1zau 2 ) ½ 2 and a~0:025. Phase diagram. The phase diagram (Fig. 3) was obtained as follows. In the deterministic limit the steady states of the system are spatially uniform. Hence the problem reduces to that of finding the steady states of a six-dimensional dynamical system. The phase space was first sampled at 50|50 discrete loci to obtain a lowresolution representation of the phase diagram. At each point all steady states (stable and unstable) were found using a multidimensional secant root-finding algorithm as described in [38], initialized with 500 different random states. For each stable state found in this way, the stability and symmetry properties were then determined. In all cases there was a unique stable steady state up to the symmetry: Either a single stable uniform state, or six equivalent stable vertex states, or six equivalent stable side states. After obtaining the low-resolution representation of the phase diagram, we used the continuity of the phase transitions in order to obtain precise phase boundary curves, by solving numerically for loci where an eigenvalue of the Jacobian matrix vanishes. Figure S1 Coarsening dynamics. Coarsening dynamics in a stochastic simulation without an orienting signal, at several time points: lt = 2.5 (A) -before amplitude saturation, lt = 10 (B), 50 (C), and 100 (D). Parameters are the same as in Fig., Figure S2 Effect of delayed application of the orienting field. Dashed lines show the average polarization in the distal direction when a bulk orienting signal is applied only from lt = 5, compared to the dynamics when the field is applied from the simulation onset (full lines). Black, red, and gray traces correspond to three different magnitudes of the applied field. These correspond, respectively, to gradients in a concentration that amount to an increase of 0.5%, 0.2%, and 0.1% from each cell to its proximal neighbor.