A reduced cell-based phase model for tissue polarity alignment through global anisotropic cues

Ordered polarity alignment of cell populations plays vital roles in biology, such as in hair follicle alignment and asymmetric cell division. Although cell polarity is uniformly oriented along a tissue axis in many tissues, its mechanism is not well understood. In this paper, we propose a theoretical framework to understand the generic dynamical properties of polarity alignment in interacting cellular units, where each cell is described by a reaction–diffusion system, and the cells further interact with one another through the contacting surfaces between them. Using a perturbation method under the assumption of weak coupling between cells, we derive a reduced model in which polarity of each cell is described by only one variable. Essential dynamical properties including the effects of cell shape, coupling heterogeneity, external signal and noise can be clarified analytically. In particular, we show that the anisotropicity of the system, such as oriented cell elongation and axial asymmetry in the coupling strength, can serve as a global cue that drives the uniform orientation of cell polarity along a certain axis. Our study bridges the gap between detailed and phenomenological models, and it is expected to facilitate the study of polarity dynamics in various nonequilibrium systems.

where φ i is the phase of cell i, which approximately describes the position of the maximum of a reaction-diffusion component, A(i) is the group of cells adjacent to cell i, η ij is the cell-to-cell direction from cell i to cell j and a ij and c ij are functions of the width d ij of the contacting surfaces (see Figs 1 and 2). In the absence of the second term, this phase model is same as the model describing spin states in ferromagnets, known as the XY model, and a special case of spatially extended phase oscillator models [16][17][18] . Such a model was employed in a previous study on PCP 10 . Our phase model is distinct from the study, as it includes novel terms representing geometric information, such as the cell shape and the relative position between neighbouring cells. As the model is easily manageable, essential dynamical properties including the effects of cell shape, external signal and noise can be clarified analytically, which have only been studied numerically in previous works using detailed models 6,8,9 . Finally, we discuss symmetry-breaking patterns in PCP [5][6][7] by using of our model. In particular, we point out that axial asymmetry in the system, such as oriented cell elongation and asymmetric distribution of coupling strength, can be a global cue for the orientation of cell polarity across the entire tissue. Our study bridges the gap between detailed and phenomenological models, it is expected to facilitate the study of polarity dynamics in various nonequilibrium systems.

Reaction-diffusion model
The entire system is composed of a population of planar cells aligned in two-dimensional space. Reaction-diffusion dynamics of each cell takes place on the one-dimensional surface, and the cells further interact with one another through the contacting surfaces between them. For simplicity, it is assumed that every cell has identical shape, which is either regular or elongated hexagonal with a perimeter of 2π, and they form hexagonal lattices, as shown in Fig. 1. Each cell obeys the equation,   Fig. 1(b), and H ij (θ i , t) vanishes if cell i does not contact cell j at θ i . As it is described later, external signals and noise may also be considered. In each cell, X i (θ i , t) is assumed to form a unimodal distribution for ε = 0; i.e., polarity is spontaneously formed. The orientation of polarity of cell i at time t is defined by the θ i value at which the first component of Here, some concrete examples are provided. For F and D , two examples are considered: (a) the real Ginzburg-Landau equation (GLE) 19 and (b) the activator-inhibitor model 20 . Both models have two variables, denoted by X i = (U i , V i ), and F and D can be written as where D 0 is set to 0.3 at which a stable unimodal distribution is obtained, and The parameter values for the latter model are taken from the reference 20 . The former is a long-wave amplitude equation, which is widely used to describe various systems near the onset of instability. The latter is a reaction-diffusion model, describing biological pattern formation 20 . In these models, given appropriate initial conditions, X i exhibits a stationary unimodal distribution for ε = 0, thus, they are suitable as dynamical models describing cell polarity. Figure 2(a,b) shows a steady profile of U i (θ i ,t) for ε = 0 numerically obtained using the activator-inhibitor model given by equation (5). As a simple example of intercellular interaction, we consider a linear coupling given by where S ij = 1 if cell i faces cell j at θ i and S ij = 0 otherwise; i.e., The coupling given by equation (6) acts as mutual inhibition between neighbouring cells through the U-component for ε > 0, causing polarity ordering, as shown in Fig. 2(c,d). Later, another type of linear coupling, i.e., θ − , is considered to show the robustness of our results.

