In vivo label-free lymphangiography of cutaneous lymphatic vessels in human burn scars using optical coherence tomography

: We present an automated, label-free method for lymphangiography of cutaneous lymphatic vessels in humans in vivo using optical coherence tomography (OCT). This method corrects for the variation in OCT signal due to the confocal function and sensitivity fall-off of a spectral-domain OCT system and utilizes a single-scattering model to compensate for A-scan signal attenuation to enable reliable thresholding of lymphatic vessels. A segment-joining algorithm is then incorporated into the method to mitigate partial-volume effects with small vessels. The lymphatic vessel images are augmented with images of the blood vessel network, acquired from the speckle decorrelation with additional weighting to differentiate blood vessels from the observed high decorrelation in lymphatic vessels. We demonstrate the method with longitudinal scans of human burn scar patients undergoing ablative fractional laser treatment, showing the visualization of the cutaneous


Introduction
Lymphatic vessels are an important part of the human circulatory and immune systems, circulating a fluid called lymph to collect excess interstitial fluid, transport lymphocytes and other white blood cells, and clear cellular debris. Several physiological and pathological processes can alter these vessels, such as the widening of the vessel lumen due to skin inflammation [1] in cutaneous wound healing and scarring; and the inadequate transport of lymph in lymphedema following cancer surgery [2]. To provide insights into such physiologies and pathologies, techniques have been developed to image the lymphatic vessels, referred to as lymphangiography. These include blue dye staining, fluorescence optical imaging, computed tomography, magnetic resonance imaging, ultrasound imaging, photoacoustic imaging and positron emission tomography [3][4][5][6]. Such approaches typically require the injection of exogenous contrast agents to label the lymphatic vessels [3], which raises issues of cost, patient safety and waiting time for in vivo imaging in humans.
Lymph is largely transparent to light in the visible and near-infrared spectrum [7], in contrast to the high degree of turbidity of the surrounding soft tissues. This characteristic of the optical properties of lymph provides a novel mechanism for imaging of the lymphatic system. In particular, it provides a useful contrast mechanism that can be exploited in conjunction with optical coherence tomography (OCT) [8]. OCT uses coherent, near-infrared light to acquire high-resolution (1-20 μm) volumetric scans of tissue microstructure. Extensions of OCT have also been developed to image other aspects of the circulatory system, specifically the blood microvasculature [7,9], such as in skin tissue [10][11][12], typically to a depth of 0.5-1 mm.
Preliminary work has explored the potential of OCT to image lymphatic vessels, most commonly based on intensity thresholding of the raw OCT signal [13]. However, imageprocessing approaches based on simple thresholding of OCT signal are subject to a range of artifacts that limit their usefulness. The OCT signal undergoes strong attenuation with depth [14] and depth-dependent variation caused by focusing optics and sensitivity fall-off, which reduces the accuracy of methods relying on a constant threshold. Alternative approaches have explored the use of adaptive thresholds [15], where the intensity values within a pre-defined window are used to calculate an optimal local threshold value, assuming a mix of both lymphatic and non-lymphatic pixels. However, the choice of the window size enforces assumptions about the vessel size, limiting the range of vessels that may be detected. More sophisticated approaches have adopted a multistage Hessian approach [16], utilizing higherorder shape information. This approach makes assumptions of the local appearance of vessels; calculating the second-order spatial derivative of the OCT signal within a small window and identifying regions containing a single long, thin region of low signal intensity, indicated by the relative sizes of eigenvalues of the Hessian matrix. Such approaches can be problematic at the intersection of vessels or in regions with multiple vessels within the computation window, where the assumptions of the local shape are no longer valid. While these methods have been primarily demonstrated on animal models, Baran et al. recently presented the first demonstration on human skin [17].
In this paper, we present a novel automated segmentation method to extract the threedimensional structure of the lymphatic vessels in OCT scans. The method corrects for variations in the OCT signal due to the confocal function and sensitivity fall-off with depth of a spectral-domain OCT system, and utilizes a single-scattering model of the OCT signal for compensation of attenuation for each A-scan, enabling reliable thresholding. A segmentjoining algorithm is then used to mitigate partial-volume effects associated with small vessels. These lymphatic images are augmented with OCT-acquired images of the surrounding blood vessel microvasculature, using speckle decorrelation with additional removal of the observed, localized decorrelation in lymphatic vessels. For the first time, we demonstrate the potential of OCT to acquire longitudinal measurements of the cutaneous lymphatic vessels. Results are presented showing the vessel networks in two burns patients undergoing ablative fractional laser treatment over a period of eighteen weeks. Longitudinal monitoring of the lymphatic system in human skin may provide a range of opportunities from new knowledge of healing processes to optimization of treatment regimes.

