Light reflectivity and interference in cone photoreceptors

: In several modes of retinal imaging, the primary means of visualizing cone photoreceptors is from reﬂected light. Understanding how such images are formed, particularly when adaptive optics techniques are used, will help to guide their interpretation. Toward this end, we used ﬁnite diﬀerence beam propagation to model reﬂections from cone photoreceptors. We investigated the formation of cone images in adaptive optics scanning laser ophthalmoscopy (AOSLO) and optical coherence tomography (AOOCT). Three cone models were tested, one made up of three segments of varying refractive index, the other two having additional boundaries at the inner/outer segment junction and outer segment tip. Images formed by the ﬁrst model did not correspond to AOOCT observations in the literature, while the latter two did. The predicted distributions of reﬂected light intensity from the latter cone models were compared to the distribution from AOSLO images, both studied with light sources of varied coherence length. The cone model with the most reﬂections at the inner/outer segment junction best ﬁt the data measured in vivo . These results show that variance in cone reﬂection can originate from light interfering from reﬂectors much more closely spaced than the outer segment length. We also show that subtracting images taken with diﬀerent coherence length sources highlights these changes in interference. Diﬀerential coherence images of cones occasionally revealed an annular reﬂection proﬁle, which modeling showed to be very sensitive to cone size and the gaps bracketing the outer segment, suggesting that such imaging may be useful for probing photoreceptor morphology.


Introduction
Adaptive optics (AO) is an important technique that enables the imaging of single photoreceptor cells in the living retina [1,2]. The enhancement of scanning laser ophthalmology (SLO) and optical coherence tomography (OCT) through AO provides cellular and subcellular detail in the study of healthy and diseased retina [3][4][5][6][7]. To better interpret what is seen in such images, it is important to know the sources of the reflected light. Having working models of light interaction with photoreceptors opens up avenues for exploration about basic retinal function and the causes of abnormalities. For example, it has been observed in AOSLO images that the reflectance of cone photoreceptors in the mosaic is variable, not only between cones, but also for individual cones over time [8][9][10] and in response to visual stimuli [11][12][13]. In AOOCT imaging it has been observed that the width of the bands corresponding to the junction between the inner and outer segment (IS/OS) and at the outer segment tips (OST) are greater than that expected from a single reflector [14][15][16], with debate about the origins of the IS/OS [17][18][19][20][21][22][23] and OST [24,25] bands. For a comprehensive account of all these phenomena, continued exploration of explanatory models is needed.
We have previously shown that a model of light propagation in the forward direction can make useful predictions about the amount of light absorbed [26]. This model used finite difference beam propagation [27][28][29] to predict how light travels down a cone. It is based on the well-established waveguiding features of cone receptors that lead to such phenomena as the Stiles-Crawford effect [30][31][32][33][34][35][36][37]. The model incorporated light polarization, unbound as well as bound light, and multiple photoreceptors with arbitrary sizes, but neglected back reflection. The modeled cone we used had three compartments (myoid, ellipsoid and outer segment) with different refractive indices.
Here we expand the technique to include reflections from refractive index discontinuities in the axial direction. We investigate two further cone models that include additional refractive index changes at the IS/OS and OST junctions. The amount of reflected light collected at a detector is calculated and the effects of interference are included.
Interference between reflecting layers in the cones has been suggested as a possible explanation for the time variance and the intrinsic response of cone intensity in AOSLO images [9,11,38,39], and it has been qualitatively shown that the reflectivity of cones is related to the coherence length of the light source [40]. In this paper we present a quantitative test of 3 cone models against AOSLO imaging data where the coherence length of an AOSLO's light source was changed. Qualitative assessments are also presented, including a differential coherence AOSLO image where only interference effects go into its formation.

