Standard Anatomical and Visual Space for the Mouse Retina: Computational Reconstruction and Transformation of Flattened Retinae with the Retistruct Package

The concept of topographic mapping is central to the understanding of the visual system at many levels, from the developmental to the computational. It is important to be able to relate different coordinate systems, e.g. maps of the visual field and maps of the retina. Retinal maps are frequently based on flat-mount preparations. These use dissection and relaxing cuts to render the quasi-spherical retina into a 2D preparation. The variable nature of relaxing cuts and associated tears limits quantitative cross-animal comparisons. We present an algorithm, “Retistruct,” that reconstructs retinal flat-mounts by mapping them into a standard, spherical retinal space. This is achieved by: stitching the marked-up cuts of the flat-mount outline; dividing the stitched outline into a mesh whose vertices then are mapped onto a curtailed sphere; and finally moving the vertices so as to minimise a physically-inspired deformation energy function. Our validation studies indicate that the algorithm can estimate the position of a point on the intact adult retina to within 8° of arc (3.6% of nasotemporal axis). The coordinates in reconstructed retinae can be transformed to visuotopic coordinates. Retistruct is used to investigate the organisation of the adult mouse visual system. We orient the retina relative to the nictitating membrane and compare this to eye muscle insertions. To align the retinotopic and visuotopic coordinate systems in the mouse, we utilised the geometry of binocular vision. In standard retinal space, the composite decussation line for the uncrossed retinal projection is located 64° away from the retinal pole. Projecting anatomically defined uncrossed retinal projections into visual space gives binocular congruence if the optical axis of the mouse eye is oriented at 64° azimuth and 22° elevation, in concordance with previous results. Moreover, using these coordinates, the dorsoventral boundary for S-opsin expressing cones closely matches the horizontal meridian.


