In vivo multifunctional optical coherence tomography at the periphery of the lungs.

Remodeling of tissue, such as airway smooth muscle (ASM) and extracellular matrix, is considered a key feature of airways disease. No clinically accepted diagnostic method is currently available to assess airway remodeling or the effect of treatment modalities such as bronchial thermoplasty in asthma, other than invasive airway biopsies. Optical coherence tomography (OCT) generates cross-sectional, near-histological images of airway segments and enables identification and quantification of airway wall layers based on light scattering properties only. In this study, we used a custom motorized OCT probe that combines standard and polarization sensitive OCT (PS-OCT) to visualize birefringent tissue in vivo in the airway wall of a patient with severe asthma in a minimally invasive manner. We used optic axis uniformity (OAxU) to highlight the presence of uniformly arranged fiber-like tissue, helping visualizing the abundance of ASM and connective tissue structures. Attenuation coefficient images of the airways are presented for the first time, showing superior architectural contrast compared to standard OCT images. A novel segmentation algorithm was developed to detect the surface of the endoscope sheath and the surface of the tissue. PS-OCT is an innovative imaging technique that holds promise to assess airway remodeling including ASM and connective tissue in a minimally invasive, real-time manner.

used PS-OCT to highlight the presence of fibrosis [16] and collagen [17]. OCT is attracting interest for the study of lung diseases due to the ability to acquire high-resolution crosssectional images during minimally invasive procedures [18].
In this paper, we demonstrate a multifunctional OCT platform that provides novel information on the airway wall using intensity, attenuation coefficient (AC), dynamic OCT (DOCT), and PS-OCT images with a distally scanning catheter in vivo, in humans. A novel automated segmentation algorithm based on the degree of polarization uniformity (DOPU) and intensity OCT images is presented. An algorithm to compensate for polarization mode dispersion (PMD) in catheter-based images improving PS-OCT imaging is described. A metric called optic axis uniformity (OAxU) is presented for highlighting fiber-like structures. The setup is demonstrated by the acquisition of images of the airways, scanning from distal to proximal, of a severe asthma patient that was scheduled for a recently-developed endoscopic severe asthma treatment called bronchial thermoplasty [19][20][21]. In this context, OCT, and in particular PS-OCT, seems to be a promising modality for minimally-invasive assessment of airway remodeling (e.g., ASM and collagen content) and the impact of bronchial thermoplasty on the airway wall structure [22].

PS-OCT imaging system
Figure 1(a) shows the diagram of the imaging system, which has been described in a previous publication by Li et al. [17]. Light from a swept-source laser with a central wavelength of 1310 nm, bandwidth of 90 nm, and a sweep rate of 50 kHz (Axsun Inc., Billerica MA, USA) is split into a sample and reference arm of a modified Mach-Zehnder interferometer. In the sample arm, a polarization delay unit (PDU) creates two depth-encoded polarization states, orthogonal in the reference frame of the lab [23]. A circulator relays the OCT light from the PDU to the catheter and finally to the polarization diversity detection module (PDDM). The transmission reference arm configuration avoids the use of a second circulator, which is known to exhibit considerable polarization distortions [24]. A fiber-based PDDM (Finisar Corp., Sunnyvale, CA, USA) recombines the sample and reference arm and splits the orthogonal polarization components into four output leads for balanced detection, synchronized by an optical A-line trigger and a built-in k-clock. Figure 1(b)-(c) present a drawing of the catheter used in this study that has a distally-scanning design, similar to the one presented by Li et al. [17]. A single mode fiber was angle-cleaved at 8° and glued to a 0.5 mm diameter custom GRIN lens with both facets polished at 8° and a working distance of 1.5 mm (GRINTECH GmbH, Jena, Germany). An alternating current micromotor, with a diameter of 1 mm, was built in-house with double wound 90 µm diameter copper wire coils around the motor housing. A micromagnet with an axle running through its longitudinal axis is encapsulated in the motor housing. The axle is held in place by two bearings that are part of the housing, allowing the magnet to rotate freely inside the housing. The alternated current running through the wires induces a variable magnetic field that actuates the magnet by aligning its poles to the field, in turn causing the rotation of the magnet-axle construct. The doubled wound coil increased the torque on the magnet compared to our earlier implementation [17] and resulted in a more reliable rotation. A 300 µm wide reflective microprism (Edmund Optics, Barrington NJ, USA) was glued on the motor axis to reflect the light to the tissue. The microprism was glued at an angle of 48° to avoid specular backreflections from the catheter sheath, which would saturate the detectors. The axial resolution of the OCT system is 12 µm in tissue (assuming a refractive index of 1.4). The fullwidth-at-half-maximum spot size of the beam exiting the catheter is 13 µm.

