Methodology 

Detection of grain boundaries (in flatmounted, rowtraced retinae):
Because the flatmounted retinae are of finite size and because grain boundaries can both initiate and terminate within the retinae (forming socalled grain boundary scars), identifying grain boundaries is not as simple as one might intuit based on experience with infinite planar crystals. Even if all Yjunctions form lines that run from the center to the periphery of the retina, if one increases the number of identical grain boundaries, the change in domain orientation associated with each grain boundary drops. Because the degree of “grain boundary”ness can vary continuously, we would like to find a way to classify a YJunction as being in a grain boundary or not being in a grain boundary.
In the flatmounted retinae, we have the positions of Yjunctions, which generate row insertions, in the traced regions. In addition, individual rows are traced (meaning that each row has a distinct curve running along it). To define whether a given Yjunction belongs to a grain boundary, we check whether the row orientation of the surrounding domain rotates appreciably in the vicinity of the Yjunction.
This method requires two arbitrary parameters:
1. the lengthscale over which we average the domain rotation (which we will call ∆r)
2. how large the domain rotation must be for a defect to be “in a grain boundary”
We describe positions on the flattened retinae via polar coordinates (where the approximate center of the retina lies at the origin). For each Yjunction, we average the domain orientation in five distinct boxes (which we will use to compute how much the row orientation changes near the YJunction):
1. The first box is centered on the Yjunction of interest (at radial coordinate r=r_0 and angular coordinate θ=θ_0). It runs from r=r_0∆r/2 to r=r_0+∆r/2 and from θ=θ_0(1/2 ∆r)/r_0 to θ=θ_0+(1/2 ∆r)/r_0 .
2. The second box runs from r=r_03 ∆r/2 to r=r_0∆r/2 and from θ=θ_0(1/2 ∆r)/r_0 to θ=θ_0+(1/2 ∆r)/r_0 .
3. The third box runs from r=r_0+∆r/2 to r=r_0+3 ∆r/2 and from θ=θ_0(1/2 ∆r)/r_0 to θ=θ_0+(1/2 ∆r)/r_0 .
4. The fourth box runs from r=r_0∆r/2 to r=r_0+∆r/2 and from θ=θ_03 ( 1/2 ∆r)/r_0 to θ=θ_0(1/2 ∆r)/r_0 .
5. The fifth box runs from r=r_0∆r/2 to r=r_0+∆r/2 and from θ=θ_0+(1/2 ∆r)/r_0 to θ=θ_0+3 ( 1/2 ∆r)/r_0 .
To calculate the orientation within each box, we first identify all row lines that run through the box. For each row line (a skeletonized line composed of individual pixels), we calculate its orientation by
1. computing the covariance matrix of the row’s pixels’ positions within the box of interest (cov; MATLAB 2016B, MathWorks, Natick, MA)
2. subsequently identifying the direction along which the row’s pixels’ positions vary the most (eig; MATLAB 2016B, MathWorks, Natick, MA)
By convention, row lines always run from the center of the retina to its periphery, so we encode each row’s orientation with a unit vector that runs along the row from center to periphery (rather than from periphery to center). Then, within the box of interest, we do a weighted average over the unit vectors (one for each row line within the box), where the weight of each row is simply the number of pixels of that row line that fall within the box. Then, we convert this weighted average of row vectors to a unit vector, which represent the row orientation within the box.
Now that we have the row orientation for each of the five boxes (unit vector (u_1 ) ̂ for box 1, …, unit vector (u_5 ) ̂ for box 5). We calculate how much the row orientation changes about the YJunction along the theta direction (∆ϕ_45=acos(dot((u_4 ) ̂,(u_5 ) ̂ ))) as well as along the r direction (∆ϕ_23=acos(dot((u_2 ) ̂,(u_3 ) ̂ ))). If the grain boundaries were purely radial, we could focus on ∆ϕ_45 and ignore ∆ϕ_23. In flatmounted retina, because of retinal deformations due to cutting and flattening, grain boundaries are frequently not wellaligned with the theta direction. Because of this, we compute: ∆ϕ_rms=√((∆ϕ_45 )^2+(∆ϕ_23 )^2 ). ∆ϕ_rms takes into account the change in row orientation in the vicinity of the YJunction about both directions.
To determine whether a Yjunction is within a grain boundary or not, we compare its ∆ϕ_rms to some cutoff (for example, ten degrees or twelve degrees or fourteen degrees). If ∆ϕ_rms is greater than the chosen cutoff, the YJunction is in a grain boundary; if not, the Yjunction is not in a grain boundary. Because by its very nature this cutoff is arbitrary, we vary this parameter (as well as ∆r) to see how the fraction of defects in grain boundaries change with these parameters. We find that the fraction of YJunctions depends very weakly on ∆r, but strongly on the cutoff in rotation of row orientation. To see if our simulations and experimental images have similar fractions of Yjunctions in grain boundaries. In addition, we fix the values of these parameters when comparing flatmounted retinae to simulations.
We note a couple possible pitfall of this method as applied to flatmounted retinae.
1. Assume that all regions of the cone mosaic have the same (anisotropic) spacings between UV cones at the time of incorporation; we expect that as the retina grows by incorporation of new cones, older regions of the retina are anisotropically deformed as the hemispheric retina dilates. Thus, using a uniform cutoff in ∆ϕ_rms (for “in grain boundary” versus “not in grain boundary”) probably overestimates the fraction of defects in grain boundaries near the center of the retina relative to the fraction of defects in grain boundaries near the periphery of the retina.
2. In regions where the row orientation is slightly warped due to deformation of the retina from the flattening process, Yjunctions may be designated as “in a grain boundary” (even if they are isolated) because of the deformation of row orientations.

