The Role of Intracellular Interactions in the Collective Polarization of Tissues and its Interplay with Cellular Geometry
Shahriar Shadkhoo Madhav Mani


Planar cell polarity (PCP), the coherent in-plane polarization of a tissue on multicellular length scales, 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 mechanisms, how the two produce the collective polarization remains an active area of investigation. Exploring a generalized reaction-diffusion model, we study the role of intracellular interactions in the large-scale spatial coherence of cell polarities, and scrutinize the stages through which polarity emerges from subcellular to tissue-wide scales, as well as their dependence on intra-cellular and inter-cellular interactions. We demonstrate that nonlocal cytoplasmic interactions are necessary and sufficient for the robust long-range polarization of tissues—even in the absence of global cues—and are also essential to the faithful detection of weak directional signals. Non-locality of intracellular interactions makes a geometric readout a possibility, namely signatures of geometrical information in tissue polarity are expected to become manifest. Therefore, we investigate the deleterious effects of geometric disorder, and determine conditions on the magnitude and spatial extension of cytoplasmic interactions that guarantee the stability of polarization. These conditions get progressively more stringent upon increasing the geometric disorder, namely narrower ranges of parameters guarantee long-range polarization. Another situation where the role of geometrical information might be evident is elongated tissues. Strikingly, our model recapitulates an observed influence of tissue elongation on the orientation of polarity. Eventually, we introduce three major classes of in silico mutants: lack of membrane proteins, lack of cytoplasmic proteins, and locally enhanced geometrical irregularities. In order to examine the predictability of our model, we adopt core-PCP as a model pathway, and interpret the model parameters and variables accordingly. Through comparing and matching the in silico and in vivo phenotypes, we interpret the effective functionalities of various core-PCP components. The success of our model in reproducing the in vivo phenotypes, helps us shed light, in particular, on the roles of the cytoplasmic proteins Prickle and Dishevelled in cell-cell communication, and make falsifiable predictions regarding the cooperation of cytoplasmic and membrane proteins in long-range polarization, faithful readout of directional cues, as well as geometrical information. 














Kavli Institute for Theoretical Physics, University of California, Santa Barbara, CA 93106, USA
Department of Engineering Sciences and Applied Mathematics, Northwestern University, Evanston, IL 60208, USA









The Role of Intracellular Interactions in the Collective Polarization of Tissues
 and its Interplay with Cellular Geometry
    Madhav Mani
    December 30, 2023
==========================================================================================================================













§ INTRODUCTION

The patterning of an organism requires the coupling of cellular states across multicellular scales. As such, the coordination of cellular processes on these longer length scales are crucial to the emergent phenotype of an organism and requires the faithful transduction of directional information across tissues. Planar cell polarity (PCP) is understood to be one of the core mechanisms responsible for such tissue-wide signaling <cit.>. At the cellular level, polarity is defined as the asymmetric localization of membrane associated proteins on the apicolateral cell junctions, which is prompted and reinforced by cytoplasmic interactions and feedback loops <cit.>. Long-range polarization arises as a consequence of juxtacrine signaling through which adjacent cells align their polarities. How these two scales are connected, and the respective roles of molecular components in establishing long-range planar polarity is the focus of our study <cit.>. These intra- and intercellular interactions are largely carried out via two PCP pathways: “core-PCP” and “Fat/Dachsous” <cit.>, each of which involves several interacting (trans-)membrane and cytoplasmic proteins. Throughout this paper, we adopt core-PCP as the reference pathway, according to which we interpret the results. However, the construction of the model is based on phenomenology, general arguments and physical assumptions, and its structure is independent of the molecular details of specific PCP pathways. 


Molecular ingredients. Generically, PCP pathways consist of membrane-bound and cytoplasmic proteins. Core-PCP pathway consists of six known proteins: the membrane proteins Frizzled (Fz) and Van Gogh (Vang), which form complexes by binding to the transmembrane protein Flamingo (Fmi) on opposite sides of cell-cell junctions, and form the asymmetric heterodimers Fz:Fmi and Fmi:Vang. In addition to the above proteins, there exist three cytoplasmic proteins, believed to mediate intracellular interactions: Disheveled (Dsh) and Diego (Dgo) bind to Fz, and Prickle (Pk) binds to Vang <cit.>. Although the presence of Dsh, Dgo, and Pk is found to be unnecessary for intercellular interactions, they facilitate the segregation of Fz and Vang to opposite sides of a cell, hence their absence impairs long-range polarization <cit.>. Furthermore, Fz:Dsh and Vang:Pk are believed to mutually suppress the activities of each other <cit.>. One of the main goals of this paper is to address the significance of such cytoplasmic proteins in stabilizing cellular polarization, as well as their interplay with cell-cell signaling in establishing large-scale polarization. 


Global cues. Although long-range polarization emerges spontaneously through cell-cell interactions, external cues are believed to be necessary for fixing the direction of polarization <cit.>. The graded distribution of regulatory factors across a tissue, i.e. morphogens <cit.>, mechanical signals <cit.>, and geometrical cues <cit.>, are speculated to provide such global orientational signals. Elongation in particular, has been observed to induce polarization. At a subcellular level, the polarization of microtubules and vesicle trafficking are also proposed to be acting as a bias to determine PCP orientation <cit.>. In the mammalian cochlea and skin, polarization is perpendicular to the elongation axis <cit.>. In mice, elongation along the medial-lateral axis has been suggested to orient the polarization along the anterior-posterior axis <cit.>. The present model aims to bring a mechanistic understanding to the potential role of cell geometry in the polarization of a patch of cells. 


Geometry and timescales. Several studies (e.g. <cit.>) have proposed underlying physical mechanisms of PCP in ordered and isotropic (i.e. non-elongated) systems. Establishment of long-range polarization during the course of development, however, can precede the formation of an ordered lattice, e.g. margin-oriented polarity in the prepupal Drosophila wing <cit.>. In particular, it is suggested that geometrical irregularities cause disruption in the polarization field <cit.>. Therefore, it is important to understand how PCP manages to propagate through disordered as well as elongated tissues. FRAP measurements of PCP proteins suggest turnover timescales to be much shorter than the timescales of cell rearrangements, justifying the study of PCP kinetics on a static tissue <cit.>.


Modeling Planar Cell Polarity. Quantitative modeling of PCP and the underlying mechanisms has been of great interest to computational biologists and biophysicists. Several classes of model have been proposed each focusing on certain aspects, from subcellular molecular circuitry in charge of single-cell polarity, to intercellular communications that give rise to propagation of polarization over large distances. The coupling between the two modules has been a key question. While individual molecular components and their roles vary among different PCP pathways, networks of these components seem to share principal functionalities. In addition to the mechanisms of interaction among different components of a PCP pathway, detection of global cues is of great importance, as the direction of polarity is eventually set by such cues. As mentioned above global cues of various kinds have been observed in experiments. While the coupling of PCP proteins to the cues of chemical origins (morphogens) is more conceivable, the readout mechanisms of other cues such as geometrical and mechanical are not easy to decipher. Therefore, an important question is how each type of these cues influence the polarity. Some models have proposed mechanism through which cells are individually polarized by gradient cues <cit.>. Others have considered scenarios where rotational symmetry breaks spontaneously and collective polarization emerges even in the absence of global cues <cit.>. The tissue polarity is then rotated in the right direction, through coupling to the global cues. The latter mechanism enjoys high sensitivity and faithful detection of global cues, and is robust against random misreading of the orientational information by individual cells. Since the origin of global cues is not always easily determined, knockdown of a certain gene does not necessarily unravel the underlying mechanism. 



A popular class of mathematical models begin with intracellular interactions, which along with cell-cell couplings—that can involve the same components—give rise to long-range alignment of tissue polarity <cit.>. Reaction-diffusion (RD) equations of different variations constitute the basis of the majority of these models. The characteristics of effective interactions between components of the PCP pathway are encoded in the RD equations. Generically, two types of protein complexes are considered that interact and localize asymmetrically on cell membranes. The interactions are assumed to activate/inhibit the like/unlike complexes, and can be of cytosolic and/or cross-membrane type. Several studied have focused on the necessary conditions, for the cytoplasmic interactions to establish long-range polarization. Well-established facts suggest that non-locality is among the essential features of these interactions <cit.>. These nonlocal interactions are mediated through diffusive cytosolic proteins or complexes thereof. It is suggested in Ref. <cit.>, that nonlocal inhibition between opposite complexes are necessary and sufficient to establish long-range polarity in ordered tissues. This is fundamentally akin to the well-known local-activation–global-inhibition mechanism, which results in the accumulation of similar components on one side and the repulsion of the other components to the opposite side <cit.>. In another study, Abley et.al. <cit.> considered various possibilities for local and nonlocal interactions between like and unlike components,  as well as different types of cell-cell couplings (i.e. direct and indirect). They concluded that intracellular interactions are crucial to the segregation of unlike complexes to the opposite compartments of cells. However, they claim that in spite of cytoplasmic partitioning, global cues are needed in order for the correlation length to exceed a few cell diameters; in the absence of cues, swirls (vortex-like) patterns of polarity appear in the steady state. 


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. Below we will see that 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. 


Outline and Results. The objectives of this paper are threefold. Starting with a generalized reaction-diffusion equation for intracellular dynamics of proteins, we address the role of intracellular interaction in establishing tissue-wide alignment of polarization in tissues with disordered and/or elongated geometries. First we demonstrate that nonlocal interactions of both kind, namely stabilizing/destabilizing, promote the cellular segregation of unlike complexes and are crucial to the global alignment of polarity. By varying the associated length scales to stabilizing/destabilizing cytoplasmic interactions, we investigate and highlight the role of stabilizing interactions, and show that non-locality of the latter enhances the correlations of the PCP field in the presence of geometrical disorder. The role of geometrical disorder is studied, and it is found that the minimum length-scale of cytoplasmic interactions that stabilizes the long-range polarity increases for larger geometrical disorder of the tissue. Furthermore, we demonstrate that in elongated tissues, nonlocal interactions stabilize the polarization axis perpendicular to that of elongation. Finally, to further signify the necessity of nonlocal interactions, and facilitate a conversation between theory and experiment, we study three classes of in silico mutants and identify phenotypic similarities with experimental observations. 




Before introducing the formalism, we shall disambiguate some terminology: “edge” and “junction”, are used interchangeably throughout the paper, depending on the context emphasizing on the geometrical/mathematical or biological aspects of the problem, respectively. Therefore, one can think of an edge as a cell-cell junction. Additionally, we use for the interactions between like-like and like-unlike complexes, either activation/inhibition or stabilizing/destabilizing interactions, respectively; again depending on the context. Finally, “defect” has been used in two different contexts: (1) defects in the polygonal network of cells, that appear as a result of geometrical quenched disorder; and (2) topological defects in the polarization field, an example of which is swirls. The two appear in separate contexts, and should not cause any ambiguity; also the latter is usually preceded by “topological".



§ MODEL AND FORMALISM

In line with the known molecular interactions in core-PCP, we introduce a set of reaction-diffusion (RD) equations that govern the binding-unbinding dynamics of transmembrane complexes. Each cell is assumed to contain a finite pool of proteins Fz and Vang, which in their active state localize on the opposite sides of cell-cell junctions, and bind to a cross-junctional Fmi-Fmi homodimer and form asymmetric complexes Fz:Fmi-Fmi:Vang. The linear densities of total (bound plus free) Fz and Vang, are denoted by f_0 and g_0, which in the absence of global cues, are assumed to be identical for all cells across the tissue. Given that the transcriptional timescales typically far exceed the kinetic timescales of protein-protein interactions, f_0 and g_0 are treated as time-independent <cit.>.



