Wavelength-dependent optical properties of melanosomes in retinal pigmented epithelium and their changes with melanin bleaching: a numerical study.

In this paper, we present the first numerical study on full metrics of wavelength-dependent optical properties of melanosomes in retinal pigmented epithelial (RPE) cells. T-matrix method was used to simulate the spheroidal shapes of mature melanosomes, and the complex refractive index was calculated by a subtractive Kramers-Kronig relation for melanin. The validity of the method was first confirmed by Mie theory, and corroborated by a comparison between visible light and near infrared (NIR) optical coherence tomography (OCT) on human retinal imaging. We also studied the changes of melanosome optical properties due to melanin bleaching by numerically reducing the absorption of melanin. This study implies a unique approach to detect melanin changes specifically in RPE by a spectroscopic contrast of optical coherence tomography.


Introduction
Melanin is an intrinsic pigment that is enriched in retinal pigmented epithelial (RPE) cells. In addition to absorbing excessive light passing though the retina, melanin also has important protective roles mediated by several biochemical mechanisms, such as scavenging reactive free radicals [1,2], quenching electronically excited states [3] and sequestering redox-active metal ions [4,5]. Those functions lead to the hypothesis that melanin bleaching (i.e. loss of absorption capability) is a key factor in the pathogenesis of age-related macular degeneration (AMD). Indeed, the melanin content in RPE drops around 50% in eyes of people aged over 70 years [4], which is correlated to the increased AMD incidence with aging. Also, epidemiologically all forms of AMD are more prevalent in the white population than more darkly pigmented races, suggesting the pathological relevance of melanin. Paradoxically, a recent clinical study showed that no significant correlation was found between melanin pigmentation and early age-related maculopathy [6]. Those seemingly conflicting literatures make the role of melanin in AMD still a controversial and inconclusive subject.
A significant part of the controversy stems from the lack of an accurate method of measuring melanin in vivo, and also the lack of understanding of the optical properties of melanin and melanosomes (the organelle in charge of synthesis, store and transport of melanin). The current clinical method of quantifying retinal melanin content is by fundus reflectance photography [6], which has two major limitations. Firstly it does not differentiate melanin/melanosomes in RPE and choroid, and secondly it does not measure melanin content directly but calculates melanin optical density indirectly by an inverse model. Both of these two limitations can have pivotal impacts in in vivo data interpretation.
On the other hand, it is still challenging to provide optical characterization of melanin and melanosomes. The challenges not only arise from the complex structures of melanin itself, but also the fact that melanin does not exist in free form in vivo, but is exclusively compacted within melanosomes. When we reviewed the literature, the measurements of absorption coefficient (μ a ) from isolated melanosomes have been documented by optical transmission microscopy [7], and later by explosive vaporization [8], and microcavitation [9]. However, little has been reported on the scattering properties of melanosomes, despite the fact that the signals collected in reflectance retinal imaging are ultimately determined by the light scattering.
In this study, we aim to provide the full metrics of optical properties of RPE melanosomes by a numerical simulation, and predict their changes with melanin bleaching. This study is based on the observation that matured melanosomes have spheroidal geometries and their sizes are binomially distributed [10]. With an assumption that melanin is uniformly distributed in the melanosomes, we adopted the T-Matrix method by Mishchenko et al. [11] to rigorously compute light scattering from randomly oriented spheroids to simulate melanosomes. The complex refractive index for melanin was calculated by a subtractive Kramers-Kronig (KK) relation. We also numerically simulated melanin bleaching by reducing μa within melanosomes and examined the changes of their optical properties. In addition, we qualitatively corroborated our simulation with the experimental data from visible light and NIR optical coherence tomography (OCT). The clinical implication is that spectroscopic OCT may provide unique advantages in quantifying melanosomes or melanin specifically in RPE in vivo.