DOPU calculation
In an interferometer, only the common polarization component of reference and sample arm interfere, meaning that the measured fringes always correspond to fully polarized light. When light propagates through birefringent or scattering media, its polarization state changes. DOPU is a metric that measures the uniformity of the reflected polarization state in a portion of space. Several algorithms have been developed to calculate DOPU [25][26][27], with one of the most commonly used being: where I is the intensity of the OCT backscattering at a certain location in the cross-sectional scan, and Q, U, and V its Stokes vector components obtained from the measured electric fields. The brackets indicate spatial averaging, over a kernel, sized 3x7 pixels (axially and azimuthally, corresponding to 21 µm and 2.63°), which has the effect of blurring the image. The DOPU was calculated for each depth-multiplexed polarization state independently, and the two results were averaged to obtain the DOPU used in the following steps. In this article, DOPU was not found to carry significant diagnostic information, but it revealed to be an ideal starting point for segmenting the inner surface of the OCT catheter. Moreover, it provided a straightforward method to identify locations associated with the presence of alveoli. These structures have peculiar characteristics, namely increased backscattering at the tissue-to-air interfaces due to non-matching refractive index and weak backscattering due to the absence of scatterers within the alveoli. The DOPU in these areas is lower because of the increased contribution of noise in the low scattering areas. As demonstrated in Fig. 2(f), the areas within the orange and the cyan rectangles have similar OCT intensity. However, the corresponding DOPU image in Fig. 2(c) shows a considerably lower DOPU for the area in the orange rectangle, due to the presence of alveoli. A threshold value of DOPU of 0.7 was empirically found to isolate the alveoli from the image. This step was particularly useful to facilitate the retrieval of PMD compensation parameters, which relies on the presence of uniformly arranged fiber-like structures and the exclusion of areas subject to depolarization.

