Light propagation and capture in cone photoreceptors

Alexander Meadway; Lawrence C. Sincich

doi:10.1364/BOE.9.005543

1. Introduction

Most of our day-to-day visual experience derives from signals that originate in cone photoreceptors. Recent experiments have shown that the functional weighting of each photoreceptor varies [1–3]. Such variation will arise, in part, from differences in synaptic strength encountered as the signals flow through retinal circuits. However, there are also optical effects that lead to differences in signal magnitude between photoreceptors. As cones vary in size and shape, the efficiency of light propagation through them will vary from cone to cone, ultimately altering the amount of light absorbed by the photopigment. Consequently it is unclear how much cone signal variation can be attributed to the biophysics of light capture versus downstream neural circuitry. Moreover, it is not known under what conditions optimal light capture can be achieved in vivo, if one is to strive for stimulation of a single cone. To help resolve these ambiguities, we developed a model of light propagation and absorption within photoreceptors to be used as a tool to study such effects.

It is well established that cone photoreceptors act as waveguides. The refractive indices of the myoid, ellipsoid, and outer segment are higher than that of the extracellular media, allowing total internal reflection. Waveguide-like behavior was first observed in photo-receptors ex vivo [4], and more recently it has been studied in vivo using adaptive optics optical coherence tomography [5]. Modal analysis of waveguide propagation has been used to understand how light is transmitted through photoreceptors [6–12], but it has weaknesses when detailed photoreceptor anatomy is considered. The typical length of a cone is about 60 μm in most areas of the human retina, except near the fovea where its length can reach 70 µm [13]. Over such short distances there may not be enough space for radiative modes (non-bound light) to exit the cell. In addition, these models only consider single photoreceptors, and if radiative modes were present, such light could be captured in adjacent cones. Modal propagation techniques may be improved by the inclusion of the radiation modes [14], but where the bound modes are discrete and few in number as in photoreceptors, the radiation modes will be continuous, so they are difficult to work with even with relatively simple conditions [15]. To incorporate the radiative modes, time domain finite difference (TDFD) methods [16] have been used to model single [17,18] and multiple human photoreceptors [19], as well as single avian photoreceptors [20], but the TDFD technique is computationally demanding due to the high spatial and temporal resolutions required for an accurate result.

In this paper we present a waveguide model of cone photoreceptors that employs a finite difference beam propagation method (FDBP) [21–23]. This method has previously been used in a limited way to simulate modal behavior as seen in optical coherence tomography images [5]. The method does not require high resolutions as needed by TDFD, so solutions are faster to compute. As calculation of modes is not necessary, propagation down arbitrarily shaped structures is easily performed, thus the effects of multiple or irregularly shaped photoreceptors can be modeled. All light is considered, both bound and radiative. The model uses a parabolically shaped inner segment, which resembles the normal morphology of cones. Absorption is included to estimate the photoresponse of the cell and polarization is taken into account. We use our model to simulate microstimulation experiments that target single cone photoreceptors with an adaptive optics scanning laser ophthalmoscope (AOSLO) [1,24,25]. We test the model against previously published results and with newly gathered psychophysical data. We also calculate the criteria for optimizing single cone stimulation with varying beam and stimulus sizes.

2. Methods

2.1 Finite difference beam propagation method

The FDBP method finds solutions to the Helmholtz wave equation with a slowly varying envelope approximation. In matrix form, this can be written as

2 i n_{0} k_{0} \frac{\partial}{\partial z} [\begin{matrix} ψ_{x} \\ ψ_{y} \end{matrix}] = [\begin{matrix} H_{x x} & H_{x y} \\ H_{y x} & H_{y y} \end{matrix}] [\begin{matrix} ψ_{x} \\ ψ_{y} \end{matrix}],

where propagation is in the z direction, i = √-1, n₀ is a representative reference refractive index of the system, k₀ is the wavenumber of the light in a vacuum and ψ_x and ψ_y are the polarized electric field components in the x and y directions. In our model, the x-y plane represents the anatomical cross section of the photoreceptor, normally orthogonal to the direction of light entry. The H_ij terms are differential operators. H_xx and H_yy are defined by

H_{x x} ψ_{x} \equiv \frac{\partial}{\partial x} (\frac{1}{n^{2}} \frac{\partial}{\partial x} (n^{2} ψ_{x})) + \frac{\partial^{2} ψ_{x}}{\partial y^{2}} + k_{0}^{2} (n^{2} - n_{0}^{2}) ψ_{x}

and

H_{y y} ψ_{y} \equiv \frac{\partial^{2} ψ_{y}}{\partial x^{2}} + \frac{\partial}{\partial y} (\frac{1}{n^{2}} \frac{\partial}{\partial y} (n^{2} ψ_{y})) + k_{0}^{2} (n^{2} - n_{0}^{2}) ψ_{y}

where n is the refractive index, which will be a function of x, y and z. We consider photoreceptors to be weakly guiding, (n_core – n_cladding) / n_core << 1, so we assume there is no cross-coupling of polarizations, allowing H_xy and H_yx to be set to zero. This is the semi-vectorial scheme.

The FDBP method represents Eq. (1) in terms of its finite differences. The field is discretized into a grid. The field at each point in this grid is related to the field at every adjacent x and y point. When the field is known in one lateral plane, FDBPM can be calculated at the next step at a distance Δz away. Thus, when the initial field and the spatial distribution of the refractive indices are known or assumed, the propagation of the field through a whole photoreceptor can be modeled step by step in 3-D.

Our implementation of FDBP follows that derived by Pedrola [26]. To avoid a large matrix problem, the alternating direction implicit method is used [27]. This method splits the step of Δz into two substeps where each considers only the x or y direction. Computing the beam propagation as a pure 3-D problem yields a matrix equation to be solved with a diagonal matrix with 7 non-zero elements in each row. This method splits the problem into two orthogonal 2-D problems, each yielding a tridiagonal system which can be solved efficiently. By separating the matrix operators into x and y component terms, Eq. (1) becomes

2 i n_{0} k_{0} \frac{\partial}{\partial z} [\begin{matrix} ψ_{x} \\ ψ_{y} \end{matrix}] = [\begin{matrix} A_{x} + A_{y} & 0 \\ 0 & B_{x} + B_{y} \end{matrix}] [\begin{matrix} ψ_{x} \\ ψ_{y} \end{matrix}] .

The operators A_i are defined by

A_{x} ψ_{x} \equiv \frac{\partial}{\partial x} (\frac{1}{n^{2}} \frac{\partial}{\partial x} (n^{2} ψ_{x})) + \frac{1}{2} k_{0}^{2} (n^{2} - n_{0}^{2}) ψ_{x}

and

A_{y} ψ_{x} \equiv \frac{\partial^{2} ψ_{x}}{\partial y^{2}} + \frac{1}{2} k_{0}^{2} (n^{2} - n_{0}^{2}) ψ_{x} .

The B_i operators are defined by

B_{x} ψ_{y} \equiv \frac{\partial^{2} ψ_{y}}{\partial x^{2}} + \frac{1}{2} k_{0}^{2} (n^{2} - n_{0}^{2}) ψ_{y}

and

B_{y} ψ_{y} \equiv \frac{\partial}{\partial y} (\frac{1}{n^{2}} \frac{\partial}{\partial y} (n^{2} ψ_{y})) + \frac{1}{2} k_{0}^{2} (n^{2} - n_{0}^{2}) ψ_{y} .

The x-polarized light may be written in discretized form, using the Crank-Nicholson scheme:

2 i n_{0} k_{0} \frac{u^{m + 1} - u^{m}}{Δ z} = (A_{x} + A_{y}) \frac{u^{m + 1} + u^{m}}{2},

replacing ψ_x with u for clarity. This equation relates the field at the (m + 1)^th step, u^m⁺¹, to the field at the m^th, u^m. Separating this equation into two substeps with an intermediate point at m + ½ gives

(\frac{4 i n_{0} k_{0}}{Δ z} - A_{y}) u^{m + ½} = (\frac{4 i n_{0} k_{0}}{Δ z} + A_{y}) u^{m}

and

(\frac{4 i n_{0} k_{0}}{Δ z} - A_{x}) u^{m + 1} = (\frac{4 i n_{0} k_{0}}{Δ z} + A_{x}) u^{m + ½} .

For the first substep, u (the x-polarized field) is discretized into a grid of N_x × N_y points, separated by the distances Δx and Δy. The field at the point (i, j, m + ½) is related to the field at (i, j, m) by

\begin{matrix} \frac{4 i n_{0} k_{0}}{Δ z} u_{i, j}^{m + ½} - \frac{u_{i, j - 1}^{m + ½} - 2 u_{i, j}^{m + ½} + u_{i, j + 1}^{m + ½}}{Δ y^{2}} - \frac{1}{2} k_{0}^{2} [{(n_{i, j}^{m + ½})}^{2} - n_{0}^{2}] u_{i, j}^{m + ½} = \\ \frac{4 i n_{0} k_{0}}{Δ z} u_{i, j}^{m} + \frac{u_{i, j - 1}^{m} - 2 u_{i, j}^{m} + u_{i, j + 1}^{m}}{Δ y^{2}} - \frac{1}{2} k_{0}^{2} [{(n_{i, j}^{m})}^{2} - n_{0}^{2}] u_{i, j}^{m} \end{matrix}

