The Fovea in Retinopathy of Prematurity

Abstract

Purpose: Because preterm birth and retinopathy of prematurity (ROP) are associated with poor visual acuity (VA) and altered foveal development, we evaluated relationships among the central retinal photoreceptors, postreceptor retinal neurons, overlying fovea, and VA in ROP.

Methods: We obtained optical coherence tomograms (OCTs) in preterm born subjects with no history of ROP (none; n = 61), ROP that resolved spontaneously without treatment (mild; n = 51), and ROP that required treatment by laser ablation of the avascular peripheral retina (severe; n = 22), as well as in term born control subjects (term; n = 111). We obtained foveal shape descriptors, measured central retinal layer thicknesses, and demarcated the anatomic parafovea using automated routines. In subsets of these subjects, we obtained OCTs eccentrically through the pupil (n = 46) to reveal the fiber layer of Henle (FLH) and obtained adaptive optics scanning light ophthalmograms (AO-SLOs) of the parafoveal cones (n = 34) and measured their spacing and distribution.

Results: Both VA and foveal depth decreased with increasing ROP severity (term, none, mild, severe). In severe subjects, foveae were broader than normal and the parafovea was significantly enlarged compared to every other group. The FLH was thinner than normal in mild (but not severe) subjects. VA was associated with foveal depth more than group. Density of parafoveal cones did not differ significantly among groups.

Conclusions: Foveal structure is associated with loss of VA in ROP. The preserved FLH in severe (relative to mild) eyes suggests treatment may help cone axon development. The significantly larger parafovea and increased outer nuclear layer (ONL) thickness in ROP hint that some developmental process affecting the photoreceptors is not arrested in ROP but rather is supranormal.

The fovea centralis is a depression near the center of the macula lutea where the retina is notably thinned. Several adaptations combine to provide the part of the retina under the fovea with the highest visual acuity (VA).¹ During development, the fovea is formed by centrifugal displacement of ganglion, Müller, bipolar, and other cells, while centripetal movement tightly packs the rod-like cone photoreceptors.^2,3 The centripetal cone packing commences prior to and completes well after the centrifugal movement of postreceptor neurons that forms the fovea.⁴ There is wide interindividual variability in both the centripetal and centrifugal processes,^5–7 with only subtle impact on ultimate VA,^8,9 in normally developing eyes. The coefficient of variance in cone density among normal eyes is greatest in the center of the fovea,^5–7 where each cone is connected to approximately three ganglion cells.¹⁰ Although the foveal area projects to only approximately two-tenths of one percent of the complete visual field (approximately 20 square degrees of >10,000 square degrees),¹¹ the foveal cones are so densely packed and the number of ganglion cells responding to those cones is so great that the fibers of the optic nerve are almost evenly split between those conveying signals from below the fovea and those conveying signals from the entire extra-foveal retina.¹² Structural abnormalities of the fovea are, therefore, of considerable interest. 

Retinopathy of prematurity (ROP) has long been associated with abnormal foveal development.^13–15 It is well-documented that the fovea in eyes with a history of ROP is shallow and that both persistent postreceptor cells and fine capillaries may be present where usually there is an avascular zone.^16,17 Reduced VA in ROP is also well-documented.^18–23 Herein, we modeled the fovea and the underlying central retinal layers in eyes with a history of preterm birth without and with ROP to obtain objective measures of foveal anatomy. We related these observations, which were based upon optical coherence tomograms (OCTs), to observations of the parafoveal cone photoreceptors displayed in adaptive optics (AO) scanning light ophthalmograms (SLOs). We also sought to determine if foveal abnormalities in ROP explain central retinal dysfunction, as represented by low VA. 

Methods

We used a Spectralis (Heidelberg Engineering, Carlsbad, CA, USA) spectral domain OCT system to capture volumetric images of the macula. In a subset of subjects, we used our previously described AO-SLO system, MAORI (Physical Sciences, Inc., Andover, MA, USA)^24,25 to capture confocal images of the parafoveal cone photoreceptors, from which we measured their density and characterized their packing geometry. Results of the OCT and SLO evaluations were compared to VA and other ocular features. 

Subjects

This study included subjects (N = 245), aged 10 to 37 years at test (mean/median = 18), without and with a history of ROP (Table). We labeled subjects “term” if they were born ≥ 37 weeks postmenstrual age (PMA) or “preterm” if they were born ≤ 32 weeks PMA. Other than preterm birth or ROP, any conditions associated with altered refractive or visual development were exclusionary; these conditions included congenital retinal detachment, glaucoma, intraocular surgery other than for the management of ROP, microphthalmia, persistent hyperplastic primary vitreous, and systemic disorders. 

Table.

 

Sample Sizes and Mean (SD) Characteristics of the Subject Population, by Group

 

