Coupling Mechanical Deformations and Planar Cell Polarity to Create Regular Patterns in the Zebrafish Retina

The orderly packing and precise arrangement of epithelial cells is essential to the functioning of many tissues, and refinement of this packing during development is a central theme in animal morphogenesis. The mechanisms that determine epithelial cell shape and position, however, remain incompletely understood. Here, we investigate these mechanisms in a striking example of planar order in a vertebrate epithelium: The periodic, almost crystalline distribution of cone photoreceptors in the adult teleost fish retina. Based on observations of the emergence of photoreceptor packing near the retinal margin, we propose a mathematical model in which ordered columns of cells form as a result of coupling between planar cell polarity (PCP) and anisotropic tissue-scale mechanical stresses. This model recapitulates many observed features of cone photoreceptor organization during retinal growth and regeneration. Consistent with the model's predictions, we report a planar-polarized distribution of Crumbs2a protein in cone photoreceptors in both unperturbed and regenerated tissue. We further show that the pattern perturbations predicted by the model to occur if the imposed stresses become isotropic closely resemble defects in the cone pattern in zebrafish lrp2 mutants, in which intraocular pressure is increased, resulting in altered mechanical stress and ocular enlargement. Evidence of interactions linking PCP, cell shape, and mechanical stresses has recently emerged in a number of systems, several of which show signs of columnar cell packing akin to that described here. Our results may hence have broader relevance for the organization of cells in epithelia. Whereas earlier models have allowed only for unidirectional influences between PCP and cell mechanics, the simple, phenomenological framework that we introduce here can encompass a broad range of bidirectional feedback interactions among planar polarity, shape, and stresses; our model thus represents a conceptual framework that can address many questions of importance to morphogenesis.

In this section we derive the explicit dynamical equations for the cell shapes and polarity protein concentrations. The dynamical equation for PCP can be rewritten from Eq. 6 in the main text as where β refers to the cell sharing the edge i with cell α, and c α,i B satisfies a symmetric equation where A and B are exchanged. λ α A is a Lagrange multiplier constraining the total number of proteins l i c α,i A within one cell. The main influence of cell shape, besides terms arising from protein number conservation, arises from the couplings J 1 and J 2 between neighboring edges: J 1 and J 2 have indeed dimensions of length and can be thought of as the typical length over which proteins interact on neighboring edges. In this description longer edges are therefore less sensitive to the concentration of their neighbors.
In addition, we assume the following form for the relaxation of the positions of the vertices and curvatures of the edges where : • e k i is the unit vector tangent to edge i at vertex k and pointing outwards from the vertex. • as described in the main text, h i is the distance between the midpoints of the curved edge i and of the straight line joining its two endpoints. α i and β i are the two cells separated by the edge i, with the assumption that h i > 0 if the curved edge bends towards cell α i .
• R i is the radius of curvature of edge i, related to h i by R i = It is clear from Eqs. ?? and ?? that, in the steady state, force balance at the vertices is satisfied and the edges are circular arcs satisfying Laplace's law, ∆P = γ R .

