*in vivo*experiments of time-domain DOT. As a result, the artifacts resulting from changes in optode coupling are reduced in the reconstructed images of the absorption coefficient, thereby demonstrating that introduction of coupling coefficients is effective in time-domain DOT. Moreover, numerical simulations, phantom experiments, and

*in vivo*studies have validated this algorithm.

## 1.

## Introduction

Recent developments in optical technologies have made it possible to provide tomographic images of oxygenated states of the human brain by applying near-infrared (NIR) spectroscopy, known as diffuse optical tomography (DOT).^{1, 2} In DOT, image reconstruction is based on the modeling of light propagation in an object, known as the photon diffusion equation (PDE). PDE describes the relationship between the optical properties in the object and the light intensities measured at the object's surface. The measured light intensities are used to solve the inverse problem based on the PDE, and the solution is the reconstructed image showing the distribution of the optical properties within the object, such as the absorption coefficient.^{3}

Because the absorption spectrum of oxy-hemoglobin is different from that of deoxy-hemoglobin, the reconstructed images of the absorption coefficient at multiple wavelengths can be converted into images of the concentrations of oxy- and deoxy-hemoglobins. Consequently, DOT images provide hemodynamic information. A promising application of DOT may be bedside monitoring of the oxygenated states of a premature infant's brain, which should be carefully monitored to prevent hypoxic-ischemic brain damage.^{4, 5, 6, 7}

However, DOT images often suffer from artifacts, which are unreasonable from a physical or physiological point of view. The artifacts are caused by a number of issues, one of which is due to the change in coupling between the optode (composed of optical fibers) and the object's surface during the measuring period. The change in coupling may easily be caused by undesirable movements of the subject; for example, movements made by premature infants. Previous studies have reported that artifacts can be removed by optode calibration, which introduces “coupling coefficients” to represent the state of optode coupling.^{8, 9, 10, 11, 12, 13}

Boas first proposed optode calibration by incorporating coupling coefficients into the algorithm for cw (continuous-wave)-domain DOT.^{8} By taking the logarithm of measured fluence rates as input data, coupling coefficients were linearly related to the measured cw data, and by employing the Rytov approximation, the reconstruction process was made into a linear inverse problem.

Schweiger also incorporated coupling coefficients into the frequency-domain DOT algorithm, where the amplitude and phase of the detected light were calibrated.^{9} They solved a nonlinear inverse problem involving the coupling coefficients and demonstrated the reduction of artifacts by showing DOT images for two-dimensional (2D)-simulation and phantom experiments.

Schmitz developed a CCD-camera and silicon-photodiode–based system and adopted a method for optode calibration with coupling coefficients.^{10} Oh formulated the DOT reconstruction algorithm with optode calibration under the Bayesian framework.^{11} Stott proposed a method to calibrate the amplitudes of the measured intensities and optode positions.^{12} Tarvainen succeeded in performing image reconstruction with optode calibration under a rotational symmetry of optode positions in frequency-domain DOT. Phantom experiments also validated their method.^{13}

Various methods of optode calibration were studied for cw- and frequency-domain DOT; however, no study has been performed on time-domain DOT. In this study, we introduce optode calibration to time-domain DOT. Changes in the amplitudes of the time-resolved measurement data are assumed to take place because of the differences found in optode coupling between the measurements obtained for rest and task states. The nonlinear algorithm for DOT image reconstruction employed in this study is based on the modified generalized pulse spectrum technique (mGPST),^{14, 15} which is one of the image reconstruction algorithms for time-domain DOT. Numerical simulations, phantom experiments, and *in vivo* experiments have validated the performance of the algorithm which incorporates coupling coefficients. As a result, the artifacts induced by changes in optode coupling are reduced in the reconstructed images of the absorption coefficient, showing that introduction of the coupling coefficients is effective in time-domain DOT and that this algorithm has been validated by numerical simulations, phantom experiments, and *in vivo* experiments.

## 2.

## Method of Image Reconstruction and Optode Calibration

## 2.1.

### Forward Model of Light Propagation in Biological issues

The radiative transfer equation (RTE) rigorously describes light propagation in biological media. This integro-partial differential equation can be solved numerically by deterministic or stochastic methods. However, solving RTE numerically may be computationally too intensive to be applied to image reconstruction. Therefore, the PDE, which is an approximation of RTE and is valid for biological objects larger than about 1 cm, is usually employed.^{1, 2, 3}

## Eq. 1

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \frac{1}{c}\frac{{\partial \phi ({\bf r},t)}}{{\partial t}} = \nabla \cdot \{ {D({\bf r})\nabla \phi ({\bf r},t)}\} - \mu _a ({\bf r})\phi ({\bf r},t) + q({\bf r},t), \end{equation}\end{document} $$\frac{1}{c}\frac{\partial \phi (\mathbf{r},t)}{\partial t}=\nabla \xb7\left\{D\left(\mathbf{r}\right)\nabla \phi (\mathbf{r},t)\right\}-{\mu}_{a}\left(\mathbf{r}\right)\phi (\mathbf{r},t)+q(\mathbf{r},t),$$*c*is the speed of light in the medium,

*ϕ*(

**,**

*r**t*) is the fluence rate at position

**and time**

*r**t*,

*D*(

**) = 1/[3**

*r**μ*

_{s}

^{′}(

**)] is the diffusion coefficient,**

