Analysis of Cone Mosaic Reflectance Properties in Healthy Eyes and in Eyes With Nonproliferative Diabetic Retinopathy Over Time

Abstract

Purpose: We investigate the reflectance properties of the cone mosaic in adaptive optics (AO) images of healthy subjects and subjects with nonproliferative diabetic retinopathy (NPDR) over time.

Methods: We acquired images of the parafoveal cone mosaic over 5 years in 12 healthy subjects and in six patients with mild NPDR. We analyzed the parameters of the cone intensity histogram distribution (mean, SD, and skewness), two metrics of the cone mosaic texture (sharpness and entropy), and two novel metrics (cone/intercone intensity and slope of the variogram). Each metric was calculated on the same four retinal locations (200 × 200 μm areas, 2° from the fovea along the four meridians) over time for each subject.

Results: The histogram distributions of cone intensities were similar between the two study groups. However, the cone/intercone intensity, slope of the variograms and entropy showed a significant difference between healthy and NPDR subjects (P = 0.036, P = 0.002, P = 0.014, respectively). All parameters, except for mean cone intensity, did not change with time in this study.

Conclusions: We observed significant differences in cone mosaic reflectance properties between healthy eyes and eyes with NPDR, in its spatial organization and in its intensity, especially between directional and nondirectional backscattering. We introduced a novel method for the study of the spatial distribution of cone reflectance, the variogram, which was able to quantify differences of the spatial dependence of cone intensities over a short range between NPDR and healthy eyes.

Diabetic retinopathy (DR) is a frequently occurring complication of diabetes mellitus and is one of the leading causes of visual impairment among adults globally.^1–3 It is a multifactorial disease apparently involving the retinal neuronal cells early in the course of disease onset and progression (i.e., neurodegenerative theory).^4–10 The neurodegenerative changes are apoptosis of several populations of retinal cells, including photoreceptors, bipolar and ganglion cells, and astrocytes. Structural and functional impairments of the neural tissue also have been supposed to contribute to the earliest alterations of the vascular structures.^11–13 

Recent clinical studies have focused on the investigation of the photoreceptor layer showing abnormalities of the cone packing density arrangement in the parafoveal region of subjects affected by diabetes, even before any sign of retinopathy was detected.^14–16 Thanks to recent advances in high-resolution retinal imaging it also has been shown how a ‘missing' cone photoreceptor as seen in an image acquired with conventional fundus illumination or AOSLO imaging systems cannot be attributed directly to the death of the cell itself, but rather to a disruption of its wave-guiding and light-reflecting ability.^17,18 

While the spatial arrangement of the cones still is the most studied property of the cone mosaic,^19–20 their light-reflecting properties recently have become the subject of an increasingly number of studies.^21–29 Nonetheless, the investigation of the light-reflecting properties of cones in retinal diseases still is limited,^30–34 even if in some clinical cases the cone reflectance has been the only apparent feature of the cones that distinguished a healthy cone mosaic from a mosaic with altered functionality.^35,36 

The study of cone reflectance deserves attention, as it potentially could lead to a deeper understanding of the cone cell physiology or pathophysiology, which cannot be inferred by their spatial distribution alone.^35–37 In this view, the development of quantitative metrics, which can be automated for extending the benefits of high-resolution retinal imaging to large populations, is expected for capturing clinically valuable information.³⁸ 

We tested a method for evaluating possible differences in the reflectance properties of the parafoveal cone mosaic in adults with DR compared to age-matched healthy subjects. Furthermore, we tracked these properties over 2 to 5 years. The scope was to carry out a preliminary evaluation of the feasibility of using adaptive optics (AO) imaging biomarkers based on cone intensity for clinical purposes. 

Methods

Subjects

All research procedures described adhered to the tenets of the Declaration of Helsinki. The protocol was approved by the local ethical committee (Azienda Sanitaria Locale Roma A, Rome, Italy) and all subjects recruited gave written informed consent after a full explanation of the procedure. 

Patients with a diagnosis of type 1 diabetes mellitus and age-matched healthy volunteers participated in this study. Inclusion criteria were at least 18 years of age and, for patients with diabetes, mild signs of nonproliferative DR (NPDR) according to the Early Treatment of Diabetic Retinopathy Study (ETDRS) severity scale,³⁹ 20/20 or better uncorrected or corrected distance visual acuity (UDVA or CDVA). Mild NPDR was defined as the presence of at least one microaneurysm and/or mild hemorrhages. Exclusion criteria were astigmatism higher than 2.50 diopters (D), the presence or a history of maculopathy (macular edema) or any other ocular disease, including lens opacity, or previous eye surgery (including laser treatments). Control subjects were recruited as healthy volunteers who had no history of systemic diseases. We will refer to the two groups as NPDR and controls, respectively. All subjects received a complete eye examination, including retinal imaging using a Spectralis SD-OCT (Heidelberg Engineering GmbH, Heidelberg, Germany) to exclude the presence of macular edema. 

