Mechanical Cell-to-Cell Interactions as a Regulator of Topological Defects in Planar Cell Polarity Patterns in Epithelial Tissues

A precise structural organization of epithelial cells is needed for the proper functioning and development of different tissues. The epithelial cell packing mechanism is associated with mechanical interactions between cells that place the tissue in a state with the lowest free energy. In addition, planar cell polarity has been recognized as another important mechanism for the epithelial organization that ensures the orientational ordering of the cells. Planar cell polarity is a consequence of an asymmetric distribution of certain transmembrane proteins, which is driven by specific intracellular signaling pathways. Mutations and other disruptions of these pathways were found to cause an impaired cytoarchitecture of the epithelium layer. Mutant cells with disrupted activities of signaling proteins basically represent topological defects in the polarization patterns of adjacent cells, which can lead to dysfunctions in tissue operation. Motivated by the fact that regional variations of bond tensions were found in the vicinity of mutant cells, we implement a computational model that combines the tissue development processes following the concept of free energy minimization and the cellular polarization. Our results reveal that a decrease in mechanical interactions between normal and mutant epithelium cells represents a conceivable regulatory mechanism that diminishes the impact of topological defects caused by mutant cells.


INTRODUCTION
Epithelial tissue is constituted by densely packed epithelium cells and has important structural and functional roles throughout embryonic development, morphogenesis, and adult life. The orderly packing and precise arrangement of epithelial cells are essential for separating biological compartments with different compositions by establishing selective permeability barriers between them (Nelson, 2009). Epithelial rearrangements are in general driven by forces that are transmitted via cytoskeletal elements and adhesion molecules within and between cells (Heisenberg and Bellaïche, 2013). At the level of single cells, the forces are modulated with changes of the subcellular actin-myosin flow, which is transmitted to other adjacent cells through cell-cell adhesion Heisenberg and Bellaïche, 2013). This force generating-transmitting motor is essential for the mechanical evolution of the multicellular system, which is, most importantly, governed by a surface energy minimization principle (Lecuit and Lenne, 2007;Yu and Fernandez-Gonzalez, 2017). Therefore, tissue pattern formation mechanisms are also associated with the development of specific packing geometries (Sugimura and Ishihara, 2013;Šiber and Ziherl, 2017). Optimized packing strategies were found essential for the functioning of different biological systems. For example, hexagonally packed hairs on the surface of the Drosophila wing enable a proper airflow during flight (Wootton, 1992). Similarly, a very regular crystalline-like distribution of cone photoreceptors in the retinal epithelium in adult zebrafish was found crucial for ensuring the optimal collection of light signal (Salbreux et al., 2012). In order to get a more precise understanding of the not well-known mechanisms of epithelial tissue formation, computational modeling approaches have recently been applied with considerable success in recapitulating properties of epithelial structure (Farhadifar et al., 2007;Fernandez-Gonzalez and Zallen, 2008;Hočevar and Ziherl, 2009;Salbreux et al., 2012;Fletcher et al., 2013;Osterfield et al., 2013;Kachalo et al., 2015;Chen et al., 2016;Ishihara et al., 2017;Ladan et al., 2019). The common point of all modeling endeavors is the implementation of the natural tendency of cells to adopt the most stable configuration with the lowest free energy (van Drongelen et al., 2018).
Nevertheless, the energy minimization principle alone cannot explain the diverse and complex organization of epithelial cytoarchitectures. In view of that, planar cell polarity (PCP) was recognized as an important feature that provides an additional ordering mechanisms needed to coordinate epithelium cell packing and orientation (Vladar et al., 2009;Sugimura and Ishihara, 2013;Butler and Wallingford, 2017;Fisher et al., 2017;Nissen et al., 2018;Axelrod, 2020). This plane polarity is a tissue-level phenomenon that is regulated by an anisotropic distribution of core PCP proteins involved in the PCP signaling pathway (Zallen, 2007;Marcinkevicius et al., 2009;Bosveld et al., 2012). Experimental studies performed on the wings of the Drosophila fly have shown that PCP of individual cells is established through the interplay between intracellular and intercellular communication accompanied by a global cue. Similarly, PCP-governed epithelium patterning mechanisms have also been reported for vertebrates (Goodrich and Strutt, 2011) and also for other types of tissues (Wang and Nathans, 2007). Moreover, it has been shown that in early developmental stages, the core PCP proteins are distributed symmetrically within the cells, whereas during development, the PCP proteins get increasingly more asymmetrically distributed, causing the cells to get more polarized (Strutt, 2001). As a result, the orientational order increases with age (Blankenship et al., 2006;Bayly and Axelrod, 2011;Markovič et al., 2014). For studying the complex spatio-temporal processes involved in PCP and the corresponding spatial arrangements, computational modeling has proven to be a very reliable approach (for review, see Axelrod and Tomlin, 2011;Wang et al., 2012), whereby some models implemented the PCP mechanisms in a phenomenological way (Burak and Shraiman, 2009;Hazelwood and Hancock, 2013;Zhu and Owen, 2013;Markovič et al., 2014;Shadkhoo and Mani, 2019), whereas other studies implemented and analyzed also biochemical details of the PCP pathways (Le Garrec et al., 2006).
Faults in PCP signaling caused by mutations in the intracellular signalization lead to an impaired organization of multicellular structures, which can impair normal physiological function and result in serious pathological conditions (Goodrich and Strutt, 2011;Butler and Wallingford, 2017;Markovič et al., 2017). Certain cells (so-called mutants) that lack the activity of particular transmembrane proteins provoke non-trivial and impaired polarity patterns by imposing a prescribed polarization in cells adjacent to the mutant cell boundary (Adler et al., 1997;Goodrich and Strutt, 2011;Hazelwood and Hancock, 2013). In particular, one type of mutant cells persuades neighboring cells to orient their polarity toward them, that is, attractive mutants (Gubb and García-Bellido, 1982;Vinson and Adler, 1987), whereas the other type causes the neighboring cells to orient their polarity away, that is, repulsive mutants (Taylor et al., 1998). In this way, mutant cells introduce cell polarity misorientations in their proximity, which can be regarded as topological defects that impair the normal functioning of the tissue (Hazelwood and Hancock, 2013;Markovič et al., 2014;Kim et al., 2016;Saw et al., 2017). Although mutants exhibit swirling patterns at the locus of the clone, a sufficiently strong global cue can retrieve the wildtype polarization (Ma et al., 2008;Shadkhoo and Mani, 2019). In addition, it has been shown that the polarity in mutants lacking graded global cues is susceptible to local geometric irregularities (Ma et al., 2008). Notably, local changes in the mechanical force balance within cells may, in turn, establish the regional variations in cell form and function that drive tissue patterning. It has been shown that the epithelial cells surrounding a mutant cell show a higher degree of variability in the length of contact lengths, which leads to a less efficient hexagonal repacking (Classen et al., 2005) and increased number of local defects (Kim et al., 2016). We argue that the remodeling of junctional components leading to variations in the strength of mechanical interactions between normal and mutant epithelium cells is a plausible regulatory mechanism that minimizes the distortion of the global orientational ordering caused by mutant cells. In the present study, we address this issue by developing a computational model that combines both features of epithelial pattern formation-the minimal free energy principle governing tissue morphogenesis and the PCP mechanisms. The combination of these two aspects enables us to analyze the impact of the interplay between localized mechanical energy balance and polarity patterns at different developmental stages and how the presence of mutant cells fabricates topological defects in epithelial tissues.