Derivation of the phase model
Our reaction-diffusion model given by equation (2) consists of N × M partial differential equations, where N and M are the numbers of cells and variables in each cell. For such a model, both analytical and numerical treatments are difficult. Therefore, we applied a perturbation method to equation (2) under the assumption of weak coupling to obtain a phase model, which consists of N ordinary differential equations and can be useful for both analytical and numerical analyses. Our method is based on the well-known phase reduction theory 13 , and it is an application of the recently developed method for oscillatory patterns reported in refs. 21,22 . Let X S (θ) be the stationary distribution of a cell in the unperturbed system (ε = 0). Because of the translational symmetry, X S (θ − θ 0 ) with any constant θ 0 is also a steady solution. The phase φ i (t) of X i (θ i , t) is defined such that X i (θ i , t) converges to X S (θ i − φ) as t → ∞ in the unperturbed system. In other words, y i (θ i , t) → 0 as t → ∞ for ε = 0, where the deviation y i (θ i , t) is defined by 0 S These eigenfunctions are assumed to form a complete orthonormal system and are normalised as m m The deviation y i can be expanded as Substituting equation (8) into equation (2), we obtain Taking the inner product with Z 0 (θ i − φ i ) and omitting O(ε 2 ), we finally obtain the phase model as . Given the functional forms of X S (θ) and Z 0 (θ), equation (17) provides a closed equation for the phases It is convenient to express Γ ij in terms of Fourier coefficients, are even, even and odd functions, respectively. By substituting these expansions into equation (18) with H ij given by equation (6), we obtain a general expression: For the regular and elongated hexagonal cell shapes shown in Fig. 1 ij . The coefficients u k and z k can be obtained for a given model.
For the GLE, the phase reduction is analytically performed. By solving Functions X S and Z 0 are shown in Fig. 3. Therefore, equation (17) with equation (22) reduces to In the GLE Z 0 is proportional to Y 0 because the linear operator is self-adjoint, i.e. = †  , in this particular model. Expressions for all other eigenfunctions are known 23 , although we only need the expressions for Z 0 and Y 0 here.
For most models, phase reduction is performed numerically by solving equation (2) for ε = 0 and its adjoint equation 22 . For the activator-inhibitor model, U S and Z U 0 ( ) are obtained, as shown in Fig. 3(b). Their Fourier coefficients are given approximately as u 0 = 0.925, u 1 = 0.397, u 2 = 0.065, z 1 = −0.180, z 2 = −0.062, and the rest of coefficients are negligibly small. We obtain Γ ij by substituting these values into equation (22). Then, our phase model becomes As shown in Fig. 4, we confirmed the accuracy of our reduction theory for both the GLE and the activatorinhibitor models by comparing the time series of the original model given by equation (2) and that of the phase model given by equation (17) with the corresponding Γ ij .
It should be noted that the phase sensitivity function Z 0 (θ) is very useful in understanding the response of the orientation of polarity to perturbation. See Fig. 3(a) as an example. If the U variable is perturbed upwards at θ = π/2, φ increases because π > Z ( /2) 0 0 , i.e., the pattern eventually shifts to the right.