For each value of i and j a system of N_x linear equations can be built:

a_{j} u_{i, j - 1}^{m + ½} + b_{j} u_{i, j}^{m + ½} + c_{j} u_{i, j + 1}^{m + ½} = q_{j}^{m},

where

\begin{array}{l} a_{j} = - \frac{1}{Δ y^{2}}, \\ b_{j} = \frac{4 i n_{0} k_{0}}{Δ z} + \frac{2}{Δ y^{2}} - \frac{1}{2} k_{0}^{2} [{(n_{i, j}^{m + ½})}^{2} - n_{0}^{2}], \\ c_{j} = - \frac{1}{Δ y^{2}}, \\ q_{j} = \frac{1}{Δ y^{2}} [u_{i, j - 1}^{m} + u_{i, j + 1}^{m}] + {\frac{4 i n_{0} k_{0}}{Δ z} - \frac{2}{Δ y^{2}} + \frac{1}{2} k_{0}^{2} [{(n_{i, j}^{m})}^{2} - n_{0}^{2}]} u_{i, j}^{m} \end{array}

which is a tridiagonal system. In matrix form this is

[\begin{matrix} b_{1} & c_{1} & 0 & 0 & . & . & 0 & 0 & 0 \\ a_{2} & b_{2} & c_{2} & 0 & . & . & 0 & 0 & 0 \\ 0 & a_{3} & b_{3} & c_{3} & . & . & 0 & 0 & 0 \\ . & . & . & . & . & . & . & . & . \\ . & . & . & . & . & . & . & . & . \\ 0 & 0 & 0 & 0 & . & . & a_{N - 1} & b_{N - 1} & c_{N - 1} \\ 0 & 0 & 0 & 0 & . & . & 0 & a_{N} & b_{N} \end{matrix}] [\begin{matrix} u_{1}^{m + ½} \\ u_{2}^{m + ½} \\ u_{3}^{m + ½} \\ . \\ . \\ u_{N - 1}^{m + ½} \\ u_{N}^{m + ½} \end{matrix}] = [\begin{matrix} q_{1} \\ q_{2} \\ q_{3} \\ . \\ . \\ q_{N - 1} \\ q_{N} \end{matrix}]

which we write as

M_{y} u_{i, j}^{m + ½} = q_{y},

and can be solved using several algorithms if u^m, n^m and n^m+^½ are known.

For the second substep for x-polarized light, from m + ½ to m + 1, the finite differences equation is

\begin{matrix} \frac{4 i n_{0} k_{0}}{Δ z} u_{i, j}^{m + 1} - \frac{T_{i - 1, j}^{m + 1} u_{i - 1, j}^{m + 1} - 2 R_{i, j}^{m + 1} u_{i, j}^{m + 1} + T_{i + 1, j}^{m + 1} u_{i + 1, j}^{m + 1}}{Δ x^{2}} - \frac{1}{2} k_{0}^{2} [{(n_{i, j}^{m + 1})}^{2} - n_{0}^{2}] u_{i, j}^{m + 1} = \\ \frac{4 i n_{0} k_{0}}{Δ z} u_{i, j}^{m + ½} + \frac{T_{i - 1, j}^{m + ½} u_{i - 1, j}^{m + ½} - 2 R_{i, j}^{m + ½} u_{i, j}^{m + ½} + T_{i + 1, j}^{m + ½} u_{i + 1, j}^{m + ½}}{Δ x^{2}} - \frac{1}{2} k_{0}^{2} [{(n_{i, j}^{m})}^{2} - n_{0}^{2}] u_{i, j}^{m + ½} \end{matrix}

where

T_{i \pm 1, j} \equiv \frac{2 n_{i \pm 1, j}^{2}}{n_{i + 1, j}^{2} + n_{i, j}^{2}}

and

R_{i, j} \equiv 2 - \frac{T_{i + 1, j} + T_{i - 1, j}}{2} .

A second tridiagonal system is created for each value of i and j;

M_{x} u_{i, j}^{m + 1} = q_{x},

with the coefficients

\begin{array}{l} a_{i} = - \frac{1}{Δ x^{2}} T_{i - 1, j}^{m + 1}, \\ b_{i} = \frac{4 i n_{0} k_{0}}{Δ z} + \frac{2}{Δ x^{2}} R_{i, j}^{} - \frac{1}{2} k_{0}^{2} [{(n_{i, j}^{m + 1})}^{2} - n_{0}^{2}], \\ c_{i} = - \frac{1}{Δ x^{2}} T_{i + 1, j}^{m + 1}, \\ q_{i} = \frac{1}{Δ x^{2}} [T_{i - 1, j}^{m + ½} u_{i, j - 1}^{m + ½} + T_{i + 1, j}^{m + ½} u_{i, j + 1}^{m + ½}] + {\frac{4 i n_{0} k_{0}}{Δ z} - \frac{2 R_{i, j}^{m + 1 / 2}}{Δ x^{2}} + \frac{1}{2} k_{0}^{2} [{(n_{i, j}^{m + 1 / 2})}^{2} - n_{0}^{2}]} u_{i, j}^{m + 1 / 2} . \end{array}

For y-polarized light a similar set of equations is built. For the first substep:

M_{y} v_{i, j}^{m + ½} = q_{y},

where v is ψ_y at a step denoted by the superscript. This has the coefficients

\begin{array}{l} a_{j} = - \frac{1}{Δ y^{2}} T_{i, j - 1}^{m + ½}, \\ b_{j} = \frac{4 i n_{0} k_{0}}{Δ z} + \frac{2}{Δ y^{2}} R_{i, j}^{m + ½} - \frac{1}{2} k_{0}^{2} [{(n_{i, j}^{m + ½})}^{2} - n_{0}^{2}], \\ c_{j} = - \frac{1}{Δ y^{2}} T_{i, j + 1}^{m + ½}, \\ q_{j} = \frac{1}{Δ y^{2}} [T_{i, j - 1}^{m} v_{i, j - 1}^{m} + T_{i, j - 1}^{m} v_{i, j + 1}^{m}] + {\frac{4 i n_{0} k_{0}}{Δ z} - \frac{2 R_{i, j}^{m}}{Δ y^{2}} + \frac{1}{2} k_{0}^{2} [{(n_{i, j}^{m})}^{2} - n_{0}^{2}]} v_{i, j}^{m}, \end{array}

where

T_{i, j \pm 1} \equiv \frac{2 n_{i, j \pm 1}^{2}}{n_{i, j}^{2} + n_{i, j \pm 1}^{2}}

and

R_{i, j} \equiv 2 - \frac{T_{i, j + 1} + T_{i, j - 1}}{2} .

For the second substep for y-polarized light,

M_{x} v_{i, j}^{m} = q_{x},

with coefficients

\begin{array}{l} a_{i} = - \frac{1}{Δ x^{2}}, \\ b_{i} = \frac{4 i n_{0} k_{0}}{Δ z} + \frac{2}{Δ x^{2}} - \frac{1}{2} k_{0}^{2} [{(n_{i, j}^{m + 1})}^{2} - n_{0}^{2}], \\ c_{i} = - \frac{1}{Δ x^{2}}, \\ q_{i} = \frac{1}{Δ x^{2}} [v_{i - 1, j}^{m + ½} + v_{i + 1, j}^{m + ½}] + {\frac{4 i n_{0} k_{0}}{Δ z} - \frac{2}{Δ x^{2}} + \frac{1}{2} k_{0}^{2} [{(n_{i, j}^{m + ½})}^{2} - n_{0}^{2}]} v_{i, j}^{m + ½} . \end{array}

Our implementation of FDBP was developed in Matlab and utilizes parallel processing to perform up to 6 propagations simultaneously. The routine takes the field at each step and uses Matlab’s backslash function to solve Eqs. (16) and (20), finding the field at each step for the x-polarization and Eqs. (22) and (26) for the y. The initial field (ψ₀) was discretized into a square grid of N × N points with Δx = Δy. To suppress reflections from the edge of the computing space an absorbing boundary condition was used. To handle this condition, before the field is propagated to the next step it is multiplied by a trapezoid window varying in value from 0 at the edges and 1. The slopes of the trapezoid chosen made a compromise between suppressing the reflections and minimizing computing space and were 0.6 μm⁻¹ for Δx = Δy = 0.1 μm.

The power of the light within a photoreceptor was calculated at each z position. This is given by the integral of the modulus squared of the field;

P (z) = \int_{0}^{r = d / 2} {| ψ (z) |}^{2} d r,

where r² = x² + y² and d is the diameter of the cone at z, as defined in the next section. We assumed that the signal produced in the psychophysical experiments (Section 2.5) reflects absorption by outer segment photopigments and is proportional to the power absorbed. To include absorption using FDBP all that is needed is a complex component in the refractive index. To calculate the absorbed power two propagations at each step are required. The field at m is propagated to m + 1 using the real and complex values for the refractive index. The power within the photoreceptor is then calculated for each propagation using Eq. (28) and the difference gives the power absorbed at that step:

P_{a b s} (z) = P_{r e a l} (z) - P_{c o m p l e x} (z) .

The total power absorbed by the photoreceptor(s) is then given by

