1. Introduction

It has been extensively demonstrated that the spectroscopic features of light measured from biological tissues contain valuable information on cellular morphology and biological functions. Recently, Spectroscopic Optical Coherence Tomography (SOCT)[1] has been developed as an important enhancement of traditional OCT, providing not only a map of reflectivity intensity but also spatially resolved spectroscopic information. As in all spectroscopy-based optical techniques for tissue characterization, the spectral contrast exploited by SOCT originates from the properties of either endogenous tissue or exogenous contrast agents. Given that SOCT is based on coherence detection approach, the mechanisms investigated by SOCT are mostly focused on chemical-dependent absorption and structure-dependent elastic scattering. Endogenous chemicals such as hemoglobin[2] and melanin[3] as well as exogenous contrast-enhancing dyes[4–7] determine the wavelength-dependent absorption profiles, while the size, shape and refractive index of cellular organelles[8, 9] as well as nanoparticle-based exogenous contrast agents[10–14] determine the scattering spectral signatures.

Due to the complex nature of biological tissue, it has always been a challenge to definitely associate the measured spectroscopic features with tissue structure and composition in the micro-nano scale. The development in high-resolution OCT system[15, 16] presents us the unique opportunity to measure spectroscopic signals from single or a small number of tissue scatterers, and has generated great research interests in correlating SCOT measurements with tissue properties in the micro-scale[17–19]. The study presented in this paper focuses on numerical simulation and phantom-based experiments of ultra-wide bandwidth SOCT measurements on dispersed scatterers with both scattering and absorption spectral profiles. These phantom studies are analogous to scenarios of imaging single cells or cellular organelles with endogenous or exogenous absorbing chemicals within turbid biological tissue, such as red-blood cells with hemoglobin or organelles enhanced with exogenous dyes. We demonstrate that it is feasible to locally extract the architecture-dependent scattering spectrum and the fluorophore-determined absorption spectrum from fluorescent microspheres via post-processing of the SOCT measurements. Furthermore, we illustrate the possibility of imaging and identifying particles labeled with dyes in a same phantom with ultra-wide bandwidth SOCT measurements.

2. Methods

2.1 Experimental setup

Our experimental setup adopts a parallel[20] common-path Fourier-domain OCT(FDOCT) with broadband thermal light source as shown in Fig. 1(a). Light from a 100W Xenon lamp (Oriel Apex, Newport) is focused by an achromatic lens L1 on an aperture with 0.8cm in diameter, which is then directly imaged onto the sample with another lens L2 and a 0.1-NA objective lens, creating a full-field illumination area with about 400µm in diameter. The effective numerical aperture of L1 and L2 is estimated as 1.8. The sample is attached to the bottom surface of a 3mm-thick quartz plate. The interference signal that is formed from light reflected by the bottom surface of the plate and light backscattered from the sample is detected by a spectrograph (Acton SP2150i) coupled to a CCD camera (Princeton instruments PIMAX-2 HB512). The bandwidth and the resolution of the spectrograph are 430–650nm and 0.4nm respectively, yielding an axial resolution of 0.8µm in water. The lateral resolution of the imaging system is of about 3µm in water as determined by the objective lens.

The post processing involves inverse Fourier transform followed by short-time Fourier transform (briefly, IFT-STFT method) from FDOCT signal. The specific data collection and processing routine are illustrated in Fig. 1(b)–(d). The spectrograph passes 1D signal in image plane through a 10µm-wide entrance slit and disperses the signal into spectra. The CCD camera takes one shot of a 2D matrix with its x-axis corresponding to a spatial position along the 1D image signal and its y-axis as the spectrum in wavelength, as shown in Fig. 1(c). Before further processing, these spectral measurements are normalized by light source spectrum and resampled for equal k-domain interval. Then IFT is performed to obtain the spatial intensity profile in the depth (z) direction. Thus a B-scan image is formed in the x-z plane of the sample as shown in Fig. 1(d). In order to obtain spectral profile from a specific depth, we performed STFT with a Gaussian-shaped window located on the corresponding z position (Fig. 1(b)). Three sub-Figs. in Fig. 1(b) demonstrate the sequence of the processing by plotting the k-domain A-scan signal, the z-domain intensity distribution after IFT, and the recovered spectra after STFT. The red Gaussian window selects our interest spatial location. Note that the z-domain data are complex numbers before and after the STFT and only the amplitude is plotted as resultant extracted spectrum. The window width in the STFT determines the trade-off between the spatial and spectral resolution, where applying a narrower Gaussian window in the spatial domain results in a wider spectral point spread function (PSF) at each sampling.

 