A total of 18 subjects were involved in the study, with a minimum of 2 years to a maximum of 5 years of separation between the baseline images and the last images. Of these subjects, 12 were healthy volunteers (mean age, 33 ± 7 years at baseline; sex, 5 M–7 F) and 6 were patients with mild NPDR (mean age, 40 ± 10 years at baseline; sex, 4 M–2 F). The diabetes duration at baseline was 16 ± 4 years. 

Image Acquisition and Processing

A flood-illuminated AO retinal camera (rtx1; Imagine Eyes, Orsay, France) was used to acquire images of the cone mosaic. All imaging sessions were conducted after dilating the pupil with one drop of 1% tropicamide. During imaging, fixation was maintained by instructing the patient to fixate on the internal target of the instrument moved by the investigator. At each retinal location, a sequence of 40 frames was acquired by illuminating a retinal area subtending 4° of visual angle in the right eye of each subject; images were acquired at several locations in the central retina covering an area of at least 5° × 4° centered on the preferred locus of fixation (coordinates x = 0° and y = 0°). All images were acquired focusing the light source at the center of the pupil, that is, inside the first four Purkinje images of the cornea, as described previously.⁴⁰ All subjects had several images taken at different time intervals, with a minimum of 1 year separation between the imaging sessions. 

Using the method described in Mariotti et al.,⁴⁰ the flat-field was calculated by averaging all unregistered frames for each subject; then, the individual frames were divided by the flat-field and registered. All image processing operations were performed using Matlab (version 9.1, R2016b; The Mathworks, Inc., Natick, MA, USA). The final images subtending 4° were montaged using the montage tool in i2k Retina Pro (DualAlign LLC, Clifton Park, NY, USA) with the aim of producing wide-field images to be used as reference to determine the position of the sample areas. 

For each eye, four areas of the cone mosaic were analyzed. The sample areas were chosen to have a size of 200 × 200 μm and to be located at 2° from the foveal center along the main retinal meridians (Fig. 1). The distance of 2° was chosen as a compromise between the ability of the instrument to resolve cones (not too close to the fovea) and the presence of rods, which at greater eccentricities start to be detectable and disrupt the cone mosaic.⁴¹ The corrected magnification factor (RMF_corr) was calculated for each eye to correct for the differences in optical magnification and, thus, retinal image size between eyes.^42–44 

Figure 1

 

Example of montage image produced with the i2K Retina software on a control subject. The montage images were used to determine the position of the center of the fovea (green dot) and of the four selected retinal locations (green boxes), which were subsequently selected on the 4° × 4° images. The distance from the fovea of the four sample areas is 2° and the size is 200 × 200 μm.

The processed images subtending 4° were registered through a two-step registration process. While the first step (coarse registration) involved normalized cross correlation and used all the structures present in the images (i.e., the cones and blood vessels), the second step refined the alignment by looking at cones only; it divided the whole images in a 3 × 3 grid and aligned them by tracking the cones that are on average the brightest over time in each grid sector. This approach permitted us to look at the same exact location of the cone mosaic over time.⁴⁰ From these registered images, the sample areas were selected at the locations determined from the montages. The sample images were normalized by dividing them with their total intensity, i.e., the sum of all the pixel intensities, and multiplying them by the total number of pixels in the area, which gave a mean pixel intensity value equal to 1.⁴⁰ No further contrast corrections were applied. 

The final data set consisted of a total of 53 observation sessions, with a total of 212 final images selected for analysis. Figure 2 shows a selection of sample images (for an overview of the study population and data set, see Supplementary Table S1). 

Figure 2

 

A selection of images used in the study. All images in this Figure were acquired at 2° temporal from the fovea at the baseline time. Scale bar: 100 μm.