Finite difference beam propagation method
The finite difference beam propagation (FDBP) method has been described in full previously [26]. Here we give only an overview. The technique is a method of solving the Helmholtz wave equation with a slowly varying envelope for given initial and boundary conditions, i.e. the initial electric field and the variation of the refractive index of the cone. The technique is 3-dimensional and uses the small angle approximation. The method is semi-vectorial, meaning that polarization is included, but cross coupling between polarizations is neglected. The method iteratively propagates a wave down the length of a photoreceptor (the z direction) by calculating what the field does in one perpendicular x-y plane based upon the previous plane using the alternating direction implicit scheme. As an example, for x-polarized light, the Helmholtz equation is where propagation is in the z direction, i = √ -1, n 0 is a representative reference refractive index of the system, k 0 is the wavenumber of the light in a vacuum and ψ x is the x-polarized electric field. In the alternating direction implicit scheme, the 3D problem is split into two orthogonal 2D problems and the operator, H, is split into two: where the A operators are defined by and Equation (2) may be discretized, using the Crank-Nicolson scheme, to relate the field at the (m + 1) th step to the field at the m th : replacing ψ x with u for clarity. The superscripts denote the step in the z direction. After separation into two substeps, Eq. (5) can be represented by systems of linear equations. For the first substep (in matrix form) M y is a tridiagonal matrix, with the coefficients and q y is a one dimensional matrix with the coefficients Each equation in the system relates a voxel of the field at u m to one at u m+½ at the lateral position i,j, which is influenced by the adjacent voxels in the y direction of u m and the refractive index. The whole system relates one row of voxels of u m to u m+½ .
The second substep is also written in matrix form; and has the coefficients and In a similar fashion to Eq. (6), each equation of this system relates u m+½ to u m+1 at the lateral position i,j, but now the influence of adjacent voxels of u m+½ are in the x-direction. T and R cause polarization dependent differences at lateral refractive index boundaries between the x and y directions. A similar result can be shown for light polarized in the y direction (see Ref. [26] or [41] for complete details).
Equations (6) and (9) can be readily solved using a number of algorithms. Our implementation was developed in Matlab and uses its backslash function to solve these and the equivalents for the y-polarization. An absorbing boundary condition is used to suppress reflections from the edge of the computing space. In this implementation only light of 842 nm was considered, as it is one of the center wavelengths used in our AOSLO imaging (see Section 2.6). At this wavelength, absorption in the outer segment is present, but insignificant enough to be neglected [42].
When light encounters a refractive index change in the axial direction, the reflection is calculated. As the calculated envelope of the field is complex, its phase is known, so the wavefront at any depth is also known. The slope of the wavefront gives the local direction of the light. This means that at a reflection boundary at any x,y coordinate, the field magnitude, polarization, and direction of the light is known, allowing the Fresnel equations to be used to calculate the reflected field, based on the refractive index difference. The reflected light is back propagated up the photoreceptor towards the external limiting membrane (ELM) which forms the entry plane of the model. The intensity at the photodetector pinhole (found in both SLO and scanning OCT systems) can be computed via the optical transfer function, so is given by the square of the convolution of the pupil function and the reflected field at z = 0. Integrating this over a pinhole gives the power that will be incident on a detector: where ψ r is the reflected field, h is the ideal amplitude pupil function, r p is the radius of the pinhole (r 2 = x 2 + y 2 ), and * denotes the convolution operation. The unreflected light is propagated forward and is calculated to keep the sum of the fields of the outputs equal to the input, so that where the subscripts −z, +z and z denote the backward reflected light, the forward propagated light, and the light incident on the boundary, respectively. The unreflected light is propagated to the next refractive index boundary, where the process is repeated. Reflections are not calculated for the back propagation, as these are at most of the order of 10 −4 in intensity, calculated using the Fresnel equations for normally incident light (see Section 2.3 for the refractive indices and boundaries). This is a negligible portion of the total reflected light. All propagations are computed with a resolution of ∆x = ∆y = 0.1 µm and ∆z = 0.2 µm.

Initial electric field
The initial electric field, ψ 0 , is calculated by Fourier transform of the field at the eye's pupil. While an ideal system would have a uniform amplitude across the beam, in practice this is difficult to achieve. Typical AO systems would have light launched into the system via an optical fiber, so have a Gaussian distribution, truncated by any optical stops, such as edges of optical elements or the eye's pupil. For calculation of ψ 0 , the eye is assumed to be 22.2 mm long, from pupil to ELM (z = 0), with a refractive index of 1.33. Values for the width of the Gaussian at the pupil and the diameter of the truncation were measured from the system used in experiments. The Gaussian had a 1/e width of 4.75 mm (1/e 2 in the intensity), truncated at 5.85 mm. For normalization, the power of the light, |ψ 0 | 2 , at z = 0 is set to 1.

Cone models
In our previous study, the cone photoreceptor was modeled as three parts with differing refractive indices. Here we call this Model 1, with a myoid (n m = 1.3605) and ellipsoid (n e = 1.394) that make up the inner segment, and an outer segment (n os = 1.39) [23]. The shape of the cone is determined by where d is the cone diameter, d os is the diameter of the outer segment and d 0 is the diameter of the inner segment at z = 0, defined at the ELM. The length of the inner segment, l is , is equal to the sum of the length of the myoid (l m = 13 µm) and ellipsoid (l e = 16 µm). The length of the outer segment (l os ) was 35 µm. All lengths were as used previously [26]. The cone is surrounded by a medium with a refractive index of n ex = 1.34 [31,32,43,44]. To obtain a reflection from the ELM, the refractive index of the first z slice of the cone was set to 1.3665, giving a boundary between the 1 st and 2 nd voxels. This value gives a reflection of an approximate magnitude seen in OCT images [45]. The representative refractive index (n 0 ) was set to 1.378. For d 0 at a particular eccentricity (ε), cone density (ρ) was calculated by based on the measurements of Song et al. [46], then converted to cone diameter using from a fit obtained by Liu et al. [47]. We used a value of 1.5 µm for d os . Model 2 has the same three segments, but is modified in two ways. The first is that the refractive index boundary between the myoid and ellipsoid is graded over a distance of 2 µm. In Model 1 the strongest reflection would be expected from the myoid-ellipsoid boundary at a depth of ∼12 µm. A reflection here is not observed in humans or other primates, but has been observed in ground squirrels [48], where it is weak compared to the IS/OS reflection. So a strong reflection here is not likely to fit previous observations well. The mitochondria in the ellipsoid that are responsible for its higher refractive index can be observed in electron microscopy images to be more disordered at the proximal end of the ellipsoid than the distal, having staggered positions (see, for example, Fig. 2 of [49], Fig. 19 of [50], or Fig. 5 of [51]). The graded refractive index is an approximation of subcellular anatomy that avoids the need to run at considerably higher enough resolution to capture mitochondrial details [23]. Phase contrast microscopy using refractive index differences to produce contrast in images also suggests that a graded refractive index here is justifiable [52]. A second change was made based on measurements of the IS/OS and OST bands in AOOCT images. It has been found that the widths of these bands are thicker than expected from a single reflector; their widths are greater than the coherence length of the source used [14][15][16]. An explanation for this was that there may be multiple discrete reflectors in these regions; by adding extracellular gaps to the cone model, this theory is tested here.
Model 3 is an extension of Model 2 that includes an additional layer of intracellular matter at the distal end of the inner segment, representing a space where mitochondria do not come into contact with the cell wall at the IS/OS (see, for example, Fig. 2 [54]). This layer is 0.2 µm thick (our ∆z) and has a refractive index equal to that of the myoid. While Model 2 tests whether or not multiple reflections are present at the IS/OS and OST, it only has 2 boundaries at each region. It is possible that there are more reflectors present, and Model 3 is a first test of increasing that complexity.