Detection of grain boundaries (in simulated cone mosaic):
The basic method for detection of grain boundaries in simulated retinae is identical to that for flatmounted (experimental) retinae with one exception:
Instead of having a semicontinuous curve, composed of pixels, that traces each row, we have row “bonds” between peaks in the simulated density. These row bonds are line segments as opposed to collections of pixels. Because of this difference, to calculate the row orientation within a given box, we find row bonds that fall entirely within the box. We compute a weighted average over row bonds, represented as vectors that run (mostly) from center of simulated retina to its periphery; each weight is equal to the length of the corresponding row bond.

Numerical solutions of anisotropic phasefield crystal model on cone:
The free energy F for an anisotropic phase field ψ:
F=∫_S〖(ψ[r+(1+∇_s^2 )^2 ]ψ+ψ^4/4〗)dr ⃗
∇_s^2 = b^2 ∂^2/〖∂x〗^2 +1/b^2 ∂^2/〖∂y〗^2 ; b>1.
(stretched along the x direction)
∇_s^2 = b^2 1/r ∂/∂r (r ∂/∂r)+1/b^2 1/r^2 ∂^2/〖∂θ〗^2 ; b>1.
(stretched along the r direction in polar coordinates)
Note that r (scalar) is the degree of undercooling.
The particledensitymodulation field ψ is evolved to minimize the free energy F while conserving the mean of the field (ψ_0=∫_S〖ψdr ⃗ 〗/∫_S〖dr ⃗ 〗):
∂ψ/∂t=∇^2 ([r+(1+∇_s^2 )^2 ]ψ+ψ^3)
For solving this equation on the cone, we first map the cone to a flat plane, which does not generate distortions because the cone is isometric to the plane. We, then, use the Laplacian for polar coordinates, respecting the periodic boundary condition of the cone. We set up in the problem in terms of the variables u, v, ψ as defined in Backofen et al. For computational efficiency, we take the Fourier transform along any direction which is periodic (e.g., along the θ direction on the cone). We use firstorder implicitexplicit methods, treating the nonlinear term in ψ explicitly. We implement all derivatives by finite differences. We use noflux boundary conditions at each nonperiodic boundary.
We evolve the system with a fixed step size in time (∆t=0.075). The computational grid is such that there are approximately 25 grid points per lattice spacing along the circumferential direction in the initial row, and approximately 10 grid points per lattice spacing along the radial direction.
We systematically vary the parameters of the phase field crystal model; note that stretching the crystal does not change the phase diagram. In short, we take a 1D cut of the phase diagram, setting ψ_0=√(r)/2 as we vary the undercooling parameter r. For this cut, we also vary the strength of the noise in the initial conditions. These three parameters are the parameters on which we do not have any quantitative handle (relative to experiments); therefore, we perform this robustness analysis on these three parameters.
Geometry for cone mosaic growth:
The retina is approximately hemispheric. A hemisphere might, thus, seem like the most obvious choice of geometry in which to test cone mosaic growth. It is important to note, however, that the retina is not a hemisphere of a fixed radius during development. Its radius increases as new retinal cells are incorporated. As the hemisphere dilates, the existing cone photoreceptor layer must be deformed. (Note that we also discuss this dilation issue in the method for the detection of grain boundaries.) The exact way in which the existing pattern deforms, and how that affects subsequent cone mosaic formation, is beyond the scope of this paper. Our aim is to choose a minimal geometry which allows us to test the phasefield crystal model’s ability to form grain boundaries.
We choose a geometry in which we can easily tune the defect density. We choose a cone frustum, which is constructed by slicing off the top of a cone with a plane that is parallel to its base. By changing the level at which we slice the cone, we tune the number of UV cones in the initial column. We choose the top level such that there are approximately two hundred initial rows, which is consistent with the number of initial rows identified in flatmounted retinae. By changing the opening angle of the cone, we can tune the number of YJunctions required to maintain constant cellcell spacing per added column. We choose an opening angle such that two row insertions are required per added column, to maintain approximately constant cell spacing. The number of YJunctions necessary to maintain constant cellcell spacing is comparable to the number of YJunctions observed in the retinae.
Initial conditions for cone mosaic growth:
At the very top level of the cone frustum, we lay down one column of cones. We add a whitenoise mask to this initial column of cones. Because we do not have any quantitative handle on the noise in cell positions at the onset of patterning in the zebrafish retina, we vary the noise strength, exploring its effect on the subsequent pattern of defects.

Citations to related material 
 Nunley, H., Nagashima, M., Martin, K., Gonzalez, A. L., Suzuki, S. C., Norton, D. A., Wong, R. O. L., Raymond, P. A., & Lubensky, D. K. (2020). Defect patterns on the curved surface of fish retinae suggest a mechanism of cone mosaic formation. PLOS Computational Biology, 16(12), e1008437. https://doi.org/10.1371/journal.pcbi.1008437
 Defect patterns on the curved surface of fish retinae suggest mechanism of cone mosaic formation Hayden Nunley, Mikiko Nagashima, Kamirah Martin, Alcides Lorenzo Gonzalez, Sachihiro C. Suzuki, Declan Norton, Rachel O. L. Wong, Pamela A. Raymond, David K. Lubensky bioRxiv 806679; doi: https://doi.org/10.1101/806679
