The role of intracellular interactions in the collective polarization of tissues and its interplay with cellular geometry

Planar cell polarity (PCP), the long-range in-plane polarization of epithelial tissues, provides directional information that guides a multitude of developmental processes at cellular and tissue levels. While it is manifest that cells utilize both intracellular and intercellular interactions, the coupling between the two modules, essential to the coordination of collective polarization, remains an active area of investigation. We propose a generalized reaction-diffusion model to study the role of intracellular interactions in the emergence of long-range polarization, and show that the nonlocality of cytoplasmic interactions, i.e. coupling of membrane proteins localized on different cell-cell junctions, is of vital importance to the faithful detection of weak directional signals, and becomes increasingly more crucial to the stability of polarization against the deleterious effects of large geometric irregularities. We demonstrate that nonlocal interactions are necessary for geometric information to become accessible to the PCP components. The prediction of the model regarding polarization in elongated tissues, is shown to be in agreement with experimental observations, where the polarity emerges perpendicular to the axis of elongation. Core PCP is adopted as a model pathway, in term of which we interpret the model parameters. To this end, we introduce three distinct classes of mutations, (I) in membrane proteins, (II) in cytoplasmic proteins, and (III) local enhancement of geometric disorder. Comparing the in silico and in vivo phenotypes, we show that our model successfully recapitulates the salient phenotypic features of these mutations. Exploring the parameter space helps us shed light on the role of cytoplasmic proteins in cell-cell communications, and make falsifiable predictions regarding the cooperation of cytoplasmic and membrane proteins in the establishment of long-range polarization.