All preterm subjects had serial fundus examinations in the neonatal intensive care nursery similar to those used in the multicenter ROP treatment trials.^26,27 Based on the maximum acute-phase ROP noted in these examinations, we binned each preterm subject into one of three groups: “none,” “mild” ROP, or “severe” ROP. Subjects in the none group never had ROP. Subjects in the mild group had clinical disease (zone II, stage 1) that resolved spontaneously without treatment. Subjects in the severe (zone II, stage 3) group had ROP that was treated by laser ablation of the avascular peripheral retina. No subject had zone I disease, a history of retinal detachment, or surgery other than laser treatment. Written, informed consent was obtained from all adult subjects (aged 18 years or older). Parents provided consent for minor subjects, and assent was obtained from the minors. This study was approved by the Boston Children's Hospital Institutional Review Board and adhered to the tenets of the Declaration of Helsinki. 

Optical Coherence Tomogram Analyses

Image Acquisition

We collected volumetric spectral domain OCTs (“C-scans”) that spanned 15 degrees to 30 degrees horizontally (nasal-temporal, x) and 5 degrees to 25 degrees vertically (inferior-superior, y). The automated layer segmentation tool included with the Spectralis OCT was used to determine 11 distinct retinal boundaries that delineated 10 layers (from proximal to distal): the nerve fiber layer (NFL), the ganglion cell layer (GCL), the inner plexiform layer (IPL), the inner nuclear layer (INL), the outer plexiform layer (OPL), the outer nuclear layer (ONL), the photoreceptor inner segments (IS) and outer segments (OS), the retinal pigment epithelium (RPE, including the OS tips), and Bruch's membrane (BM). Values personalized to each eye, provided by the Spectralis software, were used to scale “A-scan” and “B-scan” geometry to microns. 

As noted by Lujan et al.,²⁸ most of the fiber layer of Henle (FLH) is isoreflective with the ONL, rendering it mostly invisible in standard OCT imaging and, consequently, resulting in its routine inclusion in measurements of ONL thickness obtained in segmented OCTs. However, the oblique orientation of the fibers of Henle results in a directional reflectivity that can be detected by changing the incident angle of the sample light, relative to the retina; this adds optical contrast to the FLH. In a subset of subjects (n = 46), in addition to standard OCTs (beam through the center of the pupil), we collected scans with offset pupil entry positions as far nasal and as far temporal as possible to measure the thickness of the photoreceptor axons lying on the inner side of the photoreceptor cell bodies; thus, this was an approximation of the FLH thickness obtained from these “directional” OCTs (D-OCTs). 

Automatic Detection of the Fovea

We fit a sixth-order polynomial to the outermost boundary of each scan (i.e. distal edge of BM), which we subtracted from all layers to flatten and align neighboring scans in depth, z (Fig. 1A). To determine objectively the center of the fovea from these volumes, a “meta-surface” (MS) was calculated. Specifically, at each (x, y) coordinate, the z values of the inner surface of the OPL (i.e. the OPL-INL boundary) and the external limiting membrane (ELM) were subtracted from the z position of the inner surface of the NFL (i.e. the inner limiting membrane). We used the meta surface, rather than the inner limiting membrane, because some severe ROP subjects had almost no detectable pit. Then, this meta surface was fit by minimization of the root mean square error to the equation:  

\begin{eqnarray}{z_{{\rm{MS}}}}\left( {x,y} \right) &=& {a_1}\;{e^{\left( { - \frac{1}{2}\;{{\left( {\frac{{x\;-\;{x_0}}}{{{b_1}}}} \right)}^2} - \;\frac{1}{2}\;{{\left( {\frac{{y\;-\;{y_0}}}{{{b_2}}}} \right)}^2}} \right)}} \nonumber\\ &&- {a_2}\;{e^{\left( { - \frac{1}{2}\;{{\left( {\frac{{x\;-\;{x_0}}}{{{b_3}}}} \right)}^2} - \;\frac{1}{2}\;{{\left( {\frac{{y\;-\;{y_0}}}{{{b_4}}}} \right)}^2}} \right)}}\nonumber\\ && + c\;x + d\;y + {z_0}.\end{eqnarray}

(1)

 

Figure 1.

 

Layer segmentation and foveal fitting. (A) After segmentation of 11 retinal boundaries (colored lines), OCTs were flattened and aligned by the outer boundary of Bruch's membrane (red line). A meta-surface was then derived from the indicated retinal boundaries. (B) The meta-surface was fit to a difference of three-dimensional Gaussians on a plane (shown obliquely at left and en face at right) and the parafoveal extent and foveal center were determined. The B-scan closest to the calculated foveal center was used in subsequent analyses. (C) Each of the retinal layer boundaries within the foveal region of the selected B-scan were fitted to two-dimensional Gaussians.

Equation 1 is a difference of two three-dimensional Gaussians with respective amplitudes a and standard deviations b, centered on (x₀, y₀) and superimposed on a plane with slope in x of c, slope in y of d, and z-intercept of z₀. The B-scan passing closest to (x₀, y₀) was taken to be the foveal slice (Fig. 1B). We also used this equation to objectively demarcate the parafovea, a region of the macula that circumscribes the fovea. Herein, the parafovea includes that region known to change during normal retinal development.² To quantify the parafoveal area, we optimized the values of Equation 1, eliminated the contribution of the plane by setting c, d, and z₀ equal to 0, and calculated partial derivatives of the resulting function with respect to x and y. Then we calculated values of the partial derivative with respect to y as a function of the partial derivative with respect to x, and vice versa. Zero crossings in the resulting functions produced ellipses. We took the smaller of these two ellipses to represent the parafoveal area; we recorded the length of the nasal-temporal and superior-inferior axes of this ellipse (see Fig. 1B). 