Segmentation of OCT endoscopic images
We developed a robust 3D automated algorithm to precisely segment the inner and outer edges of the endoscope sheath, the motor wires, and the surface of the airway lumen using both DOPU and intensity images. We used the segmented features to enable the PMD compensation and to improve the display of the cross sectional images.
The segmentation routine used Dijkstra's algorithm [23], which finds the shortest path linking nodes in a graph by searching the minimum value among the eight nodes neighboring a node. The graph was constructed starting from a DOPU or an intensity image, in both cases a cross section presented in polar coordinates.
Due to the blurring effect of the spatial averaging kernel, the catheter sheath structure appears as a continuous feature in DOPU images, which facilitates the detection of the inner edge of the catheter sheath. To find the inner edge of the catheter sheath in the first B-scan (n = 1), a 50 pixel thick region-of-interest (ROI) of the DOPU image is set by manual input of the upper boundary chosen by visual inspection. The ROI of the nth DOPU image is obtained by selecting a region extending 15 pixels upwards and downwards from the OCT catheter inner sheath edge from the nth-1 image.
In order to segment the inner edge of the catheter sheath, we computed the y-directional gradient of the DOPU image, where y-direction is defined as the depth direction. Secondly, we apply the shortest path search algorithm on the 50-pixel ROI of the y-directional gradient of the DOPU image where the catheter sheath inner surface is expected to be. To find the shortest path, the gradient image is transformed into a graph. The algorithm adds two extra nodes with zero cost before the first and after the last A-line of each image. These additional nodes act as a source and sink of the graph, therefore avoiding the manual endpoints selection. Once the minimum cost path is computed, the path is retrieved in the image domain, effectively segmenting the feature.
Once the inner edge of the catheter sheath is segmented, the bundled motor wires (which feed current to the motor placed at the tip of the catheter) can be localized by summing the DOPU values between the inner sheath surface (blue line in Fig. 2(c)) and the top edge of the image. The maximum value of this integral gives an approximation to the motor wire center location. The wire edges are then localized by finding the leftmost and rightmost local maxima of the DOPU integral image in a window of 300 pixels around the wire center location.
The outer edge of the catheter sheath and the upper surface of the lungs are not segmented using DOPU images because of excessive blurring. Instead, intensity images are used in the second part of the algorithm to exploit the full resolution of the images. The outer edge of the catheter sheath is found by selecting an ROI in the intensity image, which extends 30 pixels below the inner sheath surface (see Fig. 2(d)). The boundary of the outer edge of the catheter sheath is detected by using Dijkstra's algorithm on the gradient of the ROI of the intensity image. To find the tissue boundary, an ROI in the intensity images is selected, starting from one pixel below the outer surface of the catheter until the end of the imaging range. This ROI is normalized between its minimum and the maximum values. Subsequently, this ROI is denoised by applying an isotropic total variation (TV) algorithm [28], with a regularization parameter set to 50% of the ROI's standard deviation. Afterward, the denoised ROI is converted into a binary image using a threshold of 0.45. From each binary image, areas which contain zeros smaller than 50 pixels are set to ones. Finally, the lung boundaries were segmented by finding the first non-zero value from the top of the binary image for each Aline (see Fig. 2(f)). Examples of segmentation from three B-scans along the pullback are shown in Fig. 2(g)-(i), from distal to proximal locations. In Fig. 2(g) the lung tissue is almost always adherent to the outer catheter sheath, with some mucus visible next to the wires. Figure 2(h) shows a frame from the middle of the pullback, where the lumen is larger than the endoscope, and therefore the tissue appears detached. Some mucus is visible next to the wires, but the tissue is nonetheless correctly segmented. Figure 2(i) displays a B-scan from the end of the pullback where the lumen is larger than the OCT imaging range, and therefore it is not visible, apart for a faint aliased image.

Algorithm
The polarizat Some publica [29,30]  The first step of the PMD compensation algorithm is finding a unique calibration region, used to determine the PMD-induced distortions caused by the imaging system and correct for them in all the B-scans of the acquired volume. Three B-scans separated in time by 100 frames were selected from the second half (more proximal part) of each C-scan, because the distal region contained more alveoli, which are depolarizing. From each of these B-scans, three consecutive 100 A-lines-wide regions were selected starting from the edge of the motor wires' shadow, yielding 9 regions. The region with the highest quality parameter (described below) was chosen as the PMD calibration region for the whole C-scan.
Analogous to Braaf et al. [30], the detected E-field E out can be modeled in Jones formalism as follows: Here, the matrix [1,0;0,-1] represents the surface reflection at the catheter outer wall (or as in Ref [30]. the sample surface reflection). The influence of the mirror rotation is causing the surface state to change periodically and induces phase wrapping of the surface state at specific locations in the B-scan. Since the surface state was calculated over three adjacent segments of 100 A-lines, the surface state surf (k) E from a segment which did not show any phase wrapping was chosen. The inverse of the obtained surf (k) E is multiplied to the field measured from the sample to eliminate the matrix out (k) J : where the reflection matrix and the matrix ' S (k) J were included in the matrix S (k). J It is noteworthy that the matrices surrounding J s are canceling only for the same location on the catheter sheath, which was approximately achieved by limiting the analysis to a segment of 100 adjacent A-lines (corresponding to 37.5° of prism rotation). The left side of Eq. (4) becomes the identity matrix when the field E out (k) corresponds to the surface of the sample, which was required by the PMD compensation algorithm [30]. The phase retardation is then extracted from the differential Mueller matrix m. Assuming that the sample is homogeneous over a small depth z, Δ the differential Mueller matrix can be calculated as the matrix logarithm of the Mueller matrix. To speed up the calculations, the matrix logarithm was computed using an approximation [32]: where I is the 4x4 identity matrix. Azzam [33] showed that the birefringence around the Q, U, and V axes is represented by a subset of the off-diagonal elements of the differential Mueller matrix in a first-order approximation: 34 where the first and second number indices indicate row and columns of the matrix m . The major advantage of this algorithm is that its output Γ is a 3D vector whose direction represents the orientation of the sample optic axis in Poincare space while its magnitude