Formalism of The Generalized Reaction-Diffusion Model
Like in the Main Text, we denote the polar complexes by F-G, namely F ≡ Fz:Fmi and G ≡ Fmi:Vang. Labeling two adjacent cells by i and j, the average concentrations of complexes F i -G j is denoted by u ij ; consistently the concentration of F i -G j is denoted by v ij = u ji . The junctions shared by cells i and j, are denoted by ij hereafter. The formation of complexes F i -G j is upregulated by the already-bound complexes of the same polarity, and downregulated by the opposite complexes, i.e. G i -F j . The key assumption in this model is that such intracellular interactions, mediated by diffusing cytoplasmic proteins (Dsh, Dgo, and Pk), are not strictly local, but instead diffuse to mediate the nonlocal interaction between the complexes localized on different junctions. The magnitude of interaction, which is proportional to the local concentration of the diffusing messenger protein, decreases exponentially with the distance between the two complexes. The concentration u ij of F i -G j complexes at point r on the perimeter of cell i is governed by the following equation: In the above equation, {k} i stands for the labels of all neighbors k of a given cell i, and i∩j dr implies integration along the edge shared by cells i and j. The factors K u (r − r ) and K d (r − r ) are the kernels of up-and downregulating interactions, and γ −1 is the dissociation timescale of bound complexes. As mentioned in the Main Text, the kernels K(r) are assumed to be of the exponential form, which is deduced from the diffusive dynamics of the cytoplasmic proteins. The concentration of cytoplasmic protein obeys a diffusion equation with diffusion constant D and degradation time τ : ∂c(r, t) ∂t = D∇ 2 c(r, t) − τ −1 c(r, t).
Assuming the degradation of C is much faster than polarization dynamics, γτ 1, it suffices to only consider the steadystate solutions of c(r, t) in Eq. (2): Here, r corresponds to the specific point on the cell membrane, where the the concentration of C is the linear superposition of the concentration of proteins diffused from all around the cell; hence integration dr . Next, λ = √ Dτ is the diffusion length of protein C before degradation. The diffusing proteins enhance/suppress the formation of like/unlike complexes. The linear dependence of cooperative terms on the concentration u(r ) originates from the fact that in Eq. (3), the source concentration of modified proteins, c 0 (r ), is assumed to be proportional to the local concentration u(r ). This dependency is a function of the binding affinities of cytoplasmic proteins and the membrane-bound proteins. Altogether, the coefficients are lumped into the phenomenological constants α and β. We shall emphasize that the above derivations are merely to provide a qualitative relationship between the degradation rate, diffusion constants, and the cytoplasmic interaction length scales. As mentioned in the Main Text we are aware of the complicated dependencies of λ on the molecular details which are not captured by the above equations.
The functional forms of the interaction kernels read: K u (r) = N −1 u exp(−|r|/λ u ) and K d (r) = N −1 d exp(−|r|/λ d ), where λ u , λ d are the characteristic length scales of up-and downregulating interactions, respectively. The prefactors N u and N d are normalization factors, to be determined shortly. Before calculating the normalization factors, we make a detour to explain an approximation we used to reduce the computational cost of simulations.
Uniform junctional distribution of proteins. The full integro-differential equation (1) couples every pairs of points on the perimeter of a cell; solving it on a 2D array of cells is computationally expensive. A simplifying approximation is to assume uniform junctional distributions of proteins, and break the integration on the cells' perimeters down to summation over the junctions. Thus the integrations are replaced by matrix products; see Eqs. (4) and (5). For notational convenience, we introduce the following: cells are labeled by the Latin letters, i, j, . . . , and comprise edges labeled by Greek letters, µ, ν, . . . . For example: i ≡ {µ, ν, . . . }, and consistently edges, that are the cell-cell junctions, are the intersections of two sets (cells): e.g. µ ≡ i ∩ j. The junctional averaging is then defined as: In the above equation, η µ and M µ are the noise and the global cue averaged over the length of edge µ. The cooperative interactions are reduced to matrix products of K u and the vector u. Matrices K u are symmetric with respect, and their elements K µν u are purely geometrical constants obtained from the following relation: A similar relation holds for downregulating interactions K µν d . In Eq. (4), coupling constants equal α µ = α/ µ , β µ = β/ µ , and for the kernels we have K µν = K νµ = µ,ν drdr K(r − r ). The self interactions, i.e. the diagonal elements equal in which x µ = µ /λ. We dropped the subscripts u and d for simplicity. The effective stochastic noise on a junction µ reads: . It is noteworthy that in reality the core proteins, for example Flamingo (Fmi) and Frizzled (Fz), are observed to be persistently localized at the puncta, the subdomains of plasma membrane; see e.g. [1].
With this background, we now calculate the normalization factors introduced above. In order to discern the net effect of interaction ranges from the effective coefficients α, β, we choose the normalization factors N u and N d , such that the self-interaction of an edge of length µ equals α, namely: α K µµ u = α, or K µµ u = 1. The same relation holds for β and K µµ d . This choice of normalization ensures that the observed behavior upon changing λ's is purely due to nonlocal edge-edge coupling, not the effective coefficients of local interactions α, β. Satisfying this condition for all edges simultaneously is not possible, except for ordered tissues. The normalization constant is thus calculated for an edge with the average length of all edges, 0 ≡ µ . Using the definition of kernels for an edge of length 0 from Eq. (5), we obtain: Here, both r and r lie on a single edge with length 0 .
Limit of strictly local cytoplasmic interactions (SLCI). In the limit of small λ/ µ → 0, we get, K µν = δ µν , where δ µν is the Kronecker delta. Thus in the SLCI limit, the equations read, Junctional and cellular polarity. Planar cell polarity can be defined at the junctional and the cellular levels. For individual junctions, polarity is defined as the difference between the concentrations u ij = [F i − G j ] and the opposite dimer, u ji , Cell polarity, on the other hand, is referred to as the asymmetric (anisotropic) distribution of membrane proteins on the cell perimeter, thus requires integration on the cell-cell junctions.
Depending on the pathway of interest, polarity can be defined either as a vector (vectorial polarity), or a nematic (axial polarity), corresponding to heterodimer and homodimer clusters, respectively. The former is a vector identified by an angle ∈ [0, 2π), whereas the latter is a traceless tensor with the angle in the ∈ [0, π) range. In our case of study, the polarization involves two distinct complexes, Fz:Fmi and Fmi:Vang, hence vectorial polarity. Below, we introduce both definitions.
(i) Vectorial Polarization: The polarization vector (also called dipole) associated with cell i is defined as: Here R i is a reference point within cell i, with respect to which the polarization is defined. We take this to be the geometrical center of mass of each cell. The magnitude of the cell polarity |P i |, and the angle φ p i ∈ [0, 2π) measured from the x-axis, can be computed using the following relations: The polarization of the tissue with N c cells, and its global order are characterized by the following quantities: (1) magnitude of average polarization, (2) average magnitude of polarization, and (3) the ratio of (1) and (2): (ii) Axial Polarization: The second definition is a measure for cell polarity when dealing with homodimers like Fmi-Fmi complexes, and determines the axis of polarization: In the above equation, φ is the polar angle of point r on the perimeter of cell i and is measured with respect to a reference axis. The magnitude of polarization equals P i = (P 2 i,1 + P 2 i,2 ) 1/2 . Its orientation is determined by the angle φ p i that satisfies cos(2φ p i ) = P i,1 /P i and sin(2φ p i ) = P i,2 /P i . In Sec.
(3), we use the same formalism for the elongation tensor.
Correlation function and correlation length. Correlation function is a measure of alignment of polarization field. In order to investigate the temporal behavior of the spatial extension of alignment, we define the equal-time correlation functions. Consider an arbitrary cell at R i , and a vector r connecting it to another cell at R j = R i + r. With no further assumption we can define a correlation function, that is dependent on the distance r = |r| and the relative angle of r and P(R i ); we call it θ r,p . The latter appears due to the vectorial nature of polarization field; there is no a priori reason for dipoles to be correlated equally in all directions. For a tissue of N c cells, the correlation function at time t reads, The above quantity calculates at each timepoint t, the average (over all cells i) of the conditional probability that the dipole at point R j = R i + r acquires the value and direction of P(R j ), given the dipole at point R i equals P(R i ).
For θ r,p = 0, π, and θ r,p = ±π/2, we get longitudinal and transverse correlations, respectively. In spite of this angular dependence, averaging the above correlation function over θ r,p ∈ [0, 2π), returns a weighted average of correlation as a function of r = |r|, which is bounded by the longitudinal and transverse correlations from above and below, respectively. In spite of the distinction between longitudinal and transverse correlation lengths, the qualitative behavior of the correlation length could be captured by the angular average of the correlation function. Further discussion in this regard can be found below in SI. Sec. (2.3). Thus, we define radial correlation function: Correlation length can be obtained from the above equation: Here, R c = 40 (cell diameter) is the largest distance between two given cells in our simulations. A perfectly correlated polarization field, is identified by ξ = R c /2. The time-dependent correlation lengths for three important cases of study are shown in the Main Text Fig. (3c2).
Geometric disorder. In 1D systems, the edge lengths are i = 0 (1+ i ), where i ∈ [− 0 , + 0 ] with uniform distribution, and i j = 2 0 δ ij /6. In 2D the level of quenched disorder is controlled by randomizing the sites of a triangular lattice, based on which a polygonal lattice is generated using Voronoi tessellation. Suppose the nodes of a triangular lattice are perturbed: i } the loci of the nodes in the ordered triangular lattice. The magnitude of the disorder term ∆ i is uniformly distributed in range [−∆ 0 , +∆ 0 ], with local correlations: ∆ i · ∆ j = ∆ 2 0 δ ij /3, where δ ij is the Kronecker delta function. For ∆ 0 = 0, the corresponding Voronoi lattice would be an ordered hexagonal lattice, hence 0 = 0. By randomly displacing the nodes of the triangular lattice, we can distort the resultant Voronoi lattice. The edge-length disorder of the Voronoi lattice i , and density of defects n d are obtained by ensemble averaging over 10000 realizations of disordered triangular lattices of size 50×50. Defects appear for disorder ∆ d 0.25, in the triangular lattice.
Physical considerations. Given the quantitative approach of this study, we find it crucial to clarify the term "long-range", used frequently throughout the paper. The Mermin-Wagner theorem states that "true long-range" ordering is prohibited in 2D systems with continuous (e.g. rotational) symmetries, except at zero stochastic noise. The long-range order is referred to as the algebraic decay of correlation functions with distance. The magnitude of noise in our system drops as 1/ √ N mol. , with N mol. the number of molecules participating in binding/unbinding reactions. Thus in the limit N mol. → ∞, long-range order is achieved. For finite N mol. , a state of quasi-long range order can potentially exist.