AOOCT simulation
We use a simplified model of OCT to simulate image generation. The powers at the detector from the model were used to create a depth reflectivity profile. As the amplitude of an OCT signal is proportional to twice the square root of the reflectance, this is square rooted, doubled, then convolved with a Gaussian with a full width at half maximum (FWHM) given by the assumed coherence length of the source. Widths and amplitudes of the bands are measured by fitting a three-component mixed Gaussian model to the result. This mimics a method used elsewhere to measure widths in actual AOOCT images [15,16]. When the extracellular gap sizes are large, this method fails as it attempts to fit one Gaussian when there are clearly more than one present, so we limited the model to gap sizes of less than 4.8 µm. The scanning of the beam is mimicked by moving the position of the initial field across the cone at the ELM.

AOSLO simulation
OCT images are formed by depth-dependent reflections interfering with a reference beam, there is no interaction between the reflections themselves in the signal. In AOSLO imaging the reflections are not separated in this way and interference between them would occur. The observed power, P, is given by the sum of all N reflections, plus the sum of all their interference terms: where γ is the degree of coherence and is a function of the optical path difference between two reflections, L (which is the geometric path length multiplied by the refractive index, n). This is related to the spectral shape of the source; for a Gaussian shape, it is given by where l c is the coherence length of the source. The optical path length is twice the distance between two reflectors, multiplied by the refractive index. As the distance between each reflectance pair is increased there is a modulation of the signal (first exponential), but this is contained within a Gaussian envelope (second exponential). In relation to z, the modulation has a period of λ/2n and the envelope has a FWHM of l c /2n, where the divisions by two come from L being twice the path length between two reflections in the z dimension due to the double pass. As a consequence, the reflectivity of a photoreceptor in AOSLO imaging is dependent on the coherence length of a source. We used Eq. (18) to compute the intensity of reflected light from cones, calculated from the propagation of light from each reflecting layer. As with the OCT simulation, the scanning of the beam is mimicked by moving the position of the initial field across the cone at the ELM. Cone intensity is the mean power at the detector of a 3 × 3 square centered on the cone. To simulate cone intensity variation over a patch of retina, we created distributions of cones with varying extracellular gap sizes. Curves of the bandwidths in the simulated OCT images versus gap size were made. These were then used to convert inverse Gaussian distributions of bandwidth (with means and standard deviations equal to those measured by Jonnal et al. [16]) to distributions of gap size at the IS/OS and OST. OCT widths were obtained for gap sizes with a spacing of 0.2 µm (∆z), but a 5 th order polynomial fit was used for determining the gap size distribution. A Monte Carlo simulation was then run to pick cones from these and the brightness of each cone calculated, using interpolated values based on the individual reflections. We ran the simulation for Models 2 and 3 (Model 1 has no extracellular gaps). The simulation was run for the two sources used in our experiments (see Section 2.6); one with low coherence (λ = 842 ± 32.5 nm, l c = 4.81 µm) and one with high coherence (λ = 840 ± 10.5 nm, l c = 14.7 µm). Coherence lengths are their values in air and were calculated by assuming a Gaussian spectral shape. The model was only run at the center wavelength, as it would take a prohibitively long time otherwise.

Comparison to in vivo images
The AOSLO instrument used for comparison has been described previously [55,56], so is only briefly outlined here. For imaging, light was channeled from a supercontinuum source (NKT Photonics, Denmark) and bandpass filtered to produce the high and low coherence sources. This light was unpolarized. A fixed bandpass filter of 842 ± 32.5 nm (for low coherence) was kept in place where the laser light was launched into the system. To image with the high coherence source, a narrower bandpass filter of 840 ± 10.5 nm was temporarily placed later in the optical path, allowing interleaved imaging with low and high coherence light. These filters were from Semrock (IDEX Health & Science, Rochester, NY). The power at the eye using the narrow filter was 160 µW. To ensure equal power for both conditions, the light intensity when using the broad filter was adjusted using an acousto-optic modulator in conjunction with neutral density filters. Retinal images of 512 × 512 pixels were acquired at 30 Hz using a raster scan. For all images the field of view was 1.2°. Wavefront aberrations were corrected by using a Shack-Hartmann wavefront sensor and a MEMS deformable mirror (Boston Micromachines, Cambridge, MA) working in a closed loop. Detection was via a photomultiplier tube (Hamamatsu, Japan), whose gain was kept constant between conditions.
Two human subjects (one 29 year-old female and one 56 year-old male) with normal vision were imaged in two locations, each ∼2°from the fovea using the low and high coherence length sources. The AOSLO system has real time image stabilization that relies on a reference image to correct for eye movement [57]. The same reference image (acquired using the low coherence source) was used for both conditions, ensuring good registration across all images. Blood vessel maps were also taken of the overlying vasculature [58] and used to create a mask to exclude any underlying cones from subsequent analysis. To quantify the light reflected from cones, each cone in the images was selected from the averaged image of the two conditions, using the brightest pixel as the reference for position. Positions where then manually adjusted on low and high coherence images separately. Cone intensity was measured in the same way as the model, using a 3 × 3 pixel square centered about the selected position.
To analyze how cone intensity changes between the sources, the Weber contrast is calculated as C = (I high -I low )/I high for each cone (actual or modeled). Because the models may not account for all variation in cone intensity, additional noise was added. Sources for variation in a localized patch of retina that are not included in these methods include variation in scattering from structure above the cones, in cone shape (e.g. cross sectional ellipticity) [59], in inner and outer segment size [26] and in photoreceptor tilt [60,61]. This added noise was applied as a random factor increasing or decreasing the cone intensity for individual cones. The factors are picked from an inverse Gaussian distribution with a median set to 1. The width of this distribution was determined by optimizing the 2-D correlation coefficient of the modeled and actual distributions (one dimension was cone intensity in the high coherence condition and the other dimension was cone contrast).