B. Curvature corrections
Our model approximates regions of the retinal epithelium as flat, planar domains. In reality, however, the retina is a roughly hemispherical sheet. In this section, we show that our approximation is justified in the framework of the vertex model provided that the retina's radius of curvature R is sufficiently large, and we verify that this is the case for our system.
The radius of curvature R can be compared to two other length scales in the model: The typical size l of a single cell and the linear size w of the entire simulation domain. In the case of the zebrafish retina, l ∼ 5µm, and the radius of the eye at 1mpf lies between 500µm and 1000µm, so that l/R 0.01. Our simulations are run on domains typically extending over about 75µm, which implies w/R 0.15. The major changes in cell packing at the retinal margin occur over an even smaller band, about 3 cells wide; if we take this distance to be the relevant measure of domain size, then w/R 0.03. Below, we show that the leading corrections to the vertex model associated with curvature are of two sorts: Corrections associated with the projection of force vectors onto the plane tangent to the sphere at a given vertex and with the differences in length associated with following the curved sphere surface are of order l/R or smaller. In addition, although the surface of the sphere is locally planar, any mapping of a finite-sized domain of the sphere's surface to a plane necessarily introduces deformations whose magnitude increases with domain size. We find that corrections to the forces on vertices due to these deformations are of order (w/R) 2 . Thus, neither correction exceeds about 2%.
To illustrate our point in more detail, we consider a vertex model representing a cell packing analogous to the model we have described, except that points are constrained to lie on a sphere of radius R. For the sake of simplicity we further assume that edges connecting two vertices are straight lines in 3D; models which include edge curvature introduce additional terms of order at most l/R. To compare this model to a planar model, we need an explicit mapping between the two. Here, we choose a mapping in which distortion is minimal near the equator and becomes progressively larger at one moves towards the poles. It is thus particularly well-suited to simulations that focus on a long strip near the equator-or, in our case, near the rim of a hemispherical retina. The coordinates of a point lying on a sphere are the polar angle θ, taken as 0 on the equator and the azimuthal angle φ (Fig. ??). These coordinates are mapped onto the Cartesian coordinates (x, y) in the plane according to the coordinate transformation (x = Rθ, y = Rφ).
The length of the straight line joining two points A (θ a , φ a ) and B (θ b , φ b ) on a sphere is In the limit l ab R, we can write θ b = θ a + ∆θ and φ b = φ a + ∆φ and expand for small ∆θ and ∆φ. We find, keeping terms up to third order, that l ab = R ∆θ 2 + ∆φ 2 cos 2 θ a − ∆φ 2 ∆θ cos θ a sin θ a , or, in terms of the coordinates for the planar representation, This expression differs from the length of the straight line joining A and B in the planar representation, l planar ab = ((x b − x a ) 2 + (y b − y a ) 2 , by a factor proportional to sin 2 xa R at lowest order. The force exerted on vertex A by the tension γ on the edge joining A and B then has components given in terms of the (x, y) planar coordinates, The first terms in brackets in the expressions for each force component are exactly the x and y components of a unit vector pointing from A to B, and thus correspond to the expression used in a planar vertex model (Eq. 4 in the main text). The additional terms represent curvature corrections that are neglected when computing forces in the planar vertex models. These corrections are of two types: • The second term in Eq. ?? introduces a correction proportional to sin 2 xa R (x a /R) 2 . This correction is due to the deformation generated by the mapping from spherical to planar coordinates, which, for our choice of mapping, induces an increasing stretching in the y direction as one goes away from the equator of the sphere. If our simulations occur in a band of width w near the retinal margin, then x a is at most w, so these corrections go like (w/R) 2 , which is about 0.02 for the values of w and R given above.
• Corrections smaller than the main term by a factor l/R-or, in fact, of order lw/R 2 for this particular mapping-are due to effects associated with the curvature at the scale of a single cell. As noted above, in the case of the zebrafish retina epithelium, l R ∼ 0.01 so that these corrections are less than 1%. One can verify that the forces associated with 2D pressures and area changes introduce similar corrections.
Taken together, this analysis shows that the planar model gives correct forces near the margin of the eye with errors of about 2% or smaller. Changes in forces of this size will introduce an equivalent error in the position of the vertices; we do not expect such a small change to modify our basic conclusions.

C. Numerical methods
To simulate the polarity protein model, we have implemented an implicit-explicit scheme (IMEX scheme, [? ]). This procedure was required due to the stiffness introduced by the logarithmic, entropic terms in F p . To represent each concentration we introduce the auxiliary variable X α,i = − log(l i c α,i ), where we have suppressed the subscript A or B. The entropic term ensures that on every edge the total quantity of protein l × c is greater than 0, so X α,i remains real. The use of the auxiliary variable allows us to keep track of very small densities that would otherwise be below the machine precision.
The dynamical equation for the concentration of polarity proteins (Eq. 6 in the main text) can be rewritten in the compact form .4Ne denotes the vector of densities of the PCP proteins A and B on the two sides of the N e edges in the packing, and U and V are two matrices of size 4N e × 4N e whose coefficients can be determined from Eq. 6 in the main text. Only concentrations within the same cell are coupled by logarithmic terms, so, provided that the densities from the same cell are grouped together in the vector (c i ), V is banded, but U is typically dense. (V would in fact be diagonal, except that the Lagrange multiplier, Eq. ??, couples the logarithm of concentrations on different edges.) The logarithmic term is the source of stiffness in the problem, and we therefore treat it implicitly in the simulations, while the first term is treated explicitly. This treatment allows us to avoid inverting the dense matrix U . Rewriting Eq. ?? for the variable X, we obtain which we simulate using the following IMEX scheme: where n indexes successive time points spaced by ∆t, and δ ij is the Kronecker delta. The matrix on the left hand side is block diagonal, with one block corresponding to one cell. The block for cell α has the form Using the Sherman-Morrison inversion formula [? ], it is straightforward to obtain the following expression for the inverse matrix which, together with Eq. ??, yields an analytical expression for X n+1 i as a function of X n i . In addition, we renormalize the total number of proteins within each cell every 100 iterations to correct for additive numerical errors.
The dynamical equations for the vertex positions and edge curvatures are integrated with an explicit Euler scheme. The whole system is simulated with an adaptative stepsize method, with an absolute tolerance abs = 0.01 on the cell shape variables and absolute and relative tolerances abs = rel = 0.001 on the polarity protein variables.

Nematic order parameter from the distribution of polarity proteins
We detail here how we obtain a 2D nematic order parameter for the distribution of polarity proteins in a simulated cell packing. This order parameter determines the length and orientation of the bars indicating cell polarity in Fig. 5 of the main text and similar supporting figures and videos. A polar vector associated to the distribution of polarity proteins within cell α P α = (P α x , P α y ) can be computed using the following expressions: where θ refers to the angle relative to the horizontal axis, with the origin taken at the cell centroid. The global order parameter within the packing can then by obtained by averaging within each cell: where the averaging is on cells α in the packing. For simulations with the parameter choice 1 = 2 however, the global polarity vanishes since head to head, tail to tail and head to tail configurations of polarity between neighbouring cells are equally favored. In that case, the global emergent order is characterized by a nematic order parameter. Introducing the notation P α = P α (cos(θ α ), sin(θ α )), the nematic order parameter is a 2D tensor with 2 independent parameters where in the last matrix, S 2 = 2 Tr(Q 2 ) = 2 Q 2 xx + Q 2 xy is the order parameter satisfying 0 < S 2 < 1, and (u x , u y ) are the coordinates of the unit vector giving the global nematic orientation.

Q4 order parameter
We detail here how we extract a fourfold symmetric order parameter Q 4 from the crosses that are fitted within each cell according to the procedure described in the text. As for the nematic order parameter, the four-fold order parameter can be viewed as a symmetric, traceless second rank tensor. (Strictly speaking, one might expect fourfold rotational order to correspond to a fourth rank tensor. In two spatial dimensions, however, all real representations of the rotation group have the same dimension, and we may choose to express the fourfold rotational order in terms of a second rank tensor, thus taking advantage of what we already know about nematic order parameters.) The matrix is defined by the following expression, where θ α c now refers to the angle of one arm of the cross fitted to cell α where Q 4 = Tr(Q 2 4 ) = (Q 4 xx ) 2 + (Q 4 xy ) 2 is the magnitude of the fourfold order parameter and (v 4 x , v 4 y ) are the coordinates of a unit vector giving the global fourfold orientation. The average is taken over different cells α. Finally, from those expressions the average angle θ 4 = arctan v 4 y v 4 x can be obtained where 2 can be taken as 1 or −1 since the definition of θ 4 is invariant by rotation of π 2 .

Positional order parameter
In order to obtain an order parameter for the positional order of the cells, we compute the Fourier transform of the density of cell centroids n(x) = cells α δ(x − x α ), with x α the centroid of cell α, and normalize by the total number of cells N : The corresponding structure factor S is given by For a crystalline solid, the function S exhibits Bragg peaks located on the reciprocal lattice of the Bravais lattice. As shown in Fig. ??, a lattice consisting of points periodically distributed along one direction and disordered along the perpendicular direction has a structure factor exhibiting peaks in one direction and straight parallel lines in the other direction. Biological images of adult retina exhibit structure factors in between those two limits, with bright peaks in one direction corresponding to the column ordering and less pronounced peaks in the perpendicular direction, corresponding to the less pronounced correlations between the positions of individual cells in different columns (Fig. ??). As the formation of an ordered column is the clearest feature of the neural retina pattern, we have chosen the amplitude of the peak with the smallest q in the direction of the column as the order parameter: where q 0 is chosen as the lowest non-zero value of q maximizing the structure factor. In practice this value always corresponded to order along the direction of retinal growth and the formation of columns. It is also this ordering that is most clearly reproduced by the simulations (Fig. ??).

II. SUPPORTING VIDEOS
• Video S1: Simulation of the progressive growth of a cone photoreceptor mosaic similar to that observed near the retinal margin in adult, wildtype fish; the cone fate is induced and the PCP pathway is activated in successive columns of cells at regular intervals, as described in the text.
• Video S2: Simulation of the relaxation of a cell packing with polarity proteins influencing edge tensions, under a global isotropic stress.
• Video S3: Simulation of the relaxation of a cell packing with polarity proteins influencing edge tensions, under a global anisotropic stress (the packing is compressed along the horizontal direction).
• Video S4: Simulation of the progressive growth of a retinal cone mosaic, but with polar rather than nematic interactions between polarity proteins.  Figure 6A). Panel J is at the level of the OLM, panel K is the subapical region (SAR) at 2.5 µm from the OLM, and panel L is the SAR at 5 µm from the OLM. Identity of cone photoreceptor subtypes is indicated by asterisks: red, green, blue, and UV (magenta), respectively. In the SAR, Crb2a has a planar polarized distribution it is expressed at higher levels at the interfaces of red, green, and blue cones within a column. The retinal surface is curved, so the left and right sides of each panel are more apical than the center (and thus show a somewhat more polarized Crb2a distribution in J and K). Only interactions involving the concentration c α,i A are represented. B) Redistribution of polarity proteins in topological transitions. Left: shrinking of an edge to a 4-fold vertex; the polarity proteins on the vanishing edge are redistributed onto neighboring edges. Right: formation of a new edge from a 4-fold vertex; a minimal amount of polarity proteins is redistributed onto the newly formed edge. C) Schematic illustrating how cells adopt a rectangular shape in the simulations. Left: with uniform tension, cells typically adopt hexagonal shapes where edges meet at vertices with 120 • angles. Right: when the polarity protein dynamics is turned on, polarity proteins accumulate on opposite edges of cells, which results in an decrease of edge tension on 2 opposite sides and an increase in tension on the other 4 sides. To satisfy force balance, the cells must then deform their shapes towards a nearly rectangular shape. When minimum and maximum tensions in a cell satisfy γmin << γmax, the cell shape is close to a rectangle and, in a perfectly regular packing, the length of the edges with polarity proteins is approximately 3l 2 , where l is the length of a side in the original hexagonal packing. Assuming each polarity protein is concentrated on one of two opposite edges in the cell, the concentration reached is about cA = 2 3l , cB = 0 or cB = 2 3l , cA = 0. Counting the polarity proteins on both sides of the edges, this results in the estimates of γmin and γmax indicated below the schematic. 0/ ab = −0.4, J1/l ab = 0.24, J2/l ab = 0.12, Tel/ ab = 0.03, other parameters as described in the text). The simulation was started with one line of cells with a uniform polarity. Following this first line, each new row of cells is specified with a random distribution of polarity protein, and no global polarity cue is given. A) Example of simulation results: the cell packing does not show rectangular order. B) Snapshots of simulations illustrating the mechanism leading to patterning defects: polarity can point up or down, leading to domain formation. Because the resulting configuration is not favorable, the polarity reorients when the next column is specified, leading to polarity misalignment and eventually to a defect in cell packing.

III. SUPPORTING FIGURES
FIG. S8: Mapping from a hemispherical to a planar surface. For a hemisphere of radius R, Cartesian coordinates x, y in the plane are related to the polar and azimuthal angles θ and φ on the hemisphere as x = Rθ and y = Rφ.