Mean-Field and Numerical Solutions
Mean-field (MF) solutions are obtained by assuming that for all junctions, e.g. shared by adjacent cells i and j, the products f ubd i g ubd j and f ubd j g ubd i are uniform across the tissue. This assumption is backed up by the following numerical analysis. We compute the root-mean squared fluctuations of f ubd i g ubd j , for all the junctions with fully randomized labels, and compare the initial and final values. Figures (Bb1) and (Bb2) show that in 2D systems for both SLCI and NLCI mechanisms, the final distribution of f ubd i g ubd j is noticeably narrower than the initial distribution. Intuitively this approximation can be justified by the fact that the linearized reaction-diffusion equations governing u and v, obey diffusion-like equations in the continuum limit (see [3] where 1D arrays were studied). Therefore, given that the total concentrations f 0 and g 0 are uniform across the tissue, the free proteins too, spread diffusively into a rather uniform state. Here we first elaborate on the MF solutions in 1D, then discuss the 2D case.

Mean-Field Solutions in 1D
In one dimension, cells are juxtaposed in an array and are separated by junctions. The proteins localize on both sides of these junctions and form dimers. A schematic of one-dimensional arrays can be seen in Fig. (Aa). The reaction-diffusion (RD) equations governing u and v complexes on a junction shared by cells i and i + 1 read: A similar equation exists for the opposite complex v i,i+1 = u i+1,i . The general RD equations introduced in Eq. (1), reduce, in 1D with strictly local interactions, to the above equation that was originally proposed by Mani, et.al. in Ref. [3], where the properties of the solutions to Eq. (16) in 1D systems are investigated elaborately. Here we briefly review the main results and show that in the steady state, the mean field (MF) solutions of ordered systems exhibit a bifurcation at a critical value of the control parameter g 0 /f 0 , i.e. the ratio of total amounts of proteins F and G. We start by deriving equations for , by adding and subtracting Eq. (16) and its counterpart for the opposite complex:  In the above equations we have: in which F 0 is the total amount of protein in a cell, assumed to be uniform across the tissue, and f 0,i is the concentration of free F available to the two junctions of cell i. A similar relation can be written for g ubd i and g 0,i . In ordered systems and within the MF approximation, all of the above variables become independent of index i. A few lines of simple algebra shows that the critical value reads: For the values of the parameters used in the simulations, i.e. α = β = 5, and γ = 1 , κ = 10, the value of g * 0 is found to be g * 0 = 0.225, which can also be seen in Fig. (2) of Main Text. Note that in the above equations, for αβ → 0, the bifurcation point diverges, implying that cooperative interactions are essential to the emergence of collective polarization. A closer look at Eqs. (17) makes it clear that the relevant parameters in the local bistability of the polarization are f ubd i g ubd i+1 , and f ubd i+1 g ubd i . When normalized by the total concentration f 0 = F 0 / 0 as the unit of concentration, one can rewrite the above parameters in terms of f 0,i g 0,i+1 and f 0,i+1 g 0,i . Therefore, for disordered systems where f 0,i and g 0,i are no longer uniform, the local bifurcation point is randomized and the singular transition turns into a smooth cross-over.
Numerical solutions, Fig. (Ac) and (Ad), suggest that in the limit of small stochastic noise and initial bias, the steady state is not guaranteed to be uniformly polarized. The initial imbalance of protein distributions is defined as p 0 = |u 0 −v 0 |, with u 0 , v 0 , the spatial averages of initial dimers' concentrations. The bias δp 0 /p 0 , is the normalized magnitude of spatial fluctuations of initial polarity; small and large biases correspond to δp 0 /p 0 1 and δp 0 /p 0 1. While in ordered systems, a moderate initial bias suffices to achieve a uniform polarization, the patterns of polarity in highly disordered systems exhibit robustness, and are largely determined by the local geometry (disorder) of the array. Therefore as we observe in 1D, the quenched disorder imposes undesirable solutions, impairing the faithful transduction of directional information through PCP signaling. As we will see below, the situation gets only more complicated in two dimensions.