Fig. 1. Schematic sketch of the system and data processing procedure. (a) Source image are focused by Lens (L1) on the aperture (AP). Lens (L2) and NA=0.1 objective lens (OL) focuses AP filtered light on sample. The beamsplitter (BS) redirects light collected by OL to spectroscope coupled by spectrograph (SG) and CCD camera. (b–d) procedure of post-processing and image formation.

Download Full Size | PDF

2.2 Simulation

In order to provide better system characterization and interpretation to the resulted OCT image, we have developed a numerical tool to simulate the formation of OCT images with wide-field imaging configuration. Light scattered from the scatterers in the sample is calculated via solving Maxwell’s Equations. For the case of spherical scatterers illuminated by plane waves, Mie theory provides rigorous calculation on both the amplitude and phase of scattered E-field E _i in the far-field. We note that the sample illumination field in our current experimental setup can be considered as a spatially incoherent homogeneous field in paraxial region, since the illumination area on sample is a direct projection of the aperture, which spatially filters the image of light source and passes only a homogeneous portion. In our setup, the full-field illumination area is about 400µm in diameter. Although the plane-wave illumination condition is not strictly satisfied here, we speculate that the scattered signal can be approximated by that generated with a plane wave illumination as long as an individual scatterer is located inside the depth of focus, which is about 55µm in our experiment setup. This speculation is supported by the comparison between the numerical and experimental results as discussed in the following section. Here we employ a Mie-theory based computational code developed by Xu et al [21] to calculated E _i(θ,ϕ) of a sphere with specific size and complex refractive index n, where θ and ϕ are the zenith and azimuth scattering angle.

Next, we calculate the E-field distribution at the pupil of the objective lens. This E-field is contributed by two components: the scattered field by the sample and the reflected field from the reference plane. The field distribution calculation is based on the geometric arrangement illustrated in Fig. 2. The origin of the coordinate system is placed at the center of the sphere and the incident wave propagates in the +z direction. The objective lens is placed at z=-l_c, where l_c is the distance between the sphere and the pupil of objective lens. The reference plane is placed at z=-l_ref. Since it is a far-field solution, the scattered light is regarded as spherical wave spreading from the microsphere as in Fig. 2(a). According to the definition of far-field field quantities, the scattered E-field component at the pupil of the objective lens can be calculated as:

where k is the wave number, r_p is the position vector of a point on the pupil of lens and |r_p| is calculated as $\mid {\mathbf{r}}_{\mathbf{p}}\mid =l_c\sqrt{1+{\tan }^2\left(\theta \right)}.$. The reference E-field component E _ref(p) can be approximated as a Gaussian beam:

which is reflected by the reference arm. Here, _p is the radial vector in the pupil plane transverse to optical axis as shown in Fig. 2(b), w(l_c) is the beam 1/e radius and R(l_c) is the curvature radius of the wavefront at the pupil. The distance from reference plane to the focus (origin) is much smaller than the focal length so that we approximate the reflected Gaussian beam backward follows illumination profile with only a phase delay. Given that the initial phase of E _i(θ,ϕ) is zero, we add the phase delay in the last exponential term corresponding to optical pathlength difference between reference and the sphere. E_const is complex constant dependent only on the distance between the reference arm and the pupil, which remains unchanged during experiment. The Gaussian profile is designed with beam waist w ₀=3µm so that the lateral resolution is kept consistent with experimental setup. The spatial incoherent illumination will just generate stable interference fringe by light reflected from one particular Gaussian beam, since it is completely phase irrelevant with every other Gaussian beams. The total field at the objective pupil can thus be calculated as E(r_p)=E _ref(p)+E _s(r_p). In the above formulation, we note that both r_p and p indicate the position of a spot within the objective pupil and are functions of θ and ϕ. They are written in different vector definitions for convenience in the formulation.

 

Fig. 2. Schematic sketch to illustrate the image formation algorithm. (a) Scattering field from microsphere located at origin. E _s is the scattered field from the sample when reaching the pupil of objective lens. (b) Reflected reference field. The reference plane is placed between the objective lens and the scatterers. The reference E-field E _ref (dots line) is approximated as a reflected Gaussian beam, backward following the illumination Gaussian profile (solid line) with a phase delay.

Download Full Size | PDF

The next step involves synthesizing the image intensity distribution after the pupil plane E-field E(r_p) is focused by the objective lens. If the scattering particle is centered in the focal plane of the objective lens, the focused E-field distribution in the image plane is calculated as:

Here, r_q is the position vector of a point in the image plane. The integration range of the zenith angle Δθ is determined by the numerical aperture (NA) of the objective lens as NA=sin(Δθ/2) and M represents the magnification of whole imaging system, which is set to be 1 through this study. This formulation is similar to the image formation algorithm proposed in [22] for synthesizing regular wide-field microscopy images. Both algorithms use the paraxial approximation, and are thus only valid for imaging systems with low or moderate NA. However, the two algorithms are different in their phase compensation terms within the integral. In our case, we assume that the lens re-adjusts the phase distribution of E(r_p) in such a manner that for a point source at any location in the focal plane, its corresponding image in the image plane is formed with complete constructive interference. In other words, our formulation also assumes that the imaging system is free of any type of aberrations.

The above calculation is repeated for a range of wavelength covering the same bandwidth of our experimental system and thus the FDOCT data set E _I(r_p,λ) is obtained numerically with a specific scattering sphere. A B-scan wavelength-dependent intensity signal can be obtained which is further resampled into k-domain and IFT is taken to ultimately synthesize the image. The same IFT-STFT method is performed to extract and analyze the spectrum. In order to incorporate the absorbent feature of fluorescent sphere, we use the fluorophore absorption spectrum provided by manufacturer to calculate the wavelength-dependent imaginary part of the refractive index n’(λ) in simulation. The absorption spectrum was first normalized to a range between 0 and 1, and further scaled by a constant τ, which is the imaginary part of refractive index at absorption spectrum peak. By comparing the simulation and experimental backscattering spectrum, we empirically set τ equal to 0.02. The linear absorption coefficient µ_a(λ) can be calculated as µ_a(λ)=2kn’(λ) and then the molar concentration of fluorophore c also can be calculated as µ_a(λ)=cε(λ) if with known molar absorption coefficient ε(λ). We estimated µ_a=5.6×10³ cm^-1 at the absorption spectrum peak 450nm for the fluorophore trapped in the 6µm microspheres.

3. Result

Our first experimental phantom is made of fluorescent microspheres embedded in transparent gels. A solution composed of agarose powder (Sigma-aldrich) and water was microwaved until boiling and the solution appears transparent after the powder is completely dissolved. Microspheres solution is added in and stirred when the gel solution cools down at about 45–50 degree C, then the gel phantom is solidified until further cooling to room temperature. Fluorescent microspheres (yellowgreen 6µm nominal diameter, standard deviation 0.146µm, Polyscences Inc.) were experimentally imaged with our FDOCT setup. Since OCT is based on coherent scattering, the fluorescence emission signal will not be present in our interference signal. Therefore, these microspheres can be regarded as scattering structures with an absorbing spectral profile depending on the fluorophore. The particle concentration is controlled to be low enough so that the microspheres are well dispersed and their scattering spectrum can be analyzed individually. The example image of one single microsphere is plotted in Fig. 3(b). In the corresponding simulation, a 6µm sphere is placed 50µm below the reference plane and the OCT image was synthesized for a 2-D longitudinal cross-section intensity image of the 3-D space as I _I(x,z:y=0) with the sphere in focus as shown in Fig. 3(a). In these two Figs, the single microsphere is presented as multiple peaks in both experimental and simulated OCT images. The patterns of the multiple peaks are nearly identical, supporting the accuracy of our simulation methodology.

We then perform STFT with different Gaussian windows to recover both fluorophore-dependent absorption and architecture-dependent scattering feature from OCT images. First, we conduct STFT with a window located at the strongest intensity peak, where z=55µm and 50.5µm in Figs. 3(a) and 3(b) respectively. The window size is selected to yield 1.2µm axial resolution, corresponding to a spectral resolution of about 35nm, and the extracted spectra are plotted in Fig. 3(c). Both experimental and simulated data show an absorption band (420–490nm) in recovered spectral profiles, matching with the standard absorption spectrum provided by manufacturer as plotted in the same Fig. Subsequently, a wider window is applied at z=62µm in Fig. 3(a) and 57µm in Fig. 3(b) to cover the second and the third peaks during STFT. The axial and spectral resolutions of the second window are about 6µm and 8nm respectively. Figure 3(d) shows the corresponding extracted oscillating spectra from both experimental and simulated data, and they fit well with backscattering spectrum calculated from a sphere of the same size by Mie theory. This extracted spectrum is dependent on the size of the sphere. Some discrepancy was observed around 500–530nm between Mie backscattering and extracted spectra and this may be caused by the ignorance of the first peak when we applied the larger window.

 

Fig. 3. Simulated (a) and experimental (b) OCT image of 6µm yellow-green fluorescent microspheres and the extracted spectra by two different spatial Gaussian windows (c–d).