To determine the region of foveal extent in the selected B-scan, we then fitted the MS cross-section for this scan to a difference of two two-dimensional Gaussians, described by  

\begin{equation}{z_{{\rm{MS}}}}\left( x \right) = {a_1}\;{e^{\left( { - \frac{1}{2}\;{{\left( {\frac{{x\;-\;{x_0}}}{{{b_1}}}} \right)}^2}} \right)}} - {a_2}\;{e^{\left( { - \frac{1}{2}\;{{\left( {\frac{{x\;-\;{x_0}}}{{{b_2}}}} \right)}^2}} \right)}} + {z_0}\end{equation}

(2)

and took the zero-crossings in the first derivative of this formula to objectively determine the foveal center (x₀) and nasal and temporal extents:  

\begin{equation}{\rm{foveal\;edge}} = {x_0} \pm \frac{{\sqrt 2 \;{b_{1\;}}{b_2}\;\sqrt {\ln \left( {\frac{{{a_2}\;{b_1}^2}}{{{a_{1\;}}{b_2}^2}}} \right)} }}{{\sqrt {{b_1} + {b_2}} \;\sqrt {{b_1} - {b_2}} }}\end{equation}

(3)

 

Evaluation of the fovea and Central Retinal Layers

Within the region bounded by Equation 3 (Fig. 1C), to provide convenient quantification of the retinal layers, we fit the 11 layer boundaries to respective, two-dimensional Gaussians of the form:  

\begin{equation}z\left( x \right) = a\;{e^{\left( { - \frac{1}{2}\;{{\left( {\frac{{x\; - \;{x_0}}}{b}} \right)}^2}} \right)}} + {z_0}\end{equation}

(4)

where, in all cases, depth, a (µm), standard deviation, b (µm), and vertical offset, z₀ (µm), were free to vary between boundaries. We performed the fits twice: in the first case, allowing the horizontal offset, x₀ (µm), to vary between boundaries and, in the second case, optimizing x₀ as a shared, ensemble parameter. When the retinal layers were thus parameterized, we gained ready access to several measures. For example: (1) with shared x₀, the intralaminar ∆ in z₀ for the boundaries can serve as an estimate of layer thickness at the rim, and the ∆ in z₀ + a can serve as an estimate of layer thickness at the trough. (2) With x₀ free to vary, the ∆ in x₀ between layers can serve as a measure of subfoveal laminar irregularity. (3) The values of a and b of the innermost boundary (the internal limiting membrane [ILM]) can provide respective measures of foveal depth and breadth.²⁹ We note that although herein we use b as the measure of breadth because it was more analytically convenient than the commonly used full-width at half-maximum (FWHM) description of the same, we also provide the FWHM, which can be readily calculated by \({\rm{FWHM\;}} = \;2\sqrt 2 \sqrt {\ln ( 2 )} b \approx 2.35b\). Furthermore, we note that the slopes of the foveal edges at one standard deviation from x₀ and at the half-maximum point can likewise be respectively calculated as \( {\pm} {\rm{}}\sqrt 2 \sqrt {\ln ( 2 )} a/( {2b} ){\rm{\;}} \approx {\rm{\;}} {\pm} {\rm{}}0.598a/b\) (at SD) and \( \pm \sqrt e a/b \approx \pm 0.607a/b\) (at half-maximum); thus, it is convenient to think of the slope of the fovea as approximately 0.6a/b, and this is the value we adopted herein. 

The adoption of corrections for pit asymmetry, including fitting to higher-order functions, has been proposed to extract improved parameters of the fovea such as depth, breadth, and slope.^30–32 Almost by definition, the addition of more degrees of freedom provides higher-order models, such as the sum of three Gaussians, with more accurate fits of the fovea than we obtain herein using a simple, single Gaussian model, but this additional accuracy comes at a considerable cost to convenience. The price is paid in several ways, including increased computational complexity during curve-optimization (especially sensitivity to starting parameters), absence of straightforward relationships between function parameters and physical structures (e.g. “FWHM ≈ 2.35b”) necessitating complex processes for converting fits to meaningful values, and decreased clarity in the observer's mind as to what, exactly, is being measured. To assess the impact of these tradeoffs, following the fits described above, we additionally fit the ILM using a sum of three Gaussians,  

\begin{equation}z\left( x \right) = \mathop \sum \limits_{i = 1}^3 {a_i}\;{e^{\left( {{{\left( {\frac{{x\;-\;{x_{{0_i}}}}}{{{b_i}}}} \right)}^2}} \right)}}\end{equation}

(5)