Mean-Field Solutions in 2D
The systems in one and two dimensions show inherently different behavior. In 1D, the proteins have only two junctions to localize. This limited number of choices and the resultant predictability are absent in two dimensions. Due to the large number of possible steady states in 2D, the initial configuration, as well as stochastic and geometrical disorders, influence the final state. We show, in this section, that nonlocal cytoplasmic interactions (NLCI) destabilize the vast majority of unpolarized fixed points in favor of the polarized ones. Furthermore, we find an optimal range of the NLCI length scale that assists with establishment of long-range alignment. In the following, we present the results for a set of parameters which lies within the polarized regime: in the units of f 0 = γ = 0 = 1, we set α = β = 5, κ = 10, and g 0 = 1.
The numerical solutions are presented below, but first let us attempt to gain some insight using analytical MF solutions for ordered tissues. Again, MF solutions imply uniform distributions of f ubd i g ubd j for all junctions {ij}, and recall that it may be justified by noting the diffusive nature of p, s dynamics, which in turn implies the diffusive dynamics of the free membrane proteins. We checked this assumption numerically. The value of f ubd i g ubd j , normalized by the mean, is plotted for all edges at the initial and final states in Fig. (Bb1) and (Bb2). In steady state, the standard deviation of the distribution is nearly independent of initial condition as well as the model parameters, and remains below 0.05, for both SLCI and NLCI regimes. Under this assumption, the RD equations take the same form as in 1D, except the total amount of proteins Fz and Vang are shared by six junctions instead of two. For ordered tissues, by virtue of sixfold symmetry, all junctions equally absorb the proteins, namely p, s are the same on all edges. Thus three edges would carry inward, and the other three carry outward dipoles.
Trivial mean-field solutions. Trivial mean-field solutions, in addition to the MF assumption discussed above, enjoy translational invariance of polarization along each of the three axes separately. The trivial MF solutions acquire two distinct configurations illustrated in Figs. (Ba1) and (Ba2), respectively. While the former represents a state with nonzero net polarization, the latter has zero net polarization. The solutions of the second type, are however destabilized when NLCI is included in the model. The net polarity in the polarized case can be calculated as follows: where p e is the magnitude of polarization of one edge calculated above, and θ is the angle between the two adjacent edges in ordered hexagons, which equals p c = 2 s 2 − 4/αβ, for θ = 2π/3.
We also consider a different situation in which edges have unequal α, β's. At this point, different values of α, and β can have various origins, that are beyond the scope of our discussion. However in the Sec. (3), we argue that unequal parameters can be a consequence of nonlocal interactions in elongated cells. Assume the three pairs of parallel edges acquire coefficients α 1,2,3 and β 1,2,3 . We consider two scenarios: (i) α 1 > α 2 = α 3 , and (ii) α 1 = α 2 > α 3 . From the results we found in the case of sixfold symmetric lattices, the onset of bifurcation is inversely proportion to s * ∼ 1/ √ α. Thus in the first scenario, axis 1 is the first one to show instability upon increasing g 0 above g * 0 . Therefore, we have f ubd g ubd = γ/(kα 1 ), where, and similarly for g ubd . For the unpolarized edges we have: Using kf ubd g ubd = γ/α 1 , we get: Defining f eff 0 and g eff 0 , , thus: and similar expressions for g eff 0 and g ubd , we get for p 1 , s 1 : From the above analysis, we learn that in systems with unequal α's and β's, the problem essentially reduces to a one dimensional systems, with effective values for the concentrations of free proteins Fz and Vang. In the second scenario where α 1 = α 2 > α 3 , the third axis remains unpolarized as the other two are effectively more absorbent, and share the total bound proteins; thus are equally polarized with fourfold symmetry. As mentioned above, one situation of interest in which the coefficients α, β acquire unequal values on different edges, is the elongated tissues. In such systems, the above two cases correspond to Figs Nontrivial mean-field solutions. The nontrivial solutions in 2D satisfy a weaker constraint; they do not possess the translational invariance of cell dipoles. The full analytic treatment of this problem is cumbersome. We briefly touch upon the solutions to provide some insight into the richness of the basin of attraction for a system in SLCI limit. The constraint of nontrivial solutions demands the junctional polarization p and the total amount of proteins s to be identical for all junctions. Thus in practice the only constraint is that three junctions are positively polarized and the other three carry negative polarizations; i.e. three outgoing and three incoming dipoles; see Fig. (Ba3). The numerous configurations satisfying this condition, outnumber the trivial solutions by a large margin. Note that the trivial solutions also satisfy the above constraint, but are very special configurations. Therefore, if the segregation is not carried out systematically across the tissue, the membrane proteins are extremely unlikely to redistribute in a polarized fashion; unless the initial distribution to begin with is nearly polarized. As mentioned previously, nonlocal cytoplasmic interactions destabilize the majority of nontrivial solutions, and drive the system towards a segregated state where the polarized trivial solution appears. for randomized edges labels (ij), in SLCI and NLCI systems, respectively. The variance of the distributions decreases significantly over time, supporting the MF approximation.