General properties of the phase model
We focus on the phase model given by equation (27) below for the following considerations. If U S (θ) and θ Z ( ) U 0 ( ) are nearly harmonic, i.e., u k and z k with k ≥ 2 are small, an approximation for equation (27) can be obtained with This is actually the case in our activator-inhibitor model: In the coupling function derived from the activatorinhibitor model given by equation (29) the first three terms are considerably larger than the other terms. Thus, the dynamical properties of the activator-inhibitor model are expected to be similar to those of the phase model in equation (27).
We consider different cell alignments and investigate the effect of cell alignments and boundaries on the existence and stability of polarity patterns. Moreover, we investigate the effect of noise and external signals. Finally, we discuss the robustness of our results. Henceforth, without the loss of generality we set ε = 1 in numerical simulations.
Straight cell alignments. We first consider two coupled cells aligned horizontally, as shown in Fig. 2(d).

Effects of cell alignments, cell shapes and heterogeneity in coupling strengths.
We investigated the dependence of polarity pattern on complex cell alignments and cell shapes, as well as the heterogeneity in coupling strengths. It should be pointed out that in equation (27), the second and third terms facilitate the phase φ i and the mean phase φ φ + 2 i j to be oriented to the cell-to-cell direction η ij , respectively. If only the first term is present in equation (27), which is the case in the XY model, there is a family of stable solutions φ i = φ * (i = 1, …, N) with arbitrary φ * values, and the realised polarity pattern is determined by the initial conditions. However, if either b ij or c ij is nonvanishing, the in-phase state even with a particular φ * value does not exist except for special networks such as a straight chain.
To obtain useful insight into the dynamical behaviour of a complicated alignment of cells, we made an approximation in the phase model using the assumption that the neighbouring cells are nearly in phase. Under the approximation that φ i = φ j for any neighbouring cells (i.e., the in-phase state), equation (27) reduces to which can be interpreted as the effective strength and the preferred direction of the net interaction of cell i, respectively. For hexagonal lattices with each cell shape being regular hexagonal, R i vanishes for cell i that does not facing boundaries of the lattice because b ij and c ij are not depending on i, j and η ij takes the values 0, 2π/n, 4π/n, …, 2(n − 1)π/n with n = 6. The same is true for square lattices with each cell shape being square. Nevertheless, for cells at the boundary, R i is non-vanishing and η i it is approximately parallel to the boundary line. Therefore, cell polarity at the boundary tends to be parallel to the boundary line and cell polarity of bulk cells is smoothly aligned to that of the neighbouring cells. As shown in Fig. 5, this prediction is confirmed using the system with winding cell alignment. When the cell shape is elongated, R i is non-vanishing even in the bulk. In this case, η i tends to orient to the direction of a contact surface with a larger width. When the number of bulk units is much more than that of boundary units, polarity orientation is dominantly depending on the cell shape. In particular, when the cell shape is uniformly elongated as shown in Fig. 1(a), stability analysis is straightforward. In this case, equation (36) reduced to  Thus, stability depends on the sign of λ. For λ > 0 (λ < 0), which is the case for for all i are stable. For λ = 0, which is the case for = π d 3 , φ i dynamics becomes neutral, and the steady state is determined by the initial conditions. To clearly demonstrate the effect of cell elongation, a two-dimensional periodic system was considered, as shown in Fig. 6. Initially, the cells were set to be regular hexagonal. Because the boundary effects are negligibly small for the periodic boundary condition under consideration, phases can be aligned with arbitrary values determined by the initial conditions. Here, a random initial condition was used, where phases were chosen from a uniform distribution within the range (−0.5, 0.5). At t = 100, phases were almost perfectly aligned at φ i ≈ 0, which is approximately the average of the initial phases. At t = 500, the cell shape was changed to be elongated. Then, the cell polarity was aligned upwards, pointing to the direction of a contact surface with a larger width, as predicted above. This polarity pattern was maintained even when the cell shape was returned to be regular hexagonal (t > 1000).
A similar result can be obtained by considering heterogeneity in coupling strength even for regular hexagonal cell shapes. We considered the condition in which coupling strength ε in our reaction-diffusion model given by equation (2) is dependent on i, j, by replacing ε with ε(1 + α ij ). Then, the corresponding phase model becomes 3 8 and = c 1 12 . By assuming an in-phase state, this equation reduces to Now we introduce an axial asymmetry such that only the surfaces along the vertical axis, which are shown as bold lines in Fig. 7, have α ij = α, and α ij = 0 for other surfaces. In this case, we further obtain Thus, the sign of α plays the exactly same role as that of λ in equation (38); the polarity pattern is aligned along the axis with stronger coupling. In Fig. 7, numerical results can be seen obtained using our phase model with coupling heterogeneity, given by equation (40). To emphasise the effects of geometry-dependent terms, we also show results obtained using equation (40) with b and c being set to zero, corresponding to the XY model, which indeed shows no response to the coupling heterogeneity. The axial asymmetry of the system under consideration affects the dynamics only in the presence of geometric-dependent terms.
Effects of noise and external signals. Phase reduction is also possible when our reaction-diffusion model includes external signals and noise, given as with Q ≥ 0 and  Φ ∈ is shown. Cells are elongated only for 500 ≤ t < 1000 during which δ = − π π t ( ) 3 1 0 ; otherwise δ = π t ( ) 3 .
In (c) snapshots are presented. Initial conditions were chosen such that no topological defects appeared (see Methods).
Thus, we obtained a probability distribution where C is the normalisation constant. As shown in Fig. 8, the probability distribution obtained numerically from the GLE, given by equations (43) and (4), is in good agreement with equation (53).

