## 2 Subjects and Methods
### Image Reconstruction
#### 2.4.3 Inverse Problem  
We solved the inverse problem by minimizing the objective function given by  
\[\psi(\mu_{a})=\frac{1}{2}\sum_{i=1}^{N_{S}}\sum_{j=1}^{N_{D}}H_{i,j}\big{[}y_{ ij}-f_{i,j}(\mu_{a},\kappa)\big{]}^{2}+\tau\sum_{k=1}^{N_{R}}(\mu_{a,k}-\mu_{a0,k}) ^{2}, \tag{14}\]  
which is the least square method for the measured result \(y_{ij}\) and the calculated result \(f_{i,j}\); \(N_{S}\) and \(N_{D}\) are the numbers of the sources and detection fibers used in measurements, respectively; \(H_{i,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), \(H_{ij}\) = 0. We used Tikhonov regularization, and \(\tau\) is the hyper parameter; \(N_{R}\) 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 is1.6 (L) \(\times\) 1.9 (W) \(\times\) 2.5 mm (H). The value of \(N_{R}\) is 8,014 in the trachea-tissue boundary grid (Figure 3b), and \(N_{R}\) is 8,992 in the homogeneous grid (Figure 3c); \(\mu_{a0}\) is the initial absorption coefficient at the start of the reconstruction.  
The measured TOF distribution, \(Y_{ij}(t)\), contains unknown constants from the measurement setup. To remove these constants, we used the normalized first-order moment, \(y_{i,j}\), as the measured result. The form is  
\[y_{i,j}=\frac{\int_{0}^{T}tY_{i,j}(t)dt}{\int_{0}^{T}Y_{i,j}(t)dt}, \tag{15}\]  
and the calculated result \(f_{ij}\) corresponding to \(y_{i,j}\) can be computed as  
\[f_{i,j}(\mu_{a},\kappa)=\frac{\int_{0}^{T}tD_{ij,j}(t)dt}{\int_{0}^{T}D_{ij,j}( t)dt}+\frac{\int_{0}^{T}t\Gamma_{ij}(\mu_{a},\kappa,t)dt}{\int_{0}^{T}\Gamma_{i,j}( \mu_{a},\kappa,t)dt}, \tag{16}\]  
where \(D_{ij}\) (\(t\)) is IRF.  
To minimize Equation (14), we used the nonlinear conjugate gradient method with the Polak-Ribiere method [26; 27], which is implemented in TOAST++ [28]. The gradient of the objective function was obtained with the adjoint method; \(\mu_{a}\) = 0.027/mm and \(\mu_{s^{\prime}}\) = 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 \(\mu_{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 \(\mu_{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 StO\({}_{2}\) (oxy-Hb/(oxy-Hb+deoxy-Hb)) was calculated.