*r*^{16}

*μ*

_{s}

^{′}(

**) is the reduced scattering coefficient,**

*r**μ*

_{a}(

**) is the absorption coefficient, and**

*r**q*(

**,**

*r**t*) is the light source in the medium.

When the impulse light source illuminates the medium from the surface at position *r*_{s}, Eq. 1 is subject to the following initial condition [Eq. 2],

## Eq. 2

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \phi ({\bm r},t) = \delta ({\bm r} - {\bm r}_s)\delta (t),\end{equation}\end{document} $$\phi (\mathit{r},t)=\delta (\mathit{r}-{\mathit{r}}_{s})\delta \left(t\right),$$*δ*is the delta function. Also, Eq. 1 is subject to the following boundary condition [Eq. 3], which expresses that no light flux is incoming from the boundary,

*r*_{b}, except the illuminating position,

## Eq. 3

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} D({\bm r}_b)\nabla _n \phi ({\bm r}_b,t) = - \frac{1}{{2A({\bm r}_b)}}\phi ({\bm r}_b,t),\end{equation}\end{document} $$D\left({\mathit{r}}_{b}\right){\nabla}_{n}\phi ({\mathit{r}}_{b},t)=-\frac{1}{2A\left({\mathit{r}}_{b}\right)}\phi ({\mathit{r}}_{b},t),$$_{n}represents the gradient normal to the surface and

*A*(

*r*_{b}) = [1+

*r*

_{d}(

*r*_{b})]/[1-

*r*

_{d}(

*r*_{b})] is the parameter depending on the internal reflection ratio,

*r*

_{d}, given as

*r*

_{d}= −1.440

*n*

_{r}

^{−2}+ 0.710

*n*

_{r}

^{−1}+ 0.668 +0.0636

*n*

_{r}, with

*n*

_{r}as the refractive index ratio between the medium and its surroundings.

^{17}

Light intensities measured at the boundary, **Γ**(*r*_{b}, *t*), are the fluxes of the fluence rate outgoing from the boundary (Fick's law) and are given by the following:

## Eq. 4

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \Gamma ({\bm r}_b,t) = - D({\bm r}_b)\nabla _n \phi ({\bm r}_b,t).\end{equation}\end{document} $$\Gamma ({\mathit{r}}_{b},t)=-D\left({\mathit{r}}_{b}\right){\nabla}_{n}\phi ({\mathit{r}}_{b},t).$$The forward problem of Eq. 1 is solved numerically by the finite element method (FEM)^{18} under the initial and boundary conditions [Eqs. 2, 3], and the measured light intensities are obtained by Eq. 4.

## 2.2.

### Inverse Problem with Optode Calibration

## 2.2.1.

#### Reconstruction using mGPST

The mGPST is employed to reconstruct the distribution of the optical properties in the medium based on the forward modeling previously described above. The mGPST uses the Laplace transform to extract features of the measured time-resolved data. The transformed “featured data,” which effectively summarizes the measured time-resolved data, is used for successive image reconstruction.

Fundamentally, the reconstruction algorithm minimizes the residual errors defined as the norm of the difference between the vectors of the measured and calculated featured data,
[TeX:]
$\| {{\bm{\tilde \Gamma} }_m - {{\bm{\tilde \Gamma}} }}\|$
$\Vert {\stackrel{\u0303}{\mathbf{\Gamma}}}_{m}-\stackrel{\u0303}{\mathbf{\Gamma}}\Vert $
, where
[TeX:]
${{\bm{\tilde \Gamma}} }_m = \int_0^\infty {{\bm\Gamma }_m } \exp (- pt)dt$
${\stackrel{\u0303}{\mathbf{\Gamma}}}_{m}={\int}_{0}^{\infty}{\mathit{\Gamma}}_{m}\mathrm{exp}(-pt)dt$
and
[TeX:]
${{\bm{\tilde \Gamma}} } = \int_0^\infty {{\bm\Gamma }} \exp (- pt)dt$
$\stackrel{\u0303}{\mathbf{\Gamma}}={\int}_{0}^{\infty}\mathit{\Gamma}\mathrm{exp}(-pt)dt$
with the constant *p* for the Laplace transform, and **Γ**
_{m} and **Γ** are the measured and calculated time-resolved data, respectively. The vectors
[TeX:]
${{\bm{\tilde \Gamma}} }_m$
${\stackrel{\u0303}{\mathbf{\Gamma}}}_{m}$
and
[TeX:]
${{\bm{\tilde \Gamma}} }$
$\stackrel{\u0303}{\mathbf{\Gamma}}$
have *M* components, corresponding to *M* measurements for the combinations of source and detector positions, and
[TeX:]
$\bm\Gamma$
$\mathit{\Gamma}$
is numerically calculated by the forward model with the given optical properties.

We can expand the calculated featured data, [TeX:] ${{\bm{\tilde \Gamma}} }$ $\stackrel{\u0303}{\mathbf{\Gamma}}$ , in a Taylor series as follows:

## Eq. 5

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} {{\bm{\tilde \Gamma}} }({\bm\mu }_t + {\bm{\delta\mu} }_t) = {{\bm{\tilde \Gamma}} }({\bm\mu }_t) + {\bf J}({\bm\mu }_t)\delta {\bm\mu }_t + \cdots, \end{equation}\end{document} $$\stackrel{\u0303}{\mathbf{\Gamma}}({\mathit{\mu}}_{t}+{\mathit{\delta}\mathbf{\mu}}_{t})=\stackrel{\u0303}{\mathbf{\Gamma}}\left({\mathit{\mu}}_{t}\right)+\mathbf{J}\left({\mathit{\mu}}_{t}\right)\delta {\mathit{\mu}}_{t}+\cdots ,$$*N*× 1 vector representing the optical properties at the discrete