Characteristic optical signature of lymphatic vessels
Lymphatic vessels are largely transparent, resulting in low OCT signal. This is shown in the en face OCT image in Fig. 1(a), acquired from a region of a cutaneous burn scar on a human subject. The characteristically low OCT signal within the lymphatic vessel lumen is typically comparable to the background noise level, as shown in a representative averaged A-scan in Fig. 1(c) of the lymphatic vessel indicated by arrow (c) in Fig. 1(a). This low signal is expected due to the high water content (94%) and low scatterer concentration (6% solids) [18] present in lymph. In contrast, blood has a much higher portion of particles (largely blood cells, >45%) [18], with these cells behaving as exceptionally strong forward scatterers. This forward scattering in blood vessels can result in a lower OCT signal than surrounding solid tissue [ Fig. 1(b) from location labelled by arrow (b) in Fig. 1(a), identified from the speckle decorrelation signal], but we have found this value to generally be notably greater than the noise floor in our OCT scans.  [19]. (b) and (c) Averaged A-scan, respectively, from a blood and lymphatic vessel region labelled by the blue arrows in (a). (d) The attenuation-compensated signal of the data in the dashed rectangle in (c) after smoothing in the cross-sectional plane to reduce speckle and intensity compensation using calibration scans from a homogenous phantom, as described in Section 3. Scale bar in (a): 500 µm.
Another characteristic that differentiates both lymphatic and blood vessels from solid tissues is speckle decorrelation with time. In blood vessels, this is due to the movement of scatterers, largely red blood cells, while in lymph vessels we largely attribute it to the low signal-to-noise ratio (SNR) in the transparent lymph. Examples from the OCT scan in Fig. 1 are shown in Fig. 2. We calculated the speckle decorrelation of co-located B-scans acquired at sequential time points, t and ' t , using a previously published method [11,20] as: where ( , , , ) i x y z t is the amplitude of the complex OCT signal at lateral location (x, y), depth location z, and time point t ; and the speckle decorrelation ( , , ) D x y z is calculated in a widow defined by M N × pixels in the cross-sectional (x-z) plane. We show both an OCT B-scan in Fig. 2(a), with magnified inset in Fig. 2(b), and the corresponding speckle decorrelation image in Fig. 2(c), with inset in Fig. 2(d). The region corresponding to the inset is delineated by a dotted rectangle in Figs. 2(a) and 2(c). The blood vessels (e.g., label B1) show high speckle decorrelation in Figs. 2(c) and 2(d). Characteristically, there is also high decorrelation at depths below the blood vessel due to the high scattering nature of blood, as described in Ref [21]. Visually, this artifact appears as axial 'tails' of decorrelation below the blood vessels, such as that shown in the inset of Fig. 2(d). Conversely, this artifact is not present in lymphatic vessels, where the speckle decorrelation is localized to the area of low SNR within the vessel lumen. Representative lymphatic vessels are marked as L1 and L2 in Fig. 2, and the distinction between blood vessel B1 and lymphatic vessel L1 is highlighted in the magnified insets in Figs. 2(b) and 2(d). This feature is implicitly incorporated into the algorithm presented in this paper to aid in the distinction between lymphatic vessels and blood vessels. In brief, by weighting the speckle decorrelation with the co-located OCT intensity value, static tissue and low SNR lymphatic vessel areas are characterized by low values of weighted speckle decorrelation, whilst areas of blood flow present a high value for the weighted speckle decorrelation.

