Imaging the Human Thyroid Using Three-Dimensional Diffuse Optical Tomography: A Preliminary Study

Featured Application: Diffuse optical tomography of the thyroid provides great promise for more accurate diagnosis of thyroid cancers. Abstract: Thyroid cancer is usually diagnosed by ultrasound imaging and fine-needle aspiration biopsy. However, diagnosis of follicular thyroid carcinomas (FTC) is difficult because FTC lacks nuclear atypia and a consensus on histological interpretation. Diffuse optical tomography (DOT) offers the potential to diagnose FTC because it can measure tumor hypoxia, while image reconstruction of the thyroid is still challenging mainly due to the complex anatomical features of the neck. In this study, we attempted to solve this issue by creating a finite element model of the human neck excluding the trachea (a void region). By reconstruction of the absorption coefficients at three wavelengths, 3D tissue oxygen saturation maps of the human thyroid are obtained for the first time by DOT.


Introduction
The incidence of thyroid cancer has been increasing in many parts of the world over the past decades, although thyroid cancer remains among the less common cancers [1,2]. The etiology of the increasing incidence is not established. An analysis of data from the Surveillance Epidemiology and End Results (SEER) database has indicated that the increase is due to improved detection and not to a rise in occurrence [3]. However, the worldwide increase cannot be explained simply by noting earlier discovery [1], and there is increasing interest in the diagnosis and treatment of thyroid cancers.
There are four main types of thyroid cancer. Papillary thyroid carcinoma (PTC) and follicular thyroid carcinoma (FTC) are categorized among the differentiated thyroid cancers that are not usually aggressive. PTC is by far the most common form and is usually slow-growing, but it can mutate into more aggressive variants. Thyroid cancer is often diagnosed after nodules in the thyroid are discovered as a lump in the neck or by routine image examination, such as an ultrasound scan performed on the carotid arteries [4].
Other than FTC, which is more aggressive than PTC, thyroid cancers can be diagnosed by ultrasound scanning and fine-needle aspiration biopsies. Since FTC lacks the nuclear atypia seen in PTC and is encapsulated like a benign tumor (i.e., follicular adenoma), diagnosis of FTC is based on the histologic demonstration of capsular and/or vascular invasion after total thyroidectomy [5]. However, it is difficult to diagnose capsular invasion because of a lack of consensus on the histologic interpretation of findings: some argue that penetration into the tumor capsule by lesional cells is indicative of malignancy, while others argue that infiltration of tumor cells into the surrounding thyroid parenchyma through the entire thickness of the capsule is capsular invasion [5,6]. Thus, noninvasive and more sensitive diagnostic approaches to FTC would be helpful.
X-ray CT (computed tomography) and MRI (magnetic resonance imaging) are common diagnostic imaging techniques. Since, however, the ability to characterize thyroid nodules and detect small carcinomas is inferior to ultrasound scanning, X-ray CT and MRI are mainly used to evaluate extrathyroidal carcinoma extensions [7]. With 18FDG PET/CT (fluorine-18 fluorodeoxyglucose positron emission tomography/computed tomography) it is possible to detect solid tumors but not thyroid cancers, as there is uptake of 18FDG in chronic thyroiditis [8]. Thyroid cancer is accompanied by functional changes, including angiogenesis [9] and its related molecular responses [10,11], as well as structural changes. Optical imaging with near-infrared light, including photoacoustic imaging (PAI), diffuse optical spectroscopy (DOS), and diffuse optical tomography (DOT), can noninvasively evaluate hemodynamic and metabolic changes associated with thyroid cancer, and offers the potential for a diagnostic approach to FTC.
To date, PAI has been applied to the detection of angiogenesis in thyroid cancer. Yang et al. demonstrated that PAI detected small blood vessels that color Doppler flow imaging (CDFI) could not, whereas it has also been reported that acoustic reflection at large impedance-mismatched boundaries at the trachea caused artifacts in PAI images [12]. This approach could be useful in diagnosing PTC but not for FTC, because FTC is not accompanied by angiogenesis [9]. It has also been suggested that two members of the matrix metalloproteinase family (MMP-2, MMP-9) are biomarkers for FTC [10]. Using an MMP-activatable photoacoustic probe, Levi et al. demonstrated that molecular PAI offers great promise for the diagnosis of FTC, but clinical application appears difficult [13]. Lindner et al. employed custom time-resolved spectroscopy and a diffuse correlation spectroscopic system to the measurements of total hemoglobin (Hb) concentrations, tissue oxygen saturation (StO2), and blood flow indexes in the thyroid and also in muscles [14]. Several studies have suggested that hypoxia contributes to the progression of thyroid cancer, especially via the hypoxia-inducible factor-1α (HIF1-α) [11,15], and a lower StO2 has the potential to be an indicator of malignancy. However, with the DOS system of Lindner et al., it is difficult to measure the thyroid separately from muscles, while it would be possible with DOT.
Diffuse optical tomography is an imaging modality that reconstructs absorption (μa) and reduced scattering (μs') coefficients from boundary measurements [16]. From the reconstructed μa, the Hb concentration and StO2 can be calculated. The μs' reflects morphological and structural changes, such as cell swelling caused by cell injuries, whereas it is unlikely that μs' is altered by capsular and/or vascular invasion by FTC. Cytochrome c oxidase (cyt. ox.) in mitochondria is also a biological chromophore that can be measured by DOT. The redox state of cyt. ox. is dependent on the intracellular oxygen concentration, and its reduction occurs under severe hypoxic conditions [17]. From this, it may be expected that DOT will be able detect FTC by measuring the StO2 and cyt. ox. in the thyroid.
To the best of our knowledge, DOT images of the human thyroid have not been reported, except for those in our simulation study [18]. This may be attributed to the fact that the human neck is a heterogeneous medium; it is composed of the trachea, the thyroid, the vertebrae, large blood vessels, muscles, adipose tissue, and others. The thyroid is located in a narrow region between the skin and trachea, and DOT image interpretation of the thyroid is not straightforward. In this study, we attempted to solve this issue by extracting the minimum volume of the human neck that is required for the image reconstruction, with the void region of the trachea excluded, and creating a finite element grid. In a clinical study, we conducted DOT measurements of a healthy human thyroid and reconstructed the spatial distribution of μa to create a 3-dimensional (3D) image of the whole of the thyroid. Then, using the μa values at three wavelengths, StO2 maps were created.

Instrumentation
A multichannel time-domain measurement instrument (TRS-80, Hamamatsu Photonics K.K., Hamamatsu, Japan) was used to measure the human neck ( Figure 1a). This instrument was developed from TRS-20SH, which has a single channel. See [19] for details of the TRS-20SH. Briefly, the TRS-80 consists of three pulsed laser diodes at wavelengths of 763, 801, and 836 nm (the pulse width is less than 100 ps for each wavelength) with a repetition rate of 5 MHz. All laser diodes are pigtailed and coupled into one optical fiber by the fiber coupler. The fiber is coupled to eight source fibers (core diameter 200 μm, NA (numerical aperture) = 0.25) via the optical switch (switching time, about 1 s), and the tip of each source fiber is attached to a measurement target (average power 200 μW). Transmitted light through the target is detected by eight fiber bundles (diameter 3 mm, NA = 0.29) and transferred to eight high-speed photomultiplier tubes (PMTs) (H7422P-50 SEL, Hamamatsu Photonics K.K., Japan) for single photon detection. The electrical output of a PMT is transmitted to the time correlated single photon counting (TCSPC) unit. It takes 1 s to obtain the time-of-flight (TOF) distribution of photons (time step size, 10 ps) in the time range from 0 to 10.24 ns for each wavelength. The TOF distributions with adequate signal-to-noise ratios (SNRs) are acquired by summation of signals for 10 s. The criterion of adequate SNR is that the maximum count of TOF distribution is greater than 2,000 and the dark count is less than 10. The average of the full width half maximum of the instrumental response function (IRF) is 390 ps. The TRS-80 is portable and available for use at the bedside. In this study, the measurement of the neck was conducted by a TRS-80 at the bedside while the subject lay on the bed.

Human Neck Measurements
The subject is a healthy woman who was informed of the aim and procedures of this study. This study was approved by the ethical committee of Hamamatsu University School of Medicine (No. 15-082). To confirm the tracheal and thyroid positions, a cervical MRI scan was performed (Discovery MR750, GE Healthcare, Waukesha, WI, USA). The MR images were also used to make an optical fiber holder in our laboratory. Since the same neck position must be reproduced in the DOT measurements, a vacuum cushion for radiotherapy was used for the MRI scan to fix the neck position.
The holder (Figure 1b,c) was made based on the 3D MR neck images by using 3Dcomputer-aided design (CAD) software (Fusion360, Autodesk, San Rafael, CA, USA) and a 3D printer (Agilista, Keyence, Osaka, Japan). Since the thyroid is a superficial organ and the minimum distance between the skin and the thyroid was about 1 cm in this subject, 1.5 cm was used as the shortest source-detector distance. A total of 15 holes in the holder were arranged at an interval of 1.5 cm and 7 source and 8 detector fibers were inserted into the holes alternately so that the whole thyroid was illuminated ( Figure 1d). As a result, 56 datasets (7 × 8) were obtained. Each fiber was placed perpendicular to the skin surface. Optical fibers were fixed firmly to the holder by the screw cap on each hole of the holder. Moreover, the empty holes were filled with black polyoxymethylene (POM) cylinders. This fiber arrangement enables us to detect light passing through the thyroid and reduce influences of any strong absorption occurring in the internal jugular vein and the common carotid artery, which are far from the optical fibers. The black natural rubber (KSNR-10054T, HIKARI CO., LTD., Osaka, Japan) inside the holder (Figure 1c) is used to absorb light at the interface between air and the skin so that the Robin boundary condition is applied to numerical modeling of light propagation. In addition, it plays the role of a cushion between the rigid holder and the neck. The holder with the optical fibers was attached to the neck of the subject using the triumphal-arch shaped support. The subject was lying down on the vacuum cushion to be able to reproduce the neck position for the MRI scan.
Prior to the neck measurements, we measured the IRFs (Figure 2b-d) in all combinations of the source and detection fibers by using an 8ch IRF module, which enabled us to obtain eight IRFs for one source fiber simultaneously. The 8ch IRF module (7.6 (D) × 10 (W) × 10 (H) cm) has one source fiber port on one side and eight detection fiber ports on the other side ( Figure 2a). A diffuser is placed inside each port. The detection fiber ports are arranged at regular intervals on the circumference of a circle (a diameter of 30 mm) of which center is opposite to the source fiber port. The space between the source fiber and the detection fibers is filled with air, and the distance between every sourcedetection pair is 7.5 cm. Light was emitted from each source fiber about every 11 s, and light transmitted through the neck was detected at all the eight detection fibers simultaneously. The TRS-80 automatically adjusted the optical attenuators placed in front of each PMT. A total measurement time was about 120 s, including the switching time of the optical switch and the time for adjusting the attenuators. Forty-nine TOF distributions with an average photon counting rate of more than 115,000 count/s and a peak count over 3,000 counts were used for image reconstruction. Figure 2b-d show the TOF distribution at S2D2. The remaining seven TOF distributions measured at seven source-detector pairs (S2-D1, S2-D6, S3-D3, S3-D8, S5-D1, S5-D6, and S7-D6) were excluded from the image reconstruction because of insufficient SNR.

Grid Excluding the Trachea
The DOT images suffer from a void region if no care is taken for the grid used for the finite element method [20,21]. To create a grid excluding the void region of the trachea, we made use of MR neck images.
The MR neck image was processed with the open-source software platform for medical image processing (3D-Slicer, version 4.10.2) [22] to generate the 3D surface geometry of the neck. First, the T2-weighted MR image was saved as the DICOM standard, which contains multiple slices and all pixels of each of the slices related to the real position and scale. Then, 3D-Slicer loaded the MRI DICOM file and contours of the neck, shoulder, trachea, and thyroid were delineated by the Segment Editor module. Finally, the Segmentations module exported the 3D surface geometry of the various objects ( Figure 3a) in a standard triangulated language (STL) file format at LPS (Left, Posterior, Superior) coordinates. In Figure 3a, the grid region (see below) is indicated in yellow. All STL files contain real-scale anatomical structures of the neck using the coordinates of the MR image. Next, 3D-CAD imported the STL files and extracted the minimum volume that is required for the image reconstruction. Since the number of photons quickly decays when near-infrared light propagates in the neck, there is no need to make the grid using the entire volume. Here, light propagation around the thyroid was numerically studied, and this rapid decay was observed [23]. The upper and lower boundaries of the volume were set as horizontal sections 10 mm above the top row and 10 mm below the bottom row of the optical fibers, respectively. The sides of the volume were sagittal sections 10 mm away from the left and right columns of the optical fibers. The posterior limit was a plane that divided the void region approximately into halves. This grid is termed the trachea-tissue boundary grid. The processed surface geometry is shown in Figure 3b. As a comparison, we made another grid by disregarding the presence of the trachea. This grid is termed the homogeneous grid (Figure 3c).
As a last step, the processed STL files were imported into the Gmsh version 4.6.0 [24] to create 3D volume grids for use with the finite element method. The Gmsh generates a very nearly uniform volume of tetrahedral grids where nodes of the surface are not necessarily identical to the nodes from the STL file. The side length of a tetrahedron was about 1 mm. In this way, uniform grids of tetrahedrons were used for the DOT. The trachea-tissue boundary grid (Figure 3b) has 53,611 nodes and 264,526 elements, and the homogeneous grid (Figure 3c) has 57,433 nodes and 291,968 elements.

Light Transport Model
The light transport model considered in this study is the time-dependent diffusion equation (TD-DE), which is the diffusion approximation of the radiative transport equation. The TD-DE form is with the Robin boundary condition which describes the change in the fluence rate φ at time t and position vector r within a domain Ω, bounded by a surface ∂Ω. Here, T is the maximum time, c = c0/n is the speed of light in the medium with refractive index n, and c0 is the speed of light in a vacuum. The refractive index of the neck is assumed to be 1.37; μa and μs' are the absorption and reduced scattering coefficients, respectively; κ (r) = {3[μa (r)+ μs' (r)]} −1 is the diffusion coefficient; ζ is for the refractive index mismatch parameter at the surface ∂Ω; qi (r, t) = δ (r−ri, t) is the ith isotropic source term at source position ri ∈ ∂Ω (i = 1,…, NS); the measured quantity is the flux Γi, j (t) at the detection position rj ∈ ∂Ω (j = 1,…, ND). The relation between φ and Γ is given by Fick's law, where ν is the outer normal vector of ∂Ω at rj.
To solve the inverse problem in Section 2.4.3, we used the normalized first-order moment of the flux To remove unknown constants in the measured TOF distribution (T = 10.24 ns). It takes a great deal of computational time to solve Equation (1), and, as a result, we obtained gi,j by the solution of the time-independent diffusion equation [25].