Download Full Size | PDF

Assisted by numerical simulations, we can provide explanation on the origination of the intensity pattern in the OCT image as shown in Fig. 3(a) and Fig. 3(b). It has been previously hypothesized that the spectral oscillation measured at the center of the microscopic image from a spherical particle is produced by the interference of waves reflected by the front and back surface of the sphere similar to signal backscattered by a slab structure of same material and with thickness identical to the diameter of the sphere[23]. And a ring pattern around the center is also found in Ref. [23], which has a spectral oscillation matching with backscattering spectrum calculated by Mie theory. In order to verify above hypothesis, we numerically synthesized the en face microscopic image as well as OCT image from a single 6µm sphere for a range of NA’s of the imaging system. The en face microscopic image is obtained by forming intensity distribution in the image plane e as I (x,y) I without adding reference E-field while the sphere is in focus. The OCT image is simulated in the same manner as in Fig. 3(a). Figure 4(a) and (d) shows the microscopic and OCT image in the same NA as objective lens used in experimental setup. We find similar pattern of multiple peaks in Fig. 4(d) with peaks #1 and #2 centering at z=50µm, corresponding to the center of the microsphere in our simulation setup. Given that the previous study in Ref. [23] used a high NA objective lens, we have also conducted simulation for NA=0.2 and 0.3 to resolve the structure with higher lateral resolution. With NA=0.2 objective lens, the en face microscopic image in Fig. 4(b) displays a center and ring pattern similar to that has been observed previously[23]. In Fig. 4(e), all peaks in OCT image keep their axial position unchanged but in the lateral direction both peak #3 and #4 split from center residing at the edge of the microsphere; while peak #1 and #2 remain at the center. In NA=0.3 case as shown in Fig. 4(f), the lateral resolution is further enhanced and four peaks are clearly separated. By carefully inspecting the positions of the peaks in Fig. 4(f), we observed that peak #1 and #2 are located at (0, 0, 47µm) and (0, 0, 54.14µm) respectively with distance of 7.14µm. Accounting the polystyrene refractive index of 1.19 comparing to water, the physical distance between peak #1 and #2 is 6µm. With prior setup in simulation of the microsphere (6µm located at z=50µm), the peak #1 and #2 in Fig. 4(f) are situated right at the front and back surface of the sphere, which confirms the slab model proposed in [23]; and the lateral position of peak #3 is 3µm, which coincides with the radium of the sphere. Comparing the microscopic and OCT image, it is apparent that the peak #3 and #4 in Fig. 4(f) correspond to the ring in Fig. 4(c) and we speculate that it is due to a complex photon pathway traveling the periphery of the sphere. The equivalent optical pathlength of peak #3 is longer than that of peak #2, which makes peak #3 appear deeper. In NA=0.1 case, the lateral resolution of objective lens is insufficient and the peak #3 in Fig. 4(f) are merged with peak #2, together creating a complex pattern (peak #2, #3 and #4) in Fig. 4(d). The fact that the extracted SOCT spectrum in Fig. 3(d) from this complex pattern has a good match with Mie backscattering verifies the previous hypothesis in Ref. [23] that the ring pattern provides dominant spectrum contribution to far-field Mie scattering quantities.

 

Fig. 4. Simulated en face microscopic (a–c) and OCT images (e–f) of 6µm microsphere for an objective lens with NA=0.1, NA=0.2 and NA=0.3 respectively.

Download Full Size | PDF