and followed the method of Breher et al.³³ to extract foveal depth and breadth. As shown in Figure 2, their method depends upon identifying three peaks in the square of the second derivative of Equation 5 and, using the x values of these points to identify foveal landmarks on the original fit (a point near the nasal edge of the rim, the center of the trough, and a point near the temporal edge of the rim, respectively), and then relating these landmarks to one another to derive foveal parameters. In brief, we calculated the breadth of the fovea by measuring the length of the line segment connecting the first and third landmarks, which (approximately) spanned the fovea. We then calculated the depth the fovea by computing the length of a second line segment, initiating at and perpendicular to the first, which connected to the second landmark at the trough of the fovea. We then explored the relationship between these values and the values obtained using our simpler model. 

Figure 2.

 

The (x and z) coordinates of the ILM (stippled black line) were fit to a sum of three Gaussians (upper smooth green line). The square of the second derivative of the fit (lower smooth green line) is a trimodal function where the x positions of the three peaks (dots), respectively, indicate a position near the nasal edge of the foveal rim, the foveal center, and a position near the temporal edge of the foveal rim. The breadth and depth of the fovea were obtained by connecting the corresponding y positions, as indicated (straight green lines). The ILM fit to a single Gaussian (Equation 4; dashed blue curve) is shown for comparison.

Adaptive Optics Ophthalmogram Analyses

Image Acquisition

We obtained images using MAORI, as previously described.³⁴ We obtained square, slightly overlapping, 1.5 degrees × 1.5 degrees, AO-SLO videos (64 frames at 1-megapixel resolution) from 4 locations surrounding each subject's preferred retinal locus (PRL) by directing them to fixate a target (a “+” symbol) offset slightly (approximately 0.25 degrees) from each corner of the imaging raster. 

To obtain a single image for analysis, we corrected for local deformations across video frames due to artifacts—such as saccades, electronics noise, “jitter” in the line trigger, and galvanometer instabilities—by applying a nonrigid registration and averaging algorithm³⁵ to each video. Then, as shown in Figure 3A, we montaged the four offset images in Photoshop CS3 (Adobe Systems Inc., San Jose, CA, USA), the last version in which interactive “photomerging” is available, to create a single, approximately 2.5 degrees × 2.5 degrees image of the fovea. We then masked regions of these images in which the cones were not clearly displayed, leaving only areas deemed to show the highest quality cone mosaic. We used a semi-automatic routine to label every cone. We then processed each dataset (image, mask, and cone centroids) using custom MATLAB (The MathWorks, Natick, MA, USA) routines that produce several parameters describing the cone mosaic. 

Figure 3.

 

AO-SLO image processing. (A) A montage of four images surrounding the fovea (blue spot) were masked for regions showing clear cones. (B) Delaunay triangulation of the cone mosaic was performed (green lines). Cones with Delaunay line segments passing through the masked region (yellow lines) were excluded from the analyses. (C) A closeup view of the mosaic and triangulation.

Analyses of the Cone Mosaic

As shown in Figures 3B and 3C, we applied Delaunay triangulation (DT) to the AO-SLO using the cone centroids as generators. DT determined each cone's neighbors and the inter-cone distances (ICDs; centroid to centroid) thereto. For every cone, we calculated the mean and standard deviation of ICD; this latter measure served as one measure of disorder.³⁶ Sidedness (S) of the DT dual graph Voronoi tessellation (VT) is, by definition, equal to the number of DT neighbors, and describes the packing pattern (e.g. hexagonal) of the cones; variability in S (i.e. the standard deviation) served as a measure of disorder.³⁶ 

Generation of “Fovea-Relative” Data

To improve our comparisons of AO-SLO data to extant anatomic data, we first transformed the DT data from units expressed in degrees to units expressed in microns, as previously described,³⁴ by calculating the ‘‘angular subtense’’ of each image⁶ based upon each subject's respective values of axial length, anterior corneal curvature, anterior chamber depth, and lens thickness.³⁷ 

In the parafovea, the cone density (D, cones·mm⁻²) versus eccentricity (ρ, mm) relationship is approximated by a power function,^38,39 of which a common formulation is  

\begin{equation}D\left( \rho \right) = c + d\;{e^{\left( { - f\;\rho } \right)}}\end{equation}

(6)

where c is density below the fovea, c − d is density in the periphery, and f describes the rate of change. Therefore, the reciprocal of density, area per cone (A, mm²·cone⁻¹), is approximated as A(ρ) = 1/D(ρ). Because ICD is essentially the diameter of each cone, and diameter is twice the radius, the relationship of ICD versus eccentricity follows  

\begin{equation}ICD\left( \rho \right) = 2\sqrt {\frac{{A\left( \rho \right)}}{\pi }} = \frac{{2\sqrt {\frac{1}{{c + d\;{e^{\left( { - f\;\rho } \right)}}}}} }}{{\sqrt \pi }}.\end{equation}

(7)

 

MAORI lacks the lateral resolution to resolve cones below the foveal center. Therefore, to estimate the position at which ρ = 0, the center of the fovea was determined by fitting ICD for every identifiable cone position (x, y) to a Lorentzian distribution of the form  

\begin{equation}ICD\left( {x,y} \right) = IC{D_0} - \frac{f}{{\left( {\frac{{{{\left( {x - {x_0}} \right)}^2}}}{{{g_1}^2}} + 1} \right)\left( {\frac{{{{\left( {y - {y_0}} \right)}^2}}}{{{g_2}^2}} + 1} \right)}}\end{equation}