In addition, for each subject one image focused at the inner retina (i.e., the blood vessel layer) was acquired, to compare the low frequency spatial changes in the image intensity between the inner retina and the photoreceptor layers. This was done to investigate the origin of the variable spatial pattern of cone intensities in AO flood illumination images of the photoreceptor mosaic. To achieve this, the images of the inner retina and photoreceptor layers were registered; therefore, for each pair of images, sample areas of 200 × 200 μm on the same retinal location of the two layers were selected. The inner retinal images then were low pass filtered using a mask on the Fourier spectrum that kept only the frequencies smaller than the frequency of cone packing, that is, 145 cycles per degree (Yellott's ring).⁴⁵ The low pass filtered sample areas of the inner retina then were subtracted from the photoreceptor layer sample areas (Fig. 3). 

Figure 3

 

Illustration of the method used for the evaluation of the contribution of the inner retina on the cone mosaic intensity. Sample areas of the cone mosaic and of the image focused on the inner retina (200 × 200 μm) are selected on the same retinal location. The image of the inner retina is low pass filtered, to exclude all frequencies that correspond to the cones (Yellott's ring) or smaller features. The low pass filtered image then is subtracted from the cone layer image.

Cone Detection

We used the same method for analysis of the cone mosaic and reflectance presented previously.⁴⁰ Here, we give a brief summary for the purpose of clarity. 

The sample areas at the same location were averaged over time for each subject. Cones were detected and segmented on the average images with the entirely automated algorithm presented by Chiu et al.⁴⁶ We used the segmentations of the cone apertures on the images averaged over time based on previous work.^26,40 Previous studies have shown that in eyes with DR a fraction of cones presents a loss of reflectivity. ^14,19 If a cone hypothetically loses its ability to reflect light (because of cone death or other causes), the use of the coordinates on the average of different times allows us to monitor the reflectance of that point on the cone mosaic even if in one or more images in the series that cone was not reflecting. In this way, we aimed at approaching the best case scenario where cones are not detected only if they are consistently dark throughout time, in which case the only way to determine if they still are structurally intact would be a localized function analysis on the individual cones.¹⁸ Therefore, the performance of the detection is increased over a single image not only on diseased but also on healthy subjects, and the skewness towards bright cones is limited as much as possible (Fig. 4). In our particular case, the averaging over time was justified by the fact that mild NPDR is not known to affect the position of the individual cones directly. The presence of more severe conditions (e.g., advanced stages of retinal dystrophies) that can cause the migration of cones to different positions presumably would not allow the use of the same averaging process. 

Figure 4

 

Visualization of image processing procedure, as performed on one sample area of a control subject. The images acquired at three different years at the same retinal location were averaged (bottom left), and the blood vessel (blue) and cones (in green) were segmented on the average image (bottom center). A detail (bottom right) shows the cone centers (green dots), the segmentation of the cone apertures (green masks) and part of the vessel segmentation (blue mask). Scale bar: 100 μm.

Where present, blood vessel shadows were excluded using the semiautomated method from previous studies⁴⁷ and only the cones detected in the areas devoid of vessels were analyzed (Fig. 4). The use of the average images for the segmentation of the blood vessels allowed the exclusion of vessel profiles that might have developed later in the imaging process, but that still would have appeared in the average images. This can be called a “conservative approach,” in the sense that cones are excluded from analysis even if they were obscured by the vessels only at one time point. If microaneurysms where present, they would have appeared as dark spots on the cone mosaic.^14,48 We did not observe such features at our designated locations, but we want to point out that they would have been equally segmented and excluded on the base of their dark appearance with the same procedure that excludes the blood vessels. 

As the area occupied by the vessels is excluded, the analysis could be performed at the exactly defined retinal locations with clearly detectable cones, without the need to move the selection areas to a more suitable location.¹⁴ 

Parameters

Cone Intensity

All sample images were normalized to have a mean intensity value of 1 and cone intensity was measured as the mean value of the pixels inside the segmentations of the cone apertures. Parameters of the cone intensity histogram, such as the mean value, standard deviation, and skewness, were calculated for each sample area. 

The cone detection algorithm detects and separates the cone apertures from the surrounding pixels by assigning weights to the intensity gradients (light-to-dark and dark-to-light) and finding the shortest path around the local maximum.⁴⁶ We introduced a parameter to estimate the ratio of light backscattered from the cone apertures to that scattered by the remaining space, which we refer to here as “intercone space.” For each image, this parameter was defined as the ratio of the mean intensity value of the pixels inside all the cone segmentations and the mean intensity value of all the pixels outside the cone segmentations (i.e., the area not covered by colored masks in Fig. 4) and gives a single numerical value for each sample area. 

In addition, we introduced a new methodology for evaluation of the spatial correlation of cone reflectance by introducing the use of variograms, which is a technique already used in geostatistics.⁴⁹ We defined here the semivariance of cone reflectances as  

\(\def\upalpha{\unicode[Times]{x3B1}}\)\(\def\upbeta{\unicode[Times]{x3B2}}\)\(\def\upgamma{\unicode[Times]{x3B3}}\)\(\def\updelta{\unicode[Times]{x3B4}}\)\(\def\upvarepsilon{\unicode[Times]{x3B5}}\)\(\def\upzeta{\unicode[Times]{x3B6}}\)\(\def\upeta{\unicode[Times]{x3B7}}\)\(\def\uptheta{\unicode[Times]{x3B8}}\)\(\def\upiota{\unicode[Times]{x3B9}}\)\(\def\upkappa{\unicode[Times]{x3BA}}\)\(\def\uplambda{\unicode[Times]{x3BB}}\)\(\def\upmu{\unicode[Times]{x3BC}}\)\(\def\upnu{\unicode[Times]{x3BD}}\)\(\def\upxi{\unicode[Times]{x3BE}}\)\(\def\upomicron{\unicode[Times]{x3BF}}\)\(\def\uppi{\unicode[Times]{x3C0}}\)\(\def\uprho{\unicode[Times]{x3C1}}\)\(\def\upsigma{\unicode[Times]{x3C3}}\)\(\def\uptau{\unicode[Times]{x3C4}}\)\(\def\upupsilon{\unicode[Times]{x3C5}}\)\(\def\upphi{\unicode[Times]{x3C6}}\)\(\def\upchi{\unicode[Times]{x3C7}}\)\(\def\uppsy{\unicode[Times]{x3C8}}\)\(\def\upomega{\unicode[Times]{x3C9}}\)\(\def\bialpha{\boldsymbol{\alpha}}\)\(\def\bibeta{\boldsymbol{\beta}}\)\(\def\bigamma{\boldsymbol{\gamma}}\)\(\def\bidelta{\boldsymbol{\delta}}\)\(\def\bivarepsilon{\boldsymbol{\varepsilon}}\)\(\def\bizeta{\boldsymbol{\zeta}}\)\(\def\bieta{\boldsymbol{\eta}}\)\(\def\bitheta{\boldsymbol{\theta}}\)\(\def\biiota{\boldsymbol{\iota}}\)\(\def\bikappa{\boldsymbol{\kappa}}\)\(\def\bilambda{\boldsymbol{\lambda}}\)\(\def\bimu{\boldsymbol{\mu}}\)\(\def\binu{\boldsymbol{\nu}}\)\(\def\bixi{\boldsymbol{\xi}}\)\(\def\biomicron{\boldsymbol{\micron}}\)\(\def\bipi{\boldsymbol{\pi}}\)\(\def\birho{\boldsymbol{\rho}}\)\(\def\bisigma{\boldsymbol{\sigma}}\)\(\def\bitau{\boldsymbol{\tau}}\)\(\def\biupsilon{\boldsymbol{\upsilon}}\)\(\def\biphi{\boldsymbol{\phi}}\)\(\def\bichi{\boldsymbol{\chi}}\)\(\def\bipsy{\boldsymbol{\psy}}\)\(\def\biomega{\boldsymbol{\omega}}\)\(\def\bupalpha{\bf{\alpha}}\)\(\def\bupbeta{\bf{\beta}}\)\(\def\bupgamma{\bf{\gamma}}\)\(\def\bupdelta{\bf{\delta}}\)\(\def\bupvarepsilon{\bf{\varepsilon}}\)\(\def\bupzeta{\bf{\zeta}}\)\(\def\bupeta{\bf{\eta}}\)\(\def\buptheta{\bf{\theta}}\)\(\def\bupiota{\bf{\iota}}\)\(\def\bupkappa{\bf{\kappa}}\)\(\def\buplambda{\bf{\lambda}}\)\(\def\bupmu{\bf{\mu}}\)\(\def\bupnu{\bf{\nu}}\)\(\def\bupxi{\bf{\xi}}\)\(\def\bupomicron{\bf{\micron}}\)\(\def\buppi{\bf{\pi}}\)\(\def\buprho{\bf{\rho}}\)\(\def\bupsigma{\bf{\sigma}}\)\(\def\buptau{\bf{\tau}}\)\(\def\bupupsilon{\bf{\upsilon}}\)\(\def\bupphi{\bf{\phi}}\)\(\def\bupchi{\bf{\chi}}\)\(\def\buppsy{\bf{\psy}}\)\(\def\bupomega{\bf{\omega}}\)\(\def\bGamma{\bf{\Gamma}}\)\(\def\bDelta{\bf{\Delta}}\)\(\def\bTheta{\bf{\Theta}}\)\(\def\bLambda{\bf{\Lambda}}\)\(\def\bXi{\bf{\Xi}}\)\(\def\bPi{\bf{\Pi}}\)\(\def\bSigma{\bf{\Sigma}}\)\(\def\bPhi{\bf{\Phi}}\)\(\def\bPsi{\bf{\Psi}}\)\(\def\bOmega{\bf{\Omega}}\)\begin{equation}\gamma \left( r \right) = {1 \over {2N}}\mathop \sum \limits_{i,\ j = 1}^N {\left( {{I_i} - {I_j}} \right)^2}\end{equation}

where Display Formula\(N\) is the number of cone pairs in one sample image separated by distance Display Formula\(r\), Display Formula\({I_i}\) and Display Formula\({I_j}\) are the intensities of two cones separated by distance Display Formula\(r\). For each image, the semivariance Display Formula\(\gamma \left( r \right)\) was calculated for a set of separation distances r from 0 to 200 μm (the size of the sample area) using bins centered at multiples of 6 μm, which is the mean cone separation at 2° from the fovea.⁴² The semivariance, Display Formula\(\gamma \left( r \right)\), then was normalized by the total variance of the cone reflectance in the image to give Display Formula\({\gamma _{norm}}\left( r \right)\), which was plotted against cone separation in micrometers. The plot of Display Formula\({\gamma _{norm}}\left( r \right)\) against the separation is called the variogram. The variogram gives a measure of the degree of spatial dependence of the cone reflectance within the sample area (Fig. 5).  

Figure 5

 

Visualization of the calculation and meaning of a variogram of intensities on a cone mosaic. On the left, a detail of a sample area with cones divided in rings of increasing radius, which correspond to increasing distances, from the cone in the center (blue cross; Scale bar: 30 μm). The squared difference of the intensities of all the cones that are separated by distance r is averaged for different distances, producing a variogram curve (on the right, the variogram calculated on the whole sample area shown on the left). In this variogram, \({\gamma _{norm}}\left( r \right)\) flattens at 1 at approximately 30 μm (vertical blue line), which means that the reflectance of cones separated by distances greater than 30 μm is not correlated (i.e., on the left, cones that are in clusters smaller than the ring highlighted in cyan have similar reflectances, while in bigger clusters the reflectance of the cones is not correlated anymore). The variogram curve is fitted linearly at short separation distances (r < 20 μm) and the slope of the linear fit is taken as a parameter.

If the variogram assumes values lower than 1, it means that the cones intensity, separated by the corresponding distance is spatial dependent; otherwise, if the variogram assumes values close to 1, the cones reflectances are not correlated. Accordingly, the changes in the shape of the curve, that is, the presence or absence of minima or maxima, give valuable information on the spatial pattern of cone reflectance (more comprehensive description of the meaning of the curve is given in the Supplementary Material). 

There are many functions based on theoretical models that have been proposed in the literature to fit experimental variograms. As the process of choosing a fitting model for experimental variograms still does not follow a universally agreed method,⁴⁹ we focused on the analysis of the portion of the curve close to the origin. This approach allowed us to retrieve information on the degree of correlation of reflectance between cones at short range distance. We performed a linear fit on the normalized semivariance for short separation distances, that is r < 20 μm for all images and all subjects (full methodology is described in the Supplementary Materials). Therefore, the slope of the linear portion of the variograms was chosen as a metric and calculated for all the sample images (Fig. 5b). 

In addition, where available, we calculated the variograms of cone intensities after subtraction of the contribution of the low pass filtered inner retina image. 

Cone Mosaic Texture

In addition to the parameters related to the intensity of the cones, we also evaluated the textural characteristics of the cone mosaic as a whole using all the pixels of the images, that is, without segmenting between cone apertures and intercone space. Two metrics, such as sharpness and entropy, commonly used to assess image quality,⁵⁰ were used for this scope. Sharpness was defined as originally proposed by Muller and Buffington,⁵¹ as  

\begin{equation}{S_{\rm{\Gamma }}} = \mathop \sum \limits_{x,y} {{\rm{\Gamma }}_{{\rm{D}}1}}\left[ {{\rm{\ }}I(x,y)} \right]\end{equation}

where  

\begin{equation}{{\rm{\Gamma }}_{D1}}\left( I \right) = {1 \over {12}}{\left( {I - \gamma \mathop I\limits^ - } \right)^\beta }\end{equation}

where Display Formula\(I\) is the intensity of the pixels in the image and Display Formula\(\mathop I\limits^ - \) is the mean pixel intensity. For the numerical parameters we used the values Display Formula\(\gamma = 0.98\) and Display Formula\(\beta = 2\), retrieved as optimized values in previous work.⁵⁰ Entropy, which is a statistical measure of randomness that can be used to characterize the texture of images, was defined as:  

\begin{equation}E = - \mathop \sum \limits_i {p_i}\ {\log _2}\left( {{p_i}} \right)\end{equation}

where Display Formula\({p_i}\) is the probability of the Display Formula\({i^{th}}\) pixel intensity value (e.g., i = 0,1… 255 for 8-bit images). The probabilities were approximated by the histogram of the image. A more comprehensive description of these texture parameters and of their application on cone mosaic images can be found in our previous work.⁵⁰  

The pixels marked as vessel shadows were excluded from the analysis, to consider only the properties of the cone mosaic. 

Statistical Analysis

To determine the effect of different factors on the parameters, we performed a linear mixed model analysis. This type of analysis was chosen to be able to deal with missing observations and to introduce random effects given by the intersubject variability.^19,42,52,53 

The parameters calculated on all the sample images were processed in two ways to determine if there was a significant interaction between the parameters and either the condition of the subject (control or NPDR), time (year), or the retinal location (2° temporal, superior, nasal, and inferior). They were analyzed either considering the four retinal locations separately for each subject or the average of the parameters over the four retinal locations for each subject at each time point. This approach allowed us to determine if it was possible to use only one global “clinically useful value” for the parafoveal cone mosaic reflectance that would allow discrimination between healthy and diabetic retinas.¹⁵ 

Data were tested for normality with the Shapiro-Wilk test. The 95% confidence intervals (95% CI) were calculated as Display Formula\( \pm 1.96\sigma /\sqrt n \), where Display Formula\(n\) is the number of observations and Display Formula\(\sigma \) the standard deviation. Statistical analysis of the results was performed with SPSS software (version 23; SPSS, Inc., Chicago, IL, USA). In both cases, all possible interactions between the factors were tested. The different mixed models for all parameters were compared using the Bayesian Information Criterion (BIC) for the penalized likelihood, and the model that minimized it was chosen as the best fitting one. 

Based on the results from previous work describing the average changes of the mean cone intensity over a period of 3 years,⁴⁰ sample size calculation was performed to detect a difference of 0.002 between the average cone intensity for the control and NPDR groups, at a significance level of 5% and a power of 98%, assuming a standard deviation of 0.001. The sample size of the study was 18 cases (allocation ratio of 2:1). 

Results

Figure 6 shows a selection of cone intensity histograms (calculated on the same sample images shown in Fig. 2). For all the parameters, the best fit linear mixed model (minimum BIC) was shown always to be the model with only the main effects, no interactions of the fixed effects and with random intercept with respect to time only (data not shown). All P values for these models are reported in Tables 1 and 2. 

Figure 6

 

Histograms of cone intensity on the selection of images shown in Figure 2. The y-axis represents the fraction of cones over the total number in the image, on the x-axis the values of cone intensity.

Table 1

 

P Values for Main Effects Condition (Between-Groups), Time (Within-Subjects), Retinal Location (Within-Subjects), and Random Intercept for the Data at All Retinal Locations Considered Separately

 

Table 2

 

P Values for Main Effects Condition (Between-Groups) and Time (Within-Subjects) and Random Intercept for the Average of the Data Between the Retinal Locations

 

Table 3 reports the estimates from the linear mixed model of the average of the parameters at the four retinal locations. For visual assessment of the differences between the groups, Figure 7 shows boxplots of the parameter estimates for the average of the locations at all time values. 

Table 3

 

Mean ± SD Values and 95% CI of the Average Between Locations Parameters as Estimated From the Linear Mixed Model for the Condition Effect Only

 

Figure 7

 

Boxplots of marginal data at all time levels for the results averaged between the four retinal locations. The plots show the median of the parameter estimates of the two study groups with the boxes representing the interquartile range, that is, the 25th and 75th percentile of the group values. Whiskers extend to the most extreme data values that are not outliers. Outliers (plus signs) are defined as values that are further than 1.5 times the interquartile range from the median.

Cone Intensity

From Table 1, we see that time and retinal location had no significant effect on the parameters related to cone intensity (P > 0.05), except for mean cone intensity (P = 0.038 for time) and cone/intercone intensity (P = 0.005 for location) considering all the values at different retinal locations. Retinal location was not a significant factor (apart from one parameter), justifying the use of the mean of the parameters at the four locations as a global estimate of the parafoveal cone reflectance properties for each subject. 

Considering the global estimates of the parameters (Table 2), time was not significant for all parameters (P > 0.05), while the condition had a significant effect for mean cone intensity (P = 0.047), cone/intercone intensity (P = 0.036), and the slope of the variograms (P = 0.002). 

The mean cone intensity was significantly higher in controls (1.0033 ± 0.0003) than in NPDR eyes (−0.0010 ± 0.0005 lower than controls, Table 3) if the four locations were averaged (P = 0.047), but not if the four locations are considered separately (P = 0.068). The other parameters describing the shape of the cone intensity distributions, standard deviation, and skewness, were not significantly different between the two study groups (P > 0.05, with separate locations and average of locations). The ratio of cone to intercone intensity was significantly higher in controls than NPDR eyes, both considering all locations (P = 0.031) or the average of locations (P = 0.036, NPDR being −0.004 ± 0.002 lower than controls 1.002 ± 0.001, Table 3). 

The linear slope of the variograms was significantly steeper in controls than in NPDR eyes, both considering the values at the four retinal locations separately (P = 0.002) and averaged (P = 0.002), with global estimates 0.039 ± 0.002 for controls and NPDR being −0.011 ± 0.003 lower than controls (Table 3). The average variogram curves for controls and NPDR eyes are shown in Figure 8. For most of the curves, the slope at the origin and shape of the curve were constant over time, confirming the lack of effect given by time (P = 0.331 for average of locations, Table 2). 

Figure 8

 

Marginal mean variogram curves for the study (red curve) and control (blue curve) groups with ±1 SD (dashed curves). The mean curves were calculated averaging all the curves for all the retinal locations of the subjects according to their group. The study group had a shallower slope than controls at the origin of the variogram curve (<20 μm), indicating a more pronounced correlation between cone intensities separated by short distances than the control group. Most of the NPDR cases showed variogram curves having a shape not leveling to a maximum value, but rather some showed a peak and then decreased again, while others showed two peaks (see Supplementary Material).

Most parameters also showed a significant contribution from the random deviation of the intercept for each individual around the respective group means (Wald Z test for random intercept; see the last column in Tables 1, 2). 

The subtraction of the inner retinal layer variations in intensity affected the shape of the histograms. In particular, after subtraction, the variograms showed a steeper slope at the origin and a leveled curve with the suppression of the peaks, which means that the low-frequency intensity variation was largely removed. Figure 9 shows two representative curves (one healthy and one NPDR), before and after the subtraction process. 

Figure 9

 

Variograms of cone intensities for two subjects (C2 and NPDR3), before and after the subtraction of the low frequency intensities as measured from the inner retina images. The subtraction process causes a flattening of the curves and a steepening of the slope at the origin. Where present, the peaks also are suppressed.

Cone Mosaic Texture

The mosaic texture metrics, sharpness, and entropy, were not influenced by retinal location or the acquisition time (P > 0.05). On the other hand, a significant difference due to the patient's condition was found; the difference between controls and NPDR was statistically significant only for entropy (P = 0.014 all locations, P = 0.14 average of locations), with NPDR being 0.16 ± 0.06 higher than controls at 5.54 ± 0.04 at the average of locations. Sharpness, on the other hand, was comparable between groups (P = 0.410 all locations, P = 0.408 average of locations) and had a significant contribution of the random intercept (P = 0.010 all locations, P = 0.007 average of locations). The random intercept term for entropy was not significant (see the Wald Z test for random intercept in Tables 1, 2). 

Discussion

The mean cone intensity was significantly different between healthy and NPDR eyes, but only if the average between the retinal locations was considered (Tables 1, 2). In addition, if all data at different locations were taken into account, there was a dependence on time, and this was the only parameter for which this was the case. This dependence on retinal location, together with the lack of significant differences in standard deviations and skewness between groups, indicated that the histograms of cone intensities are not good candidates to investigate differences in cone reflectance caused by diabetes mellitus. 

Under the assumption that the intensity of the pixels inside the segmentations is primarily directional light backscattered from the outer/inner segment (IS/OS) junction and pixel intensity outside the cone segmentations corresponds to nondirectional backscattered light,^54–56 in healthy subjects there was a significantly higher fraction of directional backscattered light than in those with NPDR, given by the higher ratio of cone/intercone intensity. A significant difference in cone intensity in areas where the cones still were detectable, although with abnormal packing density arrangement, has been found previously.^14–16,19 A possible explanation could be the alteration of the wave-guiding properties due to the pathology of the individual cones.¹⁶ The significance of the retinal location for this factor could suggest that this variation of cone intensity, even if important enough at all locations so that its average still is significantly different between groups, might not be constant across the central retina (e.g., caused by variable effect of intraretinal scattering, and so forth). Another possible explanation could be given by the rods, which are not resolved in these images but may be present sporadically in between the cones at these eccentricities. In this case, the difference in the ratio would indicate a change in the rod reflectance, or a combination of changes in rods and cones. Even if significant, the difference between the ratio in the two groups had a small numerical value (Table 3), as in the case of the mean intensity, which could make the use of this parameter for clinical applications hard to achieve. Further investigation of the ratio of cone/intercone intensity at different illumination angles could improve the understanding of how this phenomenon is correlated with cone function and if it eventually also is related to cone spatial properties, such as spatial density or packing arrangement. 

We chose to use a sample area of 200 μm side, which was bigger than what has been used in similar studies,^14,15 to have an area big enough to evaluate the reflectance property of several clusters of cones in all subjects. From the analysis of the spatial distribution of cone reflectance through the variograms, we were able to infer a significant difference between the two groups and this result was consistent between the average of the locations and for all locations considered. The NPDR curves showed a shallower slope at short distances (<20 μm), reflecting the fact that the cones appear to be clustered in bright and dark patches more than in healthy subjects. In addition, the average variogram curves in NPDR eyes showed peaks at long distances (>50 μm), which were not found in healthy subjects. These peaks corresponded to the distance at which there was a greater difference in the observed cone intensity. In NPDR eyes, it is likely that the position of the peaks reveals the size of the bright and dark intensity patches present in the images (see Supplementary Materials). The intersubject variability in the shapes would make challenging the definition of a characteristic “standard” variogram shape in an adult population at this moment. For this reason, we analyzed the initial slope of the variograms, which showed potential as a straightforward method to quantify the degree of spatial dependence of cone intensities over a short range. The clustering of cone reflectances did not show a significant dependence on the retinal location, implying that difference in these two aspects of cone reflectance between healthy and NPDR eyes might have different causes. This dependence can be caused by the physiology/pathophysiology of the photoreceptor layer itself or of other layers. Further work is needed to investigate suitable models to fit experimental data and develop a more complete procedure to compare curves from different subjects or retinal locations. In addition, more work could be done to understand how this approach can be extended to other retinal diseases. Strictly related to the multiple humped shape of the variogram found in NPDR eyes, clustering of the cone reflectivity has been qualitatively observed in DR^14,19 and inherited retinal diseases. For example, in albinism, clustering has been related to the melanin distribution in the retina;³⁰ in choroidermia, hyperreflective clusters of cones have been related to alteration in the RPE.³³ Degeneration of the RPE also has been indicated as a cause of disruption of cone reflectivity in age-related macular degeneration, possibly before the formation of drusen.³⁴ 

The inner retina has been shown to influence the low frequency pattern of cone intensity in AO flood illumination images of the cone mosaic (Fig. 9). The presence of subclinical or clinically visible abnormalities located in the inner layers of the retina (e.g., microaneurysms, microvascular abnormalities, subtle retinal edema) might further explain the greater clumping found in NPDR eyes than controls.¹⁶ A thickening of the inner retinal layers, which could not be considered as macular edema, has been found previously in diabetic eyes.^14,19 Further work is ensured to correlate the changes in the cone reflectance properties with the abnormalities in spectral domain-optical coherence tomography (SD-OCT) cross-section images of the same subject. Another possible explanation could be a high degree of alignment between neighboring cones, with clusters of cones pointing towards the pupil center appearing here brighter than clusters of cones pointing away from the center. This hypothesis could be investigated again, as in the case of nondirectional backscattered light, with observations at different illumination angles.^27,54 

In any case, even if the source of the modifications in cone reflectance was the result of shadowing artifacts and not of abnormalities in the cones, their quantification still can be used as an indicator of the disease. 

We also analyzed the texture of the parafoveal cone mosaic with two metrics originally used for the assessment of image quality. The reason behind the choice of only these two metrics, sharpness and entropy, from the cited study⁵⁰ is that these investigators found other quality metrics, such as variance, contrast, and kurtosis, to be highly correlated with sharpness, while entropy was not. In this way, we were able to explore different aspects on the images, such as the definition of bright features over the background (sharpness) and the textural properties of the image (entropy). Entropy was significantly higher in the AO images of the cone mosaic in NPDR eyes than controls, while sharpness was very similar between the two groups. A possible explanation could be that sharpness measures the definition of the cone apertures as bright dots over a darker background. In this sense, this characteristic seems not to depend on the condition of the retina. Entropy, on the other hand, relies on the histogram distribution (Fig. 6). The fact that entropy involves taking the logarithm of the histogram makes it more sensitive to extreme intensity values, and so it is able to detect a difference in images with apparently similar histograms (in a previous study, sharpness was highly correlated with the image variance, while entropy was not).⁵⁰ 

We did not observe the characteristics of cone reflectance or mosaic texture of the same areas in AO flood illuminated retinal images to change significantly with time. In a previous study,⁴⁰ we found that although the cones individually change reflectance over time in a random way, the distribution of the cone intensities maintains the same shape over time, and this result was confirmed here. This also was true for the retinas affected by NPDR, as we also observed no significant interaction between time and condition. A possible explanation for the lack of change, especially for the NPDR group, is that the time range considered was too short compared to the development time of the disease. Furthermore, the retinal location proved not to be a significant factor influencing the cone reflectance in AO flood retinal images, meaning that the cone mosaic has similar intensity and textural characteristics at different locations at the same distance from the fovea. This result confirmed the validity of using the average of the values of the four locations as one global parameter with potential to be translated to clinical studies. 

An important aspect of this work is the automation of the analysis. We used our previously described method⁴⁰ for the analysis of large retinal patches (≥2 × 2° or ≥0.61 × 0.61 mm), showing how performing the detection on time averaged images further improved the performance of the automated cone detection algorithm, which is 97% of cones correctly detected using only one image. The improvement in the automated detection performance using more images taken at different times also can be justified for the NPDR images, which can present a decrease in cone density but no major deterioration that would require manual supervision. On similar images, the investigators in a different study¹⁴ have used an automated algorithm with poorer performance⁵⁷ and the percentage of cones that had to be corrected manually was not greater than 9% also for NPDR cases, which is compatible with the performance of the same algorithm on healthy retinas.⁵⁷ 

Limitations of this study included the small sample size (especially the NPDR cases), though it had enough power to verify the hypothesis of the study, and the incomplete number of observations with time in all cases. Future work could include a more consistent number of subjects at different stages of the disease and a longer time range, which could lead to observation of change on the same retina with time. Another improvement could be achieved using complementary observation modalities, such as AO-scanning laser ophthalmoscopy (AO-SLO) and/or AO-OCT, as well as using different illumination angles, which could provide a better resolution closer to the fovea (<2°) and insights into the causes of differences in cone reflectance. Finally, analysis of the function of the retina in the selected regions also could be performed, to determine the clinical significance of the differences between healthy and NPDR eyes. 

In conclusion, we observed significant differences in cone mosaic reflectance properties between healthy eyes and eyes with mild NPDR, in its spatial organization and in its intensity, especially between directional and nondirectional backscattering. We did not observe significant changes of the parameters with time in any group. We performed a largely automated analysis of cone reflectance and introduced a novel method for the study of the spatial distribution of intensity, the variogram, which was able to quantify differences of the spatial dependence of intensity values of the cone layer at a short range between NPDR and control eyes and a tendency of cones in NPDR to appear clustered in clumps of similar intensities. 

Acknowledgments

Supported by the National Framework Program for Research and Innovation PON (grant no. 01_00110), the Italian Ministry of Health (5x1000 funding), Fondazione Roma, and the Irish Research Council (GOIPG/2013/775). 

Disclosure: L. Mariotti, None; N. Devaney, None; G. Lombardo, None; M. Lombardo, None