Correction of lymphatic vessel threshold values versus depth
To correct for variations in OCT SNR with axial depth, we have utilized a single-scattering model of the near-infrared light [14,22,23] with the A-scan expressed as: where ( ) i z is the amplitude of the complex OCT A-scan signal; and represents the single-scattering component with ρ and t μ being the initial value of the reflectance and the attenuation coefficient of the tissue, respectively. ( ) F z and ( ) S z are the confocal function and sensitivity fall-off function, representing the depth-dependent modulation of the OCT signal, respectively, due to the confocal setup of the focusing optics and the finite spectral resolution of the detector in spectral-domain OCT. This model, after removal of the effects of ( ) F z and ( ) S z , shows a linear decay in the logarithmic OCT signal, assuming sufficient homogeneity in the tissue.
As seen in Fig. 1(c), lymphatic vessels can cause local fluctuations in the OCT signal, but do not significantly alter the overall linear decay of the logarithmic signal caused by the background tissue due to their transparency. Therefore, our proposed algorithm uses this single-scattering model to estimate and compensate the OCT threshold values with depth when calculating the appropriate signal threshold value to identify the lymphatic vessels. This use of a model of the underlying optics is distinct from purely image-based approaches, such as adaptive thresholding [15], which are based on assumptions of the distribution of lymphatic and non-lymphatic pixels within a moving window. Our use of a tissue opticsbased approach has enabled us to utilize a large proportion of each A-scan in the calculation of the appropriate attenuation correction factor at each lateral (x, y) location. Figure 3(a) summarizes the data processing flow of the method. Prior to the estimation and compensation of attenuation, ( ) F z and ( ) S z are corrected using a calibration phantom. This phantom comprises polystyrene microspheres (diameter: 0.5 µm) diluted in deionized water as a homogeneous, low-scattering target. From the size and refractive index contrast of the spherical scatterers, the scattering properties are estimated from Mie theory [24]. A representative A-scan [ Fig. 3(b)] is computed by averaging intensity values in the lateral direction, and used to correct for the confocal function and sensitivity fall-off within the sample. This normalization A-scan 0 ( ) i z from the phantom is modelled using the singlescattering model presented earlier in Eq. (2) as:

Methods
where 0 ρ and 0 t μ are the initial value of the reflectance and the attenuation coefficient of the phantom. Using the approach described in Refs [22,23], the correction for the confocal function and sensitivity fall-off with depth is performed by dividing the tissue A-scan by the phantom A-scan, leading to the corrected OCT signal as: where a is the logarithm of the ratio of the proportionality in Eqs. (2) and (3). The corrected measured OCT signal, as shown in Eq. (4), is then used to calculate the attenuation coefficient at each lateral (x, y) location. We use weighted least-squares regression to fit the singlescattering model over the axial range, with the weight at each axial location being the corresponding logarithmic OCT signal. This weighting minimizes the impact of the low-SNR lymphatic vessel regions on the calculation of the signal attenuation.
To allow a single threshold to be applied across all axial depths, the algorithm then compensates the corrected OCT signal using the estimated attenuation coefficient, adding a depth-dependent compensation function to the corrected intensity, as shown in the equation:   Fig. 1(d). The lymphatic vessels are then segmented in 3-D by applying a single threshold across the entire volumetric C-scan. Partial-volume effects [25] will result in breaks in small lymphatic vessels, particularly when their diameters are comparable in scale to the imaging resolution, such as the small vessel indicated by the red arrows in Figs. 3(c) and 3(d). To partially mitigate this artifact, we have augmented the algorithm with post-processing, using a segment-joining algorithm. Lymphatic vessels are identified to a depth of 400 µm below the tissue surface for each scan, forming a 2-D projection image in the en face (x, y) plane. This depth range is chosen to avoid possible artifacts in deeper regions caused by low OCT SNR. A skeletonization algorithm is used to extract the center line of each vessel [26]. The branch points of the center lines are then identified by detecting connected center lines using a standard morphological operation implemented in Matlab (bwmorph). These are then removed to decompose the vessel network into separate segments. The two endpoints of each segment are identified, and a local estimate of vessel orientation at each endpoint is computed by finding the line of best fit to points within a small spatial window. The axial depth of each endpoint is also recorded from the 3-D volume. Nearby segments are then joined together where corresponding endpoints are within a threshold distance with similar orientation (to within an empirically chosen angular threshold) and depth. The same threshold values are maintained across all data sets. The missing pixels are linearly interpolated between matching endpoints, resulting in a single, continuous lymphatic vessel, as shown in Fig. 3(e). Finally, the area density of the resulting lymphatic vascular image is automatically quantified as the ratio of the total lymphatic vessel area (in the x, y direction) to the total transverse tissue area.
A separate data processing flow [red boxes in Fig. 3(a)] was incorporated to provide a complementary, co-located image of the blood vessels. This is based on the speckle decorrelation method using Eq. (1), as described in [20]. Briefly, co-located OCT B-scans were used to calculate the speckle decorrelation, in which the blood and lymphatic vessels show high decorrelation. As discussed earlier, blood vessels characteristically exhibit a higher SNR than lymph vessels. Thus, speckle decorrelation values were weighted by the corresponding averaged, thresholded OCT signal value, resulting in high weighted speckle decorrelation values in blood vessels and low weighted speckle decorrelation values in lymphatic vessels and static tissue. A 2-D projection image of the blood vessels was then generated by finding the maximum value of the weighted speckle decorrelation within 400 μm of the skin surface, matching the depth range for projecting the lymphatic vessels and allowing for comparison of the two networks in the same OCT scans.