(8)

where x₀ and y₀ define the foveal center (Fig. 4A). This calculated center was evaluated, “by eye,” by subjective inspection of the cone mosaic and by the central vascular pattern.⁴⁰ We did not find reason to overrule the Lorentzian model. Once the foveal center, (x₀, y₀), was so determined, the ICD data were fit to Equation 7 (Fig. 4B). 

Figure 4.

 

Mean intercone distances (ICDs) from Delaunay triangulation. (A) The ICDs (yellow dots) were fitted to a Lorentzian (gray mesh) to objectively determine the center of the retina (blue dot), where cones could not be individually resolved. (B) The eccentricity versus ICD distribution (black dots) fitted to Equation 7 (green line).

Statistical Comparisons

We performed all statistical tests in MATLAB using linear-mixed modeling (FITLME followed by ANOVA routines). Each model included, at least, a categorical test factor (i.e. fixed effect) for group and a continuous test factor for age (excepting the test of age itself). All tests other than that of age and those of AO-SLO results (which were obtained monocularly) also included a factor for eye (left versus right). If a model included any additional factors, these are described in the corresponding section of the results. We performed model estimation using the restricted maximum likelihood (REML) approach. We adopted a threshold for statistical significance of 95% certainty (α = 0.05). Where a significant effect of group was found, we performed pairwise post hoc testing (COEFTEST routine) with the threshold for statistical significance set at 99% certainty (α = 0.01); we report as “marginally significant” those cases wherein the post hoc comparison was significant in the 95% to 99% range of certainties. 

Results

Summarized clinical features of the overall subject population are given in the Table. The term subjects were significantly older than the other groups (P = 4.85·10⁻¹⁷). It was for this reason that we included the additional factor for age, throughout. Using models with group and age factors, we also found significant intergroup differences in VA (following conversion to the logarithm of the minimum angle of resolution [logMAR]; P = 5.23·10⁻⁸), spherical equivalent refraction (P = 1.29·10⁻⁹), axial length (P = 0.00226), anterior chamber depth (P = 0.00316), and lens thickness (P = 6.14·10⁻⁷). Corneal curvature did not differ significantly among the groups. Neither age nor eye (left versus right) was a significant factor in any of these analyses. 

Indeed, we note that, throughout our analyses, we never found any significant effect of eye for any tested parameter—that is, there was no systematic difference between eyes. Therefore, for simplicity, we omit further discussion of the eye factor. We furthermore describe the effect of age only when it was significant. 

Optical Coherence Tomography

Automatic detection of the fovea worked well, requiring over-riding by a human observer just over 1% of the time. In these few cases, we excluded the measurement of the parafovea (i.e. fitting Equation 1) from further analyses. Mean parafoveal extents are shown in Figure 5. We evaluated these ellipses using a model with an axis factor having nasal-temporal (horizontal) and superior-inferior (vertical) levels, and a factor for group × axis interactions. There was a significant effect of group (P = 0.0170); there were no significant differences among the term, none, and mild groups, but the severe subjects’ parafovea were slightly but significantly enlarged relative to none and mild and moderately significantly enlarged relative to term (P = 0.0294). The interaction effect was also significant (P = 0.000123), indicating that the parafoveal enlargement in severe subjects was greater along the horizontal than the vertical axis. Furthermore, we found that the parafovea was not round but significantly larger in the nasal-temporal than superior-inferior dimensions (P = 6.79·10⁻¹⁰³). 

Figure 5.

 

Mean parafoveal extent derived from fits of a meta-surface (see Fig. 1) to Equation 1.

Mean foveal layer boundaries for the standard (i.e. central pupil) OCT are plotted in Figure 6 for the 4 groups, alongside the mean thicknesses of each of the 10 layers below the fovea center and at the rim, for each group, calculated using the respective ∆s in z₀ + a and in z₀ of the layers fit to Equation 4 with shared x₀. By inspection of these layer thicknesses, we noted the well-established observation that preterm birth in general, and ROP in particular, cause inner retinal layers to be quite thick where they are normally nearly absent, with increasing ROP severity associated with increasing thickness of postreceptor laminae. Notably, layer thickness at the rim remains stable across the four groups. We evaluated layer thickness (10 levels), at the two positions (below the center and at the rim of the fovea), for significant differences, using a model with additional factors layer and position. Of course, the layers were of different thicknesses (P < 2.23·10⁻³⁰⁸) and were, overall, thinner below the fovea than at the rim (P = 1.07·10⁻²¹⁸). There was also a significant effect of group (P = 5.85·10⁻¹⁷); all pairwise post hoc group comparisons, other than none versus mild, were significant. 

Figure 6.

 