Finite Element Method
To find the optimum solution of the time-independent diffusion equation, we used the finite element method (FEM). With φ h (r, t) as the solution to Equation (1), where uk (r) (k = 1, …, NN) is the basis function and NN is the number of nodes in the grid. We introduced a vector Φi (t) = (…, Φk, i (t), …) T of nodal coefficients, and to discretize the zeroth-order and the first-order moments of TD-DE (Equation (1)), we used a Galerkin approach [25]. The forms are The individual integrals are given by The flux is expressed as Here, χk (rj) = 1 when the element τk belongs to rj ∈ ∂Ω and χk = 0 otherwise. Then, Equation where Pj = (…, Pk, j, …) T .

Inverse Problem
We solved the inverse problem by minimizing the objective function given by ∑∑ ∑ (14) which is the least square method for the measured result yi,j and the calculated result fi,j; NS and ND are the numbers of the sources and detection fibers used in measurements, respectively; Hi,j = 1 except for seven source-detection combinations that were excluded due to insufficient SNR. For the excluded pairs (i, j) = (2, 1), (2,6), (3,3), (3,8), (5,1), (5,6), and (7, 6), Hij = 0. We used Tikhonov regularization, and τ is the hyper parameter; NR is the number of points of the reconstructed grid used to solve the inverse problem. The grid consists of rectangular voxels converted from the FEM grid, and one voxel size is 1.6 (L) × 1.9 (W) × 2.5 mm (H). The value of NR is 8,014 in the trachea-tissue boundary grid ( Figure  3b), and NR is 8,992 in the homogeneous grid ( Figure 3c); μa0 is the initial absorption coefficient at the start of the reconstruction. The measured TOF distribution, Yi,j(t), contains unknown constants from the measurement setup. To remove these constants, we used the normalized first-order moment, yi,j, as the measured result. The form is where Di,j (t) is IRF.
To minimize Equation (14), we used the nonlinear conjugate gradient method with the Polak-Ribière method [26,27], which is implemented in TOAST++ [28]. The gradient of the objective function was obtained with the adjoint method; μa = 0.027/mm and μs' = 0.95/mm were employed as the initial values for all the three wavelengths (763, 801, and 836 nm), which were taken from [14]. During the iteration process, the μs' value was fixed.
The thyroid is a richly vascular organ and consists of a large number of thyroid follicles of which walls are made up of a large number of cells, mainly follicular cells. Follicles are filled with colloid (thyroglobulin). Thyroglobulin contributes to the absorption coefficients in the three wavelengths used in this study [29]. Since, however, it is difficult to estimate its contributions quantitatively, treating thyroglobulin as water, it is assumed that Hb and water are the main chromophores of the thyroid. Wavelengths of 763, 801, and 836 nm are commonly used to measure tissue Hb, and taking account of the water contribution to μa at each wavelength, Hb concentrations are calculated. The molar absorption coefficients of oxygenated Hb (oxy-Hb) and deoxygenated Hb (deoxy-Hb) were taken from [30] and the absorption coefficient of water was taken from [31]. Based on the assumption that the water concentration was 78% [32], the value of StO2 (oxy-Hb/(oxy-Hb+deoxy-Hb)) was calculated. Figure 4 shows horizontal sections of DOT images of μa at 801 nm, which is an isosbestic point for oxy-Hb and deoxy-Hb, at the 15th iteration by the trachea-tissue boundary grid (Figure 3b) at different slices. The DOT images are shown on the T2weighted MR images. Figure 4b shows the section at the middle of the isthmus of the thyroid. Figure 4a,c are sections 3 mm above and below the middle of the isthmus of the thyroid. The trachea-tissue boundary grid region is shown in Figure 4c (dashed white line). The color bar shows the reconstructed μa values. In Figure 4a-c, the outline of the thyroid is denoted by white lines. The shape of the thyroid, where the μa values are about 0.075/mm, is reconstructed in the DOT images, while the position of the reconstructed thyroid shifts slightly towards the skin. The μa values of the background tissue (muscles) were about 0.02/mm. Although the right internal jugular vein (IJV) was included in the grid, it was not reconstructed in the image because the optical fibers were arranged far from the IJV. In fact, the μa values in the IJV were not updated and the initial ones remained (0.027/mm) during the iteration process.   Figure 4. The StO2 in the background tissues, mainly muscles, was around 65%. Although the μa values at the isthmus of the thyroid were much smaller than those in the bilateral lobes, the StO2 distribution was relatively homogeneous, and StO2 of the thyroid was 82% to 92%. As a result, the isthmus of the thyroid was more clearly imaged compared to that in Figure 4.  Figure 6 shows horizontal sections of DOT images of μa at 801 nm at the 35th iteration by the homogeneous grid (Figure 3c). Although the existence of the thyroid may be read from these reconstructed images, the shape of the thyroid is not correctly reconstructed.

Discussion
This is the first study reporting the reconstruction of the spatial distribution of μa and StO2 in the whole human thyroid with DOT. The neck is composed of various organs, and the refractive index inside the neck is heterogeneous. Notably, between the tissue and the trachea, there is a strong refractive index mismatch. To date, whether the thyroid can be imaged accurately by DOT has been regarded as a seriously nontrivial question, but in this paper, we have shown that the whole shape of the thyroid can be reconstructed if the grid for the calculations excludes the trachea void region. However, the reentrance of light, which propagates through the trachea, was not taken into account in the present study. It is expected that reconstructed images will be improved if the reentrant light is taken into account [23].
In this paper, the thyroid was reconstructed at a position closer to the skin than its actual position. This can be understood from the structure of the objective function used in our numerical reconstruction. Since the diffuse fluence of near-infrared light decays exponentially with increased depth in biological tissue, measurement sensitivity decreases along the depth direction. Spatially variant regularization has been proposed to compensate for this [33,34]. However, we did not employ depth-dependent regularization, which is one reason why the position of the reconstructed thyroid shifted slightly towards the skin.
The obtained μa values in the present study were much larger than the published values that have been obtained with time-domain measurements. In a study preceding this, the μa values were estimated by fitting the TOF distribution derived from the analytical solution of the DE for a semi-infinite homogeneous medium to the measured TOF distribution [14]. However, it cannot be assumed that the human neck is homogeneous, and as the thyroid is covered by muscles, it is likely that the published μa values should be partly attributed to the μa of muscle tissue. Here, the present study indicates the possibility of selectively measuring the μa of the thyroid, and the μa values in the present paper may also be assumed to deviate from the true situation to some extent, because the reduced scattering coefficient was not updated during the iterations of the nonlinear conjugate gradient method. An accurate quantification of μa requires the simultaneous reconstruction of both μa and μs' [35], which are subjected to artifacts caused by cross-talk between the two parameters. To suppress cross-talk artifacts, the use of multiple datasets for the objective functions, such as a combination of skew and Laplace transform of the TOF distribution (Equation (3)), has been proposed [36].
The present DOT is still insufficient for the quantitative measurement of optical properties and Hb concentrations; however, the map of StO2 ( Figure 5) indicates that the thyroid is distinguishable from the background tissue (muscles): the StO2 of the thyroid was around 90%, while that of the muscles was around 65%. A high StO2 is the more favorable option for detecting cancer hypoxia, because the difference in StO2 between normal and cancerous tissue is thought to be larger than that when the StO2 of the thyroid is lower. Cytochrome c oxidase is a more sensitive indicator of cell hypoxia than StO2 [17]. However, TRS-80 with three wavelengths cannot measure the redox state of cyt. ox., and an instrument with four or more wavelengths is required.
Regarding treatment of thyroid diseases, it is important that the shape and size of the thyroid are tracked over time, which can be conducted with ultrasound scanning. Hence, several 2D segmentation algorithms for thyroid ultrasound images have been developed, and 3D segmentation is now possible [37,38]. Three-dimensional visualizations of the thyroid are also crucial for surgical intervention [39]. Recent advances in segmentation and image processing techniques have enabled 3D thyroid reconstruction from 2D ultrasound slices [40,41].
The present study has the three main limitations described above, but it has still paved the way for functional optical imaging of the human thyroid. The combination of DOT with ultrasound scanning will be complementary in the course of diagnosis and treatment of thyroid cancers.

Conclusions
This paper presents 3D absorption and StO2 images of an in vivo human thyroid by DOT. The key factor is the use of a trachea-tissue boundary grid. This study has shown that DOT can be used to selectively reconstruct the thyroid, and a natural next step will be to estimate the value of the absorption coefficient by the simultaneous reconstruction of μa and μs'. This would make DOT available for use in thyroid cancer diagnosis. It can be expected that DOT will have the potential to visualize the tissue hypoxia of thyroid cancer.  Informed Consent Statement: Informed consent was obtained from all subjects involved in the study.