Experiment
A spectral-domain OCT scanner (TELESTO II, Thorlabs Inc., NJ, USA) was used in this study to collect skin OCT scans. The scanner has lateral and axial resolutions of 13 μm and 5.5 μm (in air), respectively, with a center wavelength of 1300 nm. OCT C-scans were acquired over volumes of 6 × 3 × 3.6 mm or 6 × 4.5 × 3.6 mm (x × y × z) with a constant sampling voxel size of 5.9 × 7.5 × 3.5 µm. The smaller volume was used in some of the early scans using this protocol. Note that the quantification by area density factors out differences in the volume scanned. A pair of co-located B-scans were acquired at each lateral location to facilitate speckle decorrelation calculations.
OCT scanning of burn scar patients was approved by the Human Research Ethics Committee of Royal Perth Hospital and The University of Western Australia. Two burn scar patients (Patient 1 and 2) were recruited for this study with written consent. Both Patient 1 (age: 55 years; gender: male; cause of burns: explosion; age of scar: 6 months; scan location: abdomen) and Patient 2 (age: 23 years; gender: male; cause of burns: explosion; age of scar: 10 months; scan location: abdomen) had full-thickness burns covering the entire abdomen. They received laser treatment at Royal Perth Hospital using the same protocol as described in Ref [20]. In brief, they received three courses of treatments using an ablative fractional CO 2 laser (Ultra-Pulse ® , Lumenis Inc., CA, USA) at an interval of ~6 weeks between two sequential treatments, as shown by the green bars in Fig. 4. For each patient, a 5 × 10 cm confluent scar area was treated with a microscopic treatment zone spot size of 0.12 mm at a beam density of 5%. OCT scans were collected from the laser-treated region and an adjacent untreated region at multiple time points before, during and after the treatment, as indicated in Fig. 4. To mitigate possible motion artifact, a thin, 10 mm square metal marker with a center round hole (diameter: 5 mm) was affixed to the scar during scanning, allowing the use of an image registration algorithm [27]. A customized imaging setup described in Ref [20] was used for data acquisition with ultrasound gel applied on the scar surface to reduce the backreflection from the air-skin interface by refractive index matching. The scanning protocol described in Ref [20] was adopted to acquire co-located OCT scans at multiple time points using the physical landmarks on the scar surface as markers to locate the same region in successive scans. The longitudinally-acquired OCT scans were processed off-line and manually correlated to track changes in lymphatic vessels.