Simulation of AOSLO images
Simulated AOSLO images were created from Models 2 and 3 to enable comparison to in vivo images. The positions of cones in a simulated mosaic were generated by using a random close packing function [62], to introduce some jitter in the position of equally sized cones. The distribution that feeds the packing distances was an inverse Gaussian with a median equal to the cone spacing at the ∼2°eccentricity being modeled (see Section 2.3). Once the positions were set, the cones were given extracellular gap sizes based on the OCT width results of Model 2, described above. The 2D reflection profile of each cone was simulated within a 15 × 15 pixel array. When these arrays were placed on the mosaic, we used interpolation to handle any mismatch between the arrays and the mosaic pixel positions. Where these arrays overlapped, the values were summed together. The same cone positions and gap sizes were used in Models 2 and 3 for ease of comparison.

Differential coherence imaging
For the visualization of changes between images obtained with different coherence lengths (actual and simulated) we subtracted the low coherence image from the high to yield a differential coherence image. Although care was taken to equalize the gain and power at the eye for each condition during AOSLO imaging, it was still necessary to normalize these images to their medians before making the subtraction. Also before subtraction, the averaged cone images were registered using BUnwarpJ [63], a plug-in for ImageJ [64]. Although the image stabilization produces registered images, the correction is not perfect, having a frame to frame error on the order of a few pixels [55]; this secondary registration ensured this error was minimal. In the ideal case, the differential coherence image only shows the differences between the interference terms. In mathematical terms, the value of a pixel in such an image will be proportional to the difference in powers at the detector: where P 1 and P 2 are the powers at the detector for two sources with differing degrees of coherence (γ). Substituting in Eq. (18), only the interference terms remain:

Propagation of reflected light through modeled cone photoreceptors
To appreciate the complexity of light reflected from cones, we begin by illustrating the propagation patterns from each of the reflective surfaces in each of the three cone models, along with the simulated OCT images the models produce (Fig. 1). To ease comparisons, all cones were modeled for a 2°eccentricity site, having the same dimensions (d 0 = 6.5 µm, d os = 1.5 µm) and extracellular gaps (2.4 µm). The wavelength used was 842 nm. The initial electric field was offset laterally to highlight the waveguiding behavior of the cones. Despite the incident beam being off-center, waveguiding in the outer segment yielded a reflection pattern that was centered in the cone by the time it reached the ELM. The FWHM axial resolution for the OCT simulations was 4.2 µm, with the center-to-center spacing of the cones being 7.2 µm. There is little difference between the models in the propagation of light on its inward journey (first column of Fig. 1). The light also does not appear to be significantly affected by losses due to the reflections. This validates the assumptions in our previous study, where reflections were not included [26], and suggests that the simplest model is sufficient to study forward propagated light.
In Model 1, reflections are produced at 4 boundaries: (1) from the ELM (propagation not shown), (2) from the myoid-ellipsoid junction, (3) from the IS/OS, and (4) from the OST. The brightest reflection originates from the myoid-ellipsoid junction, and the dimmest at the IS/OS, a pattern that does not agree with typical OCT images. Model 2 has three OCT bands and matches OCT images better: the reflection from the ELM is the weakest, with the reflection from the IS/OS being the brightest. The OCT pattern for Model 3 is almost identical to Model 2, arising from more boundaries but with smaller refractive index steps giving smaller reflectivities. In all models, light in the outer segment is well guided, exhibiting essentially single mode-like propagation, with most light seemingly contained in the LP 01 mode.
Because we were aiming to estimate the range of gap sizes present in a population of cones, we first used modeled OCT signals to see how gap sizes would affect OCT band widths (Fig. 2(A)). Up to gap sizes of about 3 µm, OCT bands at the IS/OS and OST have similar widths, diverging thereafter with the IS/OS increasing in size less rapidly than the OST band. The bands at the IS/OS and OST are broader than the FWHM resolution of 4.2 µm, because signals from multiple reflectors are being combined. From these curves it is then possible to use published OCT width statistics (Fig. 2(B), statistics from [15,16], see Methods) to obtain gap size distributions for each of the models (Fig. 2(C)). These gap size distributions are used to examine interference behavior in Section 3.2.   [15,16]) were used to derive gap size distributions (C) from the functions in A, based on 100,000 draws from the distributions in B. Inner segment diameter, d 0 , was 6.36 µm, corresponding to 2°eccentricity.
Variation in gap size also leads to variation in the reflected light from cones in modeled AOSLO images. Because the gap sizes are on the order of a few wavelengths of light used for imaging, both constructive and destructive interference is predicted (Fig. 3). The interference-related modulation in cone intensity extends much further for the high coherence source than the low in both models. From these plots, it is clear that the intensity of cone photoreceptors is dependent on the gap size and the coherence length of the source. In both models, changing the coherence length has little effect on mean intensity (Model 2: low coherence, 3.04×10 −4 , high coherence, 3.04×10 −4 ; Model 3: low coherence, 1.91×10 −4 , high coherence, 1.97×10 −4 ). Model 2 leads to higher intensity variation than Model 3 for both sources (low coherence 0.27 vs. 0.18, high coherence 0.38 vs. 0.28, normalized RMS deviation).