T-matrix for simulating optical properties of randomly oriented spheroids
The T-matrix code developed by Mishchenko et al. was used to simulate the optical properties of randomly oriented spheroids in this study. The code was downloaded from Ref [11], and modified to allow updating simulation parameters iteratively at varying wavelengths. A Matlab script was written to update the geometry, size distribution, complex number of refractive index (RI) and wavelength.
In the T-Matrix method, the intensity of the scattered beam is described by a Stokes scattering matrix F(Θ) [12] where Θ (degree) is the scattering angle between the incidence beam and scattered beam within the scattering plane. The elements a 1 , a 2 , a 3 , a 4 and b 1 , b 2 denotes the eight elements in the Stokes scattering matrix, following the notations by Mishchenko et al [12]. The elements have been also commonly denoted by S 11 , S 22 , S 33 , S 44 and S 12 , S 34 in other literatures [13], and their conversion is explained in a later section. The Stokes vector of the scattering beam can then be transformed from the Stokes vector of the incident beam by the following equation: where C sca (μm 2 ) is the scattering cross section per particle, and R is the distance from the particle. The Stokes vector follows the conventional definition as a (4x1) column having the Stokes parameters I, Q, U and V, describing the light polarization [13]. I, Q, U and V is defined by, , where E ║ and E ┴ are parallel and perpendicular components of the electric field to the plane of scattering, respectively. For randomly oriented scatters and unpolarized incidence, the backscattering cross section C b (μm 2 ) follows the conventional definition of radar backscattering cross section, and is defined by the (1,1)-element of F(Θ) in Eq. (1): (4) The F matrix, the scattering (C sca ), extinction (C ext ), absorption cross section (C a ), and anisotropy factor ( g ) are readily calculated by the T-matrix code for the full metrics of optical properties.
In order to derive the μ a (μm −1 ) and μ s (μm −1 ), we scaled the cross sections to macroscopic properties by where f v being the volume fraction within a local volume V (μm 3 ). Within the melanosomes, f v is equal to 1 and V is calculated form the equivalent-sphere radius r c (μm) calculated by Tmatrix.
In order to take into account the natural size distribution of melanosomes, the T-matrix code includes various analytical functions, in particular the gamma function [12], in which r (μm) is the equivalent-sphere radius. The three other constants γ (unitless), α (unitless), and r c (μm) defines the distribution, where r c scales the mean size of the spheroids, α is the exponent of the power law term and γ determines how fast the exponential term decays with r. By taking γ = 3, r c = r 0 , and the gamma function approximates the normal distribution, defined by a Gaussian function where r 0 (μm) and σ (μm) are the mean radius and standard derivation, respectively: To average the cross sections over a size distribution, a weighted integration by the size fraction p(r) is calculated by

Spheroidal geometries of melanosomes in retinal pigmented epithelium
The melanosomes are generally in a shape of prolate spheroids, with the long axis length ranging roughly from 0.5 to 3μm. We adopted the geometric characterization of melanosomes from a transmission electron microscopy study for matured bovine RPE [10]. In this study, two shapes of spheroids were found with one being more rounded and the other being more elongated. We abbreviate these two types of spheroids as S r and S e . The subscripts r, e are short for rounded and elongated, respectively. The mean value and standard deviation of the long and short axis for S r and S e are summarized in Table 1. As the relative standard deviations are fairly consistent for both axes in each shape (S r : 0.18 for long axis, 0.14 for short axis; S e : 0.24 for long axis, 0.22 for short axis), we simplified the shapes of melanosomes by two types of spheroids with fixed aspect ratios (S r : 1.27 and S e : 3.52), and each with their axis length normally distributed, as shown in Fig. 1(c).
The T-matrix code randomizes the orientation and uses one equivalent-sphere radius r c to denote the size as in Eq. (6). The value of α is taken as α = (α l α s 2 ) 1/3 , where l and s denote the long and short axis, respectively. The values of α, γ and r c are summarized in Table 1 as well.  [10]. The mean values r 0 of the normal distribution is listed with their standard deviations σ in the parenthesis. b The abbreviation of S r and S e is for rounded spheroids and elongated spheroids, denoting the two types of shapes of melanosomes in RPE.