*N*positions in the medium and is to be reconstructed, subscript

*t*indicates the number of iterations,

**J**is the Jacobian matrix consisting of first-order partial derivatives of the featured data with respect to [TeX:] $\bm\mu$ $\mathit{\mu}$ , and δ [TeX:] $\bm\mu$ $\mathit{\mu}$ of the

*N*×1 vector is a small change in [TeX:] $\bm\mu$ $\mathit{\mu}$ . The terms higher than the first-order derivative in Eq. 5 are neglected. Now, by assuming that [TeX:] ${{\bm{\tilde \Gamma}} }$ $\stackrel{\u0303}{\mathbf{\Gamma}}$ ( [TeX:] $\bm\mu$ $\mathit{\mu}$

_{t}+ δ [TeX:] $\bm\mu$ $\mathit{\mu}$

_{t}) is equal to [TeX:] ${{\bm{\tilde \Gamma}} }_m$ ${\stackrel{\u0303}{\mathbf{\Gamma}}}_{m}$ , the updating process for true optical properties is expressed as follows. First, the following linear system equation for δ [TeX:] $\bm\mu$ $\mathit{\mu}$

_{t}is solved by the algebraic reconstruction technique,

## Eq. 6

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} {{\bm{\tilde \Gamma}} }_m - {{\bm{\tilde \Gamma}} }({\bm \mu }_t) = {\bf J}({\bm\mu }_t)\delta {\bm\mu }_t. \end{equation}\end{document} $${\stackrel{\u0303}{\mathbf{\Gamma}}}_{m}-\stackrel{\u0303}{\mathbf{\Gamma}}\left({\mathit{\mu}}_{t}\right)=\mathbf{J}\left({\mathit{\mu}}_{t}\right)\delta {\mathit{\mu}}_{t}.$$## Eq. 7

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} {\bm\mu }_{t + 1} = {\bm\mu }_t + \delta {\bm\mu }_t.\end{equation}\end{document} $${\mathit{\mu}}_{t+1}={\mathit{\mu}}_{t}+\delta {\mathit{\mu}}_{t}.$$## 2.2.2.

#### Optode calibration employing coupling coefficients

Actually, DOT reconstructs the change in the optical properties between the task and rest states. Therefore,
[TeX:]
${{\bm{\tilde \Gamma}} }_m$
${\stackrel{\u0303}{\mathbf{\Gamma}}}_{m}$
consists of the ratios of the measured intensities at the task states to those at the rest states. If the optode coupling between the optodes and the object's surface at the time of the task measurement is different from that of the rest measurement, then the change in the optode coupling will appear as artifacts in the reconstructed images. To reduce artifacts induced by a change in optode coupling, coupling coefficients are introduced into the reconstruction algorithm. To calibrate the changes in the optode coupling, the calculated featured data are multiplied by coupling coefficients^{8, 9} as,
[TeX:]
$\alpha _i \beta _j \tilde \Gamma _{ij}$
${\alpha}_{i}{\beta}_{j}{\stackrel{\u0303}{\Gamma}}_{ij}$
, where
[TeX:]
$\tilde \Gamma _{ij}$
${\stackrel{\u0303}{\Gamma}}_{ij}$
are the elements of the calculated featured data obtained at the *i*’th source and *j*’th detector positions, and *α*
_{i} and *β*
_{j} are the coupling coefficients to calibrate the change in the optode coupling at the *i*’th source and the *j*’th detector positions, respectively. Note that the featured data multiplied by the coupling coefficients are equivalent to the time-resolved data multiplied by the coupling coefficients because the Laplace transform linearly transforms the time-resolved data. No change in the optode coupling is represented by the coupling coefficients of unity, and the losses of light due to the changes in the coupling result in the coupling coefficients becoming smaller than unity.

Because it is impossible to know the values of the coupling coefficients in these experiments, they are unknowns that must be reconstructed similarly to the optical properties. The coupling coefficients and the optical properties are simultaneously optimized to minimize the residual error by iteration. By introducing the coupling coefficients to the measured data for task state, the Jacobian matrix and the vectors of the optical properties and their small changes appearing in Eq. 5 are modified as follows:

## Eq. 8

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} {\bf J} = [{\bf J}({\bm\mu }),{\bf J}(\alpha),{\bf J}(\beta)], \end{equation}\end{document} $$\mathbf{J}=\left[\mathbf{J}\right(\mathit{\mu}),\mathbf{J}(\alpha ),\mathbf{J}(\beta \left)\right],$$## Eq. 9

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \delta {\bm\mu } = [\delta {\bm\mu }^{\rm T},\delta \alpha ^{\rm T},\delta \beta ^{\rm T}]^{\rm T}, \end{equation}\end{document} $$\delta \mathit{\mu}={[\delta {\mathit{\mu}}^{\mathrm{T}},\delta {\alpha}^{\mathrm{T}},\delta {\beta}^{\mathrm{T}}]}^{\mathrm{T}},$$## Eq. 10

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} {\bm\mu } = [{\bm\mu }^{\rm T},\alpha ^{\rm T},\beta ^{\rm T}]^{\rm T}, \end{equation}\end{document} $$\mathit{\mu}={[{\mathit{\mu}}^{\mathrm{T}},{\alpha}^{\mathrm{T}},{\beta}^{\mathrm{T}}]}^{\mathrm{T}},$$*α*and