Robustness.
We discuss the robustness of our results described above against changes in our model equations. Our numerical simulations for regular hexagonal cell shapes were performed using equation (27) . We verified that these results do not change qualitatively for small changes in a ij , b ij , c ij values. Qualitative change is certainly expected when the stability condition given by equation (35) is disturbed.
We also observed that there is no qualitative difference between the phase models reduced from the GLE and the activator-inhibitor model, as shown in Fig. 4. This suggests that higher harmonics in Γ ij does not considerably affect dynamics, at least when they are small.
We can also consider different types of couplings in our reaction-diffusion model other than that in equation (6). For example, we considered ij i i j i j j which describes mutual inhibition, as well as equation (6). By assuming U S (θ) and θ Z ( ) U 0 ( ) are nearly harmonic, we obtain the following approximation . The phase model given by (55) is again a gradient system with the potential function given by The second term in equation (55) differs from that in equation (27), thus different dynamical properties may appear. However, in case of GLE given by equation (4), we obtain Thus, the corre- sponding phase model is a special case of (27) where b ij = 0. It is straightforward to confirm that the existence and stability analysis performed above does not change in this case.

Discussion and Conclusion
A theoretical framework for understanding the general dynamical properties of the alignment process of cellular polarity has been proposed. We derived the phase model as a reduced model of coupled reaction-diffusion systems and investigated polarity dynamics using the phase model. The first term of our phase model facilitates the polarity ordering between interacting cells, which is the same as the in-phase coupling in the Kuramoto model and the ferromagnetic coupling in the XY model 18 . The remaining terms include geometric information of the system, i.e., cell-to-cell directions between neighbouring cells and the width of contacting surfaces. Therefore, our phase model exhibits polarity dynamics that depends on the shape of individual cells and the alignment of cell populations, as well as heterogeneity in coupling strengths. In particular, we show that the axial asymmetry of the system facilitates the formation of globally oriented polarity patterns. The advantages of our method are substantial. Whereas our reaction-diffusion model given in equation (2) is an N-set of reaction-diffusion systems with multiple variables, this complicated system is reduced to an N-dimensional system of ordinary differential equations as given in equation (17). Using the steady state of the unperturbed reaction-diffusion system, we obtain the coupling function Γ ij , by which we can perform various analytical and numerical analyses, which are presented in this paper. Although the studied phenomena are nonlinear, our framework enables us to obtain analytical results even in the presence of noise.
Finally, we discuss the symmetry-breaking patterns of cell polarity in biological tissues [5][6][7] with the help of our model. As reviewed by Aw and Devenport 5 , although PCP aligns over long distances in skin and wing, the global cues that orient tissue polarity are not well understood. This review highlights two plausible choices. One is a factor expressed in tissue-wide gradients along the axis of polarity, supported by experimental and theoretical studies [see list of references in 5 ]. The other is mechanical tension applied to the tissues, which may act over long distances. Aw et al. recently reported that in mouse skin, an axial asymmetry in a PCP component (Celsr1, an atypical cadherin) emerges during the process of mechanical deformation along the anterior-posterior (AP) axis; i.e., the concentration of Celsr1 is higher on the junctions perpendicular to the AP axis. They demonstrated that such Celsr1 asymmetry emerges spontaneously during neighbour exchange, because the junctions perpendicular to the AP axis are persistent and there is sufficient time for Celsr1 to accumulate on those junctions, whereas other junctions are nascent. They speculated that such axial asymmetry contributes to the formation of polar asymmetry, as indeed developed in the skin. As also reviewed by Aw and Devenport 5 , in the Drosophila wing, a similar axial asymmetry in PCP components is formed during the process of mechanical deformation, and cell polarity is eventually aligned along the AP axis 6 . Our model can provide an understanding of how the axial asymmetry in the system contributes to the formation of globally aligned patterns of cell polarity. Because PCP components are essential to cell-to-cell communication for polarity alignment, the concentrations of PCP components can naturally be associated with the coupling strength; i.e., coupling is expected to be stronger on cell-to-cell junctions with higher concentration of PCP components. Therefore, we interpret the vertical axis of Fig. 7(b) as the AP axis, with which polarity is eventually aligned. This symmetry-breaking phenomenon does not emerge in the XY model, i.e., equation (27) with b ij = 0 and c ij = 0, because in that case, the model has rotational symmetry even when asymmetry in the coupling strength is considered. The symmetry-breaking phenomenon occurs in our phase model because it has geometry-dependent interaction terms originating from geometry-dependent interactions in our reaction-diffusion model. Further investigation on the robustness of global alignment against randomness in cell shapes and coupling strengths is required. For this aim, we need to extend our theoretical framework to arbitrary cell shapes, and investigations along this line are in progress.