MATERIALS AND METHODS
In this section, we describe the computational approaches used to study spatial arrangements of cells and polarity patterns in epithelial tissues. To this purpose, we combine two different models. The first model is based on the theoretical framework of Farhadifar et al. (2007), and we use it for simulating the tissue morphogenesis under different conditions. The model simulates cell proliferation and packing by taking into account cell compressibility and junctional forces arising from cortical contractility and adhesion, as specified in the continuation. For the generated epithelial structures, we then apply the second mathematical model proposed by Hazelwood and Hancock (2013) for analyzing the evolution of the PCP patterns. The PCP model provides a phenomenological description of the PCP patterns and includes the cell's ability to maintain its intracellular polarization, the ability to interact with the polarity of adjoining cells, and the ability to sense and respond to a global cue, whereby detailed molecular interactions and processes have been integrated into a simple polarity measure. For both of the models, we utilize a Monte Carlo approach to find the ground state of the corresponding energy functions, as described in more detail in the next subsections.

Biomechanical Vertex Dynamics Model
During development of an organism, the shape of its organs and tissue structures is commonly governed by deformations of epithelial tissues. A common theoretical approach to generate in silico cellular structures of the epithelium is to consider the individual cells as 2D structures with a finite number of vertices connected into a polygon representing the cell membrane. Tissue development is then governed by reconfiguring the positions of vertices leading to intermediate mechanically stable states in which the system forces are in balance. To create a genuine epithelium structure, we make use of the vertex model described in detail in Farhadifar et al. (2007). In this model, the dynamics of the system is ruled with the following formulation of the total energy of the system: The total free energy function of the system E (Eq. 1) is governed by three terms. The first term refers to the compressibility of the cell, where K α is the compressibility parameter, A α is the actual area, and A 0 the preferred area of the cell α. In all simulations, the preferred area A 0 is considered to be equal to the surface of a hexagon with side lengths equal to l = 1, which leads to a value of A 0 = l3 √ 3/2. The second term accounts the line tensions at junctions between individual cells with the parameters ij and l ij corresponding to the line-tension strengths and the lengths of the edge between vertex i and j, respectively. The last term considers cell contractility with the contractility coefficient α and the perimeter of the polygon around a cell (L α ). The contractility coefficient α reflects contractile forces generated by myosin contraction in the cell cortex (Koride et al., 2018;Staddon et al., 2018). Because actin-myosin contractility tends to reduce the perimeter of each cell and additionally affects the line tension, smaller perimeters, and smaller cell boundary lengths result in a state with lower free energy. The main parameters of Eq. (1) are schematically shown in Figure 1. In general, the following parameter values will be used in our simulations: K α = 14, ij = 1, and α = 0.85, unless stated otherwise. Such a parameter set favors a hexagonal-like packing and thereby mimics well the circumstances in realistic epithelial tissues (Farhadifar et al., 2007). We are interested in alternations of global patterns of cells caused by variations in the strength of mechanical interaction among junctions between normal and mutant cells, as an abnormal polarity of junctional complexes around the edge of a mutant propagates to neighboring cells and affects cellular packing (Goodrich and Strutt, 2011;Kim et al., 2016). To this purpose, we numerically simulate the evolution (growth) of epithelial tissue around the mutant cell. The initial configuration of the simulation is a 2D lattice of 25 hexagonal cells, which have side lengths equal to 1. The cell in the center is defined as a mutant cell. In continuation, a non-mutant cell is selected, and its area is increased in a stepwise manner until its area is doubled. This is achieved by moving the vertices of the selected cell radially outward. The increment equals 1% of the distance between the corresponding vertex and the cell's center of mass. After each increment of the cell's area, we let the whole system relax toward its mechanically stable state with the lowest free energy defined by Eq. (1). When the area of a given cell is doubled in comparison with the initial area, the cell is divided. This is achieved by generating a line with a random slope through the center of the cell, which is determined as the center of mass of all its vertices. At the two intersection points of the line with the edges of the cell, two new vertices are created and connected with a new edge. This edge separates the two cells, which have been spawned from the initially selected cell. Finally, the system relaxes toward its mechanically stable state, as described by Sugimura and Ishihara (2013). The Metropolis Monte Carlo method is used to simulate the relaxation process. In each iteration step, individual vertices are randomly displaced within a circle centered around the selected vertex with a radius of 0.04. After each displacement, the total free energy of the system is computed and compared with the overall free energy of the system before the displacement. If the new energy (E new ) is lower the old energy (E old ), the rearrangement is accepted, and the processes are repeated. If, on the other hand, the old energy is lower than the new one, the change is accepted with a probability e −(E new −E old )/N , where N equals the number of cells. If the average absolute change in the system's total free energy after 1,000 successful Monte Carlo steps is less than 10 −8 , the Metropolis Monte Carlo algorithm is stopped, and the current configuration of the system is regarded as the final state with the lowest free energy resembling a mechanically stable state. Following these steps, this process is repeated until there are 100 cells. In our simulations, open boundary conditions are considered. It should be noted that the shape of the mutant cell is allowed to adjust its shape and is involved in the energy minimization process but is never divided.
To visualize the impact of the line-tension paremeter mn on the spatial on the spatial packing of epithelium cells, we show illustrative calculations performed on a 2D grid of 25 cells, with the mutant cell located in the middle of the grid. The initial configuration of this simulation was a lattice with 25 hexagonal cells. All cells had equal side lengths l ij = 1.0. To show the robustness of the Metropolis Monte Carlo algorithm, the initial configuration has been randomly disrupted. The line-tension strengths for edges between any non-mutant cell pair were set to ij = 1.0. In Figure 2, we present the resulting epithelium structures that have relaxed to their mechanically stable state for three different values of line-tension strengths mn . We can observe that an almost perfect hexagonal grid is attained when mn = 1 (all interactions in the tissue are the same), whereas the other two scenarios ( mn = 0 and mn = 2) lead to accordingly adjusted cellular structures.
In all further simulations, where the PCP model was implemented on the top of tissue morphogenesis simulation, a lattice of 25 equal hexagonal cells was considered as our initial configuration. The cell in the center of this grid was assigned as a mutant cell. Tissue growth was then simulated as described above. We generate epithelium structures for different line-tension strengths ( mn ∈ [0, 2]). The obtained epithelial structures serve then as a basis for the study on how different scenarios of structural deviations, caused by the mutant cell, can affect the overall PCP pattern during tissue development.