P_{t o t a l} = \sum_{0}^{z = L} P_{a b s} (z),

where L is the total length of the cone. To estimate a subject’s change of perception to a changed stimulus we assume that all cones have equal weighting and sum linearly.

2.2 Modeled cone properties

FDBP modeling of the propagation of light down a cone photoreceptor requires a 3-D map of the refractive indices. To develop this map, we first consider the shape of the cone. Previous modeling simplified cone photoreceptors into one of two shapes. The first had a cylindrical myoid and outer segment, joined by a tapered ellipsoid [6,10,28]. The second had two steps, again with a cylindrical myoid and outer segment, but a cylindrical ellipsoid with a diameter between that of myoid and outer segment [17]. Our model uses an inner segment shape that diminishes parabolically and a cylindrical outer segment. This shape is qualitatively similar to the morphology of parafoveal cones and follows the suggestion that cones act as optical collectors [29].

In the model, the diameter, d, changes along its length by

d (z) = {\begin{cases} \frac{d_{o s} - d_{0}}{l_{i s}^{2}} z^{2} + d_{0} for 0 \leq z \leq l_{i s} \\ d_{o s} for z \geq l_{i s} \end{cases}

where d_os is the diameter of the outer segment, d₀ is the diameter of the inner segment at z = 0, defined at the external limiting membrane (ELM) where the anterior face of the inner segment begins, and l_is is the length of the inner segment, equal to the length of the myoid (l_m) plus the length of ellipsoid (l_e).

Exact dimensions of the three photoreceptor segments for cells situated between 1° and 5° eccentricity from the fovea are not well reported in the literature, so estimates must be made for lengths and diameters. Cone spacing as a function of eccentricity is well characterized from measurements in vivo and ex vivo [30–33]. To estimate d₀ from eccentricity, we fit a 2nd order polynomial to the data presented by Song et al. [33]. Cone density is converted to cone spacing, under the assumption of a hexagonal packing structure. It is then assumed that d₀ is a factor of 0.82 of cone spacing [7,34–36], yielding a curve of

d_{0} = - 0.049 φ^{2} + 1.22 φ + 2.68

where φ is the eccentricity in degrees. Values obtained from this curve are consistent with other results in the literature in the eccentricity range of interest. For d_os we use a value of 1.75 μm for all eccentricities, measured from Fig. 1 of Boycott et al. [37] and allowing for 15% shrinkage. Our modeled values of d₀ ranged from 4 to 8 μm. For the lengths of the inner segments, we use the perifoveal values for l_m ( = 13 μm) and l_e ( = 16 μm) from Spaide and Curcio [13], but the foveal value for l_os ( = 35 μm). These values are in reasonable agreement with optical coherence tomography measures of the inner and outer segments [38].

Fig. 1 Refractive index traits of the cone photoreceptor model. The photoreceptor shape, defined by Eq. (31), is shown in 2-D and 3-D representations. Color code depicts the different refractive indices (n) of the compartments: red = myoid, green = ellipsoid, blue = outer segment and white = surrounding media. Dimension variables are d₀ = 5 μm, d_os = 1.75 μm, l_m = 13 μm, l_e = 16 μm, and l_os = 35 μm.

Download Full Size | PDF

Where multiple cones are modeled, a hexagonal packing structure is assumed [32,39]. All cones are given the same properties. A 3-D representation of the arrangement of multiple cones is shown in Fig. 1.

Refractive indices for the myoid, ellipsoid, and outer segment were n_m = 1.3605, n_e = 1.394 and n_os = 1.39 respectively [13]. The extracellular media is given a refractive index of n_ex = 1.34 [5–7,40]. We note that no wavelength dependent values for refractive index could be found in the literature, and the values used were all based on measurements in the visible range.

The imaginary component of the refractive index accounts for the absorption of a media and its value is given by the extinction coefficient (κ). For all model runs 543 nm light was used, this being the wavelength of the light in our validation experiments (see Section 2.4). The extinction coefficient at this wavelength was calculated from an estimated absorption coefficient (α = 0.0167 μm⁻¹ at 535 nm) from macaque M cones [41]. Using κ = wαλ/4π, where w is a weighting based on the normalized spectral sensitivity [42], a value of κ = 7.22 × 10⁻⁴ was found for λ = 543 nm, a wavelength equally absorbed by both L and M cones.

If the value for the reference refractive index, n₀, is chosen poorly, errors may appear in the form of instability or numerical dissipation. For straight, single mode waveguides the effective refractive index of the fundamental mode may be chosen [23], or, if the guide is narrow, the index of the surrounding media can be used [26]. Neither of these options are suitable in this situation because of the morphology of our modeled volumes, especially with the inclusion of multiple cones. Instead, n₀ was chosen by optimization, as follows.

If absorption in the outer segment is ignored (only the real part of n_os is used and the absorbing boundary condition not used), the power at each propagation step should equal the input power, P(z) = P₀. If any errors are present, the power will deviate. Using a single cone, the electric field is propagated initially using an estimate for n₀. The reference index is then adjusted, and the propagation repeated until the power’s deviation from a value of 1 is minimized over all longitudinal distances. A value of n₀ = 1.37 was found using a cone with d₀ = 5 μm, d_os = 2 μm and λ = 534 nm and gave a peak power error of less than 1% in the range d₀ = 3 to 8 μm. This n₀ value was used in all model runs presented here.

2.3 The initial field, ψ₀

In this model ψ₀ is generated by assuming the electric field at the eye’s pupil and performing the Fourier transform to calculate the field at the focus of the system, at the ELM (z = 0). This method is advantageous as it allows easy adjustment of optical variables in the model that may be changed or measured in an experiment, such as defocus, other aberrations, or beam size.

When calculating the resolution of a system it is commonly assumed that the light across the beam has a uniform amplitude distribution. This is rarely the case in a real system. Free space lasers have a Gaussian intensity distribution profile, as does light launched into a system via an optical fiber. As such, when calculating ψ₀ here we used a Gaussian amplitude distribution that is then truncated to the beam diameter (D). The width of the Gaussian used was measured at the pupil of the adaptive optics system described in Section 2.4. The field is defined by the optics collimating the light as it exits the source or input fiber, the magnification of the various telescopes that relay the light to the eye, and any optical stops, such as mirror or lens edges. For typical adaptive optics systems with pupils ranging from 5.5 to 7.5 mm, the point spread function (PSF) from the truncated Gaussian has suppressed sidelobes, but has lower resolution, as measured by the full-width half-maximum (Fig. 2). The eye length is assumed to be 22.2 mm with a refractive index of 1.33. For normalization, the power of the light, the integral of |ψ₀|², is set to 1.

Fig. 2 Shape of the initial electric field ψ₀. Modulus of the electric field PSF calculated by Fourier transform of the pupil using a uniform beam profile (A) or with a truncated Gaussian beam (B). (C) Field profiles at y = 0 through the uniform (black line) and truncated Gaussian pupil fields (red line). For both, λ = 543 nm and D = 5.85 mm. The Gaussian function of the truncated field had a 1/e width of 4.75 mm (1/e² width for the intensity). Both fields are normalized to their peak value. The intensity full-width half-maxima were 1.6 µm for the uniform beam and 2.0 µm for the truncated Gaussian beam.

Download Full Size | PDF

2.4 Adaptive optics microstimulation

Cone-targeted microstimulation with a multi-wavelength AOSLO was used for the experiments described in Section 2.5 and is the basis for our modeling of light capture. The instrument has been described in detail previously [24,43], so its operation will be outlined only briefly here. For imaging, light from a superluminescent diode (Superlum, Ireland) centered at 842 ± 25 nm was used. Retinal images of 512 × 512 pixels were acquired at 30 Hz with a raster scan. The image size for all measurements was 1.2°. Optical aberrations were corrected by utilizing a Shack-Hartmann wavefront sensor and a MEMS deformable mirror (Boston Micromachines, Cambridge, MA) working in a closed loop. The wavefront sensor uses a portion of the 842 nm light picked off before detection. The green (543 ± 11 nm) and red (711 ± 12 nm) channels available for stimulation are separated from the spectrum of a supercontinuum laser (NKT Photonics, Denmark). In these experiments only green light was used, as 543 nm is approximately equally absorbed by human L and M cones. Reflected light from each channel may be detected and used to form an image using photomultiplier tubes (Hamamatsu, Japan). All wavelengths were corrected for longitudinal chromatic dispersion by suitable positioning of the input fiber and confocal pinhole planes.

Stimuli were controlled by acousto-optic modulators (Brimrose Corporation, Sparks, MD), temporally registered to the raster scan, and allowed 10-bit modulation of the beam intensities. We report stimulus intensities in arbitrary units from 0 to 1, where 0 represented the background light and 1 the maximum amount that could be delivered. Typically, for a 3 × 3 pixel stimulus, the background level incident on the cornea was ~2.7 log quanta, while the maximum was ~4.3 log quanta. This background level was enough to eliminate any rod contributions to the threshold measurement [24]. In order to accurately and repeatedly stimulate the same cone target, eye movements were tracked via image stabilization in real time [44,45]. During an experiment, the center of a cone was selected manually in the stabilized image for stimulus placement; this pixel defined where the geometric center of any stimulus was delivered.

Calculation of the light absorbed by a raster scanned stimulus in the model must take into account the spatial and temporal status of the beam. Since each physical stimulus point is separated in space and time, the modeled stimulus must be discretized. The discretization creates a grid of points that have a spacing chosen to be dense enough to approximate a real stimulus, but also sparse enough so that the time taken to run the model does not become prohibitive. This spacing was set to 0.9 μm, less than the Airy radius of the eye at 543 nm (1.9 μm). We then calculated the PSF centered at each stimulus point and propagated this through the computing space. Two stimulus shapes were used for experiments here; square and circular. The square stimulus is defined by the length of its edges. There was always a point that defined the center of any stimulus. This point was used to define the stimulus position relative to the center of a targeted cone. The circular stimulus is discretized in a similar way, its size is defined by its diameter. The response to a modeled stimulus is assumed to be proportional to the average power absorbed over the entire discretized stimulus.

2.5 Psychophysical evaluation of the light capture model

We used the FDBP model to make predictions of light capture for delivered stimuli in our AOSLO, and compared the predictions to results gained in two different experiments. In the first experiment, the stimulus was varied in 3 ways: (1) changing the beam size at the source of the input to the optical system, thereby altering the focused spot size, but keeping the total power at the cornea constant; (2) stimulus geometry, either circular or square; and (3) stimulus size. We expected that changes in the initial electric field lead to changes in the resulting propagation and consequently the amount of light absorbed in the outer segments. Two male subjects with normal vision were tested at 3° eccentricity near the horizontal meridian, where the inner segment diameters were expected to be ~5.9 μm using Eq. (32). The diameter of the stimulus beam was varied by an iris placed at the input of the stimulus light to the system. Two diameters were used: 5.85 mm (large) and 2.1 mm (small). To compensate the greater power of the large beam, neutral density filters of an appropriate amount were added to the beam’s path. The beam sizes chosen were constrained by the neutral filter steps sizes available in order to match powers. The power of the source was monitored by tapping 1% of its output to a photodiode using a fiber splitter before the acousto-optic modulator. This was used to ensure that the source power had stabilized before the start of the experiment and to check that the power did not fluctuate appreciably during the experiment. To ensure stimulus light was delivered to the target cones, transverse chromatic aberration (TCA) was measured and compensated [1,24,25,43]. TCA was also measured at the end of the experiment to check that there was no positional drift in the delivered stimuli.

Perceptual threshold is defined as the stimulus intensity that yielded a 0.5 probability of seeing. Thresholds were estimated by the method of constant stimuli, using 7 levels of stimulus intensity, presented pseudorandomly for each of 20 trials. To test threshold change, the large/small diameter of the iris was chosen at random for the first set of trials, the diameter was then switched for a second set of trials, and then the first condition was repeated to assess consistency. Subjects were dark adapted for at least 15 min before any trials began. Psychometric functions were generated by fitting a logistic function to the frequency of seeing data by parametric bootstrapping. Measured results were then compared to modeled results using the same parameters for the AOSLO and stimuli, quantified as a threshold ratio between the large and small beam conditions.

In the second experiment, we modeled the conditions in Harmening et al. (2014) where increment threshold measurements were used to compare the detection of a stimulus delivered to the cone center versus to the gaps between cones [24]. In that experiment the ratios of the thresholds for cones and gaps using stimuli of various sizes and at varying eccentricities were studied. We applied the FDBP model to this task. We ran the model for a stimulus size of 2.3 μm (3 pixels) over a range of 1 to 5.6° eccentricity. Because the experiment was performed on the same AOSLO system, we used the same initial conditions as in the beam size experiment.

The data in Harmening et al. (2014) suggested that the cone-gap measurement was sensitive to defocus, needing ~0.06 D to fit their model to the results. Therefore, we also varied the Zernike coefficient for defocus from −0.1 to 0.1 μm when calculating the initial field. This is equivalent to −0.08 to 0.08 D for a 5.85 mm pupil. We calculated the gap to cone threshold ratio, assuming that the response was linear to the amount of light absorbed. The best fitting defocus was found by calculating the least squares of the original data to the model data.

All psychophysical experiments followed the tenets of the Declaration of Helsinki and were approved by the Institutional Review Board for Human Use at the University of Alabama at Birmingham. Written informed consent was obtained from the participants after the nature and possible consequences of the study was explained.

3. Results

3.1 Propagation of light through cone photoreceptors

To help understand how light capture by photoreceptors depends on the initial conditions of the electric field and photoreceptor size, we generated 8 examples of propagation. Figure 3(A) shows propagation through water (n = 1.33, n₀ = 1.33), with a uniform field at the pupil, yielding the expected field profile. The result using a truncated Gaussian pupil (Fig. 3(B)) is similar, but is blurred due to the Gaussian field profile at the pupil, as this acts as a spatial filter. The next panels highlight how the model shows waveguiding, with photopigment absorption taken into account. Normal propagation is depicted in Fig. 3(C), with the beam centered on a photoreceptor, concentrating the light into the outer segment. In this example (as well as in Fig. 3(E) and 3(F)) a beating of the light can be observed. This is a display of modal behavior and it shows that the FDBP method captures this behavior. With the wavelength, refractive indices, and diameter used here, the outer segment supports two linearly polarized propagation modes, LP₀₁ and LP₁₁. As these propagate at different speeds down a guide, a beating of the field magnitude occurs as they move in and out of phase with each other [46]. The exact shape of the electric field in the outer segment is dependent on the initial field, how the inner segment guides the light to the outer, and the radiation modes. Figure 3(D) shows a displaced beam, incident between two photoreceptors, with light being guided down two cones. With a small beam (Fig. 3(E)), the focused spot is wider than in the case of the larger beam (Fig. 3(C)) and more light is guided into the outer segments of the adjacent cones. In Fig. 3(F) the photoreceptors’ initial diameter is increased from 4 to 5 μm; the result is similar to Fig. 3(C), but with clear differences, most notably a higher peak field magnitude within the outer segment but with a more confined distribution. Figure 3(G) shows light incident on the central cone at angle of 7.5°. Total internal reflection still guides the light towards the outer segment, but more leaves the photoreceptor than with the normally incident light. Notably, for the bulk of the light there is only one reflection in the inner segment before light reaches the outer segment. To show that aberrated beams can be modeled easily, we provide an example of a beam with coma (Fig. 3(H)).

Fig. 3 Examples of light propagation in cone photoreceptors. Propagation in water (n = 1.33) using a uniform field at the pupil (A) or with a truncated Gaussian (B). When light is incident centrally on a photoreceptor (C) beating of the field is evident in the outer segment. (D) Propagation when beam lands between photoreceptors (shifted 0.61d₀). (E) Diminished light intensity occurs when the beam size is set to 2.0 mm. (F) Propagation down larger photoreceptors (d₀ = 5 μm) leads to the largest peak intensity in these examples. (G) Propagation with light entering at an angle of 7.5°. (H) Example with a Zernike coefficient corresponding to coma of 0.1 μm. Modulus of the electric field for the initial fields (right) and the fields in the x-z plane for x-polarized light (left) are normalized to the peak electric field, which occurred in panel F. Light in the negative z direction was calculated by propagating the complex conjugate of ψ₀ through a refractive index of 1.33. For all examples, unless stated, beam diameter = 5.85 mm, 1/e width of Gaussian = 4.75 mm, λ = 543 nm, d₀ = 4 μm, d_os = 1.75 μm, dx = dy = 0.1 μm, dz = 0.2 μm. Values for refractive index and lengths of segments are given in the Methods.

Download Full Size | PDF

In Fig. 4(A) we compare the calculated power at each z position inside the central photoreceptor for conditions C-H of Fig. 3, using Eq. (29). We also show the amount of power absorbed along the outer segment (Fig. 4(B)), though it should be noted that this value is dependent on Δz, here equal to 0.2 μm. In none of these propagation conditions does the power decrease monotonically. This implies that light not only leaves the receptor as radiated light, but the changing shape, particularly at the junction between the inner and outer segment, allows light to re-enter. The transition between inner and outer segments can cause abrupt changes, most obvious in the case of the larger photoreceptor (Fig. 4(B), red), due to the sharp change in photoreceptor shape here. Looking at the inner segment (z < 29 μm), some noise in the curves can be seen. This is due to the discretization of the photoreceptor into voxels and may be reduced by increasing the resolution, at the expense of a longer computing time.

Fig. 4 Power propagated and absorbed along photoreceptor length. (A) Fraction of power present along length of central photoreceptor. (B) Fraction of power absorbed in outer segment of central photoreceptor. Values calculated from the propagation conditions in Fig. 3(C)-3(H), shown color coded here.

Download Full Size | PDF

In Table 1, we give numerical data for the total power absorbed in the top, center, and bottom cones shown in Fig. 3 calculated using Eq. (30), the total power absorbed by all cones, and the fraction of the total power absorbed by the central cone. Although the four other cones are included in the propagation and in the sum for the total, they are omitted individually for simplicity.

Table 1. Power absorbed in outer segments for propagations in Fig. 3.

View Table | View all tables in this article

Some of these results are straightforward. For example, if the beam covers more cones, as with the small pupil or coma cases, then the fraction of absorbed power in the center cone is reduced. Other results may be more counter-intuitive. For instance, the large cone absorbs less total light than the smaller cone because the steeper inner segment walls lead to more radiative losses, as the electric field shape is less optimal for propagation down the outer segment. The propagated field does achieve a higher single point magnitude in the large cone, but overall it has lower power than the small cone at other positions within the outer segment. Despite the large cone absorbing less light, the ratio of center cone to total absorption is similar to that of the smaller cone. This can be explained by two parallel effects: while the light coupling from inner to outer segment is poorer, this is matched by the cones being larger and more spread out, so less light is incident on and guided down the surrounding cones. Another interesting result is the power absorbed by the bottom cone in the shifted case. More light is absorbed by this cone than may be expected due to unbound light from the center cone entering its outer segment. A similar effect is seen in the angled case, where unbound light can be seen to leak into the adjacent photoreceptors. Where the light is centrally placed on a photoreceptor and directed along its axis, it is much more efficiently guided to the outer segment. Here, as in all runs, we use unpolarized light, with the orthogonal components propagated separately. No significant differences in the results between the polarizations were found.

3.2 Optimization of light capture by single cones

In the previous section it is clear that predicting light absorption in a photoreceptor engages complex relationships between cone size and beam diameter. We have not yet considered the influence of stimulus size, which will naturally have major consequences for the amount of light captured by any array of cones. Here we systematically varied these factors in the model to delimit the parameters that optimally activate a single cone with adaptive optics microstimulation.

Considering only the central cone, we found that a 4 μm cone has the highest absorption of 0.21 with a beam size of 5.25 mm and a point stimulus (defined as 0 μm, Fig. 5(A)). This is greater than a 6 μm cone which has peak absorption of 0.14 at 4.5 mm with a point stimulus (Fig. 5(D)). As the surrounding cones absorb little light with point stimuli, neither peak moves when all cones are considered (Fig. 5(C), 5(F)). Despite having peaks in slightly different places, at no point do the 6 μm cones absorb light more efficiently than the 4 μm cones. As might be expected, larger stimuli and smaller beams at the pupil complicate the picture; both cause light to be incident on the surrounding photoreceptors (Fig. 5(B), 5(E)).

Fig. 5 Power absorbed by cones for d₀ = 4 μm (A-C) and 6 μm (D-F) using a square stimulus. Absorption is based on the initial field having unit power. Grayscale values are interpolated (bicubic) from data computed at intervals of 0.25 mm for beam size and 0.9 μm for stimulus size. Absorbed power contours are in steps of 0.01 (thin lines) and 0.05 (thick lines).

Download Full Size | PDF

By dividing the power absorbed in the central cone by the sum of the central and surrounding cones, we can find the proportion of power absorbed by the central cone, which tells us what the optimal conditions are for getting the most light into a single cone. Starting with a point stimulus, the proportion of light captured in a targeted 4 μm cone reaches 0.99 with a 5.5 mm beam, while a 6 μm cone achieves 0.99 with a 4.75 mm beam (Fig. 6). From an experimentalist’s perspective, the model’s point stimulus delivers a very small amount of energy (a product of power and dwell time), and makes psychophysical testing a challenge with our microstimulation instrument. More realistic stimuli have dimensions on the order of a cone diameter or larger. If we take a 4 μm stimulus and 6 mm beam as an example, Fig. 6 shows that ~75% of the light will be absorbed by a 4 μm cone while ~98% of the light will be picked up a 6 μm cone. Thus optimal microstimulation of one cone in an array is going to be a function of stimulus size, cone size, and beam size. There are regimes shown in Fig. 6 where nearly all the light can be captured by one targeted cone. While the absorption by neighboring cones can be small, it is never non-zero and they capture more light when they are smaller and more closely packed.

Fig. 6 Fraction of total power absorbed by the central cone for d₀ = 4 µm (A) and 6 µm (B). Absorbed power contours indicated on color bar.

Download Full Size | PDF

It should now be clear that there will be an optimum beam size for the most efficient delivery of a stimulus to a cone of a given size. While the optimum delivery depends on many factors, we illustrate this optimum in Fig. 7 for a standard case. Here a point stimulus of varying beam diameter is delivered to the central cone on axis while the diameter of the inner segment is varied. By finding the peak absorption for each cone size, we can see how the optimum beam size varies. As shown by the solid red line in Fig. 7, the optimum beam size is non-linear with respect to cone diameter, but remains in the range of 4.5 to 5.5 mm. We find that there is also a secondary ridge, shown by the dashed line. This corresponds to conditions where the surrounding cones are absorbing light more efficiently, as can be seen in Fig. 5.

Fig. 7 Optimum beam size for most efficient light absorption in all 7 cones, for a range of cone sizes. Solid line follows the peak power absorption for a given cone size with a point stimulus. A secondary peak exists (dashed line) because of absorption in surrounding cones.

Download Full Size | PDF

3.3 Psychophysical validation of the light capture model

We took two approaches to assess whether predictions from the light capture model could be revealed empirically. In the first experiment, we altered the beam size but kept the beam power constant, and measured its effect on psychophysical thresholds for several stimulus configurations. Figure 8(A) shows an example stimulus (a 7.92 μm square) superimposed on the cone targeted in Subject 1. When the beam size was changed from 5.85 mm to 2.1 mm (equivalent to a full-width half-maximum PSF change from 2.0 to 4.4 µm), the psychometric threshold shifted to the left, indicating that less light was required to reach threshold for the smaller beam (Fig. 8(B)). This result is expected because the large stimulus in combination with the large focal spot places light onto surrounding cones as well as the central targeted one (for this subject at this location, smaller stimuli were not detectable). This experiment configuration is depicted as the top pair of points on the modeled threshold ratio map shown in Fig. 8(C) for square stimuli with various beam sizes. The predicted threshold change from the model for this experiment was 0.68, which was within 1 SD of the observed threshold change of 0.76 ± 0.08. All of the remaining experiments (also depicted in Fig. 8 as red-green dot pairs) had observed threshold ratios that fell within 1 SD of the predicted outcomes for different stimulus shapes and sizes (Fig. 8(C), 8(D) and Table 2). Errors were calculated by combining the 95% confidence interval of the psychometric fits in both conditions. Because all modeled values were within the measurement error, our results suggest that sensitivity to beam size and stimulus geometry changes—as predicted by the light propagation model—can be measured at the psychophysical level.

Fig. 8 Effect of beam and stimulus size on light capture. (A) Retinal site tested for increment thresholds with 7.92 μm stimulus outlined over targeted cone at 3° eccentricity in Subject 1. (B) Psychometric data and fits for two different beam sizes. Threshold is defined at the power level yielding a 0.5 probability of detection (dashed lines). Data points are scaled to the number of successfully delivered trials at each intensity (out of 20). (C) Map of predicted change in threshold ratio from a 5.85 mm pupil for a square stimulus. Superimposed on the map are the experiment conditions illustrated as red→green dot pairs, representing a shift from large to small beam sizes. The top red-green pair is for the cone tested in panels A-B. (D) Data as in C, but with a circular stimulus. Color coded threshold ratios are interpolated (bicubic) from data computed at intervals of 0.2 mm for pupil size and 0.9 μm for stimulus size. Stimulus size is the edge length of a square or diameter of a circle. Contours have steps of 0.01 normalized units.

Download Full Size | PDF

Table 2. Measured vs modeled results for change in threshold from a 5.85 mm pupil

View Table | View all tables in this article

The threshold ratio maps provide useful information about expected detection threshold values for a range of conditions. Minima (shaded in blue in Fig. 8(C), 8(D)) are found with a pupil size of 1.4 mm and a stimulus size of 7.8 μm for the square shape, and a pupil size of 1.4 mm and a stimulus size of 8.6 μm for the round shape. This corresponds to instances where stimulus light would be strongly absorbed in the surrounding cones relative to the central. This occurs more readily for square shapes because the corners of the stimuli are more likely to reach into neighboring cones, as expected. The minima are located at a stimulus size similar to the cone spacing and where the beam size is approximately half the optimum size for the center cone. Although each individual cone may absorb less light than the central photoreceptor, the cumulative effect is large. The minimum is slightly lower using the circular shape than the square because such a stimulus is naturally closer in shape to the cones themselves and less light will be captured by surrounding cones. A second region of lowered threshold is seen with a point stimulus (0 μm). This corresponds to the beam size that best couples light to the outer segment of the central cone, and is equivalent to the maxima shown in Fig. 5.

In the second experiment, we applied the model to data measured on the same AOSLO, where stimuli were directed either to the center of a cone or to a gap (Fig. 9(A) [24]). The difference in light capture between these two conditions was then examined with the model. Figure 9(B) shows the amount of light absorbed for the cone condition, with varying eccentricity and defocus. Peak absorption occurs with negative defocus (with respect to the ELM) and a small cone size, while less absorption occurs with positive defocus. Figure 9(C) is a similar plot for the gap condition. Here, however, the minimum occurs with negative defocus. Calculating the ratio of the absorption of light in the gap versus cone conditions yields a threshold ratio map (Fig. 9(D)). The model shows that this ratio varies considerably with respect to defocus, especially at large eccentricities, caused by the opposing asymmetries about zero defocus of the two conditions. By using least squares, the best fitting defocus to the gap-cone threshold ratios reported by Harmening et al. (2014) occurs with a defocus of −0.026 D (Fig. 9(E)), which is close to the defocus of −0.06 D estimated from the simpler model in that study. Although certain factors could account for this small difference in computed threshold ratios (see Discussion), the reasonable agreement suggests that the light propagation model accounts for the basic properties of absorption with adaptive optics microstimulation.

Fig. 9 Gap/cone threshold ratio model of light absorption. (A) Schematic of stimulus conditions from [24], with light either targeted on a cone or in a gap. (B) Modeled fraction of light absorbed when a 3 × 3 pixel stimulus is centered on a cone at varying eccentricities and defocus. (C) As in B, when the same stimulus is placed in a gap between cones. (D) Gap/cone threshold ratio as a function of eccentricity and defocus. (E) Comparison of modeled versus observed ratios (circles from Fig. 5(C) in [24], red line shows linear fit). Threshold ratios for different defocus values (gray lines) showed that the best fit occurred at −0.026 D.

Download Full Size | PDF

4. Discussion

The finite difference beam propagation method can be used to model waveguiding of light down a photoreceptor, a behavior that would be expected when the intracellular contents have a higher refractive index than the extracellular environment. The parabolic shape of the inner segment serves to funnel light from a larger area at the entrance of the inner segment to a smaller area at the start of the outer segment, increasing the probability of light capture. As shown in Section 3.1, the model reveals modal behavior, but yields more detail, as it includes non-bound light. For a given cone size, one can compute an optimum beam size that predicts the most efficient delivery of light into the outer segment. The existence of an optimal beam size is also predicted from traditional modal analysis of the photoreceptors, where the most efficient coupling between input light and modes is found when they are closely matched in shape [7]. Although an exhaustive test of this prediction is beyond the scope of the paper, our psychophysical results in Section 3.3 suggest the method has some validity. While we do not fully investigate the effects of angularly displaced light here, the example propagations shown in Fig. 3(C) and 3(G) (normal and angled at 7.5°) suggest a Stiles-Crawford parameter of ρ = 0.057 mm⁻². This is in line with measured estimates for 543 nm light (for example, see Fig. 1 of [47]).

In none of the model runs was polarization found to be significant. However, this is largely due to the fact that photoreceptors were isotropic about the propagation axis, i.e. circular in cross-section. Polarization is likely to become important where anisotropies arise. For example, if the inner segments have an elliptical cross-section, propagation would become polarization dependent [48]. Given the subtle irregularities seen in cone photo-receptors (for example, see Fig. 2 of Curcio et al. [35]) this may induce small differences in the efficiency of light capture with respect to polarization for individual photoreceptors, but not in bulk due to randomness in their orientations. Our model could be used to explore this phenomenon.

Figures 5 and 6 emphasize the importance of incorporating the surrounding cones in any model of single cone light capture. Even when the stimulus size is not much larger than a cone, significant amounts of light may fall on adjacent receptors and contribute to a physiological response. This is highlighted further by looking at beam sizes of less than 1.5 mm, depending on the inner segment diameter. Here a focused spot is larger than the area covered by 7 cones. For these conditions, an additional ring of 12 photoreceptors should have been included if one aimed to estimate a perceptual threshold. Had these additional photoreceptors been added to the model, the average light absorbed with a small beam at the eye’s pupil would have increased. With adaptive optics microstimulation, consideration of such an expanse of cones is not necessary because beam sizes of 6 mm are typical [1,24,25,49,50].

In order to determine how well light delivery can be confined to a single cone we estimated what parameter values would be most effective (Fig. 7). Unsurprisingly, for optimal stimulation, a small stimulus and focused spot from a large beam are required, but the model provides detail about what parameters should be used for a given condition. Even under conditions where the model shows significant amounts of light falling on the surrounding cones, insight may be provided so that results may be properly interpreted. Figure 7 shows that the optimum beam size for efficient stimulus delivery is dependent on the cone size. We chose a point stimulus here for illustrative purposes, and the optimum beam size will be dependent on factors including wavelength and stimulus size. A system with a large beam and small focal spot may not always deliver the most light to the targeted cone. It is interesting that the most efficient stimulus delivery is not linear with beam size. For d₀ up to about 6 μm, the optimum beam size decreases with increased cone size because the shape of the field departs from the optimal for guiding light into the outer segment. This is consistent with other modal coupling models [7,12]. For cone sizes greater than ~6 μm, the optimum beam size increases and then plateaus as waveguiding becomes less of a factor due to the increasing diameter of the inner segment and light simply couples more directly to the outer segment. This result implies that the exact value of optimal beam size is dependent on both the inner and outer segment diameters. We also show that the most efficient delivery is not the only case where stimulus delivery is largely confined to a single cone. For instance, in Figs. 5 and 6 where the most efficient stimulus occurs for a 6 μm cone with a 4.5 mm beam size, one can still achieve more than 90% of the captured light absorbed by the targeted cone with beam sizes > 3 mm and a stimulus size of ~6 μm. All told, it seems clear that a number of optical parameters will effect light capture and therefore lead to variation in the functional weighting of cones in any retinal circuit.

The model is clearly limited by what is known of the properties of the cones. As noted in the Methods, cone dimensions at the resolution of our model is not available in the literature at the 0-5° eccentricities investigated here. In particular, the variation of cone outer segment diameter is not known to the detail required. Not only would a larger or smaller diameter change the area in which light enters the outer segment, Eq. (31) shows that it will also change the shape, and thus propagation in the inner segment. While the variation may only be small [51], there may be differences in light absorption that may not be apparent in our modeled results. Also, there is also too little detail about refractive indices in cones. No work has been published quantifying each segment’s refractive index as a function of wavelength. For example, the refractive indices we used were measured using visible phase contrast microscopy with broadband visible light [52]. This hinders detailed investigation into the effects of wavelength on the coupling of light to photoreceptors, or the effect of dispersion on broadband stimuli. As these all factors will alter the propagation of light, more detail would be desirable for accurate modeling.

While the model shows nearly all the stimulus light may be absorbed by the targeted cone, this is idealized for a diffraction-limited system. It is worth emphasizing that, in practice, single cone stimulation is compromised further due to several factors. Errors in the real-time optical correction may broaden the PSF, and interocular scattering cannot be mitigated [53,54]. In addition, the precise focal plane in the retina is not known with adaptive optics imaging; this is generally selected by choosing the focal plane with the sharpest image, which is not necessarily at the ELM as assumed here. All these factors could be included in the model in the step where the initial electric field is calculated.

Our model does not consider reflections at the boundaries between refractive index changes, but these are not expected to impact its predictions about light capture significantly. The reflectance of the whole retina, inner limiting membrane to the posterior sclera, has been estimated to be <1%, dependent on wavelength [55]. If we apply the Fresnel equations to the boundaries in the model, assuming a normally incident, unpolarized beam, the largest reflection would be from the boundary between myoid and ellipsoid. This would have a reflectivity of 1.48 × 10⁻⁴, with Δn = 0.0335. The reflection from the inner/outer segment junction is much less, having a Δn = 0.004, giving a reflectivity of 2.06 × 10⁻⁶. Similarly, the model does not consider intracellular scattering, the largest cause of which would probably be the ellipsoidal mitochondria. These mitochondria are prolate spheroids with respect to the photoreceptor axis, and of a small size, giving conditions where strong forward scattering is to be expected [56]. While scattering by these structures may not be zero, there may only be a small effect on the propagation.

Despite these limitations, the FDBP model predictions about light capture are in accord with experimental data. In the case of the cone-gap experiment, defocus had to be included to get a reasonable fit. The amount of defocus required was small: −0.026 D is the equivalent of the focus being shifted 7.4 μm into the inner segment from the ELM. Our model revealed an asymmetry in the defocus, such that a negative shift in defocus from the ELM is not equivalent to a positive shift with equal magnitude. This underscores the importance of including direction of the light when modeling. Significant differences in the propagation can be found depending on whether the beam is converging or diverging as it is incident on the ELM. At larger eccentricities in the cone-gap modeling, it was found that the effect of defocus is opposite for the two conditions, a certain defocus will increase the amount of light collected in one, while it is diminished in the other. This makes such an experiment very sensitive to defocus. This is contrary to using a weighted Gaussian acceptance function as a model of light capture [57], like that used in Harmening et al. (2014). As this model is intensity based, there is a loss of directional information, so it may only give accurate results for stimulation with a flat wavefront, i.e. one with no defocus or aberration. The difference in defocus estimate between the two studies could arise from the fact that the focal plane relative to the ELM is unknown in such experiments, as noted above. For the beam-stimulus size experiment, no defocus was needed to produce a good fit. At 3° eccentricity there is a change in power absorbed with varying defocus, but the model predicts that the change would be similar for both the large and small beam sizes, so this experiment would be less sensitive to defocus than the cone-gap experiment. This underscores how the effects of defocus are conditional on the experiment being performed.

As the model makes good predictions for adaptive optics threshold measurements without the inclusion of scattering or reflection, it adds confidence to the extant literature that the waveguiding nature of cones is an inescapable feature for the proper characterization of cones. The Gaussian acceptance function used by others is derived from measurement, but is attributed to the fundamental waveguide mode [57], so it could be considered an incomplete waveguide model. Vohnsen [36] presented a model of rod photoreceptors that uses no waveguiding, only scattering of light by the layered pigment discs in the outer segment to explain phenomena such as the Stiles-Crawford effect. Outer segment disc scattering could be included in our model by treating the layered photosensitive discs as a fiber-Bragg grating and using a bi-directional technique to include its reflection [58]. However, it is worth noting that the high periodicity of the discs (spacing of ~15 nm [59]) suggests that any reflectivity will peak at the low end of the UV spectrum, and be weak in the visible region. Inclusion of this extra detail will not likely affect the results produced by the our model significantly, so it is probably not worth the cost of the much greater computing time that an iterative bi-directional technique with high enough resolution would require, though the problem could be reduced to 2 dimensions. Another non-waveguide model has been proposed that uses the overlap of the volume of outer segments and the incident beam [60]. This volumetric model considers the effect of multiple cells, but neglects the inner segment and any guiding effects it may have. A non-waveguiding model may be appealing in its simplicity and thus fast implementation, and may give useful results for some situations, but without waveguiding such a model will never fully describe all the interaction between light and photoreceptors.

The finite difference beam propagation method provides a useful tool to model photoreceptors and is not limited to the conditions tested here. We only modeled cone photoreceptors, but the techniques used could be applied to rods if appropriate values are used. The size and shape of the photoreceptors and initial electric field can be arbitrarily defined within the limits of the technique’s assumptions, so the method could be applied to many other conditions. For instance, the effect on light capture of a particular optical aberration could be explored by specifying a suitable ψ₀, or the impact of misshapen photoreceptors such as those occurring in retinal disease could be modeled in great detail. The model enables the investigation of many such features to study how cone signals are altered and what their impact on visual function might be.

Funding

NIH (R01-EY023581, R01-EY023591, P30-EY003039); Eyesight Foundation of Alabama.

Acknowledgements

We thank A. S. McKeown for help with psychophysical data collection and analysis.

Disclosures

The authors declare there are no conflicts of interest related to this article.

References

1. L. C. Sincich, Y. Zhang, P. Tiruveedhula, J. C. Horton, and A. Roorda, “Resolving single cone inputs to visual receptive fields,” Nat. Neurosci. 12(8), 967–969 (2009). [CrossRef] [PubMed]

2. G. D. Field, J. L. Gauthier, A. Sher, M. Greschner, T. A. Machado, L. H. Jepson, J. Shlens, D. E. Gunning, K. Mathieson, W. Dabrowski, L. Paninski, A. M. Litke, and E. J. Chichilnisky, “Functional connectivity in the retina at the resolution of photoreceptors,” Nature 467(7316), 673–677 (2010). [CrossRef] [PubMed]

3. P. H. Li, G. D. Field, M. Greschner, D. Ahn, D. E. Gunning, K. Mathieson, A. Sher, A. M. Litke, and E. J. Chichilnisky, “Retinal representation of the elementary visual signal,” Neuron 81(1), 130–139 (2014). [CrossRef] [PubMed]

4. J. M. Enoch, “Wave-guide modes in retinal receptors,” Science 133(3461), 1353–1354 (1961). [CrossRef] [PubMed]

5. Z. Liu, O. P. Kocaoglu, T. L. Turner, and D. T. Miller, “Modal content of living human cone photoreceptors,” Biomed. Opt. Express 6(9), 3378–3404 (2015). [CrossRef] [PubMed]

6. A. W. Snyder and C. Pask, “The Stiles-Crawford effect-explanation and consequences,” Vision Res. 13(6), 1115–1137 (1973). [CrossRef] [PubMed]

7. B. Vohnsen, I. Iglesias, and P. Artal, “Guided light and diffraction model of human-eye photoreceptors,” J. Opt. Soc. Am. A 22(11), 2318–2328 (2005). [CrossRef] [PubMed]

8. J.-M. Gorrand and F. C. Delori, “A model for the assessment of cone directionality,” J. Mod. Opt. 44(3), 473–491 (1997). [CrossRef]

9. A. W. Snyder, “Coupling of modes on a tapered dielectric cylinder,” IEEE Trans. Microw. Theory Tech. 18(7), 383–392 (1970). [CrossRef]

10. A. W. Snyder and M. Hamer, “The light-capture area of a photoreceptor,” Vision Res. 12(10), 1749–1753 (1972). [CrossRef] [PubMed]

11. D. G. Stavenga, “Waveguide modes and refractive index in photoreceptors of invertebrates,” Vision Res. 15(3), 323–330 (1975). [CrossRef] [PubMed]

12. B. Vohnsen, “Photoreceptor waveguides and effective retinal image quality,” J. Opt. Soc. Am. A 24(3), 597–607 (2007). [CrossRef] [PubMed]

13. R. F. Spaide and C. A. Curcio, “Anatomical correlates to the bands seen in the outer retina by optical coherence tomography: literature review and model,” Retina 31(8), 1609–1619 (2011). [CrossRef] [PubMed]

14. R. Sammut and A. W. Snyder, “Contribution of unbound modes to light absorption in visual photoreceptors,” J. Opt. Soc. Am. 64(12), 1711–1714 (1974). [CrossRef] [PubMed]

15. J. A. Besley and J. D. Love, “Supermode analysis of fibre transmission,” IEE Proc., Optoelectron. 144(6), 411–419 (1997). [CrossRef]

16. K. Yee, “Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media,” IEEE Trans. Antenn. Propag. 14(3), 302–307 (1966). [CrossRef]

17. A. M. Pozo, F. Pérez-Ocón, and J. R. Jiménez, “FDTD analysis of the light propagation in the cones of the human retina: an approach to the Stiles-Crawford effect of the first kind,” J. Opt. A, Pure Appl. Opt. 7(8), 357–363 (2005). [CrossRef]

18. M. J. Piket-May, A. Taflove, and J. B. Troy, “Electrodynamics of visible-light interactions with the vertebrate retinal rod,” Opt. Lett. 18(8), 568–570 (1993). [CrossRef] [PubMed]

19. L. Fischer, A. Zvyagin, T. Plakhotnik, and M. Vorobyev, “Numerical modeling of light propagation in a hexagonal array of dielectric cylinders,” J. Opt. Soc. Am. A 27(4), 865–872 (2010). [CrossRef] [PubMed]

20. M. B. Toomey, P. Olsson, R. Frederiksen, M. C. Cornwall, R. Oulton, A. Kelber, J. C. Corbo, N. W. Roberts, and D. Wilby, “Optics of cone photoreceptors in the chicken (Gallus gallus domesticus),” J. R. Soc. Interface 12(111), 20150591 (2015). [CrossRef] [PubMed]

21. W. P. Huang, C. L. Xu, S. T. Chu, and S. K. Chaudhuri, “A vector beam propagation method for guided-wave optics,” IEEE Photonics Technol. Lett. 3(10), 910–913 (1991). [CrossRef]

22. W. P. Huang, C. L. Xu, and S. K. Chaudhuri, “A finite-difference vector beam propagation method for three-dimensional waveguide structures,” IEEE Photonics Technol. Lett. 4(2), 148–151 (1992). [CrossRef]

23. W. P. Huang and C. L. Xu, “Simulation of three-dimensional optical waveguides by a full-vector beam propagation method,” IEEE J. Quantum Electron. 29(10), 2639–2649 (1993). [CrossRef]

24. W. M. Harmening, W. S. Tuten, A. Roorda, and L. C. Sincich, “Mapping the perceptual grain of the human retina,” J. Neurosci. 34(16), 5667–5677 (2014). [CrossRef] [PubMed]

25. K. S. Bruce, W. M. Harmening, B. R. Langston, W. S. Tuten, A. Roorda, and L. C. Sincich, “Normal perceptual sensitivity arising from weakly reflective cone photoreceptors,” Invest. Ophthalmol. Vis. Sci. 56(8), 4431–4438 (2015). [CrossRef] [PubMed]

26. G. L. Pedrola, Beam Propagation Method for Design of Optical Waveguide Devices (John Wiley & Sons, 2016).

27. J. Xiao, X. Liu, C. Cai, H. Fan, and X. Sun, “An improved three-dimensional full-vectorial finite-difference imaginary-distance beam propagation method,” Sci. China, Ser. F Inf. Sci. 49(4), 516–532 (2006).

28. B. O’Brien, “A theory of the Stiles and Crawford effect,” J. Opt. Soc. Am. 36(9), 506–509 (1946). [CrossRef] [PubMed]

29. R. Winston, “The Visual Receptor as a Light Collector,” in Vertebrate Photoreceptor Optics, J. M. Enoch and F. L. Tobey, eds. (Springer Verlag, 1981), pp. 325–336.

30. T. Zhang, P. Godara, E. R. Blanco, R. L. Griffin, X. Wang, C. A. Curcio, and Y. Zhang, “Variability in human cone topography assessed by adaptive optics scanning laser ophthalmoscopy,” Am. J. Ophthalmol. 160(2), 290–300 (2015). [CrossRef] [PubMed]

31. C. A. Curcio, K. R. Sloan Jr., O. Packer, A. E. Hendrickson, and R. E. Kalina, “Distribution of cones in human and monkey retina: individual variability and radial asymmetry,” Science 236(4801), 579–582 (1987). [CrossRef] [PubMed]

32. L. Sawides, A. de Castro, and S. A. Burns, “The organization of the cone photoreceptor mosaic measured in the living human retina,” Vision Res. 132, 34–44 (2017). [CrossRef] [PubMed]

33. H. Song, T. Y. P. Chui, Z. Zhong, A. E. Elsner, and S. A. Burns, “Variation of cone photoreceptor packing density with retinal eccentricity and age,” Invest. Ophthalmol. Vis. Sci. 52(10), 7376–7384 (2011). [CrossRef] [PubMed]

34. W. H. Miller and G. D. Bernard, “Averaging over the foveal receptor aperture curtails aliasing,” Vision Res. 23(12), 1365–1369 (1983). [CrossRef] [PubMed]

35. C. A. Curcio, K. R. Sloan, R. E. Kalina, and A. E. Hendrickson, “Human photoreceptor topography,” J. Comp. Neurol. 292(4), 497–523 (1990). [CrossRef] [PubMed]

36. B. Vohnsen, “Directional sensitivity of the retina: A layered scattering model of outer-segment photoreceptor pigments,” Biomed. Opt. Express 5(5), 1569–1587 (2014). [CrossRef] [PubMed]

37. B. B. Boycott, J. E. Dowling, and H. Kolb, “Organization of the primate retina: light microscopy,” Proc. R. Soc. London, Ser. B, Biol. Sci. 255(799), 109– 184 (1969).

38. V. J. Srinivasan, B. K. Monson, M. Wojtkowski, R. A. Bilonick, I. Gorczynska, R. Chen, J. S. Duker, J. S. Schuman, and J. G. Fujimoto, “Characterization of outer retinal morphology with high-speed, ultrahigh-resolution optical coherence tomography,” Invest. Ophthalmol. Vis. Sci. 49(4), 1571–1579 (2008). [CrossRef] [PubMed]

39. J. Hirsch and R. Hylton, “Quality of the primate photoreceptor lattice,” Vision Res. 24(4), 347–355 (1984). [CrossRef] [PubMed]

40. R. Barer, “Refractometry and interferometry of living cells,” J. Opt. Soc. Am. 47(6), 545–556 (1957). [CrossRef] [PubMed]

41. F. I. Hárosi, “Application of Fourier transform smoothing and statistical techniques to the determination of spectral parameters,” J. Gen. Physiol. 89(5), 717–743 (1987). [CrossRef] [PubMed]

42. D. A. Baylor, B. J. Nunn, and J. L. Schnapf, “The photocurrent, noise and spectral sensitivity of rods of the monkey Macaca fascicularis,” J. Physiol. 357(1), 575–607 (1984). [CrossRef] [PubMed]

43. W. M. Harmening, P. Tiruveedhula, A. Roorda, and L. C. Sincich, “Measurement and correction of transverse chromatic offsets for multi-wavelength retinal microscopy in the living eye,” Biomed. Opt. Express 3(9), 2066–2077 (2012). [CrossRef] [PubMed]

44. D. W. Arathorn, Q. Yang, C. R. Vogel, Y. Zhang, P. Tiruveedhula, and A. Roorda, “Retinally stabilized cone-targeted stimulus delivery,” Opt. Express 15(21), 13731–13744 (2007). [CrossRef] [PubMed]

45. Q. Yang, D. W. Arathorn, P. Tiruveedhula, C. R. Vogel, and A. Roorda, “Design of an integrated hardware interface for AOSLO image capture and cone-targeted stimulus delivery,” Opt. Express 18(17), 3283–3291 (2010).

46. D. Gloge, “Weakly guiding fibers,” Appl. Opt. 10(10), 2252–2258 (1971). [CrossRef] [PubMed]

47. T. J. Berendschot, J. van de Kraats, and D. van Norren, “Wavelength dependence of the Stiles-Crawford effect explained by perception of backscattered light from the choroid,” J. Opt. Soc. Am. A 18(7), 1445–1451 (2001). [CrossRef] [PubMed]

48. A. W. Snyder and W. R. Young, “Modes of optical waveguides,” J. Opt. Soc. Am. 68(3), 297–309 (1978). [CrossRef]

49. W. S. Tuten, P. Tiruveedhula, and A. Roorda, “Adaptive optics scanning laser ophthalmoscope-based microperimetry,” Optom. Vis. Sci. 89(5), 563–574 (2012). [CrossRef] [PubMed]

50. Q. Wang, W. S. Tuten, B. J. Lujan, J. Holland, P. S. Bernstein, S. D. Schwartz, J. L. Duncan, and A. Roorda, “Adaptive optics microperimetry and OCT images show preserved function and recovery of cone visibility in macular telangiectasia type 2 retinal lesions,” Invest. Ophthalmol. Vis. Sci. 56(2), 778–786 (2015). [CrossRef] [PubMed]

51. Q. V. Hoang, R. A. Linsenmeier, C. K. Chung, and C. A. Curcio, “Photoreceptor inner segments in monkey and human retina: mitochondrial density, optics, and regional variation,” Vis. Neurosci. 19(4), 395–407 (2002). [CrossRef] [PubMed]

52. R. L. Sidman, “The structure and concentration of solids in photoreceptor cells studied by refractometry and interference microscopy,” J. Biophys. Biochem. Cytol. 3(1), 15–30 (1957). [CrossRef] [PubMed]

53. T. J. T. P. van den Berg, L. Franssen, B. Kruijt, and J. E. Coppens, “History of ocular straylight measurement: A review,” Z. Med. Phys. 23(1), 6–20 (2013). [CrossRef] [PubMed]

54. L. Franssen, J. Tabernero, J. E. Coppens, and T. J. T. P. van den Berg, “Pupil size and retinal straylight in the normal eye,” Invest. Ophthalmol. Vis. Sci. 48(5), 2375–2382 (2007). [CrossRef] [PubMed]

55. T. T. J. M. Berendschot, P. J. DeLint, and D. van Norren, “Fundus reflectance–historical and present ideas,” Prog. Retin. Eye Res. 22(2), 171–200 (2003). [CrossRef] [PubMed]

56. S. Asano and G. Yamamoto, “Light scattering by a spheroidal particle,” Appl. Opt. 14(1), 29–49 (1975). [CrossRef] [PubMed]

57. D. I. A. MacLeod, D. R. Williams, and W. Makous, “A visual nonlinearity fed by single cones,” Vision Res. 32(2), 347–363 (1992). [CrossRef] [PubMed]

58. H. Rao, R. Scarmozzino, and R. M. Osgood, “Bidirectional beam propagation method for multiple dielectric interfaces,” IEEE Photonics Technol. Lett. 11(7), 830–832 (1999). [CrossRef]

59. R. H. Steinberg, S. K. Fisher, and D. H. Anderson, “Disc morphogenesis in vertebrate photoreceptors,” J. Comp. Neurol. 190(3), 501–518 (1980). [CrossRef] [PubMed]

60. B. Vohnsen, A. Carmichael, N. Sharmin, S. Qaysi, and D. Valente, “Volumetric integration model of the Stiles-Crawford effect of the first kind and its experimental verification,” J. Vis. 17(12), 18 (2017). [CrossRef] [PubMed]

Condition (Figure)	Absorbed power (normalized to input)				Ratio of center cone to total
Condition (Figure)	Top cone	Center cone	Bottom cone	Total (all cones)	Ratio of center cone to total
Centered (3C)	0.0004	0.1924	0.0004	0.1947	0.99
Shifted (3D)	0.0505	0.0600	0.0016	0.1145	0.52
Small beam (3E)	0.0052	0.1371	0.0052	0.1686	0.81
Large cone (3F)	0.0003	0.1418	0.0003	0.1435	0.99
Angled (3G)	0.0038	0.0982	0.0028	0.1070	0.92
Coma (3H)	0.0038	0.0893	0.0315	0.1314	0.68

Subject		Modeled threshold change	Measured threshold change
1	6.48 circular	0.73	0.74 ( ± 0.07)
2	6.48 circular	0.73	0.85 ( ± 0.14)
2	3.60 square	1.11	1.08 ( ± 0.11)
1	7.92 square	0.68	0.76 ( ± 0.08)

Condition (Figure)	Absorbed power (normalized to input)				Ratio of center cone to total
Condition (Figure)	Top cone	Center cone	Bottom cone	Total (all cones)	Ratio of center cone to total
Centered (3C)	0.0004	0.1924	0.0004	0.1947	0.99
Shifted (3D)	0.0505	0.0600	0.0016	0.1145	0.52
Small beam (3E)	0.0052	0.1371	0.0052	0.1686	0.81
Large cone (3F)	0.0003	0.1418	0.0003	0.1435	0.99
Angled (3G)	0.0038	0.0982	0.0028	0.1070	0.92
Coma (3H)	0.0038	0.0893	0.0315	0.1314	0.68

Subject		Modeled threshold change	Measured threshold change
1	6.48 circular	0.73	0.74 ( ± 0.07)
2	6.48 circular	0.73	0.85 ( ± 0.14)
2	3.60 square	1.11	1.08 ( ± 0.11)
1	7.92 square	0.68	0.76 ( ± 0.08)

Light propagation and capture in cone photoreceptors

Abstract

1. Introduction