Methods
Numerical implementation and visualisation. In Fig. 2, we obtained U i (θ, t) in steady state by numerically solving the equations of the reaction-diffusion model given by equations (2) and (5). The arrows indicate θ i at which U i (θ i , t) is maximal. Figure 3(a) shows the analytical results, given by equations (24) and (26), whereas Fig. 3(b) shows the numerical results obtained using the activator-inhibitor model, given by equations (2) and (5), for ε = 0 and its adjoint equation Figure 4 shows the numerical results: Solid lines were obtained from the reaction-diffusion models, given by equations (2), (4), and (5). Symbols were obtained from the phase models, given by equations (27) and (29). For the reaction-diffusion models, the phase φ i (t) was numerically determined as follows. Firstly, the first Fourier component of U i (θ,t) was calculated as Then, the phase φ i (t) was given as the solution to where C(t) ≥ 0 and φ i (t) are real. In this way, θ i = φ i (t) approximately coincides with the maximum of U i (θ i ) when U i (θ i ) is nearly harmonic because i i i i Figures 5-7 show the numerical results obtained by the phase model. Initial phases were taken randomly from a uniform distribution − − . . in Figs 5-7, respectively. For such initial conditions, no defect in the polarity patterns appears.
In Fig. 8, numerical results were obtained from the reaction-diffusion model given by equations (43) and (4) with an inclusion of additive noise θ = p t p ( , ) ( ,0) i i i (1) and external signal G i = (cos(ψ − θ i ), 0) with ψ = π. The phase was determined in the same manner as Fig. 4. The theoretical probability distribution was obtained using