| Methodology |
-
Detection of grain boundaries (in flat-mounted, row-traced retinae):
Because the flat-mounted retinae are of finite size and because grain boundaries can both initiate and terminate within the retinae (forming so-called grain boundary scars), identifying grain boundaries is not as simple as one might intuit based on experience with infinite planar crystals. Even if all Y-junctions 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 Y-Junction as being in a grain boundary or not being in a grain boundary.
In the flat-mounted retinae, we have the positions of Y-junctions, 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 Y-junction belongs to a grain boundary, we check whether the row orientation of the surrounding domain rotates appreciably in the vicinity of the Y-junction.
This method requires two arbitrary parameters:
1. the length-scale 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 Y-junction, we average the domain orientation in five distinct boxes (which we will use to compute how much the row orientation changes near the Y-Junction):
1. The first box is centered on the Y-junction 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_0-3 ∆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 θ=θ_0-3 ( 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 Y-Junction 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 flat-mounted retina, because of retinal deformations due to cutting and flattening, grain boundaries are frequently not well-aligned 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 Y-Junction about both directions.
To determine whether a Y-junction 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 Y-Junction is in a grain boundary; if not, the Y-junction 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 Y-Junctions 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 Y-junctions in grain boundaries. In addition, we fix the values of these parameters when comparing flat-mounted retinae to simulations.
We note a couple possible pitfall of this method as applied to flat-mounted 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, Y-junctions 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 flat-mounted (experimental) retinae with one exception:
Instead of having a semi-continuous 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 phase-field 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 particle-density-modulation 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 first-order implicit-explicit methods, treating the non-linear term in ψ explicitly. We implement all derivatives by finite differences. We use no-flux boundary conditions at each non-periodic 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 phase-field 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 flat-mounted retinae. By changing the opening angle of the cone, we can tune the number of Y-Junctions required to maintain constant cell-cell 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 Y-Junctions necessary to maintain constant cell-cell spacing is comparable to the number of Y-Junctions 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 white-noise 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
|
|---|