*β*are the vectors consisting of the coupling coefficients

*α*

_{i}and

*β*

_{j}, respectively,

**J**( [TeX:] $\bm\mu$ $\mathit{\mu}$ ),

**J**(

*α*), and

**J**(

*β*) are the Jacobian matrices obtained as the first derivatives of the featured data with respect to [TeX:] $\bm\mu$ $\mathit{\mu}$ ,

*α*, and

*β*, respectively, and

*δα*and

*δβ*are the vectors representing the small changes in

*α*and

*β*, respectively. The image reconstruction is conducted with the modified variables expressed by Eqs. 9, 10. The vectors [TeX:] $\bm\mu$ $\mathit{\mu}$ ,

*α*, and

*β*are the unknowns to be reconstructed.

## 3.

## Validations of Optode Calibration by Numerical Simulations

## 3.1.

### Reconstruction without Optode Calibration

Numerical simulations are carried out to validate the effectiveness of the optode calibration stated above. Numerical phantoms are 2D circular media with a radius of 40 mm, having the reduced scattering coefficient, *μ*
_{s,B}
^{′} = 1.0 mm^{−1}, and the absorption coefficient, *μ*
_{a,B} = 0.005 mm^{−1}, as the background optical properties. The refractive index is given as 1.33. The circular media at the rest state are homogeneous, and those at the task state have targets of a strongly absorbing circular region with a radius of 5 mm, *μ*
_{a,T} = 0.01 mm^{−1} (Δ*μ*
_{a} = *μ*
_{a,T} – *μ*
_{a,B} = 0.005 mm^{−1}), and *μ*
_{s,T}
^{′} = 1.0 mm^{−1}. *μ*
_{s}
^{′} is assumed not to change (*μ*
_{s,T}
^{′} = *μ*
_{s,B}
^{′}) and only *μ*
_{a} is reconstructed.

Sixteen optodes (Nos. 1 to 16) working as both the sources and detectors are placed at the medium surface with equal spacing. One of the optodes is selected as the source, and the re-emitted light of the source is detected by the other optodes. The optodes are then sequentially selected as the source. Therefore, a total of 16 ×15 time-resolved measurement data is numerically generated for one image.

The simulated measurement data are generated by solving Eq. 1 in 2D with FEM, employing the meshing of 2048 elements and 1089 nodes, whereas the FEM calculation for the image reconstruction employs the meshing of 816 elements and 473 nodes, as shown in Fig. 1. The positions and numbering of the optodes and the coordinates of *x* and *y* are also shown in Fig. 1.

For the first step, reconstructions are performed for three cases, (a), (b), and (c), without a change in the optode coupling to evaluate the performance of the mGPST algorithm in the absence of optode calibration. The coordinates of the target centers of the strongly absorbing circle are (*x, y*) = (0, 0), (20, 0), and (35, 0) for cases (a), (b), and (c), respectively.

The true images of *μ*
_{a} at the task state, the reconstructed images, and the profiles of the reconstructed *μ*
_{a} along the *x*-axis without a change in optode calibration are shown in Fig. 2. The reconstructed differences in *μ*
_{a} (Δ*μ*
_{a}) between the task and rest states are about 10% to 50% of the true values, depending on the cases. It is generally known that in DOT images, the reconstructed changes in *μ*
_{a} are generally smaller than their true values,^{11} whereas the true positions of the target are well recovered (Fig. 2). The reconstructed image for case (a) has an oscillating distribution, which typically appears in inverse problems. We found that the mGPST can reconstruct *μ*
_{a} well in these cases without a change in optode coupling.

## 3.2.

### Reconstruction with Optode Calibration

Image reconstructions of *μ*
_{a} with the incorporation of optode calibration were performed. The numerical phantoms are identical to those of cases (a), (b) and (c) described above.

Changes in optode coupling are assumed to take place at optode No. 1 located at (*x, y*) = (40, 0), when the measurements are performed for the task state. Because of the change in optode coupling, the measured light intensities are assumed to change as follows:

For case (a): When optode No. 1 is the source, the intensities measured by all other detectors decrease by 10%. When optodes other than No. 1 are the source, the intensity measured by optode No. 1 decreases by 10%, and the intensities measured by all the other optodes are not affected.

For cases (b) and (c): Similar to case (a), with the exception that the decreases in the measured intensities are 3%.

The results of the reconstructions are shown in Figs.
3, 4, 5, for cases (a), (b), and (c), respectively. Each figure provides the reconstructed *μ*
_{a} images with and without optode calibration, the true and reconstructed *μ*
_{a} distributions along the *x*-axis, and the reconstructed *α* and *β*. In case (a), when the optode calibration is not incorporated, the artifact in *μ*
_{a} is approximately 0.007 mm^{−1} and appears in the vicinity of the boundary of the circular medium as shown in Fig. 3(a-1). This artifact is stronger than the target image at the center of the medium and is caused by a decrease in the measured intensities, which can be interpreted as an increase in the absorption coefficient by the reconstruction algorithm in the absence of optode calibration. On the other hand, the image reconstructed with optode calibration shows no artifact, as seen in Fig. 3(a-2). In Fig. 3(a-3), where reconstructed distributions are represented along the *x*-axis with and without optode calibration, it was found that regardless of the incorporation of optode calibration, the reconstructed difference in *μ*
_{a} (Δ*μ*
_{a}) at the target is approximately 0.0013 mm^{−1}, which is about one-fourth of the true Δ*μ*
_{a} of 0.005 mm^{−1}.