Complex refractive index of melanosomes in RPE
T-matrix method also requires the complex refractive index (RI) n(ω) = n(ω) + ik(ω) within melanosomes. The imaginary part of RI k(ω) represents the absorption, and is calculated by where c (3e 8 m/s) is the speed of light in vacuum, ω (rad/s) is the angular frequency of the light; and μ a (µm −1 ) is the absorption coefficient of the melanosome. The wavelengthdependent value of μ a has been previously characterized as μ a = 6.49 × 10 12 λ -3.48 (cm −1 ), where λ is the wavelength expressed in nanometers [14]. The real part of the RI n(ω) is calculated by a subtractive KK relations [15], where n(ω 0 ) is the RI at a reference frequency ω 0 . P is the Cauchy principal values of the integration. The real part of the RI for melanin is generally found to be within 1.6-1.7 [16,17]. We took n(ω 0 ) = 1.65 at 650nm, and the wavelength-dependent complex RI of melanosome can be calculated as shown in Fig. 1(a), 1(b). We note that the real part of the RI is fairly consistent with variations within ± 0.0015 in the 0.5-0.9 μm wavelength range, while the imaginary part of RI decays at longer wavelengths. The sharp dip at ~350nm in Fig. 1(a) is due to the finite range of the integration, as indicated by the gray area.

Mie theory for validating T-matrix code
The open source Matlab code by Christian Maetzler [18] was used to generate the optical properties for spherical particles in order to compare the T-matrix code [19]. Mie theory calculats the scattering matrix S to describe the scattering field form spherical scatters, where E ║i and E ┴i are parallel and perpendicular components of the incident field to the plane of scattering, respectively, and E ║s and E ┴s are those of the scattered wave. The relation between incident and scattered Stocks parameters follows [13], where k is the wave number and R is the distance from the particle. A conversion from the Stocks scattering matrix to F(Θ) in Eq. (1) can be made by, The Mie theory code outputs the efficiencies for extinction, scattering, absorption, backscattering (Q ext , Q sca , Q abs , Q b ) and anisotropy factor g. The corresponding cross section per particle can then be calculated by 2 , where r (μm) is the radius of the sphere. Equation (9) was used to include the size distribution of scatters for Mie theory as well.