It has been shown that Fz can in principle bind to Fmi-Fmi homodimers and make a Fz:Fmi-Fmi complex without a Vang molecule on the other side of the junction <cit.>. Therefore a thorough analysis requires separate RD equations for Fz and Vang. Ignoring this effect for simplicity, we define the unit of junctional polarity as a Fz:Fmi-Fmi:Vang complex for now; hence treating Fz and Vang on the same footing. We discuss the effect of this asymmetry in Sec. (<ref>). For notational convenience, we denote such complexes by F-G, where F ≡ Fz:Fmi and G ≡ Fmi:Vang. At any point  on a junction shared by cells i and j, the concentrations of bound [F_i-G_j] and [G_i-F_j], are denoted by u_ij() and v_ij(), respectively. Consistent with this notation, we have u_ij() = v_ji(). The key assumption in this model is that within each cell, the formation of a dimer at a point  is nonlocally enhanced by like dimers, and its dissociation is enhanced nonlocally by opposite dimers. In short, bound Fz nonlocally stabilizes Fz and destabilizes Vang; and vice versa. This represents an positive feedback between Fz and Vang in the adjacent cells: promoting (inhibiting) Fz in the same cell indirectly brings more (less) Vang to the other side of the junction. The nonlocal interactions are diffusively mediated through the cytoplasmic proteins Dsh, Dgo, and Pk, and/or their associated feedback loops. Figure (<ref>) is an illustration of the relevant cytoplasmic interactions.