Numerical Solutions in 2D
While in relatively ordered tissues with NLCI, and in the absence of orientational cues, the polarity is determined purely by chance (stochastic noise and initial conditions), highly disordered systems show robustness against such random factors. Instead, the fixed points of polarization fields are determined collectively by the geometry of the lattice. In finite-size systems, the geometrical disorder provides a bias towards one orientation over others. This effect gets progressively more pronounced with increasing range of NLCI and/or level of the quenched disorder. External cues of sufficiently large magnitudes, however, reorient the polarity. A typical configuration of the steady states is illustrated in Main Text Fig. (3a2). We find a range of interaction length scale, 0.2 λ/ 0 0.7, for which the NLCI guarantees the long-range alignment of polarization. The directional correlation shows a peak at around λ/ 0 0.5. This range of λ also depends on the geometric disorder, which hinders the establishment of polarization, and increases the lower bound of the range. For example, with 0 = 0.6, the range changes to 0.25 λ/ 0 0.7. Before addressing the segregation, we note that for λ/ 0 0.7, all edges within a cell strongly couple to each other, rendering the system fragile against stochastic noise.

Segregation mechanisms in SLCI vs. NLCI. Through a comparison of the dynamics of SLCI and NLCI in Figs. (3b1)
and (3b2) of the Main Text, the role of cytoplasmic nonlocal interactions in cell-cell interactions becomes evident. While the average polarization P (t) remains negligible in SLCI limit, it saturates to the average magnitude Q(t) in systems with NLCI, which reflects the angular correlation of cell-cell polarities. More elaborately, during the first stage of dynamics, the nonlocal cytoplasmic interactions prepare each and every cell for later intercellular communications. The coarsening and propagation of polarization is then carried out by cell-cell interactions, which increase with the magnitude of cellular dipoles Q i (t). Here, an important question arises, regarding the interpretation of Q(t): Can Q(t) be thought of as a measure of cellular segregation of proteins? A naïve guess would be that since Q i (t) is oblivious to the direction of polarity, it only measures the magnitude of the dipoles per cell, which is a candidate for quantifying segregation. This, however, warrants a careful investigation, as one can see that Q(t) shows arguably similar behavior in both SLCI and NLCI cases. Does this rule out the lack of segregation in SLCI? One possibility is that while the average Q(t) increases rapidly like in the NLCI case, the segregation is not achieved consistently in all cells. This is in part due to initial condition. Consider a tissue with randomly distributed proteins on the membranes. The magnitude of a cell's dipole increases due to the localization of some of the free proteins on the membranes. Above the polarization instability, namely for large enough g 0 , edges with higher initial concentration of a certain protein absorb more free protein of the same kind due to the cooperative interactions. Therefore, the final polarity depends, among other factors, on the initial condition. This is a separate effect from nonlocal interactions, and is built in the nonlinearity of RD equations regardless of the length-scale of nonlocal interactions. In order to directly measure the segregation, we define partial polarities as follows: and similarly P G i is obtained by substituting v i (r) for u i (r). Perfect segregation then corresponds to P F i = −P G i in the steady state, regardless of the initial distributions. The less the level of segregation, the less the deviation from the initial condition over the time evolution: P F i ∼ P F i (t = 0) and P G i ∼ P G i (t = 0). Therefore, by comparing the final to initial partial polarities, one can easily signify the differences between NLCI and SLCI mechanisms, in terms of the cytoplasmic segregation. In order to quantify the segregation level, and compare that in the two mechanisms of SLCI and NLCI, we introduce the following quantities: 1. the spatial average of the magnitude of: S i = P F i + P G i . Note the vector sum, thus for perfect segregation we get: S i = 0. (Overlines mean spatial average over all cells at a given time). 2. the standard deviation of Q i normalized by its mean,