Interference and differential coherence effects
The dependency of the intensity of reflected light from cones on the coherence length is also seen in in vivo images. We made AOSLO images from the same 9 arcmin patch of retina using both high and low coherence sources (Figs. 4(A)-4(B)). The high coherence image has more intensity variation than the low coherence image (the variance of the normalized pixel values were 0.020 versus 0.015; the pixel values were not normally distributed). Subtracting the two images to create a differential coherence image isolates the effect of interference (Fig. 4(C)). Here it is apparent that the direction and magnitude of the intensity differential is quite variable, suggesting that gap size variations on the order of nanometers exist between individual cones. We highlight 4 cones that have large changes, but others have little or no change at all, as expected from the interference data in Fig. 3.   Fig. 4. AOSLO images of the cone photoreceptor mosaic at 1.7°inferior retinal eccentricity taken with (A) the high coherence source, and (B) the low coherence source. (C) Normalized differential coherence image (B minus A) showing the direction and magnitude of intensity changes due to interference. Four example cones are highlighted: red = cones with negative differential intensity, blue = cones with positive differential intensity.
To quantify the effect of coherence length on imaging, we examined larger areas of retina, which entailed avoiding cones that were situated beneath blood vessels, as light path distortion from overlying vessels was not built into the model. We generated a blood vessel mask for each cone mosaic site (one example shown in Fig. 5), allowing us to identify un-occluded cones for further analysis. From the 4 sites imaged, we measured the reflected light from 2,285 cones under high and low coherence conditions.
To estimate the predicted change in reflected light caused by altering the coherence length of the sources we began with a Monte Carlo simulation, drawing from the gap size distributions in Fig. 2(C) (see Methods for details). For each cone we calculated the reflectance change between high and low coherence images as a Weber contrast value C (Figs. 6(A)-6(B)). These 2-dimensional histograms of cone intensity versus contrast show the probability of a cone intensity change between the two coherence conditions, for a given intensity in the high coherence condition. Both Model 2 and 3 predict that cone intensity changes generally fall along a line with a slope close to 1. This indicates that a cone's intensity in either model is likely to move closer to the mean cone intensity (which is equal for both coherence conditions) when imaged with low coherence light, no matter what its brightness is in the high coherence condition. The predictions from Model 2 show more variance than those of Model 3, which is expected from Fig. 3 where the range of modulation from interference is shown to be higher in Model 2.
To adequately compare the differential coherence contrast data between the models and the AOSLO images, a noise factor needed to be added to the models to account for factors not present in the models (see Methods for details). The noise factor was drawn from an inverse Gaussian distribution and optimized to fit the measured AOSLO cone contrast data (Fig. 6(E)). The noise added to the models spread the data out in the direction of the intensity axis (Figs. 6(C)-6(D)). Because the noise is a factor affecting each reflection equally, the spread is proportional to the original cone intensity, but the added noise does not influence the distribution of contrast. When the two dimensions of the plots are considered, it can be seen that Model 2's distribution peaks at low cone intensity and high negative contrast, whereas Model 3 peaks closer to the mean intensity and zero contrast. Model 3 has the better correlation coefficient with the in vivo AOSLO data (R = 0.96; Model 2, R = 0.78), indicating that it has better predictive abilities with respect to changes in reflectance under differing coherence conditions. The distributions of cone reflected light intensity arising from different coherence sources reveals the main difference between Models 2 and 3 (Fig. 7). As these distributions are normalized to their mean values, a differing median is indicative of two distributions being different. Model 2 in the high coherence condition yielded the only distribution that differed from the measured data (Model 2: high coherence, p < 0.01; low coherence, p = 0.94; Model 3: high coherence, p = 0.52; low coherence, p = 0.38; Wilcoxon rank sum test). For low coherence light, the variance of both models is driven largely by the fitted noise (these distributions are much narrower than the measured data without this added noise). As only Model 3 fits both coherence distributions well, it makes the best predictions of the contrast histograms in Fig. 6. These results suggest that the IS/OS junction is more adequately described as a site with more than just two reflective surfaces.
Similar results can also been seen qualitatively when comparing simulated AOSLO to actual AOSLO images (Fig. 8). In the high coherence condition, the Model 2 image shows greater variance than both Model 3 and the actual AOSLO, which appear to be similar. Variance in the low coherence images is lower for both models, similar to the lower variance in the AOSLO image. There are excess dark regions in the AOSLO images, due to the shadowing of blood vessels. The differential coherence images from the models also show the range of positive and negative interference effects as seen in Fig. 4(C), confirming the idea that interference plays a role in the reflected light pattern from cones.  N = 100,000). Incorporating a generic noise factor in the reflectance intensity (C, D) spreads the modeled data along this axis, allowing for unaccounted sources of intensity variation. (E) Cone intensity vs. contrast distribution measured in AOSLO images(N = 2,285 cones) best matches the data from Model 3. The noise factor added to both models was optimized to produce the highest correlation coefficient with the measured data. The histogram key shows examples of how a cone's reflectance would change within each quadrant of the plot. For example, when a cone is located in the bottom left quadrant (I high < 1, C < 0), its intensity lightens, a cone in the top right quadrant (I high > 1, C > 0) darkens. Fig. 7. Comparison of modeled vs. measured probability distributions of cone reflection intensity. (A) Normalized cone intensity data in Model 2 matches the measured data in the low coherence condition (red) but not the high coherence condition (black). (B). Normalized cone intensities data from Model 3 matches both coherence sources. Both models incorporate the optimal noise factor derived for Fig. 6.  Fig. 8. Comparison of simulated and actual AOSLO images. Images with the high coherence source have consistently higher variation in cone intensities than images produced with a low coherence source. Cone positions and gap sizes are identical in all modeled images. Images are normalized to the high coherence condition. Actual images were taken at an inferior eccentricity of 1.7°. Differential coherence images (low minus high, right column) exhibit variation in reflection due to interference as well as occasional annular intensity profiles centered on cones. Differential data are represented in root mean square (RMS) units. Bottom panels show magnified views of 4 matched examples from the modeled and in vivo data (outlined in yellow and red), illustrating annuli with two types of polarity between center and annulus.