The RD equations governing the binding/unbinding dynamics read:

    du_ij()/dt   =  κ f_i^ubdg_j^ubd(1+α∑_{k}_i∫_i∩ k d'𝒦_uu(-')u_ik('))
        - γ u_ij()(1+β∑_{k}_i∫_i∩ k d' 𝒦_uv(-')v_ik('))
        + η(, t) + ℳ_ij() exp(-t/τ_0).

The notations adopted here are as follows: ∑_{k}_i is the summation over all neighbors {k} of cell i; and ∫_i∩ kd', represents integration over the junction shared by cells i,k. In the above equation, the first and second terms on the r.h.s. correspond to the formation and dissociation rates, respectively. κ is the bare rate of formation. Based on the assumption that the diffusion of unbound proteins is rapid we posit that the pool of free proteins are uniformly accessible to the perimeter of a cell. Therefore, the formation rate of u_ij() is proportional to the densities of unbound, i.e. cytoplasmic Fz in cell i, f^ubd_i, as well as that of Vang in cell j, g^ubd_j. In the second term, the dissociation rate is proportional to local concentration of the dimer itself, with the bare rate γ. The formation/dissociation processes are amplified by like/unlike dimers, respectively, through the nonlocal terms, α∑_{k}_i∫_i∩ k d'𝒦_uu(-')u_ik(') and β∑_{k}_i∫_i∩ k d' 𝒦_uv(-')v_ik('), that characterize cooperative formation and dissociation. The functional form of the kernels 𝒦_uu(-') and 𝒦_uv(-') and their coefficients α and β, are introduced below. Finally, the last term η(,t) is a stochastic Gaussian white noise: ⟨η(,t)⟩= 0, and ⟨η(,t)η(',t')⟩ = η_0^2δ(-')δ(t-t'), which arises from the molecular noise of chemical reactions and stochasticity in the upstream signaling pathways. The former, modeled as a Poisson process, is speculated to be the dominant source of noise <cit.>, with its magnitude scaling as η_0 ∼ 1/√(N_mol.), where N_mol. is the number of molecules per area of the lateral interfaces. More precisely, the number of participating molecules is N_mol.u/f_0, where u is approximately the average value of u, and f_0=1 is the unit of concentration. Using the variance of the number of reactions per unit time, that is given by the r.h.s. of Eq. (<ref>), the noise level, is estimated to be of order η_0≃ 0.01 – 0.1, for N_mol.≃ 1 – 5 × 10^3, the approximate number of Frizzled molecules in the Drosophila wing <cit.>. The last term models global cues with local magnitude ℳ_ij() and time scale τ_0, over which the corresponding gene is expressed. For simplicity, we assume that the cues only couple to one of the complexes, say F. Global cues are discussed in Sec. (<ref>B), separately. Unless mentioned otherwise, the cues are assumed to be zero. 

The densities of unbound Fz and Vang are obtained by subtracting the densities of bound proteins from the total densities.

    f_i^ubd = f_0 - 1/𝒞_i∫__i d u_i(),

and a similar relation for g_i^ubd. Here 𝒞_i is the perimeter of cell i, and u_i() represents a u-complex sitting at point  on the perimeter of cell i, i.e. (F_i-G); we dropped the index j of the neighbors, because we only care about the total bound F, not the specific cell with which this F-G complex is shared. The same notation is applied to v_i():
The cell polarity with respect to the centroid of cell i at 𝐑_i, is defined as:

    _i = 1/2 ∫__id -𝐑_i/|-𝐑_i| (u_i() - v_i()).
 
All the quantities are expressed in units of f_0 = γ = ℓ_0 = 1, where ℓ_0 is the average length of the cell-cell junctions. In the following we will see that in certain situations, an alternative, but related, definition of polarization simplifies our analyses. We define two junctional variables: the cross-junctional polarity p_ij() = u_ij() - u_ji() = u_ij() - v_ij(), and the total concentration of localized complexes, s_ij() = u_ij() + v_ij(). Given p_ij,s_ij, the cell polarity can be readily extracted; see SI. (2).


Nonlocal Interactions. The kernels 𝒦_uu(-') and 𝒦_uv(-'), identify the functional form of the interactions between like and unlike complexes, respectively, and are taken to be exponentially decaying: 𝒦_uu() = N_uu^-1exp(-||/λ_uu) and 𝒦_uv() = N_uv^-1exp(-||/λ_uv), where λ_uu, λ_uv are the characteristic length scales of u-u and u-v interactions, respectively. The prefactors N_uu  and N_uv are normalization factors, to be determined shortly. Before calculating the normalization factors, we shall make a detour to explain an approximation we used that greatly reduces the computational cost simulations. The full integro-differential equation we introduced in Eq. (<ref>) couples every pairs of points on the perimeter of a cell. Assuming uniform distributions of localized proteins on each junction, we simplify the equation to obtain effective equations for edges, where the integration over the kernels are replaced by matrix products, see Eqs. (<ref>) and (<ref>). To avoid confusions with cell indices, we use Greek letter μ to label a single edge. First, we introduce junctional concentrations:

    u_μ = ∫_μ d u_μ(),

where the integral is taken over the concentration of complex u('), on all the points ' along junction μ. Same definition applies to v_μ. With this approximation, one can recast Eq. (<ref>) into:

    du_μ/dt =   κ f_i^ubdg_j^ubd(1+α∑_ν𝒦_uu^μνu_ν) - γ u_μ(1+β∑_ν𝒦_uv^μνv_ν)
        + η_μ(t) + ℳ_μ exp(-t/τ_0).

In the above equation, η_μ and ℳ_μ are the noise and the global cue averaged over the length of edge μ. The cooperative interactions are now reduced to simple matrix products of a matrix 𝒦_uu and a vector u, wherein the matrix elements 𝒦^μν_uu are purely geometrical constants obtained using the following relation:

    𝒦_uu^μν = ∬_μ,νd_ν d_μ 𝒦_uu(_ν - _μ).

A similar relation holds for destabilizing interactions 𝒦_uv^μν. The edge-edge coupling coefficients 𝒦_uu^μν, are only a function of the cell geometries; once calculated for a given tissue, the full matrix can be used throughout the course of integrating the dynamics of protein concentrations; hence reducing computational runtime by orders of magnitude. 

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_uu and N_uv, such that the self-interaction of an edge of length ℓ_μ equals α, namely: α𝒦_uu^μμ = α, or 𝒦_uu^μμ = 1. The same relation holds for β and 𝒦_uv^μμ. This choice of normalization, by fixing the edge self-interactions, ensures that the observed behavior upon changing λ's is purely due to nonlocal edge-edge coupling and 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 a hypothetical edge with average length of all edges, ℓ_0 ≡ℓ_μ. Using the definition of kernels, for the “average” edge 𝒦_uu^avg. from Eq. (<ref>), we obtain:

    N(λ)    = ∬ d d' exp(-|' - |/λ) 
       =  2λ^2(e^-ℓ_0λ^-1 + ℓ_0λ^-1 - 1).

Here,  and ' move along the length of a single edge with length ℓ_0. We omitted subscripts “uu” and “uv” for simplicity. 


Limit of Strictly Local Cytoplasmic Interactions (SLCI). In the limit of small λ/ℓ_μ→ 0, we get, 𝒦_μν = δ_μν, where δ_μν is the Kronecker delta. Thus in the SLCI limit, the equations read,

    du_μ/dt =     f^ubd_ig^ubd_j(1+α u_μ) -γ u_μ(1+β v_μ) + η_μ(t).


Interpretation of The Model Parameters. Besides the geometrical and stochastic disorder parameter, our model consists of four independent dimensionless parameters: κ/γ , α, β , λ/ℓ_0. Here we try to make connections between the effective roles of PCP components and these parameters. First we note that the cooperative interactions represent the effective couplings of two complexes, not the interactions between Fz's or Vang's proteins separately. In other words, α integrates the stabilizing interactions of Fz by Fz, and Vang by Vang; similarly for β, i.e. destabilizing effects between Fz and Vang. The functional form of the kernels can be interpreted as interactions mediated by diffusing cytoplasmic proteins with diffusion constant D and the degradation rate τ_c^-1, such that γτ_c≪ 1; thus λ = √(Dτ_c) . The diffusion timescale of cytoplasmic proteins for D≃ 0.5 (μ^2/s), and ℓ_0≃ 5 (μ), is of the order of ≃10 (min), much shorter than polarization dynamics which occurs on timescales of a few hours. On the other hand, we know that Dsh and Pk, both promote localizing similar complexes, and suppress the opposite ones. In particular, Dsh locally promotes the localization of Fz, whereas Pk destabilizes that in a nonlocal fashion by inhibiting the membrane localization of Dsh and antagonizing Fz accumulation, hence indirectly stabilizing Vang localization <cit.>. Given that the diffusion lengths of these proteins are independent of their role (stabilizing/destabilizing), in the main text we assume λ_uu=λ_uv=λ. The cases of λ_uuλ_uv are elaborated on in SI. (4.4). Coefficients α, β, which parametrize the strength and length scale of the cooperative interactions, depend on the concentrations of bound cytoplasmic proteins Dsh and Pk, which in turn depend on the abundance of total cytoplasmic proteins, and their binding affinities with Fz and Vang; see SI. (1). In summary, both Dsh and Pk contribute to the magnitudes and length scales of the cooperative interactions α,β, and λ. In order to shed light on the principal mission of Dsh and Pk in the long-range correlation of polarization, in Sec. (<ref>) we compare the in vivo phenotypes of dsh^- and pk^- with the in silico phenotypes of α,β and λ. Finally, the formation of the polar complex Fz:Fmi-Fmi:Vang, among other factors, is contingent on the presence of Fmi and formation of Fmi-Fmi dimers. Therefore, we expect κ/γ, the ratio of formation and dissociation rates of the complexes, to be an increasing function of the binding affinities of Fz:Fmi as well as Fmi:Vang.


Correlation Function. The correlation function of polarization is defined as a measure of alignment of polarity. In order to investigate the temporal behavior of the spatial extension of the alignment, we define the equal-time correlation functions as follows. Consider an arbitrary cell at _i, and a vector  connecting it to another cell at _j = _i +. With no further assumption we can define a correlation function, that is dependent on the distance r = || and the relative angle of  and (_i); we call it θ_,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,

    S(r,θ_,p ;t) = N_c^-1∑_i(_i ;t)·(_i +  ;t).

The above quantity calculates the average conditional probability that the dipole at point _j = _i +, takes on the value and direction of (_j), should the polarity at point _i be (_i). For θ_,p = 0, π, and θ_,p = ±π, we get parallel and perpendicular correlations, respectively. In spite of this angular dependence, averaging the above correlation function over θ_,p∈ [0 , 2π), returns a weighted averaged of correlation function as a function of r=||. This function is bounded by the the longitudinal and transverse correlations from above and below, respectively. Intuitive arguments are provided in this regard, in Sec. (<ref>). Thus, we define radial correlation function as follows.

    S( ;t) = N_c^-1∑_i∫_0^2πdθ_,p/2π (_i ;t)·(_i +  ;t).

Correlation length can be obtained from the above equation:

    ξ(t) = ∫_0^R_c dr r S(r;t)/∫_0^R_c dr  S(r;t).

Here, R_c = 40 (cell diameter) is the farthest cell with which the correlation is calculated, for a given reference cell. A perfectly correlated polarization field, returns: ξ = R_c/2. 

A simpler measure for the global orientational order is 𝒪(t) = P(t)/Q(t), where P(t) = |⟨(t)⟩|, and Q(t) = ⟨ |(t)|⟩, in which ⟨∙⟩ denotes spatial average. Thus 𝒪(t) saturates to unity for perfect alignment. However, we will see below that Q(t) can be misleading in interpreting the segregation of proteins, and in not a reliable  







Geometrical Disorder. The edge lengths are ℓ_i = ℓ_0(1 + ϵ_i), where -ϵ_0≤ϵ_i ≤ + ϵ_0 with uniform distribution, and ⟨ϵ_iϵ_j⟩ = ϵ_0^2δ_ij/6. In 2D, lattice defects, i.e. non-hexagonal cells appear above a certain level of quenched disorder corresponding to ϵ_0 ≃ 0.25. In disordered cases, we use ϵ_0 ≃ 0.5 and density of defects n_d ≃ 0.6, corresponding roughly to the statistics of larval and prepupal Drosophila wing <cit.>. In order to gain more intuition on the geometrical properties of the network of cells, as a function of ϵ_0, we performed statistical analysis on synthetic tissues constructed from Voronoi tessellation of randomly positioned cell centroids. The results of these analyses are shown in SI. (3) and SI. Fig. (1).




One-dimensional Arrays of Cells. First, we briefly review the one-dimensional case as it is more amenable to analytical treatment, and captures some features of 2D systems. In particular when the sixfold symmetry of a hexagonal lattice is broken—either spontaneously or explicitly—down into a twofold (+/-) symmetry, the system behaves much like a one-dimensional array. Examples of this scenario in 2D include, partitioned cells in which the two complexes F and G are segregated to the opposite sides of the cells, and elongated cells where the sixfold rotational symmetry is broken. These cases are discussed in full details, in sections (<ref>) and (<ref>), respectively. 

In one dimension, cells are juxtaposed on a line. A pair of adjacent cells i and i+1, share a junction of length ℓ_i,i+1, which hosts membrane-bound complexes F and G. A schematic of a 1D array is depicted in Fig. (<ref>a). At the mean-field (MF) level for an ordered system in 1D, edges separate cells i and i + 1, i.e. for all edges ∀{i,i+1} : ℓ_i,i+1 = ℓ_0, we get u_i,i+1 = u and v_i,i+1 = v, hence f^ubd_i = f^ubd and g^ubd_i = g^ubd. We switch variables from (u,v) to (p,s), where p_i,i+1 = u_i,i+1 - v_i,i+1, and s_i,i+1 = u_i,i+1 + v_i,i+1. The steady-state solution, defined as ∀ i  :  p_i,i+1 = p, and ∀ i  :  s_i,i+1 = s, exhibits a bifurcation from unpolarized to polarized state, as the control parameter g_0/f_0, is increased above a critical value <cit.>. The MF polarization reads:


    s =(f_0 + g_0) - √((f_0 - g_0)^2 + 4γ/κα),


    p = ±(s^2 - 4/αβ)^1/2.


From the second equation, the bifurcation takes place at s^* = 2/√(αβ). In terms of actual control parameter g_0 we get,

    g_0^* = γ/κα/1-√(1/αβ)+√(1/αβ).
 
This result indicates the divergence of the critical value s^* (or g_0^*), for αβ→ 0, implying that the emergence of polarization requires cooperative interactions. Numerical solutions are presented in Fig. (<ref>b), for a system with constant F_0, G_0, the total number of Fz and Vang per cell. In ordered systems, one can choose F_0 , G_0 or equivalently f_0 , g_0 to be constant across the tissue. In disordered systems, the two conditions are not equivalent. Only in 1D and in order to emphasize a possible effect of disorder we assume F_0 , G_0 to be constant. As such, the concentrations available to the junctions of a cells, are no longer equal to those in the adjacent cells. In a disordered system with ℓ_i,i+1 = ℓ_0 + ϵ_i,i+1, the concentrations f_0,i = 2F_0/(ℓ_i-1,i + ℓ_i , i+1) and g_0,i = 2G_0/(ℓ_i-1,i + ℓ_i , i+1) are thus randomized. In SI. (4.1), we show that the relevant quantities for edge polarization are g_0,i+1/f_0,i and g_0,i/f_0,i+1, both of which are nonuniform as a result of length disorder. The ratio g_0,i/f_0,i which can be interpreted as “local critical point”, is thus randomized as well. Therefore in disordered systems (ϵ_0 0), the collective singular behavior of the ordered systems, at the well-defined critical point, is smeared out, and the second-order transition is replaced by a smooth cross-over (Fig. (<ref>b)).


Numerical solutions, Fig. (<ref>c) and (<ref>d), 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 is defined as δ p_0/p_0, the normalized magnitude of spatial fluctuations of initial polarity. Thus, small and large bias limits 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 are robust and largely determined by the local geometry (disorder) of the array. Therefore we observe that already in 1D, the quenched disorder imposes undesirable solutions, impairing the faithful transduction of directional information through PCP signaling. As we will see in the following, the situation gets only more complicated in two dimensions, even for ordered tissues. One of the main goals of this paper is to find mechanisms to circumvent these issues associated with two-dimensional tissues. 



§ INTRACELLULAR INTERACTIONS:
 LOCAL OR NONLOCAL?

The systems in one and two dimensions show inherently different behavior. In 1D, the proteins have only two junctions at which the can 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 a great portion 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 qualitative changes upon varying the above parameters is discussed as well. 

The numerical solutions are presented below, but first let us attempt to gain some insight using analytical MF analysis for ordered tissues. We define the MF approximation in 2D as uniform distribution of f_i^ubdg_j^ubd and f_j^ubdg_i^ubd for all junctions {ij}. The validity of this approximation can be justified by noting the diffusive nature of p, s dynamics, i.e. the amount of membrane-bound proteins (see Ref. <cit.>), which in turn is concomitant with the diffusive dynamics of the concentrations of free cytoplasmic proteins. We checked this assumption numerically. The value of f_i^ubdg_j^ubd, normalized by the mean value, is plotted for all edges for initial and final states in SI. Fig. (3). In steady state, the standard deviation of the distribution normalized by the mean, is very much independent of initial condition as well as the model parameters, and remains below ≃ 0.05, within the ranges explored in this paper, for both SLCI and NLCI regimes. Owing to the sixfold symmetry of equilateral cells, we expect the steady-state magnitudes of p, s to be the same on all edges. Thus, three edges will be carrying inward, and the other three outward dipoles; otherwise the MF criterion is violated. 



The simplest of the MF solutions, here referred to as the trivial solutions, preserve translational invariance along each of the main three axes of a hexagonal lattice separately. There exist two types of trivial MF solutions: polarized and unpolarized. The two types can be seen in Figs. (<ref>a1) and (<ref>a2): six configurations with nonzero net polarization, obtained by 60 (deg) rotations of the configuration (a1); and two unpolarized states, obtained by flipping all the dipoles in (a2). There exist other types of solutions meeting the MF criterion, in which translational invariance of cell polarity is not preserved, but f_i^ubdg_j^ubd is uniform for all edges; see Fig. (<ref>b). We call these nontrivial MF solutions. The cellular polarities are randomly oriented and long-range correlation is absent in the nontrivial MF solutions. Straightforwardly, as one can see in Fig. (<ref>b), the number of possible configurations of this type hugely outnumbers the eight trivial solutions. Analytical investigation of such solutions are beyond the scope of this study. Important intuitive arguments are provided regarding this class of solutions, and their behavior is elaborated on in SI. (4.3). We shall emphasize that the above classification of MF solutions is independent of the locality or lack thereof of the cytoplasmic interactions, and is purely based upon MF criterion, and symmetry arguments. However, the class of solutions which is realized in a system, is strongly dependent on the cytoplasmic interactions. We will see below that systems with strictly local interactions are incapable of segregating the proteins to the opposite sides of cells, whereas nonlocal interactions mediate repulsive interaction between the unlike proteins that amounts to partitioning the cells. Therefore, the nontrivial MF solutions are realized in the SLCI regime, and the polarized trivial solutions are the result of NLCI regime. Unpolarized trivial solutions Fig. (<ref>a2), although might appear by accident in SLCI systems, are destabilized by nonlocal interactions. 


Since the analysis of nontrivial solutions provides an intuitive argument as to why SLCI is insufficient to obtain long-range polarization, we highly encourage the reader to peruse SI. (4.3). The RD equations in 2D are precisely the same as those in 1D, except the pools of proteins Fz and Vang are shared between six junctions instead of two. Therefore, the MF concentrations in steady state, are identical to those of 1D case. The corresponding equations and solutions in 2D are derived in SI. (4.2). This can also be understood intuitively, by noting that in a polarized trivial MF state depicted in Fig. (<ref>a1), each cell is partitioned into a positive and a negative side. Each partition can be thought of as an “effective edge” in the 1D case. Although the lengths of these “effective edges" are different from those in 1D case, the concentrations are equal. Now, consider a hexagonal lattice; each edge carries the same |p|,s, equal to those found in 1D case. For the polarized trivial MF solution, the net cellular polarization equals p_cell=p_edge (1 - 2 cos(2π/3))=2√(s_edge^2-4/αβ) , where p_edge, s_edge are junctional polarity and sum of localized proteins, as were defined in the case of a 1D array of cells. In Sec. (<ref>), we saw that in ordered systems, minimum concentration of G above which polarized state is stable, i.e. g_0^*, increases with increasing αβ, as well as γ/κ (see Eq. (<ref>)). Using the units and parameters introduced in Sec. (<ref>), we get g_0^* ≃ 0.23, for the SLCI critical point in MF approximation and in ordered systems. For g_0 > g_0^* the polarized steady state, namely the polarized trivial MF solution, can be realized in special cases where the initial distribution of the proteins is not far from the steady state, and a small global cue assists with redistribution of proteins. Therefore, the efficacy of SLCI is strongly dependent on the initial condition. 






 §.§ Strictly Local Cytoplasmic Interactions (SLCI)

The regime of SLCI is defined as λ/ℓ_0→ 0, i.e. local cytoplasmic interactions. In Sec. (<ref>) we simplified Eq. (<ref>) in the SLCI limit. Previously we argued qualitatively, that due to the absence of “repulsive” interactions between unlike complexes in SLCI limit, the segregation is not accomplished properly. Moreover, using MF arguments we showed that there exist steady-state solutions where segregation is not enforced by cytoplasmic interactions. As such, starting from a random initial condition with no global cue, a system with SLCI will almost surely settle in a nontrivial fixed point, by only local redistributions of proteins on the same junction, which is possible through SLCI. An example of such configurations is depicted in Fig. (<ref>b). A natural system in which such solutions appear are mutants, in which mutation is induced by loss-of-function of cytoplasmic proteins. In such systems, the partial randomly oriented polarization is achieved through local activation/inhibition between the similar/opposite complexes (see Sec. (<ref>), for further discussion). 

Simulations of systems with random initial states and weak global cues testify to the lack of long-range order of polarity in SLCI regime. A generic steady state of such systems, the rose-plot of the angular distribution of dipoles, the time evolution of Q(t), P(t), their ratio 𝒪(t), as well as the correlation lengths ξ(t), are plotted in Figs. (<ref>a1), (<ref>b1), (<ref>c1), and (<ref>c2), respectively. To facilitate the comparison, Figs. (<ref>c1) and (<ref>c2) include corresponding quantities of other cases of study, that are discussed in the next two subsections. One important observation is the rapid dynamics of polarization in the case of SLCI, in particular P(t), acquiring its steady state value within a time scale of the order of t ≲ 5 γ^-1. This is exactly due to the huge basin of attraction of nontrivial fixed points in SLCI limit; any initial state is close to a fixed point, to which it is quickly attracted in the absence of cues. 





 §.§ Nonlocal Cytoplasmic Interactions (NLCI)

Nonlocal interactions, as introduced in Sec. (<ref>), are mediated by cytoplasmic proteins. Intuitively, NLCI facilitates the segregation of unlike complexes to the opposite sides of cells, by nonlocally activating the like, and inhibiting the unlike complexes formation; one can think of this is effectively “attracting” the like, and “repelling” the unlike complexes. Segregation makes the system behave more like a one dimensional lattice, by splitting each cell into two compartments. Indeed, the trivial MF solution of type (a1) in Fig. (<ref>), becomes increasingly more valid as the range of nonlocal interactions is increased up to an upper limit for λ to be found shortly. Segregations implies that adjacent edges of a cell prefer to carry the same polarities, and alternating between inward and outward is not favorable. As such, it is conceivable that configurations like the two trivial MF solutions with zero-net polarization in Fig. (<ref>, are destabilize by NLCI. The rest of this section discusses the results of our simulations. 

While in relatively ordered tissues with NLCI, and in the absence of orientational cues, the orientation of polarization is determined purely by chance, namely stochastic noise and initial conditions, highly disordered systems show robustness against such random factors, and 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. Below, where we discuss topological defects in polarity patters, we provide more evidence for this behavior. External cues of sufficiently large magnitudes, however, reorient the polarity towards the favored direction. See below for further discussion. A typical configuration of the steady states is illustrated in Fig. (<ref>a2). We find a range of interaction length scale, 0.15≲λ/ℓ_0≲ 0.7, for which the NLCI guarantees the long-range alignment of polarization, with the standard deviation of the dipoles' directions less than 30 (degrees). The directional correlation shows a peak at around λ/ℓ_0 ≃ 0.4-0.5. The aforementioned range of λ also depends on the disorder. Geometrical disorder 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. Not surprisingly, the angular correlation at a certain disorder is dependent also slightly lower in tissues with disorder on the geometrical disorder; see SI. Fig. (4). However, we believe that this is, in part, due to the simplification associated with uniform junctional distribution of proteins, in the absence of which the dipoles are have more rotational degree of freedom. Before addressing the functional range of λ, we note that for λ/ℓ_0≳ 0.7, all edges within a cell strongly couple to each other, hence preventing the segregation and making the system fragile to stochastic noise. Interestingly, this regime becomes relevant in one of the mutants discussed in Sec. (<ref>).

For the range of interest, the dynamics of Q and P, shown in Fig. (<ref>c1), imply that the spontaneous emergence (i.e. with no global cue) of collective polarization from an initially random distribution, consists of two distinct stages: (i) the segregation of PCP proteins within each cell, and saturation of the amplitude of polarity, accompanied by the formation of polarized local domains, which is followed by (ii) the subsequent coarsening and alignment of the domains across the tissue. The first and second stages are carried out mostly through intracellular and intercellular interactions, respectively. Below, we argue that the former, is indeed the key to long-range polarization. 


Segregation mechanisms in SLCI vs. NLCI. Through a comparison of the dynamics of SLCI and NLCI in Figs. (<ref>b1) and (<ref>b2), 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 we think of Q(t) as a measure of cellular segregation of proteins? A naïve guess would be that since Q_i(t) is oblivious of 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 accomplished consistently in all cells, namely some cells are highly segregated while others are not. 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:




    _i^F = ∫__id -𝐑_i/|-𝐑_i| u_i(),

and similarly _i^G is obtained by substituting v_i() for u_i(). Perfect segregation then corresponds to _i^F = - _i^G in the steady state, regardless of the initial distributions. The less the segregation strength, the less deviation from the initial condition over the time evolution: _i^F∼_i^F(t=0) and _i^G∼_i^G(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 measures: 


(a) the spatial average of the magnitude of: 𝐒_i = _i^F + _i^G. Note the vector sum, thus for perfect segregation we get: 𝐒_i = 0. (Overlines mean spatial average over all cells at a given time).

    S(t)/Q(t) = |𝐒_i(t)|/Q(t)= |_i^F(t) + _i^G(t)|/Q(t).

(b) the standard deviation of Q_i normalized by its mean, 

    δ Q(t)/Q(t)≡(|Q_i(t) - Q(t)|^2)^1/2/Q(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. (<ref>a1) and (<ref>b1), 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. (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. In SLCI, as suspected, 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. (<ref>). In (a2) and (a3), the initial (red) and final (blue) cell-by-cell distributions of S_i/Q_i for SLCI and NLCI, respectively. Again, the near-zero final values of the ratio for NLCI shows perfect segregation. In (b2) and (b3), cellular polarity normalized by spatial average at the corresponding time-point, i.e. initial and final. The width of the distribution shrinks dramatically in NLCI, whereas it remains comparable to its initial value in the SLCI case. 





Unequal interaction ranges (λ_uu≠λ_uv). Here, we discuss various cases of unequal λ_uu and λ_uv. In order to thoroughly investigate and distinguish the role of the two interactions, we run simulations on identical lattices, with identical initial distributions, for the following cases: (i) 0.1 ≲λ_uu = λ_uv≲ 0.8, (ii) λ_uu = 0.01 and 0.1≲λ_uv≲ 0.8, and (iii) 0.1≲λ_uu≲ 0.8 and λ_uv = 0.01. In summary, while the length-scale of nonlocal interactions vary from ≃ 0.1 to ≃ 0.8, local interactions are modeled by λ = 0.01. The magnitude of geometric disorder takes the values, ϵ_0 = 0 , 0.2 , 0.45 , 0.6. Borrowing the abbreviation of the celebrated local-activation–non-local-inhibition mechanism, i.e. LA-NLI, corresponding to (ii), the three regimes can be labeled as, (i) NLA-NLI, (ii) LA-NLI, and (iii) NLA-LI. In SI. (4.4) these three regimes are discussed in more detail, and the results of several simulations are illustrated in SI. Fig. (4), supporting the following conclusions. The most important findings are as follows: (1) First and foremost, as expected, regime (iii) NLA-LI is not able to stabilize long-range polarity. A localized complex, practically does the opposite of what is required for segregation, by activating similar complexes nonlocally (i.e. on other edges), and not repelling unlike complexes to the opposite side, but only inhibiting them in a small vicinity of itself. Therefore, we mostly focus on the first two regimes, namely NLA-NLI and LA-NLI: (2) While for ϵ_0≲ 0.5, the angular correlations of polarization fields are comparable in NLA-NLI and LA-NLI, for ϵ_0 > 0.5, the angular correlation arising from NLA-NLI is arguably better than LA-NLI's, suggesting the importance of nonlocal activation in highly disordered tissues. (3) For small geometrical disorder, the range of λ_uv for which the polarization is established, extends from the above to λ_uv≃ 0.85 in LA-NLI case, whereas in NLA-NLI λ_uu = λ_uv≃ 0.7. (4) For small ϵ_0 ≃ 0.4, where both LA-NLI and NLA-NLI mechanisms work equally well, NLA-NLI manages to establish the polarity faster, by up to a factor of two. This discrepancy grows with geometrical disorder; of course for ϵ_0 ≳ 0.6, LA-NLI fails to polarize the tissue, yet we observe the patterns reach the steady state more slowly compared to the NLA-NLI case. 




Stability analysis. We performed numerical stability analysis on the above cases. Starting from a state of aligned dipoles, for identical geometries and stochastic noise, we find that unlike NLA-NLI systems, in highly disordered tissues, long-range correlation is destroyed in LA-NLI systems for large enough stochastic noise. Figures (<ref>) demonstrate the results for ϵ_0 = 0.45 and ϵ_0 = 0.6. In each case the initial condition is the same for both LA-NLI and NLA-NLI cases Figs. (<ref>a1) and (<ref>b1). For NLA-NLI and LA-NLI, we use λ_uu = λ_uv = 0.5ℓ_0, and λ_uu=0.1ℓ_0, λ_uv = 0.5ℓ_0, respectively. The final distributions are shown in Figs. (<ref>a2), (<ref>a3), and (<ref>b2), (<ref>b3) for ϵ_0 = 0.45 and ϵ_0 = 0.6, respectively. By comparing the final distributions we realize that while for ϵ_0 = 0.45, the final distributions of LA-NLI and NLA-NLI remain relatively narrow around the mean values, in the case of large geometrical disorder ϵ_0 = 0.6, LA-NLI fails to establish correlated polarity. Note that even in the case of NLA-NLI, the final polarities are rotated compared to the initial condition, which is due to the disordered geometry determining the final direction of polarity; the stochastic noise is large enough to drive the polarization from its false fixed point, to the actual one dictated by geometry, yet the angular coherence of the polarization is preserved. In agreement with our simulations (not shown here), this rotation is absent in (nearly) ordered tissues. Therefore NLA-NLI seems to be necessary for long-range polarization to be stabilized in highly disordered geometries.
 

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: ℳ_ij() = ℳ_0 (x_ij - X_i); where ℳ_0 is the slope of the gradient, 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 over a certain τ_res and τ_0, requires weaker cues in systems with NLCI than those with SLCI. For persistent cues (γτ_0 ≳ 100), NLCI responds to signals as small as ℳ_0 ≃ 0.05 over γτ_res≃ 2, whereas SLCI requires at least ℳ_0 ≃ 0.5. The minimum ℳ_0 increases for smaller τ_0's. For example, in NLCI with ℳ_0≃ 0.05 a nearly persistent signal is required, whereas for larger ℳ_0 ≳ 1, even a rapid transient signal γτ_0≃ 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 geometrical 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 <cit.>, 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 (or row) of cells, with significantly larger asymmetry 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. 


Longitudinal vs. transverse correlations. A few remarks are in order regarding correlations. As discussed above, correlation function and length as introduced in Eqs. (<ref>) and (<ref>), 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 lateral correlations. This discrepancy leads to formation of (transient) vortex-like structures. The correlation lengths shown in Fig. (<ref>) do not take into account this effect, and are angular averages of the correlation length. 




Vortices and saddles. Several theoretical <cit.> and experimental studies <cit.> 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; geometry <cit.>, or genetics <cit.>. In Drosophila wing, where Fz and Vang localize distally and proximally, respectively, the coupling to global cues is believed to be dependent on the existent of Fat <cit.>. While small fat^- clones (≲ 10 cell diameters), exhibit little deviation from wild-type polarization, larger clones show swirls, 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 sector of the parameter space wherein the model parameters lie. Generally speaking, long-lived (steady) defects show up in parts of parameter space that are in between a polarized and an unpolarized sector. We observed two distinct types of steady defects: (a) As λ is increased from SLCI to NLCI regime, there exists a narrow range 0.05≲λ≲ 0.15 (for κ = 1), over which vortices and saddles appear as long-lived structures; Fig. (<ref>a). (b) Another situation that shows qualitatively similar behavior is for g_0≲ g_0^*, i.e. under-expression of one of the membrane-bound proteins, that is interpreted as a global mutation within the context of our model; Fig. (<ref>b). Such patterns are indeed observed in Vang mutants <cit.>. Local mutations of various kinds are fully discussed in Sec. (<ref>), and SI. Sec. (6). An interesting observation regarding the second type is that for a specific disordered tissue, upon cranking up g_0, from ≃ 0.2 to ≃ 0.3, the characteristics of the steady-state pattern of polarity remains very much the same, except the magnitude of polarity is increases. The swirls and branches gradually merge and align for g_0≳ 0.3, and long-range polarity is stabilized. The range of g_0 over which similar patterns are stable depends on the level of disorder. The more disordered the tissue, the wider the range and the more stable the patterns. The above behavior is independent of the initial condition and stochastic noise, implying that  the polarity pattern is fixed by the the microscopic geometry of the tissue when in regimes where the correlation lengths are of the order of a few cell diameters. Finally, we would like to make a remark on the stability of patterns. In tissues with small disorder, “steady” defects might eventually disappear over very long timescales, and for sufficiently large stochastic noise.  
In general, as briefly mentioned above as well, for NLCI, the geometrical information is read by nonlocal interactions which locally biases the dipoles. As g_0 increases, the correlation length increases and the global direction is chosen collectively. Nonetheless, the direction is determined by the bias provided by geometry, in the presence of large geometrical disorder and/or small stochastic noise. 













§ THE EFFECT OF TISSUE ELONGATION ON POLARIZATION

Elongation is suggested to be acting as a global cue in some systems, e.g. mammalian cochlea and mice medial-lateral skin <cit.>. Furthermore, it has been shown in the same study that the perpendicular polarization is not due to a naïve incorporation of length in the definition of polarization, but that the short junctions are indeed depleted of proteins. Here we show that NLCI, through increasing the strength of the cooperative self-interactions α_s,β_s, enhance the stability of F-G complexes on longer junctions. Intuitively, unbound proteins receive, on average, stronger attractive and repulsive signals from complexes localized on longer junctions. This effect results in larger α_s,β_s, thus enhanced junctional polarity. In SI. Fig. (5b), we plot the dependency of self-interactions on the length of junctions, for different λ's.

The elongation of a cell, is characterized by a traceless and symmetric nematic tensor, with the diagonal and off-diagonal elements, ±ε_i,1 and ε_i,2, respectively. The index of tissue elongation reads ℰ = N_c^-1∑_i=1^N_c(_i,1^2 + _i,2^2)^1/2; see SI. (5). In ordered tissues, there are two possible choices for elongation axis. Elongation parallel to a pair of parallel edges reduces the sixfold symmetry to a twofold associated to long junctions, and a fourfold, Fig. (<ref>a1); and vice versa if the elongation is perpendicular to a pair of edges; Figs. (<ref>a2,<ref>a3). While (a1) and (a2) are polarized perpendicularly to the axis of elongation, (a3) exhibits parallel polarization. The latter is destabilized by NLCI inhibiting localization of unlike complexes on adjacent junctions. In geometrically disordered tissues, elongation is a mixture of (a1) and (a2), both of which give rise to perpendicular polarization.



Elongated systems involve three different length scales. With L the length of long junctions, we have: (i) λ≲ℓ_0< L, (ii) ℓ_0≲λ< L and (iii) ℓ_0< L<λ. For (i) and (ii), NLCI leads to the separation of positively and negatively polarized edges. The third regime λ≳ℓ_0 , L, however, is unstable, like the regime λ≃ℓ_0 discussed in isotropic case. Time evolution of Q(t), P(t), 𝒪(t) and ξ(t) are shown in Fig. (<ref>c1) and (<ref>c2). Figure (<ref>b) illustrates the steady-state of the polarization field in an elongated tissue with ℰ = 0.4. In order to investigate the detection of this global cue, we ran simulations on identical tissues but with different values of elongation ℰ = 0 – 0.5; see SI. Fig. (5c). We find that the perpendicular polarization appears at ℰ^* ≃ 0.1, which is in a very good agreement with that in <cit.>. In order to see (a) whether the observed polarity is a collective effect or is due to single-cell geometry, and (b) that polarity is not a trivial geometrical effect, we make the scatter plots of magnitude and angle of the polarities vs. the magnitude and angle of their elongation for individual cells; see Figs. (<ref>d1) and (<ref>d2). The infinitesimal correlations between the two indicate that the polarization vector is not a local but a collective effect. 



§ LOCAL MUTATIONS AND THE ASSOCIATED PHENOTYPES

PCP mutants exhibit lack of orientational order, which is induced autonomously and/or non-autonomously by mutant clones <cit.>. The phenotypes resulting from loss- and gain-of-function of various components are commonly used to specify the roles of the corresponding proteins. Here we introduce three distinct classes of mutations within the context of our model: Absence of (I) one or both of the membrane proteins, (II) cytoplasmic proteins, in clones embedded in wild type (WT) or mutant backgrounds; and (III) enhanced geometrical irregularity in a patch of cells. 

Since the polarized junctional complexes are assumed to be Fz:Fmi-Fmi:Vang, Fmi-associated mutations cannot be tested separately; in our model fmi^- translates into double mutants fz^-Vang^-. While some studies have reported distinct phenotypes in Vang^-fz^- and fmi^- (e.g. <cit.>), others such as <cit.> have seen great resemblance between the two, namely minimal non-autonomy. This is considered as a piece of evidence in favor of bidirectional signaling hypothesis <cit.>, and lends more support to our simplifying assumption of similar roles of Fz and Vang in the primary complexes. 

As discussed in Sec. (<ref>), the roles of Dsh and Pk are captured effectively by the magnitude of cooperative interactions α,β, as well as the length-scale of the cytoplasmic interactions λ. The dependencies of these parameters on Dsh and Pk are complicated and depend on the concentrations as well as the feedback loops between the two components. In order to uncover the major roles of Dsh and Pk in core PCP pathway, we compare the in silico phenotypes induced by α , β→ 0 and λ→ 0, to the in vivo phenotypes dsh^- and pk^-. Experimental observations reveal minimal non-autonomy of dsh^- and pk^- clones in WT backgrounds <cit.>. However, the autonomous effects of the two are discernible: while dsh^- cells remain nearly unpolarized, pk^- clones seem to be almost perfectly polarized parallel to the WT background <cit.>. Comparison with our results (see SI. Fig. (6)), suggests that the mutants generated by diminished cooperative interactions look very similar to dsh^- clones, with no cellular polarization, and almost zero non-autonomy. The λ-mutants, on the other hand, resemble, to some extent, the pk^- clones, though with imperfect angular coherence within the clone (a realization of non-trivial MF solutions depicted in Fig. (<ref>b)); indeed they look more like dsh^-pk^- double mutants. As such, we predict that while both Dsh and Pk contribute to the magnitude and the length scale of nonlocal cytoplasmic interactions; Dsh is mainly in charge of the local interactions (α , β), whereas nonlocal interactions are mediated by both Dsh and Pk. This hypothesis has the following important prediction. We know that the minimum concentration of Vang required for the emergent polarization, increases as α,β decrease; Sec. (<ref>). Therefore we predict that under-expression of Vang can be—to some extent—compensated for, by over-expression of cytoplasmic proteins, mainly Dsh. This is an elegant manifestation of the collaboration between cytoplasmic and membrane proteins in establishing long-range polarization. Also note that our results suggest that λ is not a simple diffusion length of proteins, but is indeed modified as the concentrations of cytoplasmic proteins change, implying nonlinear effects, perhaps due to their interactions that assist with cytosolic diffusion.

Results of simulations and the predicted phenotypes of our model for all mutant classes are shown in SI. (6), and their schematics are tabulated in Fig. (<ref>). The first and second rows belong to type-I mutants. The left and right panels of the third row corresponds to type II and III, respectively. Comparing the phenotypes of membrane proteins with cytoplasmic ones, we notice two distinctions: (1) Non-autonomous effects of the former are stronger than those of the latter. This is consistent with the fact that cell-cell communications occur through membrane proteins, whereas cytoplasmic proteins are the carriers of intracellular interactions. (2) Both cytoplasmic mutants attract the dipoles, implying the localization of Fz at the WT-mutant boundaries in mutants lacking cytoplasmic proteins.



In order to test our model's predictions in the case of clones with high geometrical irregularities, studied in Ref. <cit.>, we simulate type III mutants both in the absence and presence of a global cue. In experiments geometrical irregularities are induced by PTEN^-, (marked by darker shades Fig. (<ref>)). Ma, et.al. <cit.>, found that while the alignment of polarization field is preserved in single mutants fat^- or PTEN^-, the angular correlation is disrupted significantly in double mutants fat^-PTEN^-, implying that geometrical disorder is an obstacle to the faithful propagation of polarization. In our simulations we used the statistics of cell area in PTEN^- patches from Ref. <cit.>, and introduced a patch with strong geometrical disorder and shrunken cellular area, that smoothly dissolves into a rather ordered background. In the absence of a global cue, the polarity field shows strong aberrations with swirl-like patterns centered at the mutant patch. Adding the global cue unwinds the swirls and alignment of WT polarization field reappears; see SI. (6). Disrupted polarity in cells with altered geometry can be understood in our model, by noting that NLCI sustains polarity, only within a certain range of λ/ℓ_0. Upon decreasing the cell size, this ratio exceeds the upper bound of the NLCI functional range, and the polarization is destabilized. 

Comparison with experimental observations <cit.> (type I), and <cit.> (type II), and <cit.> (type III) reveals qualitative similarities between the in vivo, and in silico phenotypes, suggesting that our model is capable of capturing the salient functions of different PCP components, and their coupling with global cues as well as cell geometry, all of which are concomitants of NLCI. 



§ SUMMARY AND DISCUSSION

In an attempt to understand the role of cytoplasmic interactions in PCP, and based upon the well-established facts deduced from the experimental studies, we devised a generalized reaction-diffusion model by incorporating nonlocal cytoplasmic interactions. Although we utilized the knowledge on core-PCP to construct the components of our model, no pathway-specific assumptions are made regarding the molecular details and interactions. Thus, we suspect as long as the dominant mechanism of cytoplasmic transport of proteins can be modeled by reaction-diffusion-like dynamics, our model should be able to capture and predict, at least the qualitative behavior of tissue cell polarity. We explored different scenarios of intra- and intercellular interactions, and particularly specified for a generic set of model parameters, the optimal range of cytoplasmic interactions length scales to achieve long-range polarization: 0.2≲λ/ℓ_0≲ 0.7. Investigating the cases of unequal λ_uu and λ_uv, reveals that in disordered systems with ϵ_0 ≳ 0.6, the angular correlation drops significantly compared to that with identical λ's. We further examined the response of polarization to external cues, and concluded that NLCI is essential to detecting directional signals in even moderately disordered tissues. 

A direct consequence of NLCI is the readout of cellular geometry. Of particular interest to our study is tissue elongation, as a symmetry-breaking cue. We showed that NLCI is responsible for collective stabilization of polarity perpendicular to the elongation axis. Agreement with the observed value of elongation at which the detectable perpendicular polarization appears, is suggestive of the NLCI as the dominant PCP mechanism in systems like mammalian cochlea and skin. We shall emphasize that this prediction is only valid under the following assumptions: (a) polarization is predominantly induced by reaction-diffusion processes, and (b) lattice dynamics are negligible on the native PCP kinetic timescales. Other mechanisms such as the polarization of microtubules, anisotropic stress, and relative timescales of cell division and PCP relaxation, are also known to influence the direction of polarity <cit.>. In order to examine the predictive power of our model, we studied three classes of mutants and found arguably similar phenotypes to the experimental observations, which helps with interpretation of the model parameters and prediction of the roles of Pk and Dsh in core-PCP pathway. 


Comparison with other models. Although our model is not the first semi-phenomenological approach to the problem of PCP, we believe that the features included in this model, capture a broader range of recently observed phenomena. Some prior studies (e.g. <cit.>) consider one-dimensional arrays of cells. Two-dimensional systems, however, call for a careful investigation of the cytoplasmic mechanism of segregation. Other successful models such as that studied in Ref. <cit.>, infer effective interactions from the observed response of the polarization to a combination of processes; cell elongation, cell rearrangements, and divisions. However the model does not aim at explaining the underlying cellular level mechanisms. Among the models derived from intracellular interactions, the ones put forward by Burak and Shraiman in Ref. <cit.>, and Abley et.al. <cit.>, are closely related to ours. The former investigates the role of nonlocal inhibitory interactions and demonstrates that LA-NLI is sufficient to fully drive the intracellular segregation of Fz and Vang, in ordered tissues, and for correlated polarity to emerge in the absence of external cues. However, as the authors point out, the role of nonlocal activation (i.e. NLA-NLI), as well as geometrical disorder remain to be investigated. We tried in this paper to address these questions. Apart from the modeling point of view and interesting aspects of the behavior of polarity, we found it crucial to consider the possibility of nonlocal activation from a phenomenological perspective. Given that the stabilizing and destabilizing interactions between the two complexes are carried by the same set of cytoplasmic proteins, it is a conceivable possibility for the two length scales to be comparable, i.e. λ_uu≃λ_uv. Therefore, for identical set of model parameters as well as initial conditions, we examined the behavior of tissue polarity for different values of the two length scales, as the geometrical disorder increases. Interestingly, we observe that the two seemingly unrelated factor that are missing in <cit.>, become important in relation with one another. In highly disordered tissues, nonlocal stabilizing interactions are more efficient in the cytoplasmic segregation of PCP proteins; hence NLA-NLI is a reliable mechanism to stabilize the long-range polarity, at least in disordered tissues. Further evidence was also provided through stability analysis of initially polarized disordered tissues, for systems with NLA-NLI versus LA-NLI. 

The second study by Abley et.al. <cit.>, includes the two components missing in the former paper, and finds that local polarity alignment is achieved through intracellular partitioning, accompanied by direct cell-cell coupling. However, they find that long-range polarization requires global cues, in the absence of which longitudinal and lateral coordination together give rise to steady-state swirls over domains of size of a few cell diameter. Swirling patterns have been observed repeatedly, and as discussed in Sec. (<ref>), are a consequence of unequal longitudinal and lateral correlations in the early stages <cit.>. Burak and Shraiman found that, in ordered lattices, large stochastic noise causes swirls in the absence of global cues, that disappear over long timescales, perhaps beyond relevant developmental dynamics. Although we observe steady swirls in parts of our parameter space, they appear only transiently, and long-range polarization is stabilized—even without global cues—in other sectors of parameter space. This is observed in both ordered and disordered systems. As far as the phenomenology of nonlocal activation and inhibition between complexes goes, the model proposed in <cit.> appears to be quite similar to ours. Therefore, as also mentioned in Sec. (<ref>), we suspect that although their simulations cover a relatively wide range of parameters, the discrepancy in the findings originate from the choice of parameter; e.g. distance from the critical point g_0/f_0, or relative formation and dissociation rates of the complexes κ/γ. In summary, we find stable fixed points with long-range spatial correlations in the absence of global cues, though they might require global cues as a drive to settle in the fixed point within the developmental timescales. 








With regards to the mutant phenotypes, former studies such as <cit.>, have proposed mathematical models that successfully reproduce the in vivo phenotypes. Adopting an effective phenomenological approach to determine the minimal set of criteria essential to the large-scale PCP alignment, we tried to keep the number of model parameters as few as possible, and focus on the major mechanisms at work. The successful recapitulation of the perpendicular axes of polarity and elongation, as well as the in vivo phenotypes, speaks to the predictive power of our model in identifying the roles of PCP components. We believe that in spite of adopting the core-PCP as the reference system, the prediction of our model are applicable to pathways other than core-PCP, provided that the polarity is predominantly governed by cytoplasmic RD-type equations; for instance Fat/Dachsous PCP. 


Limitations of the model.
First and foremost, we shall emphasize again that some parameters (α,β,κ,λ) represent effective quantities arising from combinations of molecular networks coupled through feedback/feedforward loops. Therefore, our model is only to account for, and investigate, the primary roles of the molecular components. There exist other quantitative models including pathway-specific details with constants inferred from experiments; e.g. <cit.>.

The second point is the assumption of uniform distribution of proteins on each junction. We know from experiments that proteins form clusters at certain loci on the junctions called puncta. However, our approximation makes the simulations faster by orders of magnitude. This assumption, however, has a downside too. Should we allow the proteins to distribute nonuniformly on the junctions, the polarity degrees of freedom vary more smoothly around the cell, and the polarization pattern shows smaller fluctuations compared to what we obtained. Therefore, we suspect that the PCP correlations could be improved, had we solved the full RD equations on the perimeters of all the cells.


Predictions and outlook. 
A list of our model's predictions is as follows. (1) The minimum concentration of Vang required for the collective polarization to appear depends inversely on the concentration of Dsh. Therefore, polarity alignment can be retained in tissues with under-expression of Vang, by over-expressing Dsh. This prediction is crucial to understanding the collaboration of membrane and cytoplasmic proteins in transmitting directional information between abutting cells; and can be readily tested by tuning the expression levels of Dsh and Vang. (2) Since NLCI is suggested to be carried by Pk (and in part by Dsh). Thus, the role of elongation as a global cue is dependent on Pk, and the knockdown of pk and dsh must invalidate the guaranteed orthogonality of polarization and the elongation axis  (3) NLCI enhances the sensitivity to the gradient cues, namely smaller magnitudes and/or shorter time scales of the gradients are required for the dipoles to reorient, in NLCI compared to SLCI systems. This discrepancy becomes more pronounced in tissues with large geometrical disorder. Furthermore, transient global cues are sufficient to reorient the polarity in disordered tissues with NLCI, provided the decaying timescale of the cue (τ_0) is larger than the intrinsic ordering timescale of polarity, i.e. that in the absence of global cue. The latter is estimated to be of the order of ∼ 10 (hrs) in Drosophila wing; see e.g. <cit.>. Transient cues can be provided through temporary expression of a gene using techniques that allow for spatio-temporal control of gene expression, such as heat-shock and light-switchable promoters <cit.>. A natural conclusion of (2) and (3) is that Pk (and Dsh) are essential to the detection of both geometrical and molecular cues. Thus, their absence impairs the readout of all such global cues. (4) The length scale of cytoplasmic interactions are dependent on the cytoplasmic protein concentrations, due to possible nonlinear effects on diffusion constants. Interpreting λ to be an increasing function of Pk (and possibly Dsh), and given that the polarization is destabilized for λ≳ 0.8, our model predicts that excess Pk destroys the PCP alignment. Interestingly, this effect was observed in a study by Cho, et.al. in Ref. <cit.>. Since the NLCI is contingent on the presence of Pk (and Dsh), predictions (2), (3) and (4) can be tested—in the absence of the corresponding global cues—by under- or over-expressions of Pk, as the representative of NLCI. (5) The direction of the polarity in disordered tissues is chosen by the geometry, is independent of initial distribution, and shows robustness against stochastic noise and small external cues. This can be tested experimentally, by comparing the response of already polarized tissues, to global cues in transverse direction, for ordered and disordered geometries, e.g. before and after the ordering in Drosophila wing. 

Finally, in spite of several insightful findings regarding the mutual interplay of PCP and tissue mechanics <cit.>, the relevant molecular and physical mechanisms are yet to be explored. Planar polarity and cell packing are known to mutually influence one another. On the other hand cell packing is highly susceptible to mechanical tension. A natural question would be, how does PCP couple to tension at the molecular level? Furthermore, given the role of microtubules in parallel polarization through biasing the transport of membrane proteins, another important question is, under what conditions does this mechanism dominate the diffusive transport, whereby polarity seems to develop perpendicularly to the axis of elongation? Our study lays the groundwork for further investigations, by uncovering one of the scenarios through which PCP couples to cell geometry.


Methods. Dynamical simulations are carried out using Runge-Kutta method of 4th order, with the time steps of 10^-3 γ^-1, on lattices of size 40 × 40 cells. For each cell, starting from a randomized distribution of F and G proteins, the concentration of proteins are evolved according to the RD equations. For each value of geometrical disorder, and fixed model parameters, simulations are run for 500 initial conditions. A white and Gaussian stochastic noise is also added to the RD equations. All points on the perimeter of a given cell interact with each other through the kernels introduced in Sec. (<ref>). Using the assumption of uniform proteins distribution along all junctions, it suffices to compute the geometrical coefficients (𝒦_μν) of junction-junction interactions by integrating the kernels along the two junctions, and for all pairs of junctions within a cell. Therefore, the integrals reduce to matrix products. The matrices are calculated for each configuration of disordered lattices. Boundary conditions are chosen to be periodic along both axes. The bulk cues are modeled as constant gradients in each cell. A “drive" term is added to the RD equation, which for each point on the perimeter of a given cell, is proportional to its distance from the centroid of the cell. The boundary cues are incorporated by polarizing a column of the boundary cells through either a permanent or a transient bias, similar to the bulk cues. The strengths of the biases decay exponentially in time, in both bulk and boundary cues.

Acknowledgement. S.S. was supported by the Gordon and Betty Moore Foundation through Grant GBMF2919.01. M.M. would like to thank the Simons Foundation MMLS program for support.






unsrt





equationsection
figuresection



Supporting Information






§ THE MODEL AND ITS INGREDIENTS

We begin by writing the full reaction-diffusion equations for the binding/unbinding dynamics of the local density of the complexes complexes [Fz:Fmi-Fmi:Vang], namely u_ij(). The corresponding equation for [Vang:Fmi-Fmi:Fz] or v_ij(), is obtained by simultaneous replacements: Fz:Fmi ↔ Fmi:Vang and u↔ v. The primary complexes Fz:Fmi and Fmi:Vang, are denoted by F and G, respectively. Their corresponding concentrations are then f and g.


    du_ij()/dt   =κ f_i^ubdg_j^ubd(1+α∑_{k}_i∫_i∩ k d'𝒦_uu(-')u_ik('))
       - γ u_ij()(1+β∑_{k}_i∫_i∩ k d' 𝒦_uv(-')v_ik('))
       + η(, t),


where γ^-1 is the timescale associated with complex dissociation. As mentioned in the Main Text, the kernels 𝒦() are assumed to be of the form of exp(-||/λ), which was motivated by the diffusive nature of the cytoplasmic proteins carrying the interactions. Although the kernels couple the concentrations of the complexes on the boundaries of the cells, the coordinate  can in principle represent any point within the cytoplasm as well as the junctions. Suppose that a cytoplasmic protein C, obeys a diffusion equation with degradation time τ_c:

    ∂ c(,t)/∂ t = D_c∇^2 c(,t) - τ_c^-1 c(,t).

Assuming the degradation and hence dynamics of C takes place on a much faster timescale than that of PCP, i.e. γτ_c≪ 1, it suffices to only consider the steady-state solutions of protein C. With D_c, the diffusion constant of the cytoplasmic proteins C, we get the following equations for C:


    D_c∇^2 c(,t) - τ_c^-1c(,t) = 0,


    c( , t) = ∫_ d' c_0(') exp(|-'|/λ_c).


Here, =0 corresponds to the specific point at which the concentration of C is measured by superimposing the concentration of proteins diffused from '. Next, λ_c = √(D_cτ_c) is the diffusion length of protein C. The diffusing proteins enhance the formation of like complexes and suppress that of unlike complexes. This nonlinear effect in turn, depends on the respective interactions and their binding affinities with the target complexes. Altogether, the coefficients are lumped into the phenomenological constants α and β.


For notational convenience, in the following paragraphs Greek letters label edges, e.g. μ≡ i∩ j. Using Eq. (<ref>) and the definitions of f^ubd_i and g^ubd_j, we can solve for the dynamics and the steady states of u_μ(). Evidently, solving the above integro-differential equations is a cumbersome task. One possible simplification is junctional averaging: u_μ = ∫_μd' u_μ(')/ℓ_μ, in terms of which we have: 


    du_μ/dt =    f^ubd_ig^ubd_j(1+α_μ∑_ν∈ i𝒦_μνu_ν)
       -γ u_μ(1+β_μ∑_ν∈ i𝒦_μνv_ν) + η_μ(t).


In the above equation, α_μ= α/ℓ_μ, β_μ= β/ℓ_μ, and for the kernels 𝒦_μν = 𝒦_νμ =∬_μ,νd d'𝒦(-'). The diagonal elements equal

    𝒦_μμ = 2ℓ^2_μx_μ^-2(e^-x_μ + x_μ - 1),

in which x_μ = ℓ_μ/λ. The effective stochastic noise on a junction μ reads: ⟨η_μ(t)η_μ(t')⟩ = η_0^2ℓ_μδ(t-t'). It is noteworthy that, as we see in Eq. (<ref>), the junctional averaging is a very simplifying approximation that reduces the integrations into matrix products. However, in reality the core proteins for example Flamingo and Frizzled in the prepupal and pupal wing of Drosophila, are observed to be persistently localized at subdomains of plasma membranes, called “puncta” <cit.>.



§ DEFINITIONS OF JUNCTIONAL AND CELLULAR POLARITY


Planar cell polarity can be defined at either junctional or cellular level. For individual junctions, polarity is defined as the difference between the concentrations of u_ij=[F_i:G_j] and the opposite dimer, u_ji, thus p_ij = u_ij - u_ji = u_ij - vij. 

Cellular polarization, on the other hand, is referred to as the asymmetric distribution of PCP proteins in the cells. Cell polarity, depending on the symmetry of binding complexes, might be either a vector quantity, called vectorial polarity, or a nematic, which is then called axial polarity. Vectorial and axial polarities are used when PCP proteins form heterodimers and homodimers, respectively. The former is a vector identified by a magnitude and angle ∈[0,2π), whereas the latter is a traceless nematic tensor with a magnitude and an angle ∈[0,π). In our case of study, the polarization involves two distinct complexes, Fz and Vang, hence vectorial polarity. However, we shall introduce both vectorial and axial cell polarities. 


(i) Vectorial Polarization. This definition is used to calculate the “dipole moments” of the distribution of bound Fz (or Vang), around each cell. The polarization vector associated with cell i is defined as:

    _i = 1/2 ∫__i d' '-𝐑_i/|'-𝐑_i| (u_i(') - v_i(')) = p_i^x x̂ + p_i^y ŷ.

Here 𝐑_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 polarization vector is identified by a magnitude |_i| and an angle θ_i∈ [0, 2π) measured from the x-axis:

    |_i| = √((p_i^x)^2 + (p_i^y)^2),     and    ϕ^p_i = tan^-1(p_i^y/p_i^x).

The polarization of the tissue with N_c cells, and its global order are characterized by the following quantities: (1) average polarization (2) average magnitude of polarization, and (3) the ratio of (1) and (2):

    (1) : = 1/N_c∑_i=1^N_c_i ,
       (2) : Q = 1/N_c∑_i=1^N_c|_i|, 
       (3) : 𝒪 = ||/Q.


(ii) Axial Polarization. The second definition is a measure for cellular polarity that is used especially when dealing with axial nematic PCP (like Fmi (Celsr) homodimers) which merely determines the axis of polarization by measuring the traceless nematic tensor of the polarity:

    𝙿_i    = 
    [  𝒫_i,1  𝒫_i,2;  𝒫_i,2 -𝒫_i,1 ],
    𝒫_i,1   = ∫_0^2πdϕ u_i() cos(2ϕ), 
    𝒫_i,2   = ∫_0^2πdϕ u_i() sin(2ϕ).

In the above equation, ϕ is the polar angle of point  on the perimeter of cell i and is measured, with respect to a reference axis. The magnitude of polarization equals 𝒫_i = (𝒫_i,1^2 + 𝒫_i,2^2)^1/2 = |𝙿_i|^1/2. Its orientation is determined by angle ϕ^p_i satisfying cos(2ϕ^p_i) = 𝒫_i,1/𝒫_i and sin(2ϕ^p_i) = 𝒫_i,2/𝒫_i. In Sec. (<ref>), we use the same general formalism for the cell shape nematic tensor. 



§ MEASURE OF QUENCHED DISORDER AND DEFECTS IN 2D SYSTEMS

In the case of 1D systems, disorder is only in the lengths of junctions, and there is no room for the change in the topology of the network. The edge lengths are ℓ_i = ℓ_0 (1 + ϵ_i), where ϵ_i∈[-ϵ_0, +ϵ_0] with uniform distribution, and ⟨ϵ_iϵ_j⟩ = ϵ_0^2δ_ij/6. In 2D the quenched disorder refers to (a) unequal edge lengths, and (b) defects defined as the non-hexagonal polygons, tiling the plane. The level of quenched disorder is controlled by randomizing the sites of a triangular lattice, based on which the polygonal lattice is generated using Voronoi tessellation. The edge lengths of the Voronoi lattice ϵ_i, and density of defects n_d are then obtained by ensemble averaging over all the realizations of the disordered triangular lattice. Perturbing the sites of a triangular lattice, we get: _i = _i^0 (1 + Δ_i), with {_i^0}. The spatial disorder term Δ_i is uniformly distributed in range [-Δ_0, +Δ_0], with local correlations: ⟨Δ_i·Δ_j⟩ =Δ_0^2δ_ij/3, where δ_ij is the Kronecker delta function. For Δ_0 = 0, the corresponding Voronoid lattice would be an ordered hexagonal lattice. By displacing randomly the sites of triangular lattice, we can distort the resultant Voronoi lattice. In order to obtain the disorder statistics of the Voronoi lattice, i.e. variations in the edge lengths ϵ_i, as well as the density of defects n_d, we average over ensemble of disordered triangular lattices. Defects with finite (i.e. nonzero) density of defects in the limit of infinitely large systems (i.e. thermodynamic limit), appear above a certain threshold of disorder Δ_d≃0.25, in the triangular lattice (see Fig. (<ref>)). The edge-length disorder in the Voronoi lattice increases linearly with Δ_0 for Δ_0≲ 0.5.





§ MEAN-FIELD AND NUMERICAL SOLUTIONS

We define the mean-field approximation in this system as constant f_i^ubdg_j^ubd for both sides of all junctions. As mentioned in the Main Text, the validity of this assumption follows from the diffusion-like dynamics of p, s. Bound proteins redistribute across the tissue until a relatively smooth state is reached, therefore the free proteins f_i^ubdg_j^ubd too, distribute uniformly. We also check this assumption numerically; see Fig. (<ref>b1) and (<ref>b2). Here we first elaborate on the MF solutions in 1D, then discuss the 2D case where the MF solutions are divided into two different classes: trivial and nontrivial. 



 §.§ Mean-Field Solutions in One Dimension

In one dimension the cells are juxtaposed in an array and are separated by junctions. The proteins localize on both sides of these junctions and form dimers. A general scheme of one-dimensional arrays can be seen in Fig. (<ref>). 
The reaction-diffusion (RD) equations governing u and v complexes on a junction shared by cells i and i+1 read:

    du_i,i+1/dt =     f^ubd_ig^ubd_i+1(1+α u_i,i+1) -γ u_i,i+1(1+β u_i+1,i).

A similar equation exists for the opposite complex v_i , i+1 = u_i+1,i. The above equations were originally proposed by Mani, et.al. in Ref. <cit.>, namely our general RD equations introduced in the Main Text, reduce to the above equations in 1D with strictly interactions. Therefore we refer the reader to Ref. <cit.>, where the properties and predictions of Eq. (<ref>) in 1D systems are investigated elaborately. Here we briefly reproduce 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 p_i,i+1 and s_i,i+1. This can be done by adding and subtracting Eq. (<ref>) and its counterpart for the opposite complex. We get:



    dp_i,i+1/dt =  (f^ubd_ig^ubd_i+1 - f^ubd_i+1g^ubd_i) + α (f^ubd_ig^ubd_i+1u_i,i+1 - f^ubd_i+1g^ubd_iu_i+1,i) -γ (u_i,i+1 - u_i+1,i),


    ds_i,i+1/dt =  (f^ubd_ig^ubd_i+1 + f^ubd_i+1g^ubd_i) + α (f^ubd_ig^ubd_i+1u_i,i+1 + f^ubd_i+1g^ubd_iu_i+1,i) -γ (u_i,i+1 + u_i+1,i) -2γβ -u_i,i+1u_i+1,i .

In the above equations we have:

    f^ubd_i = F_0 - (u_i,i+1ℓ_i,i+1 + u_i,i-1ℓ_i-1,i)/ℓ_i,i+1 + ℓ_i-1,i =  f_0,i - u_i,i+1ℓ_i,i+1 + u_i,i-1ℓ_i-1,i/ℓ_i,i+1 + ℓ_i-1,i,



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. 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:

    p^2 =    s^2 - 4/αβ    ⇒     s^* = 2/√(αβ)
        g_0^* = γ/κα/1-√(1/αβ)+√(1/αβ).

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. (<ref>), 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/ℓ_0 as the unit of concentration, one can rewrite the above parameters in terms of f_0,ig_0,i+1 and f_0,i+1g_0,i.
Therefore, as mentioned in the Main Text, 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 between unpolarized and polarized states becomes a smooth cross-over.
Below we discuss the MF solutions in ordered two dimensional systems.





 §.§ Trivial Mean-Field Solutions

Defining the trivial mean-field approximation as the translational invariance of polarization along each of the three axes separately, 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. Therefore, the solutions are identical to those of 1D case. This can also be understood intuitively, by noting that in a polarized trivial MF state, each cell is partitioned into a positive and a negative side. Each partition can be thought of as an “effective edge” in the 1D case. Although the lengths of these “effective edges" are different from those in 1D case, the concentrations are equal. Now, in the hexagonal lattice, each edge carries the same |p|,s, equal to those found in 1D case. The net polarity of each cell thus reads, p_c=p_e (1 - 2 cosθ), with p_e 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. The trivial MF solutions take two distinct configurations illustrated in Figs. (<ref>a1) and (<ref>a2), respectively. While the former represents a state with nonzero net polarization the magnitude of which is calculated above, the latter has zero net polarization. The solutions of the second type, are however destabilized when NLCI is included in the model. 

Furthermore, we 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. (<ref>), 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. Without loss of generality 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/√(α). Therefor in case (i), axis 1 is the first axis that shows instability upon increasing g_0 above g_0^*. Therefore, we have f^ubdg^ubd = γ/(kα_1), where,

    f^ubd = f_0 - s_1ℓ_1+s_2ℓ_2 + s_3ℓ_3/2(ℓ_1+ℓ_2+ℓ_3),

and similarly for g^ubd. Again, as derived above, the axes where the polarization is zero, i.e. 2, 3, we have:

    s_2,3^2+2/γβ_2,3(γ -  f^ubdg^ubdα_2,3)s_2,3-2 f^ubdg^ubd = 0.

Using the fact that kf^ubdg^ubd=γ/α_1, we get:

    s_2,3 = 1/β_2,3(1-α_2,3/α_1) + √(1/β^2_2,3(1-α_2,3/α_1)^2+2γ/α_1).

Defining f_0^eff and g_0^eff,

    f_0^eff   = f_0 - s_2ℓ_2 +s_3ℓ_3/2(ℓ_1+ℓ_2+ℓ_3), 
    
    f^ubd   = f_0^eff - s_1ℓ_1/2(ℓ_1+ℓ_2+ℓ_3),

and similar expressions for g_0^eff and g^ubd, we get for p_1,s_1:

    p_1^2 =    s_1^2 - 4/α_1β_1,
    
     s_1 = (f_0^eff +g_0^eff) -   √((f_0^eff -g_0^eff)^2+4γ/α_1).

From the above analysis, we learn that in systems with unequal α's and β's, the system essentially reduces to a one dimensional problem, with effective values for the concentrations of the proteins Fz and Vang. There remains one more interesting case in which α_1 = α_2>α_3. In this case, 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 for the edges is the elongated tissues. In such systems, the above two cases correspond to Figs. (<ref>a1), (<ref>a2), and (<ref>a3), respectively; see Sec. (<ref>).



 §.§ Nontrivial Mean-Field Solutions

The nontrivial solutions in 2D do not possess the translational invariance of polarity vectors, but only satisfy a weaker constraint. The full analytic treatment of this problem is cumbersome. We only briefly touch upon this subject to provide some intuition into how large the basin of attraction is, for a system in SLCI limit. The only assumption in nontrivial MF solutions is that for each junction shared by adjacent cells i and j, the products f_i^ubdg_j^ubd and f_j^ubdg_i^ubd are uniform across the tissue. This assumption is backed up by the following numerical analysis. We compute the r.m.s. of the fluctuations of f^ubd_i g^ubd_j, for all the junctions with fully randomized labels, and compare the initial and final values. We see in Figs. (<ref>b1) and (<ref>b2), that for both SLCI and NLCI mechanisms, the final distribution of f^ubd_ig^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 <cit.>. 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. In SLCI limit and/or for ordered lattices, by virtue of sixfold symmetry, all junctions equally absorb the proteins. Therefore above the bifurcation point, net polarization p and the total amount of proteins s, are identical for all junctions. The only constraint is that three junctions have net positive and the other three have negative polarizations; i.e. three outgoing and three incoming dipoles; see Fig. (<ref>). 



Besides trivial MF solutions depicted in Figs. (<ref>a1) and (<ref>a2), nontrivial solutions like that in (<ref>a3) constitute another class of MF solutions. In order to make it easier to picture such configurations of dipoles, imagine we start from a trivial solution of (a1), where three adjacent junctions carry outgoing dipoles and the other three carry incoming dipoles. Flipping one of the outgoing dipoles violates the MF assumption of f_i^ubdg_j^ubd= constant, unless one of the incoming dipoles flips as well. Since each dipole is shared between two cells, in order to satisfy the constraint in any finite-size system, this flipping process must continue until it forms a loop ending at the first flipped dipole. It is easy to see that the paths can be broken down into self-avoiding loops. Furthermore, all such loops preserve the net polarization of the tissue, set by the value of the control parameter g_0/f_0. However, they disrupts the uniform polarization of the trivial solutions, namely create excitations above the uniform configuration. Incidentally this is what we observe in simulations; the net polarization at the steady state is to a high accuracy independent of initial configuration. The small deviations is however understandable: the constraint of constant f^ubdg^ubd is not guaranteed to be accurately satisfied in real systems with arbitrary initial conditions and stochastic noise. Moreover, in an ensemble of the systems starting from all possible initial conditions, all nontrivial configurations preserving the constraint as well as the net polarization, are equally accessible to the system. The above simple analysis provides insight into why the system in SLCI limit is not guaranteed to reach a state with large-scale polarization, and instead NLCI mechanism stabilizes the states with nearly uniform polarization, by promoting the adjacent junctions to carry the same kind of dipoles, hence segregation. In order for the cells to meet the two criteria simultaneously, the directions of the cell polarizations must also be parallel, whereas in SLCI case, different cells can have different |_i|, hence different orientations. 




 §.§ Equal vs. Unequal Length Scales of Cytoplasmic Interactions

In the Main Text we discussed the effects of unequal nonlocal cytoplasmic interactions that appear in intra- and interspecies kernels, i.e. λ_uu and  λ_uv. Here we elaborate on that discussion and consider various cases to investigate the role of these length scales, and the behavior of angular correlation of polarity. In each case the geometry of the tissue is held fixed. The geometrical disorder is denoted by ϵ_0.
 

(1) For ϵ_0 = 0.5, find the angular correlation for different values of λ_uu = λ_uv = 0.01 – 0.8. 

(2),(3) For ϵ_0 = 0.45 and 0.6, and λ_uv = 0.5, compare the angular correlation for λ_uu= 0.01 – 0.8.

(4) For ϵ_0 = 0.5, and λ_uu = 0.5, compare the angular correlation for λ_uv= 0.01 – 0.8.




The findings are as follows:


(1) For equal length scales, we observe that the angular correlation increases with λ = λ_uu = λ_uv. The standard deviation of polarity is found to be smaller than 40 (degrees), 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 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 Fig. (<ref>a). 


(2),(3) For geometrical disorder ϵ_0 ≲0.45, and for λ_uv≃ 0.2 – 0.7 there is no detectable difference between the angular correlation of systems with different values of λ_uu≲ 0.8. Beyond this value, the correlation declines slightly, and is eventually destroyed for larger values. In order to further examine the importance λ_uu in disordered tissues, we repeat the same simulations, but for larger geometrical disorder ϵ_0=0.6. Interestingly we realize that larger disorder impedes the long-range correlation of polarity for λ_uu≲ 0.4. The correlated polarity is retained for 0.4≲λ_uu≲ 0.7 and is declined again for larger values; see Fig. (<ref>b) and (<ref>c). 


(4) Systems with nonlocal activation and local inhibition are incapable of achieving long-range polarization. We fixed ϵ_0 = 0.6 and λ_uu=0.5 and vary the inhibition length scale between λ_uv= 0.01 – 0.8. The correlation grows as λ_uv increases; see Fig. (<ref>d). An intuitive argument for why small λ_uv fails to work is as follows: for λ_uv≳ 0.8, the unlike complexes on the opposite sides of the cell inhibit each other, e.g. Fz on the right side, inhibits Vang on the left side. Thus, the Vang proteins sitting on the left side of the same cell are destabilized by Fz from across the cell. However, as far as polarity is concerned, the Vang protein on the opposite side of the cell is indeed contributing positively to the same direction of polarity as Fz. Therefore, for λ_uv≳ 0.8, the polarization, even if initially correlated over long distances, becomes highly unstable against the stochastic noise. The same logic holds true for mutual stabilizing feedback of the like complexes if λ_uu≳ 0.8, which also contributes negatively to the alignment of the polarization field. 


Comparing the charts in Fig. (<ref>), we conclude that: (1) In general geometrical disorder distorts the alignment of polarization on large length scales. (2) Equal length scales of activation and inhibition is more efficient in retaining the long-range correlation, especially in highly disordered tissues. (3) Unlike the case of equal λ_uu = λ_uv, fixing λ_uu at an intermediate value, e.g. 0.5, allows for the other length scale to exceed λ_uv≃0.7, without losing the angular correlation. Though the correlation eventually decays for larger values, λ_uv≳ 0.9.





§ ELONGATED CELL GEOMETRY

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, should we solve the full NLCI equations in elongated cells, the longer junctions will acquire larger coefficients α,β. As in the case of one dimension, searching for mean-field (MF) solutions, we assume that in 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 dimers 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 and also from the expression given in Sec. (1) of the Main Text for nonlocal interactions, it can be understood that 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. (<ref>a1), and (<ref>a2, <ref>a3) correspond to, i.e. α_1>α_2,α_3 and α_1=α_2>α_3, respectively. Using the analysis sketched in Sec. (<ref>.2), it is easy to see that in such systems the junctions with larger α,β, are more unstable towards polarization transition (remember the critical value decreases as 1/√(αβ), as α,β increase). In Fig. (<ref>a1) the polarization points toward the middle of the junction parallel to the elongation axis, whereas (a2) polarization vector passes through a vertex. Both cases of Fig. (<ref>a1) and (<ref>a2) acquire polarizations perpendicular to the elongation axis. In disordered systems (and in real systems), the cell geometry can be a combination of both types. There also exists an unstable configuration Fig. (<ref>a3) which is unpolarized. The instability is again due to nonlocal interactions which prohibits adjacent cells carrying opposite dimers; like the twofold symmetric trivial MF solutions of equilateral cells. Following the above discussion on different coefficients cooperative self-interactions α_s, and the equations derived in Sec. (<ref>.2), we calculate these coefficients as functions of L/ℓ_0, for different λ's. Fig. (<ref>b) 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 is perpendicular to the longer junctions, as depicted in Figs. (<ref>a1) and (<ref>a2).





Measure of elongation: cell nematic tensor.
Elongation is commonly parametrized using a traceless matrix ℰ. This matrix, the magnitude and the angle of elongation with respect to x-axis read:

    ℰ_i    = 
    [  ε_i,1  ε_i,2;  ε_i,2 -ε_i,1 ],
    ℰ_i    = √(ε_i,1^2 + ε_i,2^2) ,
    ϕ^e_i    = 1/2cos^-1(ε_i,1/ℰ_i).

For a tissue consisting of N_c cells, and with ε_1 = N_c^-1∑_i=1^N_cε_i,1, we get for the average elongation:

    ℰ = N_c^-1∑_i=1^N_cℰ_i ,    Φ_e = 1/2cos^-1(ε_1/ℰ).


In terms of the elongation index ℰ, we plot in Fig. (<ref>b), the angle between y-axis and the steady-state polarization vector |Φ_p-Φ_e| (in degrees), versus different values of ℰ, 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 per definition, we have |Φ_p-Φ_e|∈[0,π/2]. Of course for ℰ = 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 ℰ = 0, happens to make a small angle with y-axis (25 degrees), hence the pronounced effect of elongation. It is clear that depending on the geometry of the lattice and initial condition, the polarization can be almost perpendicular to y-axis, even for ℰ = 0. Therefore among different simulations, we chose one with a relatively large effect of elongation.



§ MUTANTS AND THE CORRESPONDING PHENOTYPES

In the Main Text we tabulated the schematics of phenotypes of three types of mutants. The schematics in the Main Text Fig. (9) 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 lack 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. (<ref>), i.e. ds^- background, non-autonomy is extended to multiple layers of cells from the clone boundaries. In 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. (<ref>), that in agreement with experiments, e.g. <cit.>, the phenotypes show almost no non-autonomy, thought distinct autonomous patterns. Lacking Dsh removes the polarity altogether, 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 geometrical irregularity. In Ref. <cit.>, Ma, et.al. suggested that clones with enhanced geometrical disorder, induced by under-expression of PTEN, are obstacles to 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. (<ref>), the system with lack of global cue shows swirly pattern of polarization. The WT polarity is retained upon adding the gradient cue.