Results
Figure 1(a) shows an OCT scan of the immature hypertrophic scar of Patient 1, acquired between the second and third sessions of laser treatment (i.e., time point D in Fig. 4). The corresponding lymphatic and blood vessel images are shown in Figs. 5(b) and 5(c), respectively, projecting all vessels encountered to a depth of 400 µm below the scar surface. The scanned area (6 × 4.5 mm) is from the center of the blue square (10 × 10 mm) shown in the scar photo of Fig. 5(a). The dark lymphatic vessels in the OCT image in Fig. 1(a) are readily identified in Fig. 5(b). A 3-D rendering of the identified lymphatic vessels is also presented in Visualization 1, showing the volumetric distribution of the network. We note that there are artifact areas of low OCT signal in Fig. 1(a), such as those indicated by the white arrows, which are caused by occlusion from the protruding adhesive tape for attaching the metal marker to the scar surface. The single-scattering, signal intensity correction algorithm presented in this paper corrected for these artifacts in Fig. 5(b). Without such a correction, these artificial low-signal regions could be erroneously marked as lymphatic vessels using more simplistic thresholding approaches. The blood vessel image, Fig. 5(c), shows a profuse vascular network, corresponding to the red appearance of the scar in Fig. 5(a). The superimposed image of the lymphatic and blood vessels in Fig. 5(d) visualizes both networks and enables their comparison. Whilst the lymphatic vessel network comprises a sparse distribution of large vessels, the blood vessels show a dense distribution of small vessels with more branching. This results in a measured area density of 11% and 49%, respectively, for the detected lymphatic and blood vessel networks.  Figure 6 presents three time points in our longitudinal study on Patient 2, imaged at time points C (between the first and second laser treatment), D (between the second and the third laser treatment) and E (~6 weeks after the final third treatment). To provide an en face OCT image indicative of the lymphatic vessels, we show a minimum projection image over a depth of 200 to 300 µm below skin surface in Figs. 6(a), 6(b) and 6(c) for the three sequential time points C, D and E, with lymphatic vessels displaying characteristic low signal. The depth range for projection was chosen to optimize the contrast of the lymphatic vessels in these images. Note that these images are illustrative, intended only to provide an indication of the appearance of lymphatic vessels in the OCT data, and to help the reader to assess the results of automated lymphatic vessel detection. Figures 6(d), 6(e) and 6(f) show the automatically detected lymphatic vessel images (animations in Visualization 2, Visualization 3, and Visualization 4), and Figs. 6(g), 6(h) and 6(i) show the combined lymphatic and blood vessel images to a depth of 400 µm below skin surface, respectively, from time points C, D and E. The arrows in Figs. 6(d), 6(e) and 6(f) mark vessel patterns common across the scans. We calculated an initial increase in the lymphatic vessel area density from 6.0% at time C to 7.3% at time D, followed by a decrease to 5.4% at time E, corresponding to relative change of, respectively, +22% and −26%. However, the blood vessel network shows a continuous regression with the area density dropping from 59% at time C to 31% at time D to 24% at time E, which agrees with results reported in [20] on the longitudinal monitoring of vasculature in the laser-treated burn scars.
The lymphatic vessel network in Patients 1 and 2 was imaged and quantified by our method at multiple time points throughout the laser treatment regime. Longitudinal changes in the area density values for lymphatic vessels are plotted in Fig. 7, with the time of treatment marked by the green bars on the horizontal axis. We note an acute increase in lymphatic area density after the first laser ablation, from 7.8% to 9.6% and 2.8% to 4% (time point A to B) for Patient 1 and 2 at 3 days after the first treatment. The area density then shows a consistent increase as subsequent treatments were performed, followed by a regression after the final treatment.