Annular patterns from differential interference in cones
The models revealed a new interference related phenomenon that also seems to be present in in vivo AOSLO images when studied closely. For a subset of cones, the differential coherence image appears as an annular reflection, consisting of a dark spot surrounded by a lighter ring, or a light spot surrounded by a darker ring. These annuli were more prevalent in Model 3 than Model 2. Annuli were found in AOSLO images in both subjects at both locations, but not in as great numbers as seen in the modeled image. Higher magnification views are shown in Fig. 8, along with examples from the models at the same scale, which have comparable sizes. Because cone morphology varies slightly on a local scale, we examined how this variation could influence the reflection patterns observed.
To show that the visibility of the annulus in the modeled images is dependent on the gap sizes, we varied the IS/OS and OST gaps of one cone systematically (Fig. 9). The difference images show that only changes of fractions of a wavelength in the gap sizes are needed to obliterate the annulus, as would be expected from an interference effect.
We next investigated the effect of cone diameters. We ran Model 3 for cones with fixed extracellular gap sizes, but varied the inner segment and outer segment diameters (d 0 and d os ). Varying these cone dimensions changes the reflection intensity in both high and low coherence sources (Fig. 10). Intensity increases monotonically with outer segment diameter, but peaks at a diameter of 5 µm for the inner segment. In the differential coherence image, it is clear that presence of an annulus is also size dependent, with larger diameters producing annuli ( Fig. 10(C)). For these propagations, the light was unpolarized. Fig. 9. Cone reflection intensity as a function of extracellular gap size. A cone that exhibited an annulus with Model 3 in Fig. 8 was chosen as a starting condition (IS/OS gap = 1.28 µm, OST gap = 2.27 µm). The gap sizes were varied independently in steps of 0.05 µm. Variation of cone intensity is more obvious in the high coherence condition, and is greatest when the IS/OS gap is varied. The visibility of the annulus in the differential coherence images is dependent on both the IS/OS and OST gap sizes. Cone images are normalized to the maximum value among coherence conditions; differential images are normalized to the maximum deviation.
When polarized light is modeled, the reflection profiles are elongated in the x and y directions (Fig. 11(A)) and can become bifurcated with larger cone dimensions. Polarization effects are not seen in the AOSLO images because the source was unpolarized. This behavior is reminiscent of multimodal propagation, a combination of the LP 01 and LP 11 modes [32], but likely deviates from them because radiation modes are included. We dissected the origin of the annuli by examining the profile of the power returned from each layer for two outer segment diameters. When the outer segment diameter is small, all reflections had approximately Gaussian profiles with similar widths, but varied in amplitude (Fig. 11(B)). As the outer segment diameter is increased, all reflected intensities increase, but those from the IS/OS have a wider profile and can exhibit multiple peaks. In addition, the intensity of the OST reflections increases more than those of the IS/OS, so its influence on the result is increased. Thus, annuli originate from two factors; (1) the differences in reflected intensity that are gap size dependent, where the interference pattern will change when the coherence length is changed, and (2) the relative width of those reflections, dependent on the cone dimensions. To illustrate these two factors, we depict the reflection patterns from the IS/OS and OST separately, for two outer segment diameters. The smaller diameter does not produce an annular differential coherence image (Fig. 11(C)), while the larger one does because the intensity patterns from the two gaps are not only of opposite magnitude but are also of different size. Since gap sizes and cone dimensions are likely to vary in vivo, this will lead to different types of differential coherence images for each cone-whether bright or dark, circular or annular-as discerned in the AOSLO images (Fig. 8). Example cone reflectance from a high coherence (A) and low coherence source (B) for a cone with an IS/OS gap = 1.40 µm and OST gap = 1.17 µm. Differential coherence images (C) show that the presence of an annulus is sensitive to segment diameters. Varying the gap sizes will yield different reflectance profiles, including cases with negative difference images (see Fig. 9).