S(t)
The former characterizes the asymmetry of protein distributions, and the latter measures the consistency in segregation among the cells. We plot the above quantities as functions of time in Figs. (Ca1) and (Cb1), respectively. In (a1), while the ratio approaches zero in steady state for the system with NLCI, it remains finite ( 0.5) in the SLCI case, clearly showing the lack of segregation. In (b1) we see that in SLCI, the normalized standard deviation drops slightly from 0.5 at t = 0 to 0.4 in steady state, whereas in the NLCI case, it drops to nearly 0. The latter implies that segregation is fully achieved in all cells, and the individual cell polarities are very much close to the average polarity. Therefore, as suspected, in SLCI only the average value of Q grows, whereas cells are not coherently polarized across the tissue. In order to show explicitly the cell-by-cell distributions of initial and final values of S i and Q i , and compare SLCI and NLCI cases, we plotted these quantities in Fig. (C). In (a2) and (a3), the initial (red) and final (blue) cell-by-cell distributions of S i /Q i are illustrated for SLCI and NLCI, respectively. Again, the nearly zero final values of the ratio for NLCI shows perfect segregation. The cellular polarity normalized by spatial average at the corresponding time-point, i.e. initial and final, are shown in (b2) and (b3). The width of the distribution shrinks dramatically in NLCI, whereas it remains comparable to its initial value in the SLCI case.
Longitudinal vs. transverse correlations. As discussed above, and according to Eqs. (14) and (15), correlation function and length are dependent on the relative angle of the reference dipole and the connecting position vector, and at least in our case is stronger in the longitudinal compared to transverse directions. Intuitively, the polarity of a given cell points towards the edges with higher concentrations of localized proteins. These edges are shared with neighbors that are located rather longitudinally with respect to the axis of polarity. The polarities of these neighbors too, are influenced by the shared edges. Therefore longitudinal correlations are stronger than the transverse correlations. This discrepancy leads to formation of transient vortex-like structures, but the steady-state behavior of the correlation length remains qualitatively the same. We ignore this effect and compute the effective correlation length by averaging over all angles.
Vortices and saddles. Several theoretical [4,5,6] and experimental studies [7,8,9] have observed swirls and saddles as different forms of the so-called "topological defects". Such defects appear either as steady or transient patterns. Steady defects can be an indication of mutations of various origins; geometric [9], or genetic [10,11,4,12,13,14,11]. In  Table (1) of the Main Text. (b) A system with g 0 = 0.25 (i.e. close to polarization threshold), and "wild-type" interaction length scale λ/ 0 = 0.5. Both systems exhibit swirling and crossing patterns that appear as long-lived steady patterns. We picked one ordered and one disordered tissue. However, in both cases of small λ and small g 0 , defects appear irrespective of the geometric disorder.
Drosophila wing, where Fz and Vang localize distally and proximally, respectively, the coupling to global cues is believed to be dependent on the presence of Fat [10,9]. While small fat − clones ( 10 cell diameters), exhibit little deviation from wild-type polarization, swirling patterns appear in larger clones, where the polarity is aligned over clusters of (roughly) 10 cells; also implying that Fz feedback loops are left intact in fat − patches. Therefore, the propagation of polarization across neighboring cells is carried out through Fz feedback loop, and the global alignment is achieved through coupling to the cues. Our model predicts both transient and steady swirls, depending on the values of the model parameters. We observe two distinct types of steady defects: (a) As λ is increased from SLCI to NLCI regime, there exists a narrow range 0.05 λ 0.15, over which vortices and saddles appear as long-lived patterns; Fig. (Da). (b) Another situation where swirls appear is when g 0 g * 0 , i.e. under-expression of one of the membrane proteins, say Vang, that is interpreted as a global mutation within the context of our model; Fig. (Db). Such patterns are indeed observed in Vang mutants [7]. Local mutations of various kinds are fully discussed in Main Text and below in SI. Sec. (4).
An interesting observation regarding the second type, i.e. g 0 g * 0 , is that for each specific realization of the geometrically disordered tissues, the characteristics of the steady-state pattern of polarity (e.g. positions of swirls) are left almost intact, as g 0 increases from 0.2 to 0.3. The magnitudes of dipoles, however, increase. Consequently, the swirls and branches gradually merge and align for g 0 0.3, and long-range polarity emerges. This behavior, in disordered tissues, is independent of the initial condition and stochastic noise, implying that the patterns of polarization are fixed by the the cell geometry of the tissue in regimes where the correlation lengths are of the order of a few cell diameters. In general, and as briefly mentioned above, in NLCI the geometrical information is detected through nonlocal interactions. The coupling to geometry provides local biases for cell dipoles. As g 0 increases, the correlation length increases and the global direction is chosen collectively.
Equal vs. unequal length scales of cytoplasmic interactions. In the Main Text we discussed the effects of unequal nonlocal cytoplasmic interactions that appear in upregulating and downregulating kernels, i.e. λ u and λ d . Here we elaborate on that discussion and consider various cases to investigate the role of these length scales, and the behavior of the angular correlation of polarity. In each case the geometry of the tissue is held fixed. The geometric disorder is denoted by 0 .
(1) For 0 = 0.5, find the angular correlation for different values of λ u = λ d = 0.01 -0.8. (1) For equal length scales, we observe that the angular correlation increases with λ. The standard deviation of polarity is found to be smaller than 30 • , for 0.2 λ 0.7, and is maximized at around λ 0.4 -0.5. For λ 0.7 the stochastic noise destroys the long-range polarization drastically. The reason is that the stabilizing and destabilizing forces compete at all points on the perimeter of each cell, namely all points are strongly coupled to one another. Thus the segregation becomes progressively more unstable as λ is increased; see (2),(3) For geometric disorder 0 0.45, and for 0.2 λ d 0.7 there is no detectable difference between the angular correlation of systems with different values of λ u 0.8. Beyond this value, the correlation declines slightly, and is eventually destroyed for larger values. In order to further examine the importance λ u in disordered tissues, we repeat the same simulations, but for larger geometric disorder 0 = 0.6. Interestingly we realize that larger disorder impedes the Directional cues. We consider two types of cues, bulk and boundary signals, each of which may be persistent or transient. Bulk cues couple to the F complex across the entire tissue, whereas boundary cues couple only at the boundaries. For bulk cues in, say +x direction, we use the gradient cues of constant slope in each cell: M ij (r) = M 0 (x ij − X i ); where M 0 is the slope of the cue, x ij is the x-coordinate of the points on junction (ij), and X i is that of the centroid of cell i. We simulated the response of the polarization field and observed that NLCI significantly enhances the sensitivity of the polarization field to such global cues. Before proceeding, we shall mention that there exist two time scales in this analysis: the response time scale of polarization field τ res , and the persistence time scale of the cue τ 0 . The results in a nutshell, are as follows: (1) Reorientation of dipoles for fixed τ res and τ 0 requires weaker cues in systems with NLCI than those with SLCI. For persistent cues (γτ m 100), NLCI responds to signals as small as M 0 0.05 over γτ res 2, whereas SLCI requires at least M 0 0.5. The minimum M 0 increases for smaller τ m 's. For example, in NLCI with M 0 0.05 a nearly persistent signal is required, whereas for larger cues M 0 1, even a rapid transient signal γτ m 1 is sufficient to rotate the dipoles over the same time scale. (2) In accord with (1), and due to the small correlation length in SLCI systems, the detection of a cue in these systems happens over exceedingly larger time scales, compared to NLCI with the same magnitude. (3) In the case of boundary cues (i.e. a column of polarized cells), in nearly ordered ( 0 0.2) tissues with zero stochastic noise, SLCI suffice to detect the signal. Presence of geometric disorder and/or stochastic noise, however, necessitates NLCI for the dipoles to align with the cue. Given that the onset of PCP alignment precede the geometrical ordering of the tissue [2], NLCI seems to be the key to the detection of directional cues. Finally, an interesting observation is that, (4) NLCI systems appear to detect sufficiently large initial boundary signals. Initial boundary signal is implemented by polarizing a column of cells, with significantly larger asymmetry of protein distributions compared to the bulk cells. This implies that a temporary boundary signal would in principle be able to rotate the dipoles, should the cytoplasmic interactions be nonlocal.