Foveal layer boundary heights, relative to the outer boundary of Bruch's membrane, in term (blue), none (green), mild (yellow), and severe (red) subjects. At the left, the mean height of each layer boundary is plotted as a function of eccentricity. At the right, mean layer height at the fovea and rim is plotted as determined by fit of the layers to Equation 4 (with shared x₀ and the ∆s in z₀ + a and z₀ and taken as measures of layer thickness at the fovea and rim, respectively). From the inner retina (light shades) to the outer retina (dark shades), the demarcated layers are the NFL, GCL, IPL, INL, OPL, ONL, IS, OS, RPE, and BM.

The thickness of the subfoveal retina, and the depth and breadth of the overlying fovea, are plotted for each group in Figures 7A and 7B. The sum of thickness and depth is quite constant across groups, but preterm birth and then increasing severity of ROP were associated with increasing central retinal thickness and decreasing foveal depth. In respective analyses, we detected significant effects of group in the foveal depth (P = 1.83·10⁻³⁴) and foveal breadth (P = 3.06·10⁻⁶) measurements. Consistent with the above analysis of layer thicknesses, all pairwise post hoc group comparisons, other than none versus mild, were significant in the depth data. However, our pairwise intergroup comparisons of breadth detected no significant differences among the term, none, and mild groups, but found that all three were significantly narrower than the severe group. We also evaluated foveal slope, which we took as a ratio of depth to breadth (0.6a/b). In accordance with the depth and breadth data from which it was derived, there was a significant effect of group and all pairwise post hoc tests were significant, except none versus mild (which was moderately significant, P = 0.0307). In addition, we evaluated the variability in x₀ in the fits of Equation 4, wherein it was free to vary, for the seven most proximal layer boundaries (we excluded more distal boundaries since these are frequently flat or nearly so). Variability increased with prematurity and increasing ROP severity such that there was a significant effect of group (P = 0.000382) and, in pairwise post hoc testing, variability in the severe group significantly exceeded that in all others; the term versus mild difference was moderately significant (P = 0.0106). 

Figure 7.

 

Mean (± standard deviation [SD]) parameters of the fovea and subfoveal retina. (A) Thickness of the subfoveal retina and foveal depth. All pairwise intergroup differences in both parameters were significant or marginally significant (none versus mild). (B) Foveal breadth. Brackets indicate significant intergroup differences. (C) Foveal slope. All intergroup comparisons other than that bracketed (i.e. none versus mild) were significant. (D) Subfoveal laminar irregularity. Brackets indicate significant intergroup differences (none versus mild was moderately significant).

Inter-relations among depth, breadth, and VA are plotted in Figure 8. We evaluated whether either the depth or breadth of the fovea were significant predictors of VA (logMAR), beyond their association with group, using a model with depth and breadth factors. The model fit the data well (r = 0.869), but the only statistically significant effect was depth (P = 0.000595), which had a coefficient of ‒1.08 logMAR·mm⁻¹; stated simply, every additional millimeter of foveal depth is associated with approximately a log unit improvement in VA. This suggests that, once accounting for changes to the depth of the pit, additional effects of ROP on acuity are minor, and that the above-reported significant effect of group on VA is, in fact, mostly a reflection of the fact that determinants of pit depth also mediate VA which, in turn, is correlated with group. 

Figure 8.

 

Relations among foveal depth, foveal breadth, and visual acuity in term (blue circles), none (green triangles), mild (yellow diamonds), and severe (red squares) subjects. The shaded regions within the stippled lines demark the 95% prediction interval for term. (A) The black line is an orthogonal regression between depth and breadth (there is no relationship). (B) Lines (colors for each group as for symbols) are results of linear modeling. There was no significant association between breadth and visual acuity. (C) Depth was significantly associated with visual acuity.

Figure 9 plots the breadth and depth values obtained using a more involved procedure³³ (based upon fitting the ILM to the sum of three Gaussians) as a function of the depth and breadth values obtained directly from the parameters of a single Gaussian (in this case, FWHM for breadth). The correlation between the two methods was good for breadth (r = 0.764) and especially for depth (r = 0.992), but neither relationship fell directly on the diagonal. That is, both methods are measuring the same things, but in subtly different ways. 

Figure 9.

 

Relations of foveal depth and breadth assessed using the single Gaussian model adopted herein (abscissa) and the sum of three Gaussians method of Breher et al.³³ (ordinate) for term (blue circles), none (green triangles), mild (yellow diamonds), and severe (red squares) subjects. Solid lines have unity slope (on the primary axes) and zero intercept; stippled lines are orthogonal regressions through the data.