The coupling coefficients were simultaneously estimated in the reconstruction process, as plotted in Fig. 3(a-4). The coupling coefficients *α*
_{1} and *β*
_{1} for optode No. 1 are reconstructed as approximately 0.90. These values were almost identical to the 10% decrease observed in the source and measured intensities of optode No. 1 provided in the data generation. All other coupling coefficients were found to be approximately 1.01, which suggests that the optode coupling did not change. The optode calibration is shown to function well to recover the decrease in the source and measured intensities.

For cases (b) and (c), reconstruction without optode calibration shows a decrease in the source and measured intensities as reflected by the changes seen in *μ*
_{a}, which appear as artifacts as shown in Figs. 4(b-1) and 5(c-1), respectively. The optode calibration removes artifacts as shown in Figs. 4(b-2) and 5(c-2). Similar to case (a), we observed comparable changes in the reconstructed *μ*
_{a} distributions along the *x*-axis with optode calibration, as shown in Figs. 4(b-3) and 5(c-3).

The estimated coupling coefficients are also plotted in Figs. 4(b-4) and 5(c-4). Because the observed decreases in the source and measured intensities in cases (b) and (c) were 3%, which is smaller than the 10% observed in case (a), the coupling coefficients reconstructed in cases (b) and (c) are closer to unity. In case (b), *α*
_{1} and *β*
_{1} for optode No. 1 are 0.975, which is very close to 0.970, reflecting the 3% decrease in the source and measured intensities of optode No. 1. *α*
_{i} and *β*
_{j} for optodes other than No. 1 are very close to 1.00. In case (c), however, *α*
_{1} and *β*
_{1} were 0.921, which overestimates the 3% decreases in the source and measured intensities of optode No. 1, whereas *α*
_{i} and *β*
_{j} for optodes other than No. 1 were very close to 1.00. In case (c), it is difficult for the algorithm to discriminate between the artifact and the target, which coincidently appears in the same position. In this case, the reconstructed maximum *μ*
_{a} values of the target with and without optode calibration were 6.77 × 10^{−3} and 8.24 × 10^{−3} mm^{−1}, respectively, as shown in Fig. 5(c-3), whereas that without the change in optode coupling is 7.60 × 10^{−3} mm^{−1}, as shown in Fig. 2. The artifact is overestimated, and the coupling coefficients, *α*
_{1} and *β*
_{1}, become smaller compared to those observed in case (b). A similar result has been shown by Schweiger
^{9} When the target is within the vicinity of the surface, calibration should be used with caution.

## 3.3.

### Validation of Optode Calibration by Phantom Experiments

Phantom experiments using a time-resolved measurement system were conducted to validate optode calibration. The time-resolved measurement system consisted of 16 optodes composed of fiber bundles, ultra-short pulse lasers and time-correlated single photon counting units. Each optode worked as both the source and the detector. The ultrashort pulse lasers generated light pulses having wavelengths of 759 and 834 nm, duration of approximately 100 ps, a mean power of 0.25 mW, and a repetition rate of 5 MHz. The light pulses illuminated the measured object through one of the optodes while all of the optodes collected light that was re-emitted from the object. It took approximately 120 s to acquire data for one DOT image. The details of the system may be found in the literature.^{19}

The measured object, the phantom, was a cylinder made of polyacetal resin covered by a silicon rubber, as shown in Fig. 6. Covering the solid resin with the soft rubber simulated deformable human skin. We have confirmed by preliminary experiments that the soft surface causes the change in the coupling between the optode and the measured object more significantly than the hard surface.^{20} The phantom was also covered by a spongy black rubber with a thickness of about 10 mm to prevent stray light from being detected. The optodes were in contact with the phantom surface through holes in the spongy black rubber. These holes had a diameter slightly larger than those of the optodes.

The polyacetal cylinder had a radius of 35 mm and the optical properties of *μ*
_{a} = 0.00066 mm^{−1} and *μ*
_{s}
^{′} = 0.863 mm^{−1}. The silicone rubber covering the cylinder had a thickness of 5 mm, *μ*
_{a} = 0.092 mm^{−1}, and *μ*
_{s}
^{′} = 0.64 mm^{−1}.

The time-resolved measurements were made at the rest state having homogeneous optical properties of the background, and at the task state, having a target of a small cylinder with *μ*
_{a} higher than the background. At the rest state measurement, all the optodes were attached to the phantom surface. After the measurement of the rest state, the cylindrical region with a radius of 10 mm and a center at (*x, y*) = (20, 0) in the phantom was replaced by an Intralipid solution with *μ*
_{a} = 0.0026 mm^{−1} and *μ*
_{s}
^{′} = 1.05 mm^{−1} for a regional change in *μ*
_{a} at the task state. At the same time, optode No. 1 was moved 2 mm away from the phantom surface to induce changes in the optode coupling at the task state.