Elongated Cell Geometry
Measure of elongation. Elongation is commonly parametrized using a traceless matrix E. This matrix, the magnitude and the angle of elongation with respect to x-axis can be calculated as follows: For a tissue consisting of N c cells, and with ε 1 = N −1 c Nc i=1 ε i,1 , we get for the average elongation: In terms of the elongation index E, we plot in Fig. (Fb), the angle between y-axis and the steady-state polarization vector |Φ p − Φ e | (in degrees), versus different values of E, for a system with the same initial condition and same initial lattice (i.e. before stretching), but elongated along y-axis. Here, Φ p is defined as the angle of the axis between polarization and the x-axis, and over the domain of Φ p ∈ [0, π). Since Φ e = π/2 for elongation in the y direction, we have |Φ p −Φ e | ∈ [0, π/2]. Of course for E = 0, the axis of elongation is not defined, yet we measure it with respect to y-axis. It is important to note that this figure is only an example where the polarization for E = 0, happens to make a small angle with yaxis (25 • ). It is clear that depending on the geometry of the lattice and the initial condition, the polarization can be perpendicular to y-axis, even for E = 0. Therefore among different simulations, we chose one with a relatively large effect of elongation, so that the onset of perpendicular polarity becomes easily distinguishable.  Illustrations of type I, II, and III mutants. The layout of the table is the same as that in the Main Text, and the red arrows show the directions of distortion with respect to the wild-type polarity. This table facilitates a more detailed comparison of the autonomous effects that were absent in Fig. (6) of the Main Text. In particular note the differences between the polarities within the putative dsh − and pk − clones, that were induced by small (α, β) and small λ, respectively.
Elongation parameter E, measures the deviation from equilateral hexagons where E = 0. In ordered tissues, there are two possible choices for elongation axis: parallel to a pair of parallel edges, Fig. (Fa1); and perpendicular to a pair of edges; Figs. (Fa2,Fa3). While (a1) and (a2) are polarized perpendicular to the axis of elongation, (a3) exhibits parallel polarization. The latter, however, is destabilized by NLCI inhibiting localization of unlike complexes on adjacent junctions. In disordered tissues, elongation is a mixture of (a1) and (a2), both of which give rise to perpendicular polarization.
In anisotropic cells, the effective values of cooperative interactions are larger for the long junctions. Here, without deriving explicit expressions for the effective parameters in elongated system, we only argue that the longer junctions acquire larger coefficients α, β. As in the case of one dimension, searching for mean-field (MF) solutions, we assume that in the steady state the concentrations of bound proteins are translationally invariant along the three main axes of the lattice. In a self-consistent approach, the effective α, β are dependent on the concentrations of dimers on other edges. Therefore, assuming the system has reached its steady state, we can write the cooperative interactions as functions of the concentrations of proteins on all the edges, weighted by the geometrical factors originating from nonlocal interactions. For small ranges of NLCI, λ 0.7, the effect of other edges are negligible and only the self-interaction of each edge is to be taken into account. Intuitively this is because the self-interaction is a monotonically increasing function of the edge length. Therefore longer edges with NLCI, have larger values of α, β. Upon increasing λ, the interactions between all pairs of edges increase. However the qualitative behavior of effective α, β's for different edges does not change. With this in mind, and using the symmetry arguments between junctions with equal lengths, one can see that the cartoons in Figs. (Fa1), and (Fa2, Fa3) correspond to, i.e. α 1 > α 2 , α 3 and α 1 = α 2 > α 3 , respectively. Using the analysis sketched in Sec. (2.2), it is easy to see that in such systems the junctions with larger α, β, are more unstable towards polarization transition (remember the critical value is inversely related to α, β). Following the above discussion on different coefficients of cooperative self-interactions α s , and the equations derived in Sec. (2.2), we calculate these coefficients as functions of L/ 0 , for different λ's. Fig. (Fb) shows the numerical results for 0.2 ≤ λ/ 0 ≤ 1. In order to discern the effect of the elongation from that of nonlocal interactions, we normalize all the curves by the coefficients calculated for L = 0 , such that α s (1) = 1 for all λ's. The coefficient of self-interaction is a monotonically increasing function of the length of junction; hence polarity emerges perpendicular to the longer junctions, as depicted in Figs. (Fa1) and (Fa2).