Phenomenological Model of Planar Cell Polarity
In order to mathematically describe the asymmetric distribution of the six core PCP proteins, we use a phenomenological model introduced by Hazelwood and Hancock (2013), which implements the three-tiered mechanism for regulation of PCP (Tree et al., 2002). The model is based on the so-called Ginzburg-Landau method, commonly used to study ordering processes in condensed matter physics. As such, the model describes the system via a free energy function, which has the following form: During the simulation, we update the polarization of individual cells, which is given by the vector m i , where the index i indicates that the polarity vector is assigned to the ith cell. The total free energy of the systems F tot is composed of three parts. The first part, F c , accounts for the ability of individual cells to maintain a certain degree of polarity, thereby mimicking the asymmetric distribution of proximal and distal proteins. Here, the preferred inherent polarization magnitude of an individual cell is given by m A,i . The second part, F e , captures intercellular physical interactions among adjacent cells, n i , which try to impose a common direction of polarization. This reflects the maximization of intermolecular forces between cells due to ligand binding, as deviations from a uniform alignment lead to an increase in free energy. In this term, we additionally take into account that the orientation order among cells i and j, where j is an adjacent cell of cell i and is governed by the length of the junction between them, l ij . Accordingly, K ij is defined as follows: where K represents the average interaction strength among adjacent cells and l i is the average junction length of the ith cell. The latter is computed by averaging the edge length of the polygon around the ith cell. Lastly, the third part, F f , couples the polarization of individual cells with a global cue, G. The contribution of the three individual terms to the total free energy F tot is weighted with parameters B, K, and C, which reflect the contribution of intrinsic polarization, intercellular interactions, and interactions with the global cue, respectively. The main advantage of this model is that it is equipped with a small number of biologically relevant parameters thereby facilitating a qualitative assessment and comparison with experimental observations (Vu et al., 2019). We consider two typical sets of parameters B, K, and C corresponding to the two different stages in the development of epithelium (Markovič et al., 2014). In the early developmental stage, the polarization shows a weak tendency of ordering and can be modeled by setting B = K = C. In a further developed stage, the global PCP patterns are found to be practically complete aligned with the direction of the global cue G (Bayly and Axelrod, 2011;Hazelwood and Hancock, 2013). In this case, we set the model parameters to B = K and C = 10 B, indicating a stronger contribution of global alignment to the total Frontiers in Materials | www.frontiersin.org free energy of the system. We use this parameter set because in early developmental stages on a tissue level no uniform polarization of orientation is established. The latter embraces a lower contribution of the global cue signal compared with the terms B and K (Strutt, 2001). At the final stage of PCP development, the polarization direction of cells is very uniformly oriented in the same direction, which highlights that the energy term F f is contributing more to the total free energy F tot than the other two terms F c and F e . This in turn is considered by using larger values for the parameter C with regard to the parameters B and K (Strutt, 2001;Hazelwood and Hancock, 2013;Markovič et al., 2014). Therefore, the main two parameters, which determine if local or global ordering will be preferred, are the parameters K and C. If C < K, local interactions are considered to contribute more to the total free energy function compared with the contribution of global ordering. Consequently, local orientation is dominating the evolution of the system. On the other hand, if K < C, the contribution of the global ordering term is prescribed a higher weight in the energy function, and a more uniform polarization pattern is formed. The parameter, B, highlights the tendency of a cell to maintain a predefined inherent polarization. Lastly, it should be noted that only the relative relations of parameters B, K, and C are important, as they proportionally contribute to the free energy term F tot and thereby determine which energy term is the main driving force in the system. In all calculations, the values of parameters B and K equal to 1.0, and only C was varied. In order to graphically represent how well an individual cell is aligned along with the global cue, we use a color code system as presented in Figure 3.
In all simulations, mutant cells are located in the middle of the tissue. We consider only attractive mutant cells, which means that the polarization of the neighboring cells always tends to be aligned radially inward with a magnitude equal to the cells preferred magnitude of polarization m A,i . This mimics the lack of Fz in mutant cells, which induces neighboring cells to point their polarization toward the mutant (Goodrich and Strutt, 2011;Hazelwood and Hancock, 2013;Markovič et al., 2014). For all other cells, the initial orientation θ i and the magnitude of polarization | m i | are prescribed randomly. The random number interval used for the two variables is θ i ∈ (−π/2, −π/2) and | m i | ∈ (0.5, 1.5). The preferred polarization magnitude of individual cells is set to m A,i = 1. A generic Metropolis Monte Carlo simulation is then used to find the stationary state of the system, signified by the global minimum of the total free energy function given in Eq. (2). In each Monte Carlo step, one cell is randomly selected, and a new magnitude and orientation of polarization for the chosen cell are randomly generated. If the reformed configuration decreases the total free energy of the system F tot , the change in m i and θ i is accepted. In order to avoid local minima, we implemented the Metropolis algorithm, which enables us to explore the energy landscape of the system and accept a given modification with a probability proportional to the corresponding Boltzmann weight. In turn, this allows us to move up or down on the energy landscape and to escape local minima in the space of arrangements (Frenkel, 2002). To enhance performance and prevent a bias in preferring only small or only large displacements, we used an adaptive step size for angle adjustments during this energy minimization protocol.
This process is repeated until the absolute average change in F tot over 100 successful Monte Carlo steps exceed the values of 10 −8 . Finally, we quantify the degree of global alignment of the tissue polarization with the order parameter r defined as, where N equals the number of cells and θ i is the angle of polarization of the ith cell relative to the symmetry breaking direction, specified by the global cue vector, G (Markovič et al., 2014).