Discussion and conclusion
The characteristic transparency of lymph to near-infrared light allows the automated detection of lymphatic vessels in OCT scans without the need for exogenous contrast agents. This transparency distinguishes the lymphatic vessels as regions with markedly low SNR in the OCT signal. However, the segmentation of these vessels in the OCT scans is confounded by the variation in OCT signal due to the increasing signal attenuation with depth. In this study, we present a method to segment the lymphatic vessels by thresholding the attenuationcompensated signal using the single-scattering model, aided by the correction of the confocal function of the focusing optics and sensitivity fall-off of the OCT scanner. In contrast, the method in [17] compensates the attenuation at each pixel as a function of the summed intensities of all pixels within the A-scan located at greater depths, which has been shown to be advantageous for heterogeneous tissue, but lacks a mechanism to account for the confocal function and sensitivity fall-off. Whilst sensitivity and spatial resolution of the OCT scanner will limit the ability to detect very small lymphatic vessels, our approach utilizes consistent threshold values across all data sets, reducing variation between measurements in a longitudinal setting. In addition, weighted speckle decorrelation images were calculated, showing blood vessels whilst suppressing lymphatic vessels due to the low SNR in lymph. We believe that this forms a useful approach to visualize the lymphatic and blood vessel networks simultaneously, allowing characterization of both vascular networks.
The demonstration of our method on burn scar patients undergoing ablative fractional CO 2 laser treatment provides a proof-of-concept investigation for longitudinal assessment of changes in cutaneous tissue. The resulting lymphatic vessel images enable the visualization of the individual vessels with diameters of ~30-150 µm and the interconnection of the vessels to form a network. Previous studies have reported the diameters of the cutaneous lymphatic capillaries (37-76 µm) [28,29], which are broadly in agreement with our results. We hypothesize that the larger vessels identified in our images may be dilated in response to inflammation [1] attributed to the wound healing process [30]. The laser treatment in this study produced an array of localized regions of ablation to promote remodeling of the scar tissue, which induces inflammation additional to the underlying wound healing process. Our results suggest an increase in lymphatic vessel area density as multiple laser treatments were performed, and regression after the cessation of treatment. This finding is supported by the literature, as lymphatic vessels have been reported to have an enlarged diameter in inflamed skin [1].
Methods, including ours, for segmenting lymphatic vessels from OCT scans largely rely on the low signal in vessel lumens due to their transparency [7,[15][16][17]. When applied to various biological tissues, the reliability of this group of methods might be impacted by other structures with low backscattering. An example of a confounding structure is hair follicles in normal skin tissue, typically extending from the skin surface to the deep regions in an OCT scan with gradually increasing size with depth. They can be erroneously identified as lymphatic vessels due to their low OCT signal within the high-backscattering background tissue. Further image processing could be added to distinguish such artifacts from vessels using the shape and directionality of the structures, for example, by 3-D tracing. Alternatively, the use of OCT in combination with contrast agents provides another possible method to mitigate such artifacts [31]. In addition, the OCT signal of blood vessels in other types of tissue might be lower than that observed in this study of cutaneous tissue. This can lead to the misidentification of a blood vessel as a lymphatic vessel by intensity-based methods. Additional differentiation of vessels would be required, potentially using the characteristic axial 'tails' of decorrelated speckle found below blood vessels, as noted in Ref [21].
A limitation of our proposed approach is that it does not distinguish between occluded or collapsed lymphatic vessels, and areas that are devoid of lymphatic vasculature. Thus, the approach is unable to quantify lymphatic changes where the vessels are not patent. Instead, it provides a measure of the overall density of patent lymphatic vessels, via the density of low SNR pixels in the OCT scan. We note that our use of a post-processing algorithm to join vessel segments that have been separated because of partial-volume effects may erroneously remove the natural breaks in some lymphatic vessels. As this correction is most commonly applied to small vessels, the impact on estimates of vessel density is expected to be minimal. Other limitations of this study include the use of skin contact during OCT scanning, where the OCT probe was placed on and attached to the scar surface to reduce motion artifact. The resulting contacting pressure on the scar tissue could collapse some lymphatic vessels, which are known to have thinner walls than blood vessels [32]. To mitigate variation from these possible artifacts, we utilized the same operator to acquire all longitudinal scans, with the OCT probe mounted on an articulated arm and held gently on the skin. Future approaches could include pressure measures to help maintain a uniform pressure across multiple scan sessions. We note that lymphatic vessel imaging methods based on OCT have only been validated on animal models. A validation study in human skin by comparing OCT against other imaging techniques using contrast agents will be needed to further characterize the accuracy of OCT-based imaging of lymphatic vessels and the errors in the measurement of vessel density. In addition, the incorporation of the measurement of lymph flow velocity [33] into the imaging of the vessel morphology would expand the utility of our method.
We note that the distribution of lymphatic vessels is heterogeneous. In clinical usage, we anticipate that robust estimation of vessel density would require multiple acquisitions distributed over the area of interest to allow for a statistical estimation of the vessel density.
In conclusion, we have presented a novel method based on OCT to image cutaneous lymphatic vessels in humans. This method removes the depth-dependent modulations in the OCT signal due to the confocal function and sensitivity fall-off and compensates for the signal attenuation with depth. This correction of the OCT signal allows a single threshold value to be applied across the entire C-scan for the detection of lymphatic vessels, followed by a segment-joining algorithm to mitigate partial-volume effects. Augmented vascular images were generated, simultaneously showing both lymphatic and blood vessel networks. We have demonstrated this protocol in a longitudinal study of human burn scars undergoing ablative fractional laser treatment. Our integrated approach for the assessment of lymphatic and blood vessels holds promise for longitudinal monitoring of pathologies and their response to treatments.