Mutants and The Corresponding Phenotypes
In the Main Text we presented the schematics of phenotypes of three types of mutants. The schematics in the Main Text Fig. (9) only show the non-autonomous effects of each clone. Here we show the actual polarization field obtained from the simulations, in which the autonomous effects within the clones are also illustrated.
Type-I Mutants: Lack of membrane proteins. Type I lacks one or both of the membrane proteins, namely f 0 , g 0 . The non-autonomy is evident in all cases. In the first row of Fig. (G), non-autonomy is extended to multiple layers of cells from the clone boundaries. In the second row where the background lacks Vang as well, the non-autonomy is mostly limited to a single layer of surrounding cells.
Type-II Mutants: Lack of cytoplasmic interactions. The role of cytoplasmic proteins is deduced by comparing the in silico and in vivo phenotypes. We concluded that Dsh is mostly responsible for the local cooperative interactions with some contributions to the nonlocal part. Pk on the other hand is mainly involved in nonlocal interactions, through cytoplasmic diffusion. It can be seen in bottom left panel of Fig. (G), that in agreement with experiments, e.g. [4], the phenotypes show almost no non-autonomy. Their autonomous effects, however, are distinct. Lacking Dsh removes the cell polarity within the clone, whereas absence of Pk has only minor effects on the polarity of the clone.
Type-III Mutants: Disordered geometry. We tested the predictions of our model in the case of geometric irregularities. In Ref. [9], Ma, et.al. suggested that clones with enhanced geometric disorder-induced by under-expression of PTENobstruct the propagation of polarization field. This effect appears in tissues where global cues are absent. Polarity was observed to be retained when the global cue was added. As can be seen in the two right panels of the last row in Fig. (G), the system with lack of global cue shows swirling pattern of polarization. The wild-type polarity, however, reappears upon adding the gradient cue.