RESULTS
We investigate the impact of the strength of cell-cell adhesion, mn , between the mutant cell and its neighbors on the global polarity pattern. To this purpose, we generated several epithelial tissue architectures for different values of line-tension strengths around the mutant cell ( mn ∈ [0, 2]), as described in Biomechanical Vertex Dynamics Model. The interaction strength between all other cells was set to ij = 1. It should be noted that studying the role of the intercellular interaction strengths between normal and mutant cells and how they affect cellular arrangement was motivated by previous endeavors, which have suggested that intercellular contacts around mutants are modified (Goodrich and Strutt, 2011). Then, we used each of the generated epithelial structures to evaluate the impact of the defect incited by the mutant cell in the center of the global polarity pattern. Results presented in Figure 4 show six characteristic stationary polarity patterns for three values of mn and two different values of the global cue strength C. Different degrees of the global cue C correspond to different epithelium developmental stages. In the early developmental stages, the global polarization order is known to be lower (C = 1) than in the late stage, where the polarity on the global scale is much more pronounced (C = 10) (Markovič et al., 2014). It can be observed that the extent of the topological defect is lower for higher values of C. Modifications of the second parameter, mn , correspond to variations in strength of mechanical interactions between normal and mutant epithelium cells. For mn = 0, the cell-to-cell interactions do not contribute to the energy function [see Eq. (1)], and hence the constraint in the length of the cell boundary is not pronounced (see Figure 2). On the other hand, higher degrees of mechanical interactions (for example, mn = 2) provoke the opposite effect, thereby forcing the neighbors of the mutant cell to reduce their boundary lengths. Notably, it can be observed than for higher line-tension strength mn values, the disruptions of the global polarity patterns are more pronounced than for lower values of mn . Apparently, a proper adjustment of the intercellular interactions along the junctional regions around the mutants can serve as a regulatory mechanism that reduces distortions of the orientational disorder provoked by the mutant cells.
To quantify this visual assessment, we compute the average order parameter r as a function of the line-tension strength mn , for two different values of the global cue strength (C = 1 and C = 10). The results are presented in Figure 5A and were obtained by averaging over 100 realizations of the functional PCP model for individual line-tension strengths (see Biomechanical Vertex Dynamics Model) in order to ensure statistical accuracy. As expected, the overall degree of global alignment is higher for higher values of the global cue C, which is in agreement with previous studies (Axelrod and Tomlin, 2011;Hazelwood and Hancock, 2013;Markovič et al., 2014;Butler and Wallingford, 2017). Most importantly, for both values of C, the order parameter r is a decreasing function of the line-tension strength mn . This result confirms our abovementioned predictions, because reducing the costs for intercellular interactions on cell boundaries among cells that surround a mutant cell does indeed ensures a better global alignment of PCP.
To further investigate the effect of altered mechanical interactions around the mutant cell, we calculate the standard deviation in the lengths of the cell boundaries of the mutant's neighbors l mn . In Figure 5B, we show this standard deviation values std(l mn ) as a function of the line tension mn . The results point out a non-trivial behavior. Minimal deviations of the lengths l mn are attained when the cell-to-cell interaction in the whole system are homogeneous ( mn = ij = 1). It should be noted that the simulated tissue growth process is governed by a stochastic step. As a result, the final packing geometries are characterized by different cell shapes, despite the fact that the initial configuration consisted of regular hexagons. Most importantly, such heterogeneous cellular architectures are also found in realistic epithelial tissues (Farhadifar et al., 2007). For that reason, the final value of std(l mn ) is not 0 even in the most homogeneous scenario, where mn = ij = 1. Either The orientations of cell polarities are color coded as introduced in Figure 3. In all cases, the global cue is oriented in the positive direction of the x-axis: G= (1, 0). positive or negative variations of mn lead to higher variability in the boundary lengths. Notably, differences in intercellular interactions around mutants in comparison with other epithelial cells have also been observed experimentally (Classen et al., 2005), which supports our ideas on mechanical regulations of impaired polarity patterns.
By decreasing the line-tension strength mn , the neighbors of the mutant cell tend to magnify the lengths of cell boundaries among them and reduce the boundary lengths with other cells in their vicinity. This not only gives rise to the variations in lengths but also causes the formation of a kind of structural barrier, which tends to isolate the mutant cell from its environment (see Figures 2, 4). As a result, the impact of the mutant via PCP interactions on the whole epithelial tissue decreases [see Eq. (6)], and hence, the global polarity patterns become less impaired. On the other hand, an increase of the line tension strength mn causes the opposite effect, which in turn even more pronounces the topological defect provoked by the mutant cell.