Quantifying optic axis uniformity
The optic axis (OA) orientation of a birefringence sample can be determined with PS-OCT. The round-trip nature of the measurement restricts the measurement to the QU plane of the Poincare sphere. However, due to the birefringence of the fibers, the measured sample optic axes appear on a tilted plane on the Poincare sphere. Rotating this plane back on the QU plane introduces a π ambiguity because two equivalent rotations bring the optic axes plane to the QU-plane, implying that generally only the relative OA orientation is shown [34]. Tissue that is arranged in a fiber-like structure, e.g., collagen and muscle, typically exhibits form birefringence. If the fibers are oriented in the same direction, the OA will show consistent orientation throughout the sample. In this paper, we used a simple metric that evaluates the uniformity of the sample OA [35], previously described in a similar form by Yamanari et al. [36]. The OA orientation was determined from the direction of the real three-dimensional vector Γ extracted as in 2.6. By evaluating the spatial uniformity of the normalized vectors over a small portion of the sample, the optic axis uniformity (OAxU) was extracted as: where , η ν and μ are the components of the vector Γ along the Poincare space axes and γ is the length of the vector , Γ while the angular brackets represent spatial averaging over a small volume. Using this definition restricts the values of OAxU between 0 and 1. Determining the OA orientation is particularly useful for highlighting areas of a sample with uniformly arranged fiber-like structures, as pointed out by Adams et al. [11], so we suggest that OAxU could be an intuitive metric to achieve this goal. Moreover, catheter-based PS-OCT systems require to account for the varying polarization state impinging on the tissue throughout the catheter rotation, which complicates retrieval of the OA orientation [12,37,38]. By choosing an appropriately sized averaging kernel, a constant sampling polarization state can be assumed, circumventing the management of the rotating surface polarization state. In this case, OAxU was calculated over a two-dimensional spatial kernel sized 4x18pixels, corresponding to 26 µm in depth and 7° along the azimuthal direction.

DOCT based on adjacent A-line phase difference
DOCT is a popular method to highlight the presence of blood flow in OCT images. DOCT measures the phase difference between interferograms acquired from the same location at different points in time. Scatterers that have moved in that period create a phase shift of the reflected field from that location. This method requires high phase stability to measure flow in small blood vessels. In the case of endoscopic OCT, the preceding pivoting point and nonuniform rotational distortions (NURD) introduce artificial axial and lateral motion, in turn leading to phase jitter. These effects are strong in proximal scanning catheters, but there are nonetheless examples of DOCT with such devices [13,39]. To this end, distally scanning catheters are more suitable, and examples of angiographic images produced with micromotorbased catheters have been shown [40,41].
To extract DOCT images, we corrected for the phase jitter in the axial direction of the images. In our algorithm, the first step was combining the four measured fields into one by adding the complex fields [42]. In the second step, each A-line z ρ

E E
compensating for bulk motion [43,44]. Spatial averaging of the resulting fields was obtained by convoluting the complex fields with a Gaussian kernel sized twice the optical resolution in both dimensions. In formula: This process guarantees a low noise floor in the phase-difference image. Finally, the square of the phase component z, Δ ρ ϕ of the signal was Gaussian filtered with a kernel size of twice the resolution. A threshold of 0.63 rad was chosen to display the phase difference data as a binary image overlapped with the PS-OCT images.