For the 46 subjects who underwent both the standard (i.e. central pupil entry) as well as the directional (i.e. offset pupil entry) OCT, between-group comparisons of the mean positions of the retinal laminae are shown in Figure 10. In addition to the expected presence of proximal, postreceptor retinal layers in the preterm subjects, the central retinal ONL was also systematically thicker with increasing ROP severity (cyan shaded rectangle in top right panel of Fig. 10, replotted with data from all subjects in Fig. 11). We evaluated the change in layer thicknesses at the rim of the fovea produced by the D-OCT using a model with layer and pupil entry position factors, as well as factors for the interactions of group, layer, and pupil entry. Of course, there was still a highly significant effect of layer (P < 2.23·10⁻³⁰⁸) but now there was no significant effect of group, implying that the persistent consequences of prematurity and ROP on retinal layer thickness are mainly in the subfoveal retina. However, because the D-OCT reveals the FLH, we would expect that ONL and OPL thickness measurements would, at minimum, differ using this method from measurements made using standard OCT and, indeed, we found a significant layer × pupil entry interaction (P = 2.57·10⁻¹⁴); other boundaries shifted only trivially between standard and D-OCT, indicating that it was the enlarged FLH that drove this effect. If the FLH differs among the groups tested—the expected result if ROP caused changes in the routing of the cone axons to their bipolar cell targets—there would also be a group × layer × pupil entry interaction, and there was (P = 2.33·10⁻¹¹). Notably, inspection of the layers revealed that it was the mild subjects, not the severe subjects, that differed most from term. In other words, D-OCT altered measurements of OPL thickness, revealing the FLH, but the only group with a significantly altered FLH was the mild group. 

Figure 10.

 

Various comparisons (i.e. overlays) of mean retinal boundaries in the term (blue), none (green), mild (yellow), and severe (red) groups obtained using standard (i.e. central pupil; solid lines) and directional (i.e. eccentric pupil; dashed lines) OCTs. The yellow highlighted region in each bottom panel represents the fiber layer of Henle, as revealed by directional OCT. The cyan highlighted region in the top right panel is partially replotted in Figure 11 at increased scale and including all subjects.

Figure 11.

 

Mean height of the boundaries of the ONL (left panel, replotted from Figure 9 with additional data), and their difference, ONL thickness (right panel), in term (blue lines), none (green lines), mild (yellow lines), and severe (red lines) subjects.

Additionally, there was a slight but significant effect of age (P = 0.0364), which translated to an approximately 5.4 µm predicted increase in foveal rim thickness per decade. 

Scanning Light Ophthalmoscopy

We evaluated the sidedness (S) and spacing (ICD) data from each experimental group (Fig. 12) using a model with eccentricity as a continuous factor. Sidedness did not differ significantly among groups: There were no group differences in either the standard deviation in the number of sides nor the proportion of six-sided VT cells. Neither was there a group difference in the standard deviation of ICD. An additional analysis of cone density (D) versus eccentricity (Fig. 12) likewise found no significant differences between groups, although there was, of course, a highly significant effect of eccentricity (P < 2.23·10⁻³⁰⁸). Thus, the arrangement of the parafoveal cones did not seem to differ much between groups. 

Figure 12.

 

Cone density, measured using MAORI, plotted as a function of eccentricity in term (blue lines), none (green line), mild (yellow lines), and severe (red lines); as indicated, the term, mild, and severe lines at higher eccentricities are replotted from Ramamirtham et al.³⁴ (none of the subjects were not reported therein). Gray data are mean cone densities from histology, replotted from Curcio et al.’s⁵ Figure 6 and the cyan and purple shaded regions are normal ranges from respective adaptive optics studies by Song et al.⁴⁷ and Zhang et al.⁹ Whereas severe cone density was the lowest overall, this difference was not significant. The preterm subjects’ data appeared to be exiting the normal range at the lowest eccentricities.

Discussion

The processes that lead to the development of the fovea include movement of all retinal cell types, including both centripetal packing of photoreceptors and centrifugal spreading of inner retinal neurons and glia. Despite the high variability in both the centripetal and centrifugal processes, it is well-recognized that preterm birth in general, and ROP in particular, alter these processes, leading to a shallow fovea. In particular, the distance from the trough of the fovea to the rim has been reported to be reduced in preterm retinae because there is increased thickness—and, indeed, presence—of inner retinal layers. Because of this, it has been argued that ROP arrests the development of the fovea.¹⁴ It may be that the foveal avascular zone (FAZ) is involved in the formation of the fovea.⁴⁰ Certainly, the fact that ROP—a disease clinically characterized by abnormalities of the retinal vasculature—also impacts foveal development is consistent with this notion. Presumably, ROP was active in these eyes at the preterm ages when the fovea was still forming and disrupted the normal centrifugal movement of postreceptor cells. Whether this failure is concomitant or merely coincident to the vascular abnormalities remains unknown. 

In a large population of term and preterm subjects, we automated parameterization of central retinal features. We applied this approach to eyes without or with a history of ROP, including ROP so severe as to require treatment. Our results are in good agreement with those of Sjöstrand et al.,¹⁴ who performed detailed, manual, landmark-based measurements in a small number of subjects similarly stratified. For instance, our results confirmed that the fovea is shallow in preterm groups, and that the severity of antecedent ROP is associated with the degree of foveal diminishment. We derived a meta-surface, to which we fit a difference of three-dimensional Gaussians, to objectively define the parafovea. In the most severe cases of ROP (i.e. treated), the parafoveal area was significantly larger. We also found that the thickness of the central ONL systematically increased in that it was thinnest in term, thicker in none, yet thicker in mild, and thickest in severe. These results hint that the centripetal movement of foveal cones, upon which foveal cone packing is thought to depend in normal development, is not arrested in ROP but is enhanced. Indeed, we have shown that there are marked differences between term-born eyes and eyes with ROP in the growth (or stretching) of the posterior pole at the ages when the fovea is remodeling.⁴¹ Of course, other explanations, such as differences in intraocular pressure or FAZ development are also possible and, given the results of the cone counting (below), perhaps more likely. 