Results
T-matrix calculation is an effective approach to calculate the scattering properties from spheroidal geometries and various orientations. The code provided by Mishchenko et al. has been widely used in a broad range of scientific regimes including the field of biomedical optics [20]. The code with single orientation has been experimentally validated [21] and successfully applied to characterize the elastic light scattering from cell nuclei [22,23].
Adapting the similar code with random orientation of spheroids, we aim to characterize the optical properties of melanosomes in RPE. In order to ensure the code runs properly, we first compared the T-matrix calculation with Mie theory on spherical scatters. Figure 2 shows the F-matrix output by two methods on polystyrene microspheres (mean diameter d = 1 μm, n = 1.59) in aqueous medium (n = 1.33) at λ = 0.5μm. The F-matrix in Mie theory was converted from the Mie scattering matrix by Eq. (14). Two scenarios were simulated where the microspheres were either monodisperse or the size are normally distributed. In the monodisperse case, the diameter range was limited within 1 ± 1e −6 μm; while in the normally distributed case, the diameter followed a normal distribution with the standard deviation equal to 0.25μm. In each scenario, we further compared the wavelength-dependent optical properties by two methods in Fig. 3. For all the above comparisons, the results from T-matrix showed excellent agreement with Mie theory, with deviation less than 6.4%. The averaged deviation between two methods were summarized in Table 2.  Having numerically validated the method, we next simulated the optical properties of melanosomes in RPE. The S r and S e spheroidal shapes were separately calculated with their size distribution, and their results were summed by their volume fractions (S r :S e = 0.82:0.18) [10]. We first examined the conventional absorption and scattering coefficients, μ a and μ s as shown in Fig. 4. The spectrum of μ a showed good agreement with the literature values in Ref.
14. We note that μ s is at least one order of magnitude higher than μ a despite the strong absorption of melanin. This is due to the high RI of melanosomes contracting to the aqueous medium. Fig. 4. T-matrix simulation of (a) Wavelength-dependent absorption coefficient μ a and (b) scattering coefficient μ s from melanosomes in RPE. In (a), the μ a from the reference [14] (red solid curve) and T-matrix calculation (blue dash curve) were plotted together for comparison.
The spectra of C b , C ext , C sca , and g of melanosomes are shown in Fig. 5. We found that C b , increases in visible light range, peaks around 700nm, and decreases in longer wavelength. The value of C b at 800nm is about ~25% higher than that at 500nm. The spectra of C ext and C sca monotonically decrease with respect to wavelength, reducing by ~50% from 500nm to 800nm. g also exhibits monotonic decay at longer wavelengths. Next we simulated the melanin bleaching, by numerically scaling µ a by a constant in Eq. (10) (i.e. 1 for no bleaching, 0.5 for 50% bleaching, and 1e −4 for 99.99% bleaching). C b shows the largest bleaching-induced contrast in visible light range, increasing by ~11% and 33% at 500nm with 50% and 99.99% bleaching, respectively. In comparison, C ext and C sca have only increased by ~1% and 2% with 99.99% bleaching, respectively. The complete data set for no bleaching is included in Appendix.
In order to provide some experimental support for our simulation, at least qualitatively, we compared in vivo human retinal imaging by optical coherence tomography (OCT) operated in the visible light (530-620nm) and NIR (~850nm) range from a normal eye ( Fig.  6(a), 6(b)). The data is from Yi et al's previous study on human retinal imaging using visible light OCT [24]. Because the OCT signal is determined by C b for a thin layer of RPE (~20μm), the relative change between two wavelength ranges can be correlated to the result in Fig. 5(a). As shown in Fig. 6(c), fifty A-lines from three regions were averaged to generate the depth profiles. After subtracting the noise floor, the depth profiles were normalized to the intensity level at the IS/OS peak. We can see that in NIR OCT, the signal from RPE is significantly higher than that in visible light OCT, by around 27%, 48%, and 28% in three regions. This contrast corroborates with our simulation that C b is higher in NIR at ~800nm than visible range. A similar spectroscopic contrast was also observed in rodents [25].

