OCT-based label-free in vivo lymphangiography within human skin and areola

Due to the limitations of current imaging techniques, visualization of lymphatic capillaries within tissue in vivo has been challenging. Here, we present a label-free high resolution optical coherence tomography (OCT) based lymphangiography (OLAG) within human skin in vivo. OLAG enables rapid (~seconds) mapping of lymphatic networks, along with blood vessel networks, over 8 mm x 8 mm of human skin and 5 mm x 5 mm of human areola. Moreover, lymphatic system’s response to inflammation within human skin is monitored throughout an acne lesion development over 7 days. The demonstrated results promise OLAG as a revolutionary tool in the clinical research and treatment of patients with pathologic conditions such as cancer, diabetes, and autoimmune diseases.


Supplementary Figure 2 | Swept-source OCT system
Portable bench top Swept-source OCT (SS-OCT) system with an MEMS-tunable vertical cavity surfaceemitting laser (VCSEL) source engine is used. This VCSEL source (Thorlabs Inc.) is able to sweep the wavelength of the laser across a broad spectral range near 1310 nm at a fixed repetition rate of 100 kHz, giving a long coherence length (>50 mm) with ~12 µm axial resolution in tissue and 2 mm image range in human skin. The spectral bandwidth of laser output was measured ~67 nm at -3 dB. A Mach-Zehnder interferometer (MZI) built in the laser module provided real-time optical clocking, from which the MZI signal was detected by a dual-balanced detector connected to an external clock input of a 12 bit 500 MS/s data acquisition (DAQ) card. The measured sensitivity was 105 dB, almost constant up to depth of 4.25 mm in air. The sample arm consists of a probe imaging head (including galvanometric scanners) and a 5X objective lens (LSM03, Thorlabs Inc.), giving 22 µm lateral resolution.

Supplementary Figure 3 | Mouse ear imaging using proposed OCT methods
Imaging results of simultaneous optical microangiography and lymphangiography in the mouse ear after 5-minute depilation. Fig. 3a shows enface maximum projection view of the 3D OMAG image, Fig. 3b show sMIP of the 3D OCT lymphangiography. Fig. 3c is the overlay of Fig. 3a-b. Fig. 3d-f are the crosssectional views of corresponding data sets from a location delineated with the orange dashed line. Fig.  3d is the OMAG image, Fig. 3e is the OCT structure image where lymphatic vessels are shown with green dashed circles, and Fig. 3f is the overlay of Fig 3d-e. The lymphangiography image acquired using our proposed method matches with the earlier reported results in literature 5,6 , which are validated using blue dye injection. The field of view of images is 2 × 2 mm 2 .

Supplementary Methods for OCT based Lymphangiography
Accurate detection of lymphatic vessels within human skin requires OCT images with good contrast. Non-uniform intensity distribution of the OCT images due to light attenuation and heterogeneous tissue properties make it difficult to extract the lymphatic vessels from the surrounding tissue with high scattering. In order to mitigate this problem, OCT structural images are processed in axial (A-scan or Bscan based processing) and lateral directions (enface based processing) and the contrast is enhanced using the algorithms described below. Please see Supplementary Fig. 1 for the proposed processing protocol.

Attenuation compensation in A-scan
In the standard SS-OCT system, the output I of the interference between a light beam reflected from a reference mirror and that from a sample can be described by Choma et al. 7 [ ] = 1 2 [ ]( + + 2√ cos (2 ∆ )) where, k is the optical wavenumber, ρ is the detector responsivity, S[km] is the source spectral density, ∆x is the pathlength difference between the reference and the sample arms, and RR and RS are the reflectivity of the reference and the sample arm, respectively. The acquired signal ideally has values at M evenly spaced wavenumbers, so k= {k1,k2,…, km, …, kM}. The cross correlation term [ ] is derived after removing the constant offset terms from Eq. 1: Through a discrete inverse Fourier transform of Eq. 2, OCT A-scan, D, as a function of depth, z can be obtained: However, when light penetrates the tissue, it attenuates through the absorption (µa) and scattering (µs). This causes change in the acquired OCT A-scan signal by modifying Eq. 3 in its continuous form: where α is the backscattering coefficient, and D(µ) = µa +µs. When D(µ) is large at a certain depth, such as at a blood vessel, it shadows the remaining depth in that A-scan, DA. This effect contributes to the poor contrast between the tissue and the lymphatic vessels in the deeper sections. To compensate for the attenuation in the OCT images, this decay term can be removed by using the following operation 2, 8 : substituting Eq. 5 into Eq. 4 and discretising: The backscattering coefficient α can be adjusted depending upon tissue characteristics, and F is the final depth in the A-scan over which most of the light is attenuated. Note that this method does not require extracting attenuation coefficients µa and µs. Nevertheless, D[z] in Eq. 6 is the compensated reflection that allows the display of uniform images. This method assumes that the most of the light is attenuated within the recorded imaging depth range, and the backscattering coefficient is constant.

Histogram equalization of the en face cross sections
Due to the non-uniform mineral oil distribution topically applied onto the skin and the non-flat surface of the skin, the OCT image intensity is also non-uniform in enface cross sections. This situation is more severe when a large field of view is used such as in our case. Hence, additional OCT image enhancement step is needed before mapping the lymphatic vessels.
Histogram of an image is the probability distribution function of the occurrence of one specific grey level in that image. Histogram equalization is a contrast enhancement technique that increases the global contrast of images when the information content is represented by close contrast values. Histogram equalization method maps input intensities into an output image where its probability distribution function is uniform and its cumulative distribution function histogram is the same as the input image. This allows for stretching the contrast of the image to be distributed uniformly in the dynamic range 9 . However, it may produce unrealistic images and amplify background noise whilst decreasing the usable data. Histogram equalization has been generalized to local histogram equalization which divides the input image into smaller regions and applies histogram equalization algorithm to create a uniform distribution in each region 10 . Unfortunately, it creates noise amplification in homogenous regions.
To solve this problem, high peaks in the histograms, which correspond to homogenous regions in the image, are clipped and probability distribution function is non-uniformly reconstituted with Rayleigh distribution to enhance the foreground and supress the noise using equation below 4 .
= + √ 2 2 ln ( Here, y is the output intensity, x is the input intensity, ymin is the low bound, β is the Rayleigh parameter and Pinput(x) is the cumulative probability of the input image. A higher β value will result in more contrast enhancement in the image while at the same time increasing saturation and noise levels. This technique transforms non-uniform, low-contrast OCT images into a uniform, high-contrast images which become useful in the mapping of the dark areas (supposedly lymphatic vessels) in the OCT images to a final enface image.

Enface mapping of lymphatic vessels
In order to segment the dermis region in the OCT images, we first detect the first layer of epidermis by extracting the first peak of each A-line. Then, we skip the epidermal layer starting from this first peak until the dermal-epidermal junction is reached by using a pre-determined number of pixels. Finally, we save the pre-determined pixels correspond to the papillary and reticular dermal layers in each A-scan. Repeating this procedure for every A-line flattens the curvature in the initial volume, and the interested volume can be automatically segmented with accuracy without the use of an additional treatment.
After segmenting dermis, we sort the pixels in each A-scan in ascending order, based on their intensities. Then we average the first N number of pixels with minimum intensities before mapping on to an enface 2D plane. In summary, Zen_face [x, y] is calculated as: Here, Zsorted[x,y,i] is the value of the i th pixel in ascending order and M is the number of averaged pixels. We call this method sorted minimum intensity projection (sMIP). Final enface image is produced in inverted colors. The above procedures are illustrated in supplementary Fig. 1.