Furthermore, we found a significant association between VA and pit depth that seemed to account for the association with group, implying that the determinants of pit depth, rather than ROP severity alone, determine VA. We also found that the parafoveal retina becomes slightly thicker with age; the rate in our model, 5.4 µm·decade⁻¹, is in excellent agreement with the findings of a multicenter study of retinal thickness and age.⁴² 

We detected a slight but significant thinning in the FLH, relative to term, in the mild subjects. Sjöstrand et al. made the case that the increase in ONL thickness in ROP is consequent to decreased centrifugal displacement of developing bipolar cells, which, in turn, decreases the centrifugal pull on the axons of the central cones, rather than an increase in cone density.⁴³ This, they argue, would leave the cone axons and cell bodies oriented relatively vertically and lead to an increase in ONL thickness. However, the fact that ONL thickness increased monotonically with increasing severity of ROP whereas FLH thickness did not, stands as a counterpoint. Possibly, treatment of the peripheral retina normalized the development of the FLH in severe subjects, which may also be related to the larger parafovea that we identified exclusively in the severe group. 

In our AO-SLOs of the cone photoreceptor mosaic in the parafovea, we found no difference in the density of the cone photoreceptors. This stands in contrast to our previous AO-SLO observations of the cones at greater eccentricities, although it is in agreement with AO-OCT cone counts at the same locations.³⁴ Also in disagreement with our more eccentric AO-SLO data, we did not find greater irregularity in cone packing, as measured by the standard deviation in the number of sides or in intercone distances. Notably, in these more central images of the cone mosaic, there are few to no rods, whereas in our earlier, more eccentric measurements, rods were abundant and, in most images, likely outnumbered the cones; this may have made counting the cones difficult. That said, a full complement of cone photoreceptors is consistent with the normal FLH thickness that we detected in the severe subjects using OCT. Thus, we suspect that the number and distribution of photoreceptors in the ROP retina is approximately normal, although their (peripheral) organization may not be. Alternatively, because the parafovea is larger and the ONL thickest in the severe group, it is also possible that cones may have moved from further eccentricities to pack the fovea, thus leaving them sparser in the near periphery. 

With our small sample of preterm subjects, it is hard to say anything with confidence, but we are inclined to favor the hypothesis that our earlier low peripheral cone counts occurred because counting peripheral cones was difficult (as opposed to the cones being absent) for two reasons. First (as noted above), the AO-OCT cone counts at the same locations did not recapitulate the AO-SLO results. Second, our results from the rat, whose retina is similar to that of peripheral human retina, and in whom we use oxygen-induced retinopathy to model ROP,⁴⁴ do not show evidence of ROP-induced photoreceptor cell loss.⁴⁵ 

In relating our AO-SLO results to our OCT findings, we note (again) that there was no significant change from normal cone density in our small cohort of preterm subjects. Indeed, although we unfortunately lacked the ability to resolve the most central cones, the preterm cone densities in general and the severe cone densities in particular were, if anything, trending lower. This observation stands in contrast to the thicker subfoveal ONL measurements that we obtained. Thus, if some process of centripetal packing is, in fact, enhanced by preterm birth and ROP, that process is not a straightforward one. Certainly, the posterior pole is growing differently than normal in ROP eyes,⁴¹ so this requires further investigation. There are innumerable measures of cone packing available,³⁶ and sundry issues in cone identification and imaging,⁴⁶ and so it is likely that some combination of increased sample size, alternate measures of cone packing, better cone centroid detection, and improved AO correction, will yet lead to detection of subtle differences in term and preterm subfoveal cone arrangement that explain our finding. 

Finally, we note that there are some differences in key parameters obtained using our simpler, single Gaussian model and a more complex model which uses a sum of three Gaussians followed by some elegant mathematical “gymnastics,” but, overall, the simpler model seems to work just fine. Absolute differences in breadth, for instance, are probably mostly driven by the fact that the three-Gaussian model is measuring neither FWHM nor rim-to-rim (see Fig. 2), but something that our simple model cannot directly capture. Especially if one is mainly interested in measuring foveal depth which, at least in the eyes of former preterms seems to be the most systematic foveal abnormality, there appears little need for a more complicated analysis; any advantages it offers are probably negated by the easier acquisition and interpretation supplied by the use of a simpler model. 

Acknowledgments

Supported by Boston Children's Hospital Ophthalmology Foundation (J.D.A.), Department of Defense W81XWH-14-1-0592 (J.D.A.), Massachusetts Lions Eye Research Fund (R.M.H.), and National Institutes of Health EY010597 (A.B.F.). 

Disclosure: J.D. Akula, None; I.A. Arellano, None; E.A. Swanson, None; T.L. Favazza, None; T.S. Bowe, None; R.J. Munro, None; R.D. Ferguson, Physical Sciences, Inc. (E, P); R.M. Hansen, None; A. Moskowitz, None; A.B. Fulton, None