Supplemental Discussion
Assessment of the method As far as we are aware, a computational algorithm for reconstructing retinal flat-mounts into standard space and then transforming these retinal coordinates into visuotopic coordinates has never been described before. Our validation studies suggest that the method is able to estimate the original location of a point on a flattened retina to within 8 • of arc, which is equivalent to 3.6% of the nasotemporal axis. The ability to reconstruct retinae allows for the analysis of populations of labels that would otherwise be split across the cuts or tears of the flat-mount preparation. It also allows retinae from different eyes to be mapped onto a standard retina, enabling cross-eye comparisons to be made. The ability to map images into standard retinal space is a useful method of visualisation, as are the various map projections from standard retinal space. As part of the software package that implements the algorithm, we have provided routines to analyse point data on the sphere.
It is important to understand the limitations of the method and how they might be mitigated. Firstly, the algorithm can only be as good as the experimental error of dissection and fixing fed into it. There is a great deal of variability in dissected retinae: in some it is clear where the vertices of cuts and tears lie and in others is much less so; thus error may be added when the cuts and tears are marked up. The deformation measure gives an indication of the quality of a reconstruction and was higher (worse) for younger retinae ( Figure 3B), which are more delicate and harder to dissect cleanly. A high deformation measure indicates that a reconstruction should be treated with caution. In retinae where the locations of the tears and cuts are less clear, trying an alternative marking up can reduce the deformation.
Secondly, the physical model used by the algorithm contains simplifying assumptions. The retinal tissue is modelled as a mesh of masses connected by springs. This basic model is used in modelling deformable surfaces in computer animation and clothes design [1][2][3] and has been used in models of soft tissue [4]. This model should ensure that neighbourhood relations are preserved -nodes close in the 2D structure should still be close in the 3D structure -and incorporates the notion of an elastic material. However, it does not incorporate an explicit representation of tensile and shear forces.
A more physically-principled method of modelling deformation of objects is offered by the finite element method (FEM), in which the stress and strain of each triangular element in the mesh is derived as a function of the deformation of the points of the element [5]. The stress and strain matrices represent shearing within the material, and materials with varying degrees of compressibility (as quantified by Poisson's ratio) can be modelled. The FEM is used widely to describe properties of soft tissue in simulations of surgery [6]. The FEM could be applied fairly straightforwardly to the forward problem of how a retina with known cuts deforms during flattening. However, it is not possible to apply the FEM directly to the problem of mapping the flattened retina onto the intact retina. This would require knowledge of the stresses and strains in the flattened retina, whereas the only information available is the outline of the retina and the location of the optic disc.
Finally, and less crucially, the meshing algorithm is sensitive to neighbouring points on the outline being very close together. This leads to a large number of very small triangles being constructed, which can lead to a tangle of many flipped triangles in the energy minimisation phase and, ultimately, small areas of high strain in the reconstructed outline. Suitable prepossessing of outlines to remove very short edges can eliminate this problem.

Possible extensions
In terms of the reconstruction algorithm, it might be possible to automate the mark-up procedure. This would be straightforward for cleanly dissected retinae, but would be difficult for retinae with irregular cuts and tears. Furthermore marking up takes at most a few minutes using the interface supplied with the software. The other piece of information required by the algorithm is the rim angle, which is currently determined by the age of the tissue (Table 1). It would be possible to infer this automatically ( Figure 3D) but we decided that the extra complexity did not warrant doing this.

Animals
All animals used in this study are on a C57BL/6J background (Charles River, UK). All procedures were performed in accordance with the European Communities Council Directive of 24 November 1986 (86/609/EEC) under the Animals (Scientific Procedures) Act 1986. For surgery, adult animals were given an injection of Vetergesic (0.1 mg/kg) and anaesthetised either with Isoflurane (induction at 4% and maintenance at 2%) or ketamine/xylazine (3 ml/kg of stock solution: 33 mg/ml, Ketamine and 3 mg/ml, Xylazine). Animals were then secured with ear-bars in a custom-made frame and maintained at 37 • C throughout surgery. For injections into the dLGN, a skull flap was drilled and either unilateral or bilateral cortical aspiration was performed to reveal the thalamus. Subsequently, either fluorescein-conjugated dextran (Fluoro-Emerald) or tetramethylrhodamine-conjugated dextran (Fluoro-Ruby) (5% w/v; 10,000MW; D1820, D1817; Invitrogen) was pressure-injected through glass pipettes with a 120 µm internal diameter and tip diameters 30-100 µm. Multiple injections were made into the thalamus; total volumes into one thalamus ranged between 50 nl and 250 nl. For injections into the superior colliculus, a skull flap was similarly drilled and unilateral cortical aspiration was performed to reveal the entire superior colliculus. Subsequently, 5-10 nl of Red and Green latex microspheres (RetroBeads; Lumafluor) were pressureinjected into the upper 200-300 µm of the superior colliculus using similar pipettes as above, but with 30-50 µm tip diameters. Following surgery, the cavity from aspiration was filled with absorbable gelatine sponge, the skull-flap replaced and the skin sutured. Animals were allowed to recover and left for 24-48 hours after which they were terminally anaesthetised by intraperitoneal injection of pentobarbital sodium (300 mg/kg) and perfused transcardially with phosphate buffered saline at 4 • C followed by paraformaldehyde (4% w/v) at 4 • C.

Retinal dissections
With the eye still in the head, a microknife cut was made in nasal cornea at the level of the middle of the nictitating membrane. This allowed a subsequent orienting cut to be made in nasal retina. The eye was removed and dissected using a Sylgard-filled Petri dish and insect-pins. The front of the eye was removed just anterior to the corneo-scleral junction and the lens removed. Fine scissors were used to make a long cut through the sclera, choroid and retina from the nasal orienting cut towards the optic disc. Shorter relaxing cuts were made at approximately 90 o or 120 o intervals, depending on the age of the animal. The retina was freed from the choroid and cut free at the optic disc. It was then transferred to a Poly-L-Lysine coated microscope slide, mounted ganglion cell layer up and cover-slipped using Fluoromount (Sigma). Retinae were imaged using fluorescence microscopy.
To investigate the orientation of the eye, the insertions of the superior, inferior and lateral rectus muscles were marked onto the retina by passing a needle coated with Fluoro-Emerald or Fluoro-Ruby (5% w/v) through the sclera into the retina at the centre of the muscle-insertion prior to dissecting the retina. The staining for shortwave opsin was done using a polyclonal antibody against S-opsin (1:200; AB5407; Millipore) with a Fluorescein-conjugated secondary antibody (1:1000; AP123F; Millipore).

Reconstruction algorithm
The reconstruction algorithm described here has been implemented in R [7] and can be downloaded as an R package from http://retistruct.r-forge.r-project.org/. It has been tested on GNU/Linux, but should also work in MacOS and Windows. The user manual, available from the same site, contains details of the two main data formats Retistruct can process. These are either in the form of coordinates of data points and retinal outline from an in-house camera-lucida setup, or in the form of bitmap images with an outline marked up in ImageJ ROI format [8]. In this paper, retinae were imaged using an Axiophot2 microscope (Zeiss) and processed to find cell locations using ImageJ. Retinal ganglion cell locations in retinae used for Figures 1, 2 and 5, however, were plotted using an in-house camera-lucida setup. For the latter, retinae were either sampled fully or partially (one in nine 150 µm square boxes).
The steps of the reconstruction algorithm are illustrated in Figure 1. The raw data ( Figure 1A) consists of the sequence of points making up the outline, sets of data points and sequences of connected points defining landmarks, such as the optic disc. The reconstruction process then proceeds as follows: 1. The location of one of the nasal position and all the incisions and tears in the outline are marked up by an expert ( Figure 1B) using the custom-built graphical user interface. Each tear is defined by a central point, referred to as the apex, and two outer points, referred to as vertices. Tears within tears or cuts can be marked up.
2. The retinal outline is triangulated using the conforming Delaunay triangulation algorithm provided by the Triangle package [9] such that there are at least 500 triangles in the outline (grey lines in Figure 1C).
3. The tears and cuts are stitched automatically (cyan lines in Figure 1D). To do this, firstly the length of each side of every tear is computed. The fractional distance of each point in each tear is then defined as the distance along the tear of that point from the apex divided by the total length of that side of the tear. For each point on one side of a tear a corresponding point is inserted at the same fractional distance along the opposing side. Tears within tears are dealt with by excluding the child tear when computing the fractional distance. At the end of the procedure there is a set of correspondences between two or, in the case of the vertices of a child tear, three points.

4.
There is then an extra round of triangulation to incorporate the points that have been inserted into cuts and tears.
5. The points within each set of correspondences are merged and allocated positions on the 2D surface.
6. The triangulation points are then projected onto a sphere curtailed at latitude of ϕ 0 ( Figure 1D). Note that whilst we present the coordinates of the standard retina in terms of colatitude and longitude, all internal calculations use latitude and longitude. Colatitude is converted to latitude by subtracting 90 • . The latitude is estimated on the basis of measurements from intact retinae of animals of the same age as the retina under reconstruction. The radius R of the sphere is determined by the area of the flattened retina and ϕ 0 . Points on the rim of flattened retina are fixed to the rim of the curtailed sphere.
7. The optimal projection onto the curtailed sphere ( Figure 1E) is inferred by shifting the locations of the vertices on the spherical surface so as to minimise the energy measure, E, which comprises the sum of normalised squared differences between lengths of corresponding connections in flattened and spherical retina (e 2 L , see Equation 1) and a penalty term that prevents the triangles from flipping over: where L i and l i are lengths of corresponding edges i ∈ C in the flattened and spherical retina respectively, L is the mean length of an edge, |C| is the number of edges, f is a penalty function to be defined below, and A j and a j are signed areas of corresponding triangles j ∈ T in the flattened and spherical retina and A is the mean of A j . The parameter β controls the relative contribution of the area penalty, and ν controls how much larger triangles are penalised than smaller ones; they are set as described below. The signed area a j is positive for triangles in correct orientation, but negative for flipped triangles. The penalty function f is a piecewise, smooth function that increases with negative arguments: The parameter x 0 is set at 0.5. There is thus no penalty unless triangles have been squashed to less than 50% of their size in the flattened retina. The length of each edge l i is computed from the formula for the central angle between its vertices: where ϕ 1 and ϕ 2 are the latitudes of the vertices and λ 1 and λ 2 are the longitudes. The derivatives of E with respect to ϕ i and λ i are computed and used to minimise E. Optimisation proceeds by: (i) turning off the area penalty β = 0 and using the FIRE structural relaxation algorithm [10]; (ii) starting from this configuration, setting β = 8 and ν = 1 and again using FIRE to minimise the energy, which removes a number of the flipped triangles, dealing preferentially with the biggest; and (iii) minimising using the BFGS quasi-Newton method (as implemented in the R optim function) with β = 8 and ν = 0.5 to refine the optimisation. This procedure was arrived at after trials on a corpus of over 200 marked-up retinae. The procedure of reconstructing a retina takes around 1 minute (range 20 seconds to 11 minutes) on a Dell Optiplex 990 with an i7-2600 CPU running Scientific Linux 5.

Projections used to display of reconstructed data
The data on the sphere can be projected onto a polar or azimuthal equidistant plot in the plane with coordinates (ρ, λ) where the radius ρ = π/2 + ϕ. Alternatively the data can be plotted in an azimuthal equal area (Lambert) projection [11], in which area is preserved in the sense that equal areas on the sphere project to equal areas on the plane of projection, achieved by setting ρ = 2(1 + sin ϕ). Finally the data can be plotted in an azimuthal conformal (Wulff) projection in which angles are preserved; this is achieved by setting ρ = tan (π/4 + ϕ/2). For data transformed into visuotopic space (see later), we used the sinusoidal orthographic projections. The sinusoidal projection, which projects the entire globe onto the plane and preserves area, is given by x = (λ − λ 0 ) cos ϕ, y = ϕ, where λ 0 is the longitude at the centre of the projection [12]. The orthographic projection, which gives a perspective view of one side of the globe, is given by x = cos ϕ sin(λ − λ 0 ), y = cos ϕ 0 sin ϕ − sin ϕ 0 cos ϕ cos(λ − λ 0 ), where the projection is centred on (ϕ 0, λ 0 ) [13].
By convention, the polar plots are viewed as though the animal is facing towards the observer. This means that when plotting a retina from a right eye, the nasal pole on the right and N,

Analysis of reconstructed data
Mean location of data points The Karcher mean of a set of points on the sphere [14,15] is defined as the point (ϕ, λ) that has a minimal sum of squared distances to the set of points (ϕ i , λ i ) : where the function l is defined in Equation 5.
Kernel regression estimates In kernel regression, we have values y i at each data point (ϕ i , λ i ). We adapt the Nadaraya-Watson estimator [18,19] to the Fisher density and the spherical coordinates: The concentration parameter κ is set by minimising the summed squared error: whereŷ i (ϕ i , λ i ; κ) is the kernel regression estimate of y i calculated from all but the ith data point.

Contour analysis
The kernel estimates can be represented using a contour plot. This is produced by the standard method [11] of laying a 100 by 100 Cartesian grid over the azimuthal area-preserving (Lambert) projection, finding the spherical coordinates corresponding of each point on the grid, and then computing the density estimate at each point. A standard contouring algorithm is then used to plot the contours. There are a number of methods for determining contour heights [11]. We chose to define them in terms of the amount of probability mass that they exclude; e.g. the 5% contour excludes 5% of the probability mass. The area contained within each contour can be determined simply by finding the area within the contour in the area-preserving (Lambert) projection [11].

Mapping of retina onto visual space
In previous work of the mapping of visual space on the superior colliculus [20], visual space has been described in terms of azimuth θ and elevation α with respect to the long axis of the mouse. Visual space is mapped onto the retina via the optical system of the eye. This mapping depends on the angle that a ray makes with the optic axis, which we define as being oriented at azimuth θ 0 and elevation α 0 . In order to map the eye onto visual space, we consider a reference frame in which the z-axis is vertical, the x-axis is pointing forwards along the long axis of the mouse and the y-axis is pointing rightwards, as viewed by an observer facing the animal. Notionally, the optic axis of the eye starts off vertical, and a sequence of rotations is undertaken in order to move the eye to its correct orientation. First the projection of rays from the retina to visual space is determined using the approximate optics defined in the main text. Second, the coordinates of points are converted to Cartesian coordinates. Third, the coordinates are rotated about the x-axis by −90 + α 0 degrees. This ensures that the optic axis is oriented at the elevation α above the horizontal, but the direction of its azimuth will be −90 • . The points on the eye are then rotated about the original z-axis by 90 + θ 0 degrees giving the optic axis an azimuth of θ 0 , where positive values are to the right, as viewed from a point on the positive x-axis. Finally the Cartesian coordinates are transformed back to spherical coordinates in the original reference frame.