Discussion
We extended the use of the FDBP method to investigate the reflectivity of cone photoreceptor models and the potential effects of interference within cones at the IS/OS and OST junctions. Two of the models predicted a sensitivity of cone reflection intensity to the coherence length of the light source, which we confirmed through in vivo AOSLO imaging.
Our first result concerned how each of the three models behaved with respect to generating OCT signals. Model 1, made of 3 conjoined segments, did not make good predictions. It had the strongest reflection in the middle of the inner segment, where none is observed in human AOOCT imaging [15,16,45,65,66]. Consequently, it was not investigated further. Models 2 and 3 incorporated extracellular gaps and made good predictions of AOOCT imaging, specifically the broadened reflections from the IS/OS and OST that would be expected from multiple closely-apposed reflectors [14][15][16]. When applied to AOSLO imaging, Models 2 and 3 both showed that a variation in the extracellular gap sizes causes a variation in the cone reflection intensity. Model 3 in particular provided a better fit to the in vivo data, showing the strongest relationship between coherence length and intensity (Figs. 6 and 7). This suggests that the light returned from cones does not originate from simple single reflectors or thick media scatterers, as neither of these factors is subject to interference effects. In the case of single reflectors (as in Model 1), the distances are too great between IS/OS and OST for interference to occur with the sources used here, which are typical of AOSLO systems. With thick media scatterers, any constructive interference that occurs between light from different depths is likely to be cancelled out by destructive interference from other depths, assuming the scatterers are equally distributed. It is unclear how many reflectors are involved, however increasing the number of reflectors at the IS/OS or at the OST in Model 3 is likely to only modestly improve the fit to the in vivo imaging data, but this has not been tested here.
The only difference between Models 2 and 3 is the additional intracellular layer at the distal end of the ellipsoid, which is based on anatomical observations [51,53,54]. This reduces the effects of interference at the IS/OS. Maximum interference occurs when the path length between two reflectors is equal to Mλ/2, where M is an integer. Minima occur when that distance is λ(M/4 + 1/2). With only two reflections and a varied gap size, Model 2 can achieve either limit. Such limits are not achievable by Model 3 because the additional gap was of fixed dimension and did not match either interference condition. Thus the interference effect in Model 3 will always be reduced. The result is less variation in the cone intensity, giving a better match to the empirical measurements (Figs. 6 and 7). A fuller elaboration would be to vary Model 3's intracellular gap, but this would require propagations to run at higher resolutions. Such a model would still produce less variation than Model 2, but there may be more cones displaying extreme hypo-or hyper-reflectivity.
Interference effects do not explain all the differences in reflected cone intensity that are observed in AOSLO images (e.g. compare Figs. 6(B) and 6(D)). For the models to be fit to the data, additional intensity variation was required. We assume that this variation accounts for retinal factors that were not adjusted in the model that may affect cone intensity; namely, (1) scattering from overlying retinal tissue, (2) cone cross-sectional ellipticity [59], (3) segment lengths [49,67], (4) cone direction [60,61], and (5) variation in inner and outer segment diameters [49]. In fact, Fig. 10 shows that cone intensity is quite sensitive to the cone segment diameter. If these variations were explored rigorously in future studies, there may be less need to add the inverse Gaussian noise that we utilized to fit the data. However, there is almost no quantitative data with sufficient resolution for these five factors that could enable more comprehensive modeling.
One earlier study examined the variation of cone reflection intensity as a function of the coherence length of the source in AOSLO images [40]. Their model was like our Model 1, having two primary reflectors at the IS/OS and OST, and relied on interference between them. When modeled with a coherent source, their cone images revealed axial interference, as well as lateral interference for the most closely-packed cones in the fovea. In contrast, with a low coherence source, their model exhibited little interference for cones outside the fovea, as expected for such a model (see their Fig. 5(c)). This outcome of their model is at odds with their AOSLO image obtained with a 680 nm low coherent source, which shows appreciable variation in cone intensity (their Fig. 6, right panel), similar to the images we acquired. The coherence length of their 680 nm superluminescent diode source was not given, but is likely to be within the range we used (5-15 µm). Because their model did not consider reflectors spaced on the order of a few microns, it could not have predicted interference with such short coherence length sources. Considered along with our findings, both sets of results suggest that intensity variation in cones in AOSLO images can arise from interference on at least two scales: one that includes interference between the IS/OS and OST when sources with coherence length 60 µm are used, and one solely between multiple interfaces at the IS/OS and OST junctions when sources 20 µm are employed. Lateral interference remains likely to be prominent for foveal cones, but was not explored in our model or imaging.
Using phase sensitive OCT imaging, it has been observed that the outer segment changes in optical length in response to stimuli in animals [68,69] and in humans [70]. In mouse rods this appears to be due to an osmotic mechanism that leads to swelling as well as refractive index changes [71]. A similar mechanism may operate in human cones. This swelling can be small yet produce a large effect in AOSLO images. Optical path length changes of only λ/4 can cause a cone's brightness to go from a peak to a trough (Fig. 3). Cone brightness changes in response to a visual stimulus have been observed to be oscillatory and occur when the coherence length of the imaging light source is much shorter than 20 µm [11,12,38,39]. Interference effects between close reflectors can explain this intrinsic response: osmotic swelling could cause a path length change between reflectors, modulating the power of the returned light. If modulation is observed with short coherence lengths, the reflectors must be in close proximity. Where modulation is observed using long coherence light to image, interferences between the IS/OS and OST will contribute to the signal [11].
Although cones occasionally appear to be poorly reflective in the AOSLO images (Figs. 3, 6 and 7), the interference mechanism does not explain hyporeflectivity in all situations. Cones with persistently low reflectivity are present in normal retinas, yet retain normal light sensitivity [72]. These dysflective cones are also found in a variety of disease states, having very low reflectivity in AOSLO images but still retaining visual function [18,73,74]. Similarly, photoreceptors are sometimes invisible in OCT, which the models presented here do not explain: Models 2 and 3 always have bright reflections from the IS/OS and OST junctions. These dysflective photoreceptors may be imaged by other modalities such directional OCT [73,74] or split detection imaging [75,76], which suggests that waveguided reflection is disrupted in certain cases, likely from distorted cone shapes or misalignment of their axes. As the FDBP method can be used with arbitrarily shaped structures, these atypical photoreceptors may be investigated using the same technique.
Differential coherence imaging ideally shows only interference effects and emphasizes how cone intensity changes with coherence length. If there was no interference at the IS/OS or OST junctions it is unlikely that differential coherence images would show cone-related structure with the coherence lengths used here. The images in this study were not taken simultaneously and photoreceptors do have a temporal variation, but this largely has a period of 2.5 to 3 hours [11], much longer than the interval between images taken at the same location here (∼1-2 minutes). Significant variation of cone brightness on this timescale has not been reported, but variation on the scale of seconds has been seen, and is described as rare [8]. It has been reported from AOOCT imaging that some cones may have additional reflectors in the outer segment that are not permanent features [77,78], and that the OST band yields a reduced intensity during disc shedding [79]. Both effects could lead to temporal changes between images on short time scales.
The annular profiles that appeared in the differential coherence images (Fig. 8) may provide a method of measuring the cone dimensions in vivo. We explored how the shape of the annuli contains information about the inner and outer segment diameters in a separable way, keeping in mind that annuli will also be dependent on gap sizes. However, it would require imaging in animals to be able to histologically confirm cone dimensions with sufficient resolution. Our model would provide a useful guide for how to optimize such experiments.
The annuli were not as prevalent in AOSLO images as predicted by the models. There are several reasons why this may be the case. First, real-time image stabilization is subject to errors on the order of a few pixels [55], which will diminish the annulus signal when present through the averaging of individual image frames. Second, the models did not include scattering (such as from mitochondria or outer segment discs [80]), which may cause the light reflected from the OST to be more spread out after propagation through the inner segment, resulting in a broader profile. Third, as noted above, the models did not include lateral cone interference. Effects from adjacent cones may mask the annuli; consistent with this, where annuli were observed in vivo, the spacing between cones tended to be greater than average locally. Fourth, Model 3 shows that the presence of an annulus is quite dependent on photoreceptor dimensions (Fig. 10). Differences of only 0.1 µm in the outer segment diameter may determine whether a cone displays an annulus in the difference image. Our simulation in Fig. 8 used a single outer segment diameter, which may not have been correct for the patch of retina imaged in that subject. Unfortunately, measurements from histology of outer segment diameters near the fovea lack sufficient resolution (e.g. Polyak reports measures with a precision of only 0.5 µm between the fovea and 5°eccentricity [81]), and individual variation on the scale of 100 nm will frustrate exact fitting until direct histological confirmation can be performed. The inner segment diameter is better characterized, but there is large variation between cones and subjects [49]. Fifth, the models had idealized flat reflective surfaces: tilted or curved shapes are likely to produce different reflection profiles. Sixth, the refractive indices of the intracellular compartments are not well characterized. The values we used were all measured using white light [82], while imaging is typically performed in the near infrared. All of these factors deserve more considered treatment in order to improve the usefulness of differential coherence imaging.
While not empirically tested here, the modeling suggests that the annulus is in part due to multimodal behavior, which is highlighted by consideration of polarized light (Fig. 11). Such behavior of reflected light in cones has been predicted [32,37,83,84] and observed in microscopy [30,83], AOSLO images [85,86], and AOOCT images [43,87,88]. The model presented here also makes such predictions, with only the IS/OS reflections exhibiting multimodal behavior. This is consistently observed in AOOCT imaging, where the IS/OS band can generate more than one mode, while the OST band typically does not [43,87,88]. If a model only considers bound modes, it will not make this prediction. Multiple modes may exist in a waveguide when the V number is > 2.405 [89], and is dependent on core diameter, wavelength and refractive indices of the core and cladding. Where light enters the cone near the ELM, this condition is met for larger inner segment diameters, where the LP 01 and LP 11 modes may exist with a wavelength of 840 nm. Close to the IS/OS, only the LP 01 is allowed due to its narrowed diameter and despite the higher refractive index of the ellipsoid. Thus, the tapered shape acts as a filter for the higher order modes, which become unbound and radiative [34]. So when light is reflected from the IS/OS it exists only in the fundamental mode, and this will be the only light returned. Due to the short distance traveled axially by light in the photoreceptor, there may not be enough length for unbound light to radiate away, and this may still exist inside the cone. The model presented here does consider the radiative light and it is this that gives rise to the multimodal behavior ( Fig. 11(B)). As the radiative light is reflected, it may excite bound modes that may exist closer to the ELM. Distinct from the inner segment, light in the outer segment is well guided and nearly all contained in the LP 01 mode which is the only one possible when d os < 1.76 µm (Fig. 1), because the radiative modes have the space to leave that portion of the cone. This explains why multimode behavior is not seen in the reflections from the OST, because the double pass of the whole cone acts as a mode filter [90,91].
The excitation of the LP 11 mode in the inner segment causes reflections from the IS/OS to be broad compared to those from the OST (Fig. 11). This has been observed in AOOCT images. Liu et al. [43] used second moment analysis to measure the lateral size of the IS/OS and OST reflections. They found that the lateral sizes diverge past an eccentricity of 3°. Using Eqs. (11) and (12), this eccentricity corresponds to d 0 = 6.5 µm, approximately the threshold size for seeing annuli when d os < 1.8 µm (Fig. 10). As noted in our modeling results, the differential sizes of the IS/OS and OST reflections produce the annuli, so the modeled size threshold is confirmed by AOOCT reflections diverging at the same eccentricity.
The FDBP method has been expanded to include reflections and calculate interference. When the cone model incorporates close, multiple reflections it makes predictions about variation of cone reflectivity in AOSLO images and it corresponds to AOOCT imaging. These predications may also explain phenomena in the wider literature, such as the intrinsic optical response or the recently observed "splitting" of the IS/OS reflection after a stimulus [92]. The model also possibly revealed a new technique to obtain information about photoreceptor morphology in vivo. The model is not specific to the conditions presented here. The shape and size of the photoreceptor is arbitrarily defined, allowing different shapes to be investigated. On the detection side, different configurations such as split detection could be implemented by changing the pinhole shape. By expanding the modeling repertoire it should be possible to gain new information about both normal and diseased photoreceptors.