DISCUSSION
In this study, we have investigated how mutant cells affect the structural and functional organization of epithelial tissues. Because these mutant cells directly affect the orientation of polarization (Hazelwood and Hancock, 2013) and the structural properties (Classen et al., 2005) of its adjacent cells, the global properties of the developed epithelium depend on both effects. Therefore, we have combined two modeling approaches that independently one from another describe the mechanisms of mechanical development (Farhadifar et al., 2007) and the establishment of PCP in epithelial tissues (Hazelwood and Hancock, 2013). This enables us to study the interplay between mechanical interactions and the corresponding established global PCP patterns in growing cellular tissues, which is established by cell growth and cell division around a mutant cell. Special attention is given to the mutant cellinduced irregularities in polarized cytoarchitectures. Motivated by previous experimental findings, we have systematically investigated mechanical interactions at cell borders between healthy and mutant cells by varying the line-tension strength between these cells. Namely, previous studies have implied that intercellular contacts between normal and mutant cells are modified, thereby leading to local variations and structural disorder (Classen et al., 2005;Vladar et al., 2009;Salbreux et al., 2012;Heisenberg and Bellaïche, 2013;Kim et al., 2016;Butler and Wallingford, 2017). More precisely, mutant clones disturb the packing of adjacent wild-type cells, which are not able to assemble efficient boundaries with new neighbors and cannot regularize their packing geometry (Classen et al., 2005). In zebrafish mutants, adherens junctions are destabilized, and the apicobasal polarity of the retinal epithelium is destroyed (Salbreux et al., 2012). Removal of adherens junctions has also been shown to result in a substantial rounding of the cells caused by interfacial tension release (Brückner and Janshoff, 2018). It should be noted that even though the interplay between adherens junctions and PCP is incompletely understood, recent findings suggest that they are intimately linked (Tatin et al., 2013;Stocker and Chenn, 2015). In addition, in epithelial cell populations in the Drosophila wing, it has been shown that cells lacking Fz or Ds function cause trichomes proximal to the mutant to invert their polarity (Goodrich and Strutt, 2011). The functional PCP model (Hazelwood and Hancock, 2013) enabled us to investigate how a mutant cell affects global polarization in an early and final developmental stage (see Figure 4). At an earlier developmental stage, local interactions prevail (Strutt, 2001). In turn, this causes the defect in polarization initiated by the mutant cell to affect a wider area compared with the final developmental stage, where the global cue is more strongly expressed. The effect of the strength of the global cue has been shown to affect the extent of such swirling polarization patterns in such a way that increasing global cue values lower the degree of swirling (Ma et al., 2008;Burak and Shraiman, 2009;Shadkhoo and Mani, 2019).
Our findings are in agreement with the abovementioned phenomena. For example, by changing the relative contribution of the global cue strength in the PCP model, we mimic the effect of a stronger or weaker contribution of the global cue on the global polarization pattern (Ma et al., 2008;Burak and Shraiman, 2009;Shadkhoo and Mani, 2019). By adjusting the line-tension strength between the edges of a mutant and wildtype cell, we can further consider how changes in mechanical cellto-cell interactions affect the cytoarchitecture of the numerically evolved epithelium and the corresponding PCP patterns. We have computed and analyzed several epithelial structures for different sets of model parameters (see Figure 4). Our results revealed that the line-tension strength between mutant and wildtype cells does not only affect the evolved global structure but additionally influences the established global PCP pattern. This link is predominately a consequence of the increasing number of wild-type cells around a mutant cell, with higher nm values. The latter is caused by the fact that with increasing values of nm , wild-type cells tend to minimize the cell boundary length common to the mutant cell and maximize the cell boundary length with adjacent wild-type cells. In turn, this leads to an accumulation of wild-type cells surrounding a mutant cell. The neighbors of those cells mediate the improper polarity pattern to a larger proportion of the tissue, causing a lower degree of global ordering. This tendency can be observed for weak and strong global cues strengths (see Figure 5A). Regarding the structural features of cells in the vicinity of a mutant cell, our results show that the variability in the lengths of the cell sites around the mutant cell increases as the line-tension strength is either decreased or increased from the reference value of nm = 1 (see Figure 5B). This finding is in agreement with the experimental observations, which indicate that the cells surrounding a mutant cell are characterized by higher variability in their cell site lengths than are the other cells in the tissue (Classen et al., 2005). Previous reports have also suggested that modifying local strain cues can accelerate the reorientation of cells with respect to some static strains (Rens and Merks, 2017). This is a promising result indicating possibilities for regulating defects in epithelial structures and might have important therapeutic applications. It can be imagined that the defects could be regulated by tuning the line tension strength between the cells. In particular, a pharmacologically induced decrease in the strength of mechanical interactions between normal and mutant epithelium cells could act as an efficient regulatory mechanism for diminishing the topological defects caused by mutant cells. Because the current understanding of the interplay between the mechanical properties of the epithelium structure and the inherent polarity of the single cells is still limited, additional experimental as well as theoretical studies are needed. Our phenomenological framework allows only for a qualitative assessment of the processes occurring during tissue morphogenesis and the development of PCP, whereas the details about the molecular processes and interactions, which drive polarity transmission or mediate the orientation with the global cue, remained unaddressed. Most importantly, as also suggested by others (Minc and Piel, 2012), we need to integrate the mechanical properties and the PCP of the epithelium as two developmental driving processes that mutually affect each other, leading to asymmetric growth and shape-dependent cellular orientations.

DATA AVAILABILITY STATEMENT
All datasets generated for this study are included in the article/supplementary material.

AUTHOR CONTRIBUTIONS
RM, MG, and MM designed the study and wrote the manuscript. RM performed numerical simulations and analyses. All authors contributed to the article and approved the submitted version.

FUNDING
This work was supported by the Slovenian Research Agency (research core funding, nos. P1-0055 and P3-0396, as well as research project no. I0-0029).