Modeling Photo-Bleaching Kinetics to Create High Resolution Maps of Rod Rhodopsin in the Human Retina

We introduce and describe a novel non-invasive in-vivo method for mapping local rod rhodopsin distribution in the human retina over a 30-degree field. Our approach is based on analyzing the brightening of detected lipofuscin autofluorescence within small pixel clusters in registered imaging sequences taken with a commercial 488nm confocal scanning laser ophthalmoscope (cSLO) over a 1 minute period. We modeled the kinetics of rhodopsin bleaching by applying variational optimization techniques from applied mathematics. The physical model and the numerical analysis with its implementation are outlined in detail. This new technique enables the creation of spatial maps of the retinal rhodopsin and retinal pigment epithelium (RPE) bisretinoid distribution with an ≈ 50μm resolution.


Introduction
Currently, many retinal diseases are clinically evaluated by subjective examination of retinal images [1] and lower resolution visual field testing [2]. Improved quantitative analysis in clinics should increase understanding of mechanisms of early retinal disease progression by mapping local changes in bleachable rhodopsin within the rod outer segments and bisretinoid levels within the underlying RPE. Localized reduction in bleachable rhodopsin may be an early indicator of local visual cycle defects and future rod loss, whereas local increases in bisretinoid fluorescence within the RPE has been associated with RPE stress and early sign of loss of function [3][4][5][6]. Thus, reliable clinical noninvasive maps of bleachable rhodopsin distribution in conjunction with RPE autofluorescence may provide a means to assess early intervention and disease prevention strategies. In the present work we develop the methodology to quantify rod rhodopsin using standard clinical cSLO autofluorescence imaging and analysis enabled by state-of-the-art mathematical techniques. The paper is addressed to both, applied mathematicians and vision scientists. The latter can familiarize themselves with computational details that can significantly improve quantitative results. Applied mathematicians may better comprehend the biomedical aspects enabling future design of customized mathematical tools beyond rhodopsin measurements. Our efforts aim at improving the understanding between both communities to enable further synergies. Acronyms that may not be well-known to one or the other community are explained in Table 1. Readers who prefer to avoid the mathematical details can skip Sections 3 and 4.
The retina is a multi-layer neural tissue, uniquely suited for noninvasive optical imaging due to the evolutionary design of the eye and its ocular media. Located just above the choroid, the RPE is a single cell layer that nourishes overlying photoreceptors. Noninvasive autofluorescence imaging of bisretinoids in the RPE and of absorbing molecules in the overlying retina offers the possibility to sensitively monitor early changes in retinal function and early pathophysiology.
Localized rod photoreceptor losses have been observed in post mortem histology of agerelated maculopathy, and the spatial rod distribution during normal aging has been characterized in humans [7][8][9][10]. Rod rhodopsin bleaching has been widely observed [5,[11][12][13][14][15][16][17][18], and degraded dark adaptation during aging and disease progression has been attributed to reduction in visual cycle regeneration of cis-retinal and rod rhodopsin. Thus, reduced bleachable rhodopsin and pigmentary changes may locally reflect early RPE dysfunction [19][20][21][22][23][24][25]. Although the biophysical model of dark adaptation and rod rhodopsin bleaching was introduced in [5,11,15,18,26], attempts of rod rhodopsin quantification in humans from retinal bleaching using standard clinical instruments has only recently been reported [27,28] and not yet reduced to a routine clinical method. Bleaching kinetics are dependent on the actual retinal irradiance of the rhodopsin within the rod outer segments and therefore slowed by overlying retinal chromophores, such as macular pigments and hemoglobin, and by lens pigments. Lens absorption of blue laser light uniformly reduces retinal irradiance over the entire field of view in the cSLO. On the other hand, hemoglobin within retinal vessels effectively masks both the bleaching of underlying rhodopsin and RPE autofluorescence, which prevents analysis of pixels beneath visible retinal vessels. The macular pigments, lutein and the related carotenoid zeaxanthin, are concentrated within the photoreceptor nerve fibers of the fovea and reduce irradiance at 488nm at both the rod outer segment and RPE levels to an increasing degree as one approaches the center of the fovea [29][30][31]. Therefore, current bleaching and regeneration models must be extended to incorporate macular pigment contributions, and in Section 2 we shall combine the underlying biophysical models. Melanin within the RPE is more uniformly distributed [32] and is behind the photoreceptors, so that it reduces only the irradiance at the level of RPE lipofuscin and not the rod rhodopsin. Consequently, variations of melanin in the RPE and choroid within the central macula only reduce the autofluorescence amplitudes and not the rate of rhodopsin bleaching. To integrate measured autofluorescence images into the developed model, we shall use state-of-the-art variational analysis techniques from applied mathematics. The rhodopsin distribution maps are computed by optimization procedures outlined in Section 3. As a postprocessing step described in Section 4, we detect retinal vessels through image analysis and refine mathematical image inpainting methods to derive a spatial rod rhodopsin map of the human retina, in which retinal vessels are removed. Section 5 is a brief summary of our numerical approach studied in the Sections 3 and 4. We present examples and compare our rhodopsin maps with the typical distribution of rod photoreceptors in Section 6. Results are discussed in Section 7 and conclusion are given in Section 8.

Biophysical models
This section is dedicated to present a biophysical foundation of rod rhodopsin quantification from standard imaging devices. After highlighting the use of standard clinical instruments to measure retinal autofluorescence, we shall recall the rhodopsin bleaching and regeneration model in [5], refine the macular pigment model from [33], and combine both models. We shall also describe the underlying clinical protocol.

Retinal autofluorescence imaging
As light penetrates the retina, it is largely unscattered and only locally absorbed by retinal chromophores (first, hemoglobin within large retinal vessels, then macular pigments largely in the photoreceptor axons, then the unbleached opsins within the photoreceptor outer segments) before reaching the lipofuscin granules in the RPE. The emitted fluorescence from unoxidized lipofuscin bisretinoids can be excited by a range of visible wavelengths (blue to yellow) to yield long-wavelength emission in the red, which is detected noninvasively by the fluorescence imaging camera. Let Λ and λ be the excitation and emission wavelengths, respectively, and AF(Λ, λ) denotes the measured autofluorescence. The Beer-Lambert law for the double-path yields where D(Λ) and D(λ) are the optical densities of the underlying tissue at wavelengths Λ and λ, respectively, F(Λ, λ) is the fluorescence efficiency of lipofuscin, and I(Λ) is the radiant power of the excitation light. We use two different imaging techniques: First, to measure macular pigment, we developed autofluorescence imaging of the human retina at varying emission and excitation wavelengths by modifying standard fundus cameras [34]. Secondly, to subsequently quantify rod rhodopsin, we record images at a commercial cSLO camera (Heidelberg Retinal Angiograph 4.0) that delivers an average of % 3μW/mm 2 at 488nm by rapidly scanning a small laser beam (10μm) over the retina. The intensity of excitation in the cSLO is > 100-fold less than in the fundus camera so that each cSLO image is acquired over a 100ms scan with an incremental rhodopsin bleaching (% 1%). After % 40s of cSLO imaging, the rhodopsin is completely bleached (> 98%), and its initial attenuation of the excitation light at 488nm is removed. To map the rod rhodopsin from the magnitude of the brightening of the autofluorescence in registered movies, we shall study the quantitative models of bleaching and macular pigment absorption in more detail.

Rod rhodopsin models
2.2.1 Bleaching model. In a dark-adapted retina the rod opsin is predominantly in its bound form with 11-cis retinal, i.e., rhodopsin. Rhodopsin strongly absorbs green-blue (% 500nm) light, which photo-isomerizes its bond cis-retinal to all-trans-retinal inducing protein conformational changes and subsequent photoreception signaling. These various downstream forms lose the strong green-blue absorbance, hence rhodopsin bleaching refers to its conversion by light to these non-absorbing forms [5,15,18].
If we start with a dark-adapted retina, then the intensity of the local autofluorescence, AF, within a small circle of pixels should increase as the rhodopsin bleaches (becomes transparent to the 488nm exciting laser light). The unbleached rhodopsin in the photoreceptor layer initially attenuates the lipofuscin excitation by % 50% outside the central fovea, and the excitation light reaching the RPE then progressively increases as the overlying rhodopsin is locally bleached. Thus, we need to incorporate the time course in Eq (1), so that we have AFðL; l; tÞ ¼ IðLÞFðL; lÞe ÀðDðL;tÞþDðl;tÞÞ ; where we assume a constant radiation power I(Λ). The optical density changes over time, and contributing chromophores are macular pigments in the fovea and rhodopsin in the perifoveal region. There is, however, a region within the fovea, where both are present at a significant density. The macular pigment contribution was not accounted for in [5,15,18]. Here, we explicitly model the signal attenuation caused by macular pigment absorption to obtain rod rhodopsin maps covering the entire macula. Therefore, we obtain DðL; tÞ þ Dðl; tÞ ¼ D MP ðLÞ þ D MP ðlÞ þ D rh ðL; tÞ þ D rh ðl; tÞ; where D rh and D MP denote the optical density of present rhodopsin and macular pigments, respectively. If R(t) denotes the fraction of rhodopsin remaining at time t, then we obtain D rh ðL; tÞ þ D rh ðl; tÞ ¼ ðD rh ðL; 0Þ þ D rh ðl; 0ÞÞRðtÞ; where R(0) = 1 corresponds to the dark adapted retina and D rh (Λ,0)+D rh (λ, 0) is the double path optical density of rhodopsin present at time t = 0. If % rh (0) denotes the physical density of present rhodopsin at time t = 0 and its molar extinction coefficient is k rh , then we have D rh ðL; 0Þ þ D rh ðl; 0Þ ¼ % rh ð0Þðk rh ðLÞ þ k rh ðlÞÞ; which yields AFðL; l; tÞ ¼ IðLÞFðL; lÞe ÀðD MP ðLÞþD MP ðlÞþ% rh ð0Þðk rh ðLÞþk rh ðlÞÞRðtÞÞ : In order to determine % rh (0), we still need to specify R(t).
2.2.2 Regeneration model. We use the model developed in [5,15] to specify the fraction of rhodopsin remaining at time t due to the rates of rhodopsin photoactivation and of its regeneration by the visual cycle. The visual cycle is a complex process that transports the released all-trans-retinoid to the underlying retinal pigment epithelium, there converts it back to 11-cisretinal and returns it to the rod opsin. The normal time constant for visual cycle regeneration of rhodopsin increases with age from % 700sec to % 1000sec and is much greater than the rate of rhodopsin bleaching in autofluorescence imaging with 488nm cSLOs.
To specify R(t), the fraction of rhodopsin remaining at time t, we need to model the regeneration process, for which we started with Michaelis-Menten kinetics of the production of cis-retinal from cis-retinol. The change in the fraction of unbound opsin Ops(t) = 1 − R(t) is proportional to the 11-cis retinal concentration that binds to opsin, where S is 11-cis retinol, E is 11-cis RDH (RDH5), and P is 11-cis retinal. The fractional change of unbound opsin satisfies @ @t OpsðtÞ ¼ ÀkPðtÞOpsðtÞ: Thus, when rhodopsin is bleached by a steady light of illuminance I, its regeneration resembles a first-order reaction 1ÀRðtÞ K ; where K is the time constant. Therefore, the fractional change of the rhodopsin concentration satisfies where L is a "bleaching constant" that corresponds to the reciprocal of "photosensitivity" and has been measured by retinal densitometry [35,36]. Eq (3) has the analytical solution In cSLO measurements, the retinal illuminance I is much bigger than L K , and we evaluate R only at t much smaller than K. Therefore, Eq (4) reduces to RðtÞ % e À I L t ; a simplification that was already pointed out in [5,15].

Macular pigment model
To improve accuracy in rod rhodopsin mapping, we quantify macular pigments from multispectral autofluorescence image sets. In the past, we developed noninvasive multi-spectral autofluorescence imaging of the human retina by adding selected interference filter sets to standard fundus cameras, cf. [34,[37][38][39][40][41][42][43]. By exciting the fluorescent lipofuscin granules within the RPE, the Beer-Lambert model for the double path penetration enables us to effectively measure the spatial macular pigment distribution from a set of multi-spectral autofluorescence images. Quantitative measurements of macular pigments based on two-wavelength autofluorescence images have been introduced by Delori et al. in [33,44]. To more robustly analyze the spatial macular pigment distribution, we introduce a multiple-wavelength model that enables more effective self-consistency tests.
Let AF f (Λ, λ) and AF p (Λ, λ) be the autofluorescence measured at the fovea and the perifovea, respectively. While AF f depends on the specific location within the fovea, the term AF p is often replaced by a circular average at 6 degrees [33]. We denote the optical density of the foveal and perifoveal tissue by D f and D p , respectively. Let F f and F p be the fluorescence efficiencies of lipofuscin in the foveal and perifoveal regions. According to Eq (1), the foveal and perifoveal autofluorescence are given by AF f ðL; lÞ ¼ IðLÞF f ðL; lÞe ÀðD f ðLÞþD f ðlÞÞ ; AF p ðL; lÞ ¼ IðLÞF p ðL; lÞe ÀðD p ðLÞþD p ðlÞÞ : We can assume that the fluorophore at the fovea has the same composition as that at the perifovea (constant shape of its spectrum over fðL j ; l j Þg n j¼1 ), and that foveal-perifoveal differences in absorption by other pigments (retinal blood, visual pigments, and RPE melanin) are negligible. We also neglect the low macular pigment concentration in the perifoveal region, so that the optical density of macular pigment D MP at 460nm is the difference [33,44]. We use the relative extinction coefficient k MP of macular pigment, scaled to k MP (460nm) = 1, such that D MP (λ) = k MP (λ)D MP (460nm), and we obtain log AF p ðL; lÞ AF f ðL; lÞ We choose n pairs of excitation and emission wavelengths fðL j ; l j Þg n j¼1 and weights fo j g n j¼1 such that P n j¼1 o j ¼ 0: We apply the above equations for each wavelength pair (Λ j , λ j ), multiply by ω j , and add them up to obtain Since the foveal-perifoveal differences in fluorophore composition are negligible, the ratio F p ðL j ;l j Þ F f ðL j ;l j Þ does not depend on j. Therefore, we can determine D MP (460nm) by If we only choose two excitation wavelengths Λ 1 = 480nm and Λ 2 = 520nm with the weights 1 and −1, respectively, then we obtain the formula originally proposed in [33]. Our formula (5) enables several tests for self-consistency by removing or adding wavelength pairs and by changing the weights.

Combined model
To measure rod rhodopsin in the human retina, we first quantify the macular pigment density D MP (460) from multi-spectral autofluorescence image sets according to formula (5). The rod rhodopsin density can then be quantified from cSLO autofluorescence movies using the bleaching model. The macular pigment density is a parameter of the bleaching model whose exact specification increases the accuracy. By combining Eq (2) with the simplification of Eq (4), we derive the complete bleaching model Note that each flash of the fundus camera in macular pigment measurements fully bleaches rhodopsin so that our subsequent multi-spectral images are not affected by the bleaching kinetics and hence do not depend on the rod rhodopsin distribution.

Bleaching protocol
For recording rhodopsin bleaching cSLO movies, we follow a simple clinical protocol in which the human subject wears vermillion sunglasses (rod-protecting) while waiting and the photographer performs focus adjustment using the infrared reflection imaging (nonbleaching) in the cSLO. A 488nm excited autofluorescence movie (% 8 frames per sec, 1 minute long) is started that is recorded from the start with blinks every 10s, which refresh the tear film layer on the cornea. The average photon flux bleaches the rod rhodopsin after > 25s. We record the cSLO movie until steady-state rhodopsin bleaching.

Mathematical analysis of rhodopsin quantification
After clinical autofluorescence measurements are recorded supposedly satisfying Eq (6), we must solve for % rh (0) at each pixel location. This requires advanced mathematical tools outlined in detail next.

Numerical rhodopsin extraction
The recorded cSLO movie yields the autofluorescence AF(Λ, λ, t) over time in Eq (6). For notational convenience, we shall denote the model parameters by The extinction coefficients k MP and k rh are derived from the literature. We find the rhodopsin Mapping Rod Rhodopsin in the Human Retina density % rh (0) by numerically determining the parameters α, β, and γ in where T 0 % 0 is the time of initial exposure to light. Note that we have bleaching curves f(t), t 2 [T 0 , T], associated to AF(Λ, λ, t), for each pixel location (x, y). Computed parameters are optimized within the declared model, but may differ from the actual physiological parameters α, β, γ by a small error margin, or differ in certain areas due to physiological phenomena that are not present in the model. So, to avoid confusion, we will denote the numerically recovered parameters byâðx; yÞ % aðx; yÞ,bðx; yÞ % bðx; yÞ, andĉðx; yÞ % gðx; yÞ. Given a set of autofluorescence measurements f(t) satisfying the model for AF(Λ, λ, t), the numerical scheme approximates the bleaching process in the form of the double-exponential law in Eq (8). To avoid the use of double-exponents, we consider a "mean-square-log" deviation between f(t) and AF(Λ, λ, t), and numerically determine the three model parametersâ,b, andĉ that minimize E. We want to point out that this energy is different from the one used in our earlier preliminary work [38,39], where we directly enforced f(t) % ae −ce −bt . In the present approach Eq (9), we would like f ðtÞ ae Àce Àbt to be close to 1, which implies that ln f ðtÞ The remainder of the mathematical section is dedicated to develop and validate an iterative algorithm that determines the global minimizer of E in Eq (9). We shall see that the Hessian matrix of E is positive definite throughout, so that zeros of the first partial derivatives not only yield local minimizers but indeed the global minimizer. Therefore, the iterative scheme will indeed converge towards the global minimizer of Eq (9).

Implementation details
To find the actual minimizers of E, we set p = ln(a) and use a gradient descent algorithm. The latter algorithm can equivalently be derived from solving the Euler-Lagrange equations associated to Eq (9), with the use of an artificial time marching [45]. Indeed, the gradient descent of E can be formulated as a discretization of the following system where a, b and c are now supposed to be functions depending on some artificial time parameter τ. More precisely, the steady-states of Eq (11) satisfy Eq (10) and yield the minimizer of E in Eq (9), so that we can setâ cf. [45]. To actually compute the steady-states lim τ ! 1 a(τ), lim τ ! 1 b(τ), and lim τ ! 1 c(τ), we use a so-called semi-explicit finite difference scheme [46] that leads to ð ln ðf ðtÞÞ À p nþ1 þ c n e Àb n t Þdt; ð ln ðf ðtÞÞ À p n þ c n e Àb n t Þc n te Àb n t dt; ð ln ðf ðtÞÞ À p n þ c nþ1 e Àb n t Þe Àb n t dt; where τ p , τ b , and τ c are the individual step widths and initial values p 0 , b 0 , and c 0 that still need to be determined. As usual, the integral shall be replaced with the finite sum over the measured time points t 1 , . . ., t m , i.e., only ff ðt i Þg m i¼1 are available, and the index n+1 indicates the updated value of the solution at the next step, the index n corresponds to the "current" parameter values that were updated at the previous step. The second equation is explicit, i.e., we use b n on the right-hand side. The first and third equations can be implicit, i.e, the index n+1 is used on the right-hand side, because we are still able to solve the equations by ð ln ðf ðtÞÞ À p n þ c n e Àb n t Þte Àb n t dt; ð14Þ According to Eq (12), the minimizers of Eq (9) are computed from the iterations (13), (14), We shall verify the legitimacy of the above iterative scheme in Section 3.3. Before, we still need to discuss the choice of the initial values p 0 , b 0 , and c 0 for the respective model parameters. To determine their order of magnitude, we utilize a spatially averaged bleaching curve, so that the fitting is numerically stable, and we then use the resulting steady state solutions as initial values for consequent pixelwise processing. The fitting routine is, in fact, applied to the modified image stack, in which each pixel is averaged over its surrounding 8 × 8 pixels. The latter effectively reduces noise, including local distortions due to the eye movements and the mutual image registration.
We also see the importance to determine the macular pigment density a-priori because it influences the spatial distribution of the initial value p 0 = ln(a 0 ) that corresponds to α in Eq (7). In other words, the macular pigment measurements enable us to use more spatially accurate initial values of p 0 yielding faster convergence of the numerical algorithms, so that overall accuracy is improved.
From the biophysical point of view it appears reasonable to assume that β, being the quotient between the steady-light illuminance I and the bleaching constant L, is constant throughout the image. Thus, we only minimize over p and c, and we shall support this simplification in Section 6.1.

Mathematical validation of the energy minimization
If we restrict the values of parameters to a closed bounded set of the form a 2 [0, A], b 2 [0, B], c 2 [0, C], then E(a, b, c) in Eq (9) attains its nonnegative minimum. Since we keep the parameter b fixed, we only consider the first and third equations in Eq (11), which are According to a generalization of the Cauchy-Picard-Lindelöf Theorem any initial condition p (0) = p 0 and c(0) = c 0 yields a unique solution of this system of ordinary differential equations for all time provided that all second order partial derivatives of the right hand sides are uniformly bounded. Indeed, the second order partial derivatives in this case are so that the boundedness condition holds. Thus, the system (17) has solutions p(τ) and c(τ) defined for all artificial time τ. We now would like to claim that they converge towards a steady-state as the artificial time τ tends to infinity, and thatâ :¼ lim t!1 e pðtÞ and c :¼ lim t!1 cðtÞ, consistently with Eq (12), yield the minimizer of E.
To be precise, if we fix b, then Eq (17) are necessary conditions a minimizer of E must satisfy. To show that the parameters found by solving the system (17) indeed yield the minimum of E, we need to check that the Hessian matrix of E is positive definite. The latter, together with Eq (17), is a sufficient condition and equivalent to Since @ 2 p E ¼ 2ðT À T 0 Þ > 0 is always satisfied, we only need to check on the determinant. Its definition and the Cauchy-Schwartz inequality applied to the two functions e −bt and 1 yield det The inequality in these computations becomes an equality only if the two functions e −bt and 1 are linearly dependent, which would require b = 0. We fix b > 0 anyway, so that the Hessian of E is positive definite. Therefore, there is only one critical point of the energy E in Eq (9), which must be the global minimum. According to Eq (16), the gradient descent (13) and (15) yields this minimizer as n tends to infinity.

Postprocessing: inpainting retinal blood vessels
Blood in visible retinal vessels is a dominating absorber of excitation wavelengths used to excite RPE autofluorescence, so that pixels containing these blood vessels do not detect lipofuscin from the underlying RPE. The detected fluorescence signal in these pixels represents largely dark noise, which does not exhibit the bleaching behavior described by the model in Section 2.4. Further inaccuracies in pixels at retinal vessel edges are caused by the micro-movements of the eye during recording and image registration imprecision. As postprocessing, we shall remove retinal blood vessels from the rhodopsin map requiring retinal vessel detection and sophisticated mathematical techniques for image inpainting. The latter means that we mask the blood vessels in the rhodopsin maps and treat those pixels as occlusions (or missing pixels) in the context of an inpainting problem, in which surrounding pixels "diffuse" into occlusion areas. We shall recover the unknown intensity values in missing areas by using a technique that was shown to be efficient in restoring a wide variety of natural images [47].

Detecting blood vessels
Retinal blood vessel detection has already been widely addressed in the literature [48][49][50][51][52][53][54][55][56], and, in principle, we could choose one of the standard tools, see also [57,58] for related approaches. Nevertheless, we shall describe two semi-automated methods producing a mask of the pixels to be inpainted that directly evolves from the rhodopsin model itself.
1. Let f(x, y, t) denote our image stack, where (x, y) keeps track of the pixel location. For each fixed pair (x, y), the numerical minimization algorithm provides the outputsâðx; yÞ,bðx; yÞ, cðx; yÞ, as well as the number of iterations N(x, y) required for the gradient descent to converge. The "image" or "matrix" N contains edge information that enables tracing contours including the edges of the blood vessels. Thus, the speed of convergence of the numerical minimization scheme provides a tool to detect the retinal blood vessels.
2. The second method relies on the error between measured intensities and reconstructed autofluorescence from the computed model parameters. In other words, we aim to detect the pixel locations (x, y) where the cSLO measurements f(x, y, t) strongly deviate from the proposed bleaching model Eq (8). For instance, we can check the magnitude of Eðâ;b;ĉÞ.
Here, instead of the squared logarithmic deviation in Eðâ;b;ĉÞ, we use the total L 1 deviation at each pixel location, i.e, DAFðx; yÞ ¼ X m i¼1 jf ðx; y; t i Þ Àâðx; yÞe Àĉðx;yÞe Àbðx;yÞt i j; where t 1 , . . ., t m denote the time points associated to the frames in the recorded cSLO movie. Thresholding of the values between varying percentiles enables the identification of pixels, where the model does not fit the actual measurements, hence detecting the blood vessels. We denote those pixel locations by O.

Inpainting technique
Image inpainting is an active mathematical research field and several techniques have been proposed in the literature [59][60][61][62][63][64][65][66][67][68][69]. Here, we need image inpainting to remove the blood vessels as occlusions from the rhodopsin distribution map. Due to its robustness against measurement noise, we shall use a refined variational method developed in [37,47] for binary images, briefly described in the following. Since the cSLO camera output is 8-bit grayscale, we use bit-wise processing (work with 8 binary images separately) as outlined in [47], which further suppresses distortions effectively. Mathematically, we identify the images with binary functions g: and refer to [47] for the detailed motivation, see also [65]. Briefly, the first term forces the minimizerû to be close to g on [0, 1] 2 \O, so that only minimal changes in intensities outside blood vessels are allowed. The second and third term refer to the underlying Ginzburg-Landau model for u, meaning regularity conditions on u and the way g can be extended to O. For the reader familiar with such techniques, the term j u j 2 B denotes the Besov 1-2-2 norm computed using a wavelet expansion of u, and the details are described in [47]. Such variational methods based on optimization commonly yield much better results than more elementary interpolation techniques.
The pixels of retinal vessels and hence the missing area O are not perfectly well-defined due to measurement inaccuracies or micro-movements of the eye during the recordings and resulting image registration artifacts. To compensate for this, we replace the first term in Eq (21)  where the mask χ O : [0, 1] 2 ! [0, 1] smoothly changes from 0 to 1 at the "boundary" of O in [0, 1] 2 , thus gradually assigning a lower weight to the forcing u(x, y) % g(x, y) close to the unknown area. The final rhodopsin map then is a result of a two-step procedure. At the first step, three versions of the map were recovered using inpainting with slightly varying parameters of , μ in Eq (21) and smoothening of the mask χ O in Eq (22)-meaning the information near the vessels was given different levels of priority each time. At the second step, the average of the resulting maps was subject to adaptive binary wavelet thresholding, which eliminates further pixel artifacts, see [47] for details about the technique.

Summary of the numerical scheme
We first compute the macular pigment map. Second, we use the knowledge of the spatial macular pigment distribution to specify the initial values of the numerical fitting scheme when computing the rhodopsin map. Third, we detect retinal blood vessels to remove them from such maps and use image inpainting as outlined in the previous section yielding the final rod rhodopsin maps.

Results
We derived the methodology to quantify the rod rhodopsin distribution from cSLO movies that show autofluorescence brightening to a steady-state level as the rhodopsin bleaches. The present results section addresses three points. First, we shall compute the spatial variation ofb when minimizing E over all three parameters a, b, c, motivating the paradigm to keep b spatially constant. Second, we show the spatial rhodopsin distribution maps derived from numerical optimization, vessel detection and inpainting. Third, the characteristics of the computed rhodopsin maps are compared to the rod distributions derived in [7-9].

Numerical algorithm
The parameter β is a fundamental biophysical constant of rhodopsin-photon absorption and thus is presumed to be spatially invariant in our analysis, so that we only optimize over p and c [18,35]. We shall validate this approach by observing that the optimization over all 3 parameters leads to similar results. It should be mentioned that computed parametersâ,b, andĉ turn out positive without explicit enforcement. When solving Eq (10) for all 3 parameters by applying Eqs (13), (14), and (15), we see only minor variations of b(x, y) away from the blood vessels, the fovea and the optic disc. More precisely, we computed an average of the intensity curves over a large annulus-like area (5-8 degrees) away from the fovea and ran the fitting routine for the averaged experimental curve. A representative histogram of the recovered values of b, cf. Fig 1, shows strong concentration around one value, which we can then assign to b 0 .

Inpainted rhodopsin maps
The main objective of the present methodology paper is to develop the entire cycle of rod rhodopsin quantification, ranging from the biophysical model, the clinical measurements, and all the way to the mathematical optimization procedures and image analysis tools. Indeed, after retinal images are acquired and registered, building the rod rhodopsin maps from cSLO measurements is a 5-step procedure: 1. macular pigment: We spatially quantified macular pigment from multi-spectral autofluorescence image sets refined through self-consistency tests in    Temporal sequences of the intensity values (blue) and the corresponding fit (red) through minimizing E in Eq (9). Due to the low signal to noise ratio in cSLO measurements, there are large variations, but we can still recognize the overall trend. By averaging the cSLO movie over an angulus of 3 degree width, the noise and hence the variation are suppressed and the fit becomes tighter. (a) typical temporal sequence with its fitted curve at one pixel (b) temporal sequence from averaging intensities over an annulus (5-8 degrees) and fitted curve. The parameter b 0 is determined as the maximal count in the underlying histogram as described in Section 6.1, and a 0 is spatially modified according to the macular pigment maps.

preliminary rhodopsin map:
We then apply the gradient descent to perform the minimization on 8 × 8 pixel patches to reduce noise. Parameters for pixels near the fovea or near the major blood vessels show atypical behavior and the fitting changes qualitatively, Figs 4, 5. When ignoring such areas, plotting the parameterĉ yields a preliminary rhodopsin map in Fig 6a.

inpainting:
We inpaint the blood vessels to obtain the final rod rhodopsin map in Fig 7. The deviation between the preliminary rhodopsin map and those after image inpainting, assessed away from the mask, are below 5% in the amplitude, and < 1% in the mean-square error.

Rhodopsin maps compared to rod photoreceptor topology
We compare our rhodopsin distribution maps with typical rod topology as characterized in [7][8][9]. Note that we compare the distribution meaning comparison of actual densities in a qualitative but not quite quantitative fashion. Thus, local rhodopsin changes can be observed but the overall amplitude of our maps is arbitrary to some extent, so that we suppress units in our figures. There are no rods nor rhodopsin present in the fovea, and rhodopsin density increases with distance to the foveal center, Fig 7c. As observed in rod photoreceptor cells [8], the rhodopsin density increases most rapidly along the superior vertical meridian and less rapidly along the nazal horizontal meridian, see Fig 7d.

Discussion
Our analysis of cSLO lipofuscin AF image sequences in modestly (rod-protecting sun glasses and dim room lights) dark-adapted subjects is based on a physiological model of rhodopsin bleaching kinetics [18] to provide a rapid means for spatially mapping of the rhodopsin The gradient descent minimization scheme led to the optimized parametersâ andĉ, where the macular pigment distribution was used to refine the initial value a 0 . Blood vessels show up as artifacts in the resulting rhodopsin maps and must still be removed. We detect blood vessels through a scheme that evolves from the computational minimization of E related to the rhodopsin model. Spatial regions where the numerical algorithm converges significantly slower than in other image parts shall be identified as blood vessels. The number of required iterations leads to an image which can simply be thresholded to identify retinal blood vessels. (a) spatial map of parameterĉ optimal within the model. Brighter means increase in rhodopsin density, but brightest spots occur as registration artifacts along blood vessels. Therefore, the grayscale colormap is suppressed. (b) extracted mask of blood vessels derived from thresholding the number of iterations needed for convergence in the optimization scheme. Simple pixel growth in the mask would lead to complete coverage of the vessels. distribution within rod photoreceptors over a 30 degree field of view. Our procedure for mathematical optimization of these maps was guided by an interdisciplinary integration of the clinical image acquisition process, the biophysical models, and the mathematical analysis together with the image processing. The numerical scheme is modular, so that our methods for parameter fitting, vessel detection and image inpainting could be replaced with other tools, more familiar to any given clinical team. Nevertheless, our proposed methods use state-of-the-art in applied mathematics to provide improved numerical stability and reproducibility. Both models for macular pigment absorption and rhodopsin bleaching with regeneration used in our analysis, have been validated in the literature. Recently Morgan and Pugh [27] have published a method for measuring time-dependent rhodopsin bleaching based on reflected light using a cSLO device. Their method requires spatial averaging over specifically illuminated 4.8-degree retinal fields, which represents an additional clinical processing step. Our current method and analysis is advantageous in that it obtains reliable maps of bleachable rhodopsin and RPE autofluorescence levels simultaneously at relatively high (% 50μm) resolution for a large (30 degree) retinal field from a brief (1 minute) clinical imaging sequence. Melanin granules within the RPE undergo changes over the human lifespan, increasingly being converted by lysomal processing into lipofuscin-coated melanolipofuscin granules. They are uniformly dispersed among lipofuscin granules compared to the more apical distribution To derive the final rhodopsin map, we first detect the blood vessels. Instead of a binary vessel mask, we use a gradual interface at vessel borders resulting in a smooth vessel mask as proposed in Eq (22). Inpainting based on Eq (21) was used to remove the retinal blood vessels from the final rhodopsin map leading to a smooth rod rhodopsin map. We show the computed parameter γ, which is the sum of emission and excitation rhodopsin absorbance. Since the emission rhodopsin absorbance (590nm-600nm) is neglible, γ indeed is the rhodopsin absorbance at 488nm. The center of the fovea lacks rhodopsin whose density increases when moving apart from the center. Consistently with the rod distribution described in [8], rod rhodopsin increases most rapidly along the superior vertical meridian and increases least rapidly along the nasal horizontal meridian. Although we do not see a connected hot spot of highest rod rhodopsin density, we observe larger and more connected areas of highest rhodopsin density in the superior retina than in the inferior retina, again, consistent with the rod distribution in [8]. (a) cropped rhodopsin map from Fig 6a to be inpainted (b) smooth mask χ used in the spatial consistency forcing term (c) the final output showing the distribution of the rhodopsin optical density γ, i.e., the rhodopsin absorbance at 488nm. Units are suppressed to indicate that we compute its distribution rather than absolute quantitative optical densities as we are not able to validate overall amplitudes due to background in our image sets. (d) horizontal and vertical rhodopsin profiles, averaged over small stripes. of RPE melanin observed at earlier ages, cf. [32,70]. Such changes increase the local excitation light levels reaching fluorescent bisretinoids. Where melanin becomes substantially oxidized its visible absorption is reduced leading to increases in detected local autofluorescence. However, the effect of melanin changes with age or pathology should not affect the local rod irradiance and therefore neither the local kinetics nor fractional amplitude seen in our analysis of the cSLO AF imaging sequences. Thus, it seems reasonable to explicitly model the macular pigment signal attenuation but neglect the melanin contribution. The physiological variables α(x, y), β(x, y), and γ(x, y) appear to be consistent as we optimize the model parameters using the steady state of a system of ODEs. Our model and parameter optimization provide plausible maps of rod rhodopsin distributions in normal individuals using a clinically simple 1min noninvasive imaging sequence. The fitting curves enable us to detect local rhodopsin changes up to a resolution of 50μm, which is competitive to other standard clinical instrument, and higher resolution is usually only achieved by more specialized instruments. Our rod rhodopsin maps are consistent with characteristics of typical rod distribution observed in the literature [7,9], which appears to qualitatively validate our approach. Rod density increases most rapidly along the superior vertical meridian and increases least rapidly along the nasal horizontal meridian, just as in our rod rhodopsin map. We observe larger and more connected areas of highest rhodopsin density in the superior retina than in the inferior retina, again, consistent with the rod distribution in [8].
We must point out though that we claim spatial mapping of rhodopsin distribution that provides qualitative spatial rhodopsin density maps. Image background in standard fundus and cSLO cameras prevents us at this point in time from claiming quantitative spatial mapping of the actual rod rhodopsin density. Image background needs to be removed before our actual numerical scheme becomes active leading to amplitude differences that depend on the method and amount of background removal. Nonetheless, our claim of enabling detection of spatial variations in rod rhodopsin density appears valid and may be sufficient to identify clinically relevant early pathological changes. Once background removal is solved satisfactorily as an additional module in our image analysis pipeline, our proposed method would immediately enable spacial mapping of rod rhodopsin density in a quantitative fashion.
Given the recent hypotheses that increased bisretinoid levels in the RPE cause increased stress and dysfunction in Stargardt's macular dystrophy and age-related maculopathy, our ability to simultaneously map local changes in bleachable rhodopsin along with changes in bisretinoid (A2E) fluorescence may improve evaluation of these conditions. Hence, it may help better assessing early interventions to prevent progression of highly localized pathologies. By computing simpler fractional change maps between early and late cSLO bleaching image sequences, we have observed decreased bleachable rhodopsin at hyperfluorescent lesion edges in a central Stargardt's lesion as well as decreasing bleachable rhodopsin as one approaches the hypofluorescent centers of each reticular druse. Using the more integrated analysis presented here, could make such clinical observations more precise and lead to better understanding the relationships between increased bisretinoids, RPE dysfunction, and local visual cycle deficiencies that led to rod dysfunction and loss.

Conclusions
We developed noninvasive multispectral retinal fluorescence imaging by adding filter sets to standard fundus cameras enabling us to compute spatial macular pigment maps. The macular pigment maps are then combined with 1-minute-long cSLO autofluorescence imaging sequences starting at a moderately dark-adapted state to create additional maps of rod rhodopsin and RPE bisretinoid fluorescence using a biophysical bleaching and regeneration model and mathematical optimization techniques. To foster future collaborations and synergies, we outlined both the biomedical and the mathematical approach in detail. Retinal blood vessels, blocking the 488nm laser beam, are removed from the rhodopsin maps by image inpainting. We computed spatial rod rhodopsin distribution maps that match the characteristics of rod distribution observed in the literature. Further efforts are required to derive quantitative rod rhodopsin density maps from our proposed approach. Our limitation to qualitative density maps is due to the necessity of background removal in our image sets, whose methodology we have not been able to validate up to now. Our modular image analysis approach is well-suited to incorporate additional analysis steps, so that a validated background removal methodology would enable spacial mapping of rod rhodopsin density in a quantitative fashion.
We derived so far a non-invasive clinical methodology to map local rod rhodopsin distribution at % 50μm resolution using widely available clinical imaging devices. Our approach appears well-suited to identify localized changes in the rod photoreceptors and relate them to local changes within the underlying RPE bisretinoid fluorescence levels. Although this method appears robust in younger eyes and in elderly pseudophakes, we need further optimization of our methods for those elderly patients, in which age-related lens changes combine with spatially varying diffuse reflectance of the retina creating non-uniform backgrounds in our image sets. Additional refinement and validation may enable us in future studies to correlate local rhodopsin reductions and autofluorescence levels within early microscopic lesions in order to better understand their natural progression and their responses to benign disease prevention strategies that are most likely to be effective in such early disease states.