Our second phantom is consisted of a mixture of two types microspheres (6µm yellowgreen, 3µm polystyrene) embedded in 1% intralipid gel solution, which provide similar scattering coefficient as biological tissue (g~0.8, µ’_s~10cm^-1)[24]. The phantom preparation followed the same procedure as described in above section except that two types of microsphere and intralipid solutions are added in heated gel solution. The goal of this phantom study is to test the feasibility of locating and identifying contrast enhancing dyes when they are concentrated in micro-meter sized scatters embedded in a turbid medium. Here, we use the bandwidth of 450–550nm, covering the ascending section from yellowgreen fluorophore absorption spectrum. Based on our previous analysis, we next develop an algorithm to discriminate the fluorescent spheres apart from non-absorbing ones. This algorithm includes a sequence of following steps. (i). We locate the image intensity peak positions by setting a threshold as a percentage of the maximum value; (ii). We place a Gaussian window set with axial resolution of 2µm on selected peak positions and perform STFT to recover spectra; (iii). We apply least-square-fit (LSF)[17] and attribute the extracted spectra to linear combinations of standard absorbent spectra from two microsphere types. For the polystyrene sphere, we use a flat line as standard. Then we have two linear coefficient maps representing the spectral contributions from two types of spheres separately. The algorithm is illustrated in Fig. 5. 3D OCT volume is obtained as a stack of x-z B scan images by scanning the sample alone y direction. Another two 3D dataset representing linear coefficients for yellowgreen and polystyrene microspheres are then calculated by our algorithm. We use 3D-YG and 3D-PS short for these two 3D dataset in later description. As shown in Fig. 5(a), each slice shows one original B scan OCT images and the 3D objects are reconstructed in light gray. 3D-YG is then superimposed in the volume in green color representing yellowgreen microspheres positions. The en face projection in x-y plane of all the spheres is plotted on the bottom of the volume, which is re-plotted in Fig. 5(b) for later side-by-side comparison. To verify our algorithm, an en face fluorescence image from yellowgreen microspheres was captured prior to OCT scanning and is plotted in Fig. 5(c). Figure 5(d) and (e) are the en face OCT projection of 3D-YG and 3D-PS. In Fig. 5(c), two yellowgreen spheres pointed by blue arrow are identified in Fig. 5(d) and not present in Fig. 5(e). And inside the red circle there are in fact one yellowgreen and two polystyrene spheres spaced closely, judging from Fig. 5(c) and (d). Our algorithm is also able to differentiate the yellowgreen one from the other two as shown in Fig. 5(d) and Fig. (e).

 

Fig. 5. (a). (Media 1) 3D OCT image of the scattering phantom containing yellowgreen and polystyrene microspheres; volume=125*85*100µm in x-y-z. See the animation of scanning and 360 degree view in Movie1.mpg (3.5M). (b) En face projection images on x-y plane of original 3D OCT, (d) 3D-YG and (e) 3D-POLY data. (c) En face fluorescence image at the same position taken prior to OCT scanning. Bar=20µm.

Download Full Size | PDF

Figure 6 shows a close look at one x-z B scan OCT image and demonstrated the effectiveness of our algorithm. Figure 6(a) is the original B scan OCT image taken from the yellow line in Fig. 5(b) and we developed a false-color coded map in Fig. 6(b) in HSV color space to identify the positions of two different kinds of microspheres. We encoded Hue with a resulting matrix of subtracting yellowgreen coefficient map from polystyrene one and encoded Saturation by the result of summing two coefficient maps. And the Value is encoded by the original OCT intensity map as Fig. 6(a). The color-coded image shows that there is one yellowgreen sphere at about x=140µm and three polystyrene ones at around x=105µm, 125µm and 150µm. The discrimination of yellowgreen microsphere can be verified by the fluorescence image in Fig. 5(c). Thanks to our color-coded Figs, we gathered more information based on their spectral behavior and managed to differentiate spheres with their characteristic absorbing feature.

 

Fig. 6. (a). Intensity based OCT image of the scattering phantom containing fluorescent and non-absorbing microsphere; (b). False-color coded SOCT image for discrimination; Gray scale colorbar is designated for the dB intensity and Hue colorbar encodes the result of subtracting YG coefficient map from polystyrene one.

Download Full Size | PDF

4. Discussion and conclusion

Biological tissue scatterers are usually heterogeneous, nonspherical, and organized in a complex manner. This work represents a starting point for understanding SOCT data measured from localized scattering and absorbing structures. It demonstrates that in order to interpret SOCT measurements in details for even relatively simple scatterers, a thorough understanding of the full-wave optical response of the scatterer is required. Numerical techniques that solve the Maxwell’s equation for arbitrary geometry are useful for investing and understanding OCT images of a more complex structure.

Several algorithms have been developed to perform time-frequency analysis of SOCT data [1, 7, 9, 25–28]. This work adopts a straight-forward IFT-STFT method to extract spectra from a specified spatial window. Our investigation shows that in order to extract absorption spectral profile from localized structures, it is important to use a small window size in the IFT-STFT algorithm. Spectra selected by a bigger-sized window are very likely corrupted by complex scattering spectra, originating not only from the structural properties of the scatterers themselves, but also from coherent multiple scattering effects which create speckle in image and high frequency disturbance on recovered spectra.

In conclusion, we have demonstrated the feasibility of characterizing the absorption and scattering spectra of micro-scale structures in a turbid medium using a SOCT system. We have also presented an analysis using numerical OCT simulations based on full-wave solutions of the Maxwell’s Equations to elucidate the origin of the multiple peaks in the OCT image for a single microsphere. Finally we have demonstrated the possibility of identifying contrast agents concentrated in micron-sized scale in a 3D SOCT image.