Discussion
In summary, we presented a rigorous simulation using T-matrix method to characterize the wavelength-dependent optical properties of melanosomes in RPE, considering their spheroidal geometries, size distributions and complex RI. We also simulated the melanin bleaching by numerically reducing the absorption of melanosomes, and examined the changes of their optical properties. Experimental comparison between vis-OCT and NIR OCT showed that the backscattering signal from RPE is significantly lower in visible light range than NIR, which corroborates our simulation showing a similar spectroscopic trend in Cb.
There are several influential parameters in our simulation that require further discussion. In our calculation of complex RI of melanin, we used 1.65 for melanosomes at 650nm in KK relation which creates a large RI contrast to the aqueous medium. This high contrast resulted in a high μ s that is at least one order of magnitude higher than μ a , despite melanin itself being a strong absorber. In vivo experimental evidence of melanosomes' high RI was provided by confocal reflectance microscopic imaging of human skin [26]. Melanosomes in pigmented skin generated high brightness over other cellular components compared with unpigmented skin, suggesting the strong backscattering efficiency presumably resulted from the high RI contrast. The high RI of melanin was also validated in vitro using natural melanin particles from Sepia officinalis [16,17]. In those experiments, integrating spheres were used to measure μ a and μ s of the melanin suspension, and the RI of melanin was predicted using the particle size characterization and concentration. In the above literatures, the RI of melanin is generally found around 1.6-1.7 in visible and NIR range. Recently, a relevant study on the RI of melanin in bird feather showed the value of RI being 1.7-1.8 in visible light range [27]. Although in this case the melanin is from a different species, the high melanin RI is in a general agreement with previous literatures. We should also note that the calculation of KK relation has errors in both ends of the spectrum, due to the finite range of the integration as shown in Fig. 1(a). In the future, multiple subtractive KK analysis may be applied to reduce the uncertainty of the calculation over the entire spectrum with more experimental reference points. In our current study, the validity appears to hold over the spectral range larger than 400nm, as suggested in Fig. 4(a).
Another important assumption is that the melanin bleaching is not affecting the real part of RI. Given that the μs of melanosomes is largely due to the RI contrast to the medium, this assumption indicates that the melanin bleaching would not induce drastic μ s changes. This is confirmed by Fig. 4(b) and 4(c) showing that melanin bleaching induces moderate changes of μ s , even under 99.99% loss of absorption. Understanding the physical process of melanin bleaching is the key to this assumption. Watt et al. showed that the completely solubilized melanin in dimethyl formamide has almost identical absorption spectrum as the polymerized melanin in water, indicating that the mechanisms of the absorption is from the melanin protomolecules at sub-nanometer scale, instead of high-order polymerization [28]. Littrell et al. studied the structures of melanin and revealed a sheets-like structures composed of stacks of protomolecules [29]. When chemical bleaching was introduced, the structural changes starts at the length scale of ~10Å by reducing the sheet distance. Interestingly, they found that the reduction in absorbance during bleaching is not uniform as would be the case if the density of absorbing melanin particles were simply reduced without any structural change. Thus, they suggested that there is structural reconfiguration of melanin protomolecules upon bleaching, possibly reversible, while the microscopic mass density and thus the real part RI might not change significantly. In any case, the process of melanin bleaching awaits more investigation, and this assumption ultimately would need further experimental validation.
We should note that our study does not consider the ultra-structures of melanin within the melanosomes. Melanin itself is a complex macromolecule consisting of polymers of eumelanin and pheomelanin. The higher-order structures of intact melanin within the melanosome is extremely difficult to predict, and also beyond the capability of T-matrix method. Instead, we assumed that melanin are uniformly distributed within melanosomes and provide the RI contrast to cytoplasm. We also did not consider the possible variation of melanin contents among melanosomes. However, given that melanin bleaching is globally present on retina, a measurement over a large field of view could still provide an ensemble averaging that reflects the overall melanin bleaching.
The implication of the simulation and experimental evidence is that the melanin bleaching in RPE could potentially be detected by spectroscopic OCT. Given the high RI of melanin and their abundance in RPE, the signal detected in OCT may be contributed predominantly by melanosomes. Thus, spectroscopic analysis across visible and NIR ranges in OCT presents a unique approach to detect melanin loss and melanosome changes specifically in RPE, excluding any confounding factors from choroidal melanin. The depth resolving capability to investigate only RPE layer is advantageous to other fundus photography based imaging methods, which may offer better diagnostic values in AMD and a fresh look at the pathological roles of melanin.
In summary, we provided a rigorous simulation using T-matrix method to characterize the full metrics of wavelength-dependent optical properties of melanosomes in RPE. This study provides insight to specifically assess melanosomes in RPE using spectroscopic OCT. We discussed the limitation of several key assumptions in our simulation. Future experimental measurements addressing those limitations are on-going.

Appendix
We included the data of wavelength dependent optical properties of melanosomes in RPE with no bleaching in Appendix Table 3 as plotted in Fig. 5. Table 4 and 5 list the same metrics from S r and S e spheroidal shapes, separately.