AC calculation
Images based on the depth-resolved AC of the tissue were produced to improve the interpretation of images. Since the optical energy is not fully dissipated in the tissue and we did not incorporate correction for beam shape and focus position [45][46][47], the AC images were used for qualitative evaluation only. The calculation was based on the algorithm developed by Vermeer et al. [48], which uses chromatic dispersion corrected and roll-off compensated intensity OCT images in linear scale. According to Eq. (18) in Ref [48], the AC µ i of each depth location i is estimated by dividing the intensity i I of each depth pixel I by the sum of the intensities of the pixels underneath: where ∆z is the distance between two consecutive depth pixels and the factor 2 accounts for the double pass of the light into the tissue. Finally, the AC image was produced by taking the base-10 logarithm of the data and displaying them in an inverted greyscale ranging from 10 -1.8 to 10 2.7 mm −1 (45 dB range).

Imaging protocol patient criteria
This severe asthma patient participated in the TASMA trial (trial.gov number NCT02225392). The Ethics committee of the hospital approved the protocol. During bronchoscopy, OCT imaging was performed in two basal segments of the right lower lobe. In short, after advancing the tip of a therapeutic bronchoscope until the segmental airways, the OCT-catheter was inserted through the working channel with an inner diameter of 2.8 mm (Olympus, Tokyo, Japan). After that, the OCT catheter was advanced to the pleura and subsequently retracted by a manual slow-pull technique for image acquisition. The slow pullback was repeated once over the same length of approximately 10 cm to create two image data sets of the same area.

Examples of histology of the human airway wall
Human airway histological slides were obtained from patients that were scheduled for lobectomy because of non-small cell lung cancer (Ethics approved; NL51605.018.14). Exvivo non-diseased airways were dissected, sectioned and analyzed. Human airway histology sections are shown here as a comparison for the (PS-)OCT images presented in the next section. 5 µm-thick histological sections were stained with two different stains: hematoxylin and eosin (H&E) and ASM-specific desmin. Figure 4(a)-(b) show an example of each type of staining. The main histological features of human airways are distinguishable in H&E staining ( Fig. 4(a)), such as epithelium (cyan arrow), basement membrane (grey arrow), lamina propria (pink arrow), and smooth muscle (blue arrow). The red arrow in Fig. 4(a) indicates a piece of cartilage. Desmin staining locates ASM, which appears as brown-colored parts of Fig. 4(b). The main histological features can also be identified in OCT cross-sections acquired in vivo from a human bronchiole, as shown by Fig. 4(c)-(d), which present a traditional OCT intensity image and the AC image of the same frame.

Results
We acquired t of a human lu the potential OAxU image OCT imaging the PS-OCT t are shown in c

Compari
The AC ima different repr 4 the OCT es display local attenuation or scattering density, these artifacts disappear, as indicated by the light blue arrows in Fig. 5(d). The tissue in the superficial areas indicated by the same arrows resemble folded epithelial layers, typical of relaxed bronchioles, as shown in histological slides such as the ones of Fig. 4. It is noteworthy that these structures are almost indistinguishable in the intensity OCT image (Fig. 5(c)). These examples of folded epithelial layers are apparent in Visualization 1 from frame 490 to frame 560 in the top half of the image.

PS-OCT and DOCT in vivo imaging of pulmonary airways
The retrieved PS-OCT show several features that are not visible in the OCT or AC images. Some of these features are particularly striking in the OAxU images, which highlight the presence of fibers oriented consistently in a small portion of tissue. Histological crosssectional images of airway walls reveal the presence of several types of tissue with fiber-like properties, such as ASM and connective tissue layering the cartilage (see Fig. 4). ASM is expected to show high OAxU because of its consistently-oriented fibers. Cartilage is surrounded by a layer of connective tissue, the perichondrium, which suggests that its presence will appear in retardation images as a ring of high birefringence.
In this section, we present several exemplary frames from the OCT volumes, which are available as supplementary material [Visualization 1 and Visualization 2]. Figure 6 and Fig. 7 show two frames of the volume obtained in vivo from location RB8. The blue arrows in Fig. 6(c),(d) indicate a layer of tissue with strong birefringence and high OA uniformity, located about 200 µm from the lumen surface. The layer probably corresponds to ASM. The blue arrows in Fig. 6(a),(b) point at the same area in the intensity and AC images, but delineating the same structures is significantly more difficult. However, the ASM layer seems to show a reduced scattering coefficient and to be surrounded by two thin highly scattering layers, visible in Fig. 6(b). The alveoli at the black arrow location in Fig. 6(a),(b) appear as low-scattering areas. The orange arrows point at large blood vessels, clearly visible from the strong DOCT signal, which appears in red in Fig. 6(c),(d).  Fig. 7(c),(d). Fig. 7(b) Fig. 8(b)). associated ge.  folds are often only visible in the AC images. Structures in deep locations are more prominent in AC images, which show alveoli and cartilage with higher contrast compared to traditional OCT images, suggesting that AC images might be a better method to show tissue structure in OCT images of the airways. Moreover, corrected AC images might offer an additional foundation for discriminating different types of tissue based on their attenuation properties.
The PS-OCT images highlight the presence of several layers of tissue in the airway wall, which are not easily distinguishable in the structural OCT images. A metric that highlights the presence of uniformly oriented birefringent tissue in a small portion of tissue, OAxU, was used for the first time to delineate layers of ASM and connective tissue in the airway lumen, which can be analyzed in more detail in the birefringence images. As already pointed by Adams et al. [11], identifying the OA orientations provides a more robust way of highlighting the presence of ASM, compared to birefringence images. OAxU simplifies this goal by circumventing the need of complex algorithms that compensate for the rotation of the motor (or of the catheter in case of proximal scanning setups) to find the absolute OA orientation, which was the objective of three recent publications [12,37,38]. Since tissue remodeling affects the ASM in the airway wall, OAxU might be a straightforward method to assess the effect of bronchial thermoplasty. Nevertheless, the OAxU images require further investigation to improve the interpretation of this feature-rich modality. Local phase retardation images add information by providing a measurement related to the density of the fiber-like structures in the image.
The DOCT signal reveals large blood vessels close to the airway wall. The NURD of the OCT catheter hampers the detection sensitivity of slow blood flow with DOCT. Good quality angiography of the human airways has been shown exploiting blood autofluorescence [15,49]. If angiography is hypothesized to have diagnostic value for asthma, fluorescence detection could be integrated into the OCT catheter using double-clad fibers as shown in [50].
Precise segmentation of the inner and outer surfaces of the catheter is crucial for performing the PMD and axial bulk motion compensation algorithms presented in this work, therefore improving the quality and reliability of the PS-OCT and DOCT images. Moreover, the segmentation allows removing parts of the image that could be distracting, such as the area within the catheter sheath and the area of tissue obstructed by the wires.
In conclusion, the setup and the algorithms demonstrated in this article further support the suggestion that multifunctional OCT, and in particular PS-OCT, is an adequate tool in terms of specificity, resolution, and field-of-view for the study of tissue remodeling in human airways. PS-OCT imaging might qualify as a high-resolution real-time imaging technique to assess airway wall tissue remodeling in chronic airways diseases, including asthma.

Contributions
FF built the PS-OCT system and designed the OCT catheter. FF, PIB, AWG and JTA acquired the data. DvI designed and built the custom micromotor and assembled the catheter. VD developed and implemented the segmentation algorithm and wrote the AC code. JW, MGOG, and FF adapted the PMD compensation algorithm. JW, FF, MGOG, and JFdB devised the OAxU and Mueller based algorithms. JAMD performed the first in vivo animal and human measurements to assess the catheter safety. JFdB supervised the project and acquired the funding. FF wrote the manuscript with contributions from JW, VD, MGOG, AWMG, PIB, and JFdB. All authors revised the manuscript.