The measured time-resolved data for two states of the phantom were transformed into featured data, and their ratios were used as input data for image reconstruction of the change in *μ*
_{a}. The distributions of *μ*
_{s}
^{′} were known and not reconstructed for simplicity. The reconstructions were conducted with and without optode calibration.

The results of the reconstructions are shown in Fig. 7. The reconstruction without optode calibration results in an image with a strong artifact near the moved optode No. 1, as shown in Fig. 7a. Weak and unexpected artifacts were also observed at other positions. On the other hand, image reconstruction with optode calibration successfully removed most of the artifacts. The target was clearly reconstructed with *μ*
_{a} at approximately 50% of the true value. The plot of *μ*
_{a} along the *x*-axis shown in Fig. 7c indicates that the artifact is reduced by approximately 70% of the reconstructed image without optode calibration. Figure 7d shows the reconstructed values of *α* and *β* and the values of *α*
_{1} and *β*
_{1} for optode No. 1, which largely deviate from unity to compensate for the change in optode coupling. The values of *α* for optode No. 3 and No. 16 also deviate from unity and are close to *α*
_{1} for optode No. 1, which may be due to the interference between optode No. 1 and the optodes close to No. 1. However, the reason for the deviations of *α*
_{3} and *α*
_{16} from unity is not clear because *α*
_{2} for optode No. 2 is not affected. Some of the reconstructed coupling coefficients are larger than unity; this may indicate that the changes in the coupling increased the measured light intensities for some unknown reasons. One possible explanation of the coupling coefficients being larger than unity may be that light was reflected at the interface between the polyacethal resin and the silicon rubber. In the light propagation model, the reflection of light at the interface was not incorporated.

## 4.

## Calibrated DOT Image of the Absorption Coefficient from *In Vivo* Measurement Data

To further study the effectiveness of the optode calibration, the method stated above was applied to *in vivo* measurement data to reconstruct the image of *μ*
_{a}. The *in vivo* measurement data was provided from *in vivo* experiments of time-domain DOT performed on the head of a premature infant at 29 days after birth. The gestational age of the infant was 27 weeks, and his birth weight was 658 g.

The time-resolved measurement system was used again, and 16 optodes were attached to the surface of the head by the use of a circular optode holder. The circumference length of the subject head was 23.3 cm at the measurement plane. The contacts of the optodes to the surface of the skin were controlled by springs to keep the contact pressure constant during measurement. However, because the subject was not anesthetized, there was the possibility that the optode coupling had changed because of the subtle movement of the subject. An artificial respirator assisted the subject's respiration. The measurement data was acquired for the following two conditions: condition I: normocapnic condition with a respiratory rate of 25 breaths/min and a maximum inspiratory pressure of 1.6 kPa above the atmospheric pressure, and condition II: hypocapnic condition with a low arterial partial pressure of CO_{2} (PaCO_{2}) and a respiratory rate of 40 breaths/min. The maximum inspiratory pressure was the same as that of condition I. We anticipated that this hypocapnic condition would induce a slightly lowered blood volume in the brain, leading to a slight decease in *μ*
_{a} in the brain. The measurement under condition II was conducted about 10 min after the measurement under condition I. Written informed consent was obtained from the parents of the infant. This experiment was conducted upon the approval of the Kagawa University Ethics Committee, the National Institute of Advanced Industrial Science & Technology (AIST) Ethical Committee on Medical Engineering Experiments and the AIST Tsukuba Center East Office Human Engineering Experiment Committee.

Two-dimensional images of *μ*
_{a} at a wavelength of 760 nm were reconstructed with and without optode calibration. The ratios of the featured data obtained at condition II to those obtained at condition I were used for the reconstructions as well as the phantom experiments described in Sec.
3.3. The initial guess of *μ*
_{a} was given as 0.0142 mm^{−1} uniformly in the head. The reduced scattering coefficient *μ*
_{s}
^{′} was fixed as 0.6 mm^{−1}. The initial guess of *μ*
_{a} and fixed value of *μ*
_{s}
^{′} were provided by previous studies described in the literature.^{21}

Figure 8a shows the residual errors
[TeX:]
$||{{\bm{\tilde \Gamma}} }_m \, - {{\bm{\tilde \Gamma}} }\,({\bm\mu }_t)||$
$\left|\right|{\stackrel{\u0303}{\mathbf{\Gamma}}}_{m}\phantom{\rule{0.16em}{0ex}}-\stackrel{\u0303}{\mathbf{\Gamma}}\phantom{\rule{0.16em}{0ex}}\left({\mathit{\mu}}_{t}\right)\left|\right|$
as a function of the number of iterations. In the reconstruction with optode calibration, the residual error converged faster and became smaller compared to the reconstruction without optode calibration. The reconstructed coupling coefficients are plotted in Fig. 8b. The reconstructed *μ*
_{a} images are shown in Figs. 8c and 8d without and with optode calibration, respectively. The dots along the boundary indicate the optode positions. When comparing the reconstructed *μ*
_{a} images with and without optode calibration, it is clearly apparent that the artifacts that appear within the vicinity of the boundary were remarkably reduced by optode calibration. In particular, there were very small values of *μ*
_{a} found in the left frontal lobe and larger values of *μ*
_{a} observed in the left occipital lobe in the image reconstructed without optode calibration. These artifacts were reduced by optode calibration, and as a result, the decrease in *μ*
_{a} from the assumed initial value of 0.0142 mm^{−1} was clearly seen in the center of the head, as anticipated by lowering the PaCO_{2} in condition II. Thus, the optode calibration functioned effectively from the perspective of the *μ*
_{a} images. Further studies will be needed to show the effectiveness of optode calibration. This may be done by physiologically checking the reconstructed images of the blood volume and oxygenation states, as these images are obtained from the *μ*
_{a} images at two wavelengths.

## 5.

## Discussion

Optode calibration has previously been demonstrated to be effective in eliminating artifacts caused by a change in optode coupling in cw- and frequency-domain DOT. This study has applied optode calibration to time-resolved DOT using the modified GPST algorithm. Numerical simulations, phantom experiments, and *in vivo* experiments have validated its effectiveness. The changes in optode coupling that occur during measurements often cause artifacts in the reconstructed images appearing within the vicinities of the optodes, as shown in Figs.
2, 3, 4, 5. This may be explained by the following: when an optode initially is attached correctly to the object surface at the time of rest measurement and if the optode is separated slightly from the object's surface at the time of task state measurement, light intensities (amplitudes) measured at the task state are recorded smaller than those when the optode is kept attached correctly to the object's surface. The algorithm used may misunderstand that the attenuated light intensities are attributed to strong absorbers close to the optodes. As a result, artifacts with large absorption coefficients appear within the vicinities of the optodes. Conversely, when there are initially small gaps between the optodes and the object's surface at the time of rest state measurement and if the optodes are then attached correctly to the object's surface at the time of task state measurement because of the movement made by the object, light intensities measured at the task state will be recorded larger than those when the small gaps are kept even at the time of task state measurement. In these cases, artifacts with small absorption coefficients appear within the vicinities of the optodes. Both cases may take place simultaneously, particularly during *in vivo* experiments in which the subjects can easily move. The coupling coefficients *α* and *β* can correct for changes in optode coupling for both cases, as stated above.

To reduce the artifacts incurred by changes in optode coupling, the featured data are multiplied by the coupling coefficients to compensate for the light intensities and are incorporated as additional unknowns in the image reconstruction algorithm. The optimization of the coupling coefficients is conducted simultaneously with the image reconstruction. A possible disadvantage of adding unknowns is that it can aggravate the ill-posed nature of the inverse problem. However, the number of the coupling coefficients, which total 32 in this study, is much smaller than the number of the reconstructed absorption coefficients in the object, which is the same as the number of FEM nodes (473). Thus, the influence of adding coupling coefficients to the performance of the algorithm seems insignificant.

Another idea to reduce artifacts by changes in the optode coupling is to not use the optode measurement data, which reveals changes in the optode coupling. However, it is difficult to determine which optodes are associated with the changes in coupling. By using the coupling coefficient, one can fully utilize important information that is included in the measured data.

Now we can discuss what happens when more than one optode suffers from changes in coupling. During *in vivo* experiments, multiple optodes suffer from changes in coupling, and as a result, multiple artifacts appear in the image without optode calibration, as shown in Fig. 8c. These artifacts appear separately within the vicinities of the optodes in principle, but if two optodes with poor coupling existed next to each other, the two artifacts close to the two optodes may be combined into one large artifact. In addition, the artifacts are reduced by optode calibration to some extent, as shown in Fig. 8d, and the values of *α* and *β* are randomly distributed, as shown in Fig. 8b. Therefore, artifacts that are caused by changes in optode coupling, which appear separately or combined within the vicinities of the optodes and multiple artifacts, are dramatically reduced by the method of optode calibration.

Because the propagation probability of light in a medium between two optodes is highest at the positions close to the surface around the optode, as shown by the well-known banana shape, the measurement sensitivity of absorbing targets is greatest close to the surface around the optodes. Therefore, optode calibration may also be acting as a method to correct the measured data for the effects of uneven sensitivity in the medium. The reconstructed image without optode calibration is accompanied by wavy artifacts that have 8 periods in an angular direction with 16 optodes [Fig. 2a]. These artifacts, seen in Fig. 2a, disappear in Fig. 3(a-2), where the image has been reconstructed with optode calibration. Therefore, we hypothesized that optode calibration may correct for uneven sensitivity in the medium. However, it is doubtful that the depth-dependent sensitivity is corrected by simple multiplication of *α* and *β* to the time-resolved measurement data. Thus, the coupling coefficients *α* and *β* are assumed to be constant and have no role in modifying the time-resolved profiles, in which the early and late photons propagate in the shallow and deep layers in the medium, respectively. Also, the inversion process often produces wavy images in nature. Therefore, the wavy artifacts in Fig. 2a may coincidently have 8 periods, which is half the number of optodes. The possibility of correcting for the uneven sensitivity in the medium should be further investigated in the future.

## 6.

## Conclusions

Optode calibration eliminates artifacts by a change in optode coupling that is applied to the time-resolved DOT algorithm based on the modified GPST, and its effectiveness is validated by numerical simulation, phantom experiments, and *in vivo* experiments. Some numerical simulations of the time-domain DOT with optode calibration have shown that artifacts can be effectively reduced and that the targets represented by a change in the absorption coefficient are localized at the correct positions. Experiments of time-domain DOT using the phantoms with a soft surface have also validated the performance of the reconstruction using optode calibration incorporated into the modified GPST algorithm. The reconstructed coupling coefficients of the optode are subject to a change in the optode coupling between the rest and task states, which accurately reflects the change in optode coupling, for example, decreases in the source and measured light intensities between the rest and task states. However, when the target in the object is located close to the object's surface, we found that optode calibration overestimates the change in optode coupling and should be used with caution. *In vivo* experiments with optode calibration performed on the head of a premature infant functioned effectively to reduce the number of artifacts in the *μ*
_{a} image. However, further studies are needed to show the effectiveness of optode calibration by checking the reconstructed images of blood volume and oxygenation states from the physiological point of view. It should also be noted that this algorithm has been validated by numerical simulations, phantom experiments, and *in vivo* experiments.

## Acknowledgments

This work was partly supported by JSPS KAKENHI [Grant-in-Aid for Scientific Research (B)], (No. 17360095) “Investigation of oxygen transportation in biological tissues by reflection-type diffuse optical tomography” (2005 to 2006) and (No. 17390307) “Study of blood volume and oxygenation in brains of newborn infants using NIR tomographic imaging” (2005 to 2006).

## References

**,” Phys. Med. Biol., 50 R1 –R43 (2005). https://doi.org/10.1088/0031-9155/50/4/R01 Google Scholar**

*Recent advances in diffuse optical imaging***,” IEEE Signal Process Mag., 18 (6), 57 –75 (2001). https://doi.org/10.1109/79.962278 Google Scholar**

*Imaging the body with diffuse optical tomography***,” Inverse Problems, 15 R41 –R93 (1999). https://doi.org/10.1088/0266-5611/15/2/022 Google Scholar**

*Optical tomography in medical imaging***,” Scand. J. Clin. Lab. Invest. Suppl., 188 9 –17 (1987). Google Scholar**

*Cerebral monitoring in newborn infants by magnetic resonance and near infrared spectroscopy***,” Phys. Med. Biol., 47 4155 –4166 (2002). https://doi.org/10.1088/0031-9155/47/23/303 Google Scholar**

*Three-dimensional optical tomography of the premature infant brain***,” Dev. Sci., 5 (3), 371 –380 (2002). https://doi.org/10.1111/1467-7687.00376 Google Scholar**

*Basic principles of optical imaging and application to the study of infant development***,” Neuroimage, 30 521 –528 (2006). https://doi.org/10.1016/j.neuroimage.2005.08.059 Google Scholar**

*Three-dimensional whole-head optical tomography of passive motor evoked responses in the neonate***,” Opt. Express, 8 (5), 263 –270 (2001). https://doi.org/10.1364/OE.8.000263 Google Scholar**

*Simultaneous imaging and optode calibration with diffuse optical tomography***,” Appl. Opt., 46 (14), 2742 –2756 (2007). https://doi.org/10.1364/AO.46.002743 Google Scholar**

*Image reconstruction in optical tomography in the presence of coupling errors***,” Appl. Opt., 39 (34), 6466 –6486 (2000). https://doi.org/10.1364/AO.39.006466 Google Scholar**

*Instrumentation and calibration protocol for imaging dynamic features in dense-scattering media by optical tomography***,” J. Opt. Soc. Am. A, 19 (10), 1983 –1993 (2002). https://doi.org/10.1364/JOSAA.19.001983 Google Scholar**

*Source-detector calibration in three-dimensional Bayesian optical diffusion tomography***,” Appl. Opt., 42 (16), 3154 –3162 (2003). https://doi.org/10.1364/AO.42.003154 Google Scholar**

*Optode positional calibration in diffuse optical tomography***,” Appl. Opt., 44 (10), 1879 –1888 (2005). https://doi.org/10.1364/AO.44.001879 Google Scholar**

*Computational calibration method for optical tomography***,” IEICE Trans. Inf. Syst., E85-D (1), 133 –142 (2002). Google Scholar**

*Time-resolved diffuse optical tomography using a modified generalized pulse spectrum technique***,” J. Biomed. Opt., 12 (6), 062107 (2007). https://doi.org/10.1117/1.2815724 Google Scholar**

*Time-resolved diffuse optical tomography and its application to**in vitro*and*in vivo*imaging**,” Phys. Rev. E, 50 (5), 3634 –3640 (1994). https://doi.org/10.1103/PhysRevE.50.3634 Google Scholar**

*Diffusion approximation for a dissipative random medium and the applications***,” Appl. Opt., 22 (16), 2456 –2462 (1983). https://doi.org/10.1364/AO.22.002456 Google Scholar**

*Scattering and absorption of turbid materials determined from reflection measurements. 1: Theory***,” J. Math. Imaging Vision, 3 263 –283 (1993). https://doi.org/10.1007/BF01248356 Google Scholar**

*Application of the finite-element method for the forward and inverse models in optical tomography***,” Rev. Sci. Instrum., 70 (9), 3595 –3602 (1999). https://doi.org/10.1063/1.1149965 Google Scholar**

*Multichannel time-resolved optical tomographic imaging system***,” 139 –140 (2009). Google Scholar**

*Effect of the movement of probes on the images of diffuse optical tomography***,” Pediatr. Res., 58 568 –573 (2005). https://doi.org/10.1203/01.PDR.0000175638.98041.0E Google Scholar**

*Developmental changes of optical properties in neonates determined by near-infrared time-resolved spectroscopy*