Time-Domain Functional Diffuse Optical Tomography System Based on Fiber-Free Silicon

Based on recent developments in both single-photon detectors and timing electronic circuits, we designed a compact and cost effective time-domain diffuse optical tomography system operated at 1 Hz acquisition rate, based on eight silicon photomultipliers and an 8-channel time-to-digital converter. The compact detectors are directly hosted on the probe in a circular arrangement around a single light injection fiber, so to maximize light harvesting. Tomography is achieved exploiting the depth sensitivity that is encoded in the arrival time of detected photons. The system performances were evaluated on simulations to assess possible the limitations arising from the use of a single injection point, and then on phantoms and in vivo to prove the eligibility of these technologies for diffuse optical tomography.


Introduction
Diffuse Optical Tomography (DOT) in the near infrared spectrum is a valuable tool for many applications, like breast cancer analysis [1][2][3][4][5], functional imaging [6,7], small-animals imaging [8], and non-destructive assessment of food quality [9,10].In fact, the main chromophores that are present in biological tissues feature a low absorption coefficient (µ a < 0.1 mm −1 ) and a high reduced scattering coefficient (µ s ≈ 1 mm −1 ) to radiation within the so-called therapeutic window (600-1100 nm), thus permitting light to propagate for long pathlengths and to be detected when re-emitted outside of the tissue.
DOT relies on sampling the external boundary of the medium with a distribution of sources and detectors [11,12].Depending on the light modulation, different approaches to DOT can be employed: continuous-wave (CW), frequency-domain (FD), and time-domain (TD).
Due to the lower cost and complexity, DOT systems requiring a large number of sources and detectors have been implemented mainly using CW technologies.Important achievements were obtained, for instance, in the field of neuromonitoring [6], breast imaging [13,14], tomography of osteorticular diseases [15], and peripheral arterial disease [16].CW approaches have permitted an easier exploration of biomedical applications, but convey a lower amount of information as compared to FD or TD techniques, in terms, for instance, of coupling between absorption and reduced scattering coefficients or strong dependence on superficial effects, movement artifacts, and optical coupling.This requires compensating lower information content with a larger number of optodes, exploiting for instance multiple source-detector distances to discriminate scalp-skull systemic fluctuations from true brain activation [17,18].The explosion of source-detector pairs can be challenging in view of the fascinating evolution toward a fiberless, whole-scalp diffuse optical tomography systems [19].In general, the possibility to collect multiple information from the same optode is an interesting challenge to increase quantitation in DOT, reduce artefacts and confounding effects, and improve robustness in clinical settings or in natural/social situations [20].
The TD approach, consisting of the injection of picoseconds pulses and time-domain detection, can indeed expand the information content carried by the same optode.TD techniques have been demonstrated to be potentially sensitive and capable to localize deep perturbations (depth > 3 cm) when sources and detectors are placed on the same surface of the sample (reflectance configuration) [21].The unique feature of TD approaches is the possibility to increase the depth sensitivity for increased photon traveling time.Thus, for the same source-detector pair, different measurements with different spatial sensitivity can be collected.Consequently, tomographic reconstruction can be achieved in the TD domain with a lower number of optodes as compared, for instance, to CW.
Due to the complexity of the system and the cost of sources, detectors, and timing electronics, the TD approach has not been widely used in the DOT field.In fact, only a few systems in literature are capable of a full TD-DOT [22][23][24].These systems are quite sophisticated and expensive, being based on parallel chains that are based on time-correlated single-photon counting (TCSPC) boards and Photomultiplier Tubes (PMTs).Moreover, they are cumbersome because every PMT needs its own cooling and power supply system.Further, light harvesting through fibers bundles adds to the complexity of the instruments.Another limitation of TD systems is the low measurement speed since the single-photon counting acquisition is dominated by Poisson statistics, thus requiring a large-area, low noise, and high-quantum efficiency time-domain detectors.All of these aspects have limited the number of in vivo and real-time applications.
During the last few years, Silicon Photomultipliers (SiPMs) have been proposed and successfully validated as novel compact solid-state detectors for diffuse optics [21,25].Both CW and TD diffuse optics techniques are beginning to take advantage of SiPMs, since they hold the best features of PMTs (e.g., large active area) and those of solid-state detectors (high photon detection efficiency, low cost and complexity, low bias voltage, small size, robustness, and insensitiveness to magnetic fields) [26][27][28][29].When considering, in particular, TD diffuse optics, recent studies demonstrated impressive improvements in single-photon timing resolution (down to <60 ps Full-Width at Half Maximum, FWHM) [30] and the possibility to fabricate small electronic circuits that are able to efficiently control the detector when used in contact with the tissue under investigation [31].In this way, the collection capability of diffused light is maximized thanks to the exploitation of the whole numerical aperture of the detector.This technology, already allows for the easy fabrication of miniaturized time-domain single-photon detectors at a very low cost (<$100).
Interesting technologies for TD diffuse optics are represented by counters [32] and Time-to-Digital Converters (TDCs) [33].However, only the latter devices are becoming commercially available at a cost lower than $1000 per timing channel, with a high count rate, and in principle they can be profitably employed for TD diffuse optics in general, provided that proper solutions are devised to cope with the high non-linearity of their time scale (e.g., implementation of procedures for correction of the non-linearity using a reference measurement with well-known shape).
The combination of these technologies is progressively decreasing the gap between CW and TD instruments in terms of size and cost, thus possibly leading to a widespread diffusion of the TD technique thanks to its inherent advantages.This breakthrough will occur when also pulsed laser sources will undergo similar technological improvements and scaling.To some extent, this development on sources is already started as well, with promising proof-of-concept results [21,34,35].
In this paper, we present the first compact TD DOT system that is capable of providing functional tomographies at 1 Hz rate exploiting eight fiber-free detection channels based on SiPMs and an 8-channel TDC.The aim of this work is to demonstrate that these technologies are nowadays eligible for the design of small multichannel diffuse optical tomography instruments, thus paving the way to a new class of systems featuring substantial improvement in terms of dimension and complexity.
The use of SiPMs for TD DOT was already proposed in a previous work [36].Yet, in that case: (i) SiPMs were still embedded into a table-top module [37]; (ii) the detection system was still based on optical fibers; (iii) measurements were performed in a scanning mode using a 2-channel system resulting in relatively long measurement time (some minutes); and, (iv) the timing electronic circuits were traditional TCSPC boards.Vice versa, here, we explore a simple approach, with eight SiPMs that are directly embedded onto the probe in a circular arrangement and a single source in the center, permitting continuous data acquisition at 1 Hz.The adoption of a single injection point is clearly a limitation, which is partially compensated by the in-depth scanning capability that is provided by TD acquisition.Indeed, this scheme was chosen to demonstrate the feasibility of a basic tomographic setup, focusing on the detection chain.In any case, the system can be easily expanded by adding an optical fiber switch to take advantage of the multiple injection points to improve tomographic reconstructions or embedding other laser sources for spectroscopic purposes.
In the following, we first describe the system architecture.Then, we present simulated results of the expected performance of the system.Further, we also report the measurements that are carried out on solid phantoms with a switchable inclusion.Finally, we show a proof-of-concept of in vivo brain functional imaging.

Setup
Figure 1 (left) shows a block schematic of the designed instrument.The core part is the probe, which embeds the eight SiPM channels.Each channel is a custom-made optode based on a SiPM, whose design, characterization, and preliminary in vivo use were already presented in Ref. [31].To ensure uniform performances among channels, for this work we selected detectors (S13081-050CS, Hamamatsu Photonics, Japan; active area of 1.3 × 1.3 mm 2 ; for other specifications see [38]) that showed similar breakdown voltage (i.e., average value of 52.2 V, with a maximum dispersion around ±100 mV among the various devices).A detailed characterization of the employed SiPMs is reported in [30].
The overall electronics that is needed to drive the detector is a small custom-designed Printed Circuit Board (PCB, dimension: 30 × 6 mm) that is hosted inside a small plastic cap.As it is described in Ref. [31], the PCB embeds: (i) a passive network for the sensing of the avalanche signal generated inside the SiPM due to a single-photon absorption, (ii) a radio-frequency Monolithic Microwave Integrated Circuit (MMIC) amplifier that provides an amplified copy of the avalanche signal to the probe output, (iii) a capacitor network to stabilize the bias of both the detector and the MMIC, and (iv) micro coaxial connectors to interface the optode with the system.Optodes are held in place by means of a drilled neoprene plate, following a circular geometry where a source optical fiber is placed in the center and detectors are placed in a circular arrangement around the source (with source-detector separation-i.e., the radius of the circumference-of 3 cm) without using optical fibers.The choice of this radius has been piloted by simulations to obtain the best compromise between probed volume, spatial resolution, and signal-to-noise ratio of in vivo measurements.As a laser source, we employed a 689 nm wavelength pulsed diode laser (PDL 800, Picoquant GmbH, Berlin, Germany) working at a repetition frequency of 40 MHz and featuring an average output power of about 1 mW.The light is injected into the sample by means of a 100 µm core optical fiber and the intensity is adjusted by a Neutral Density (indicated by n.d. in Figure 1) variable attenuator.Eight custom-designed signal converters PCBs translate the analog pulses of the SiPM probes (about 50 mV amplitude) into well detectable Low Voltage Transistor-Transistor Logic (LVTTL) signals.To reconstruct the Distribution of Time of Flights (DTOF) of photons, each converter output is connected to one input channel of a TDC (SC-TDC-1000/08S Surface Concept GmbH, Germany).The Nuclear Instrumentation Module standard laser synchronism is provided to the TDC using an additional signal converter to obtain a LVTTL signal.In this way, the eight channels simultaneously reconstruct the DTOF curves that are emitted by the sample where detectors are placed.The TDC is linked, thanks to a Universal Serial Bus 3.0 plug, to a personal computer, whose task it is to store and display the acquired data.Two benchtop power supply units provide all of the needed polarization voltages to the system: an 8 V voltage for the avalanche amplification network hosted on the probe and a 58 V supply to bias the detectors 5.8 V above their breakdown voltage.Finally, a small fixed-voltage enclosed power supply module provides the ±5 V supplies that are required by the LVTTL translator boards.Figure 1 (right) shows a picture of the overall system.It is worth appreciating the extremely compact dimensions of the eight-channel system (for further detail about the size of the building blocks see Table 1), which can of course be further squeezed by designing custom power supply modules instead of using benchtop devices.
The use of such TDC required software elaborations in order to overcome some intrinsic limitations, in particular: (i) high differential non-linearity (DNL) of the time scale (up to 100% of the least significant bit), and (ii) maximum acceptable rate for the synchronism pulse of 7 MHz, while the laser pulse rate is 40 MHz.
volume, spatial resolution, and signal-to-noise ratio of in vivo measurements.As a laser source, we employed a 689 nm wavelength pulsed diode laser (PDL 800, Picoquant GmbH, Berlin, Germany) working at a repetition frequency of 40 MHz and featuring an average output power of about 1 mW.The light is injected into the sample by means of a 100 µm core optical fiber and the intensity is adjusted by a Neutral Density (indicated by n.d. in Figure 1) variable attenuator.Eight customdesigned signal converters PCBs translate the analog pulses of the SiPM probes (about 50 mV amplitude) into well detectable Low Voltage Transistor-Transistor Logic (LVTTL) signals.To reconstruct the Distribution of Time of Flights (DTOF) of photons, each converter output is connected to one input channel of a TDC (SC-TDC-1000/08S Surface Concept GmbH, Germany).The Nuclear Instrumentation Module standard laser synchronism is provided to the TDC using an additional signal converter to obtain a LVTTL signal.In this way, the eight channels simultaneously reconstruct the DTOF curves that are emitted by the sample where detectors are placed.The TDC is linked, thanks to a Universal Serial Bus 3.0 plug, to a personal computer, whose task it is to store and display the acquired data.Two benchtop power supply units provide all of the needed polarization voltages to the system: an 8 V voltage for the avalanche amplification network hosted on the probe and a 58 V supply to bias the detectors 5.8 V above their breakdown voltage.Finally, a small fixed-voltage enclosed power supply module provides the ±5 V supplies that are required by the LVTTL translator boards.Figure 1 (right) shows a picture of the overall system.It is worth appreciating the extremely compact dimensions of the eight-channel system (for further detail about the size of the building blocks see Table 1), which can of course be further squeezed by designing custom power supply modules instead of using benchtop devices.
The use of such TDC required software elaborations in order to overcome some intrinsic limitations, in particular: (i) high differential non-linearity (DNL) of the time scale (up to 100% of the least significant bit), and (ii) maximum acceptable rate for the synchronism pulse of 7 MHz, while the laser pulse rate is 40 MHz.When the detectors are not exposed to the laser pulsed light, it is possible to acquire the background noise alone, which is expected to be completely flat over the TDC time scale, not being correlated with the laser synchronism.Thus, to overcome the issue of the high DNL, a dark count rate acquisition is used to calculate proper coefficients (one for each time bin), which allow us to linearize the TDC time scale.This procedure has to be done only once since the DNL is very stable over different days.As it concerns the maximum synchronism rate, we divide the 40 MHz synchronism of the laser by a factor of 16 thanks to a native start divider function that is embedded into the TDC module, obtaining 16 replicas of the DTOF.Since each replica conveys only 1/16 of the total number of detected photons, we sum all of the replicas by refolding the TDC time scale in order to have a single DTOF with all of the detected photons.Since the number of bins contained into a single replica of the DTOF is not integer (i.e., actual laser pulse rate is not a perfect multiple of the TDC bin size), before applying the refolding operation we resampled the TDC time scale by increasing the number of bins by a factor of 1000.Hence, the refolding operation that is applied on the re-binned timescale also allows for the reduction of the actual bin size in the final reconstructed waveform.
The system is controlled by a home-made software, written in ANSI C, and compiled in the LabWindows ® environment.The software controls all of the operations that are related to the measurements (e.g., automatic detection of optimal laser attenuation to cope with the single-photon statistics, data acquisition, real-time visualization, data save).This is a general purpose program to control time-resolved photonics systems, which has been evolving in the course of the last 25 years and is applied to more than 10 different instruments in various application fields (e.g., optical mammography, brain imaging, broadband diffuse spectroscopy, and time-resolved fluorescence spectroscopy).Although written in C language, it follows an object-oriented scheme, with new features being added as new objects of abstract classes.In particular, the new TDC functionalities were added as new acquisition board permitting a seamless operation as compared to standard TCSPC boards.TDC data are pre-processed for refolding over replicas and to compensate for non-linearity, with very little additional computation time (<20 ms), allowing us almost real-time reconstructions of a single curve.
Figure 1 (top left) shows the Instrument Response Function (IRF) of one channel with a dynamic range of about three decades and a FWHM of 270 ps.The reported IRF is representative of all eight channels, which all show a similar IRF.
For the key performances of the system, the reader can refer to the paper by R. Re et al.
[31] where probe-hosted detectors have been extensively characterized for diffuse optics application using well-established protocols for performance assessment, such as MEDPHOT and BIP ones [39,40].

Tomographic Reconstruction
The detected TD curves are binned in 10 temporal gates of about 220 ps width, and the center of gravity of every IRF is used as temporal reference.The inverse problem is solved under the Born Approximation assuming a linear relation between the purely absorption perturbation and the relative contrast between the measurement obtained with and without the perturbation [12]: where J is the Jacobian matrix, ∆µ a is the vector of the absorption perturbation for every voxel, and y m is the vector of the time-domain relative contrasts.The Jacobian matrix is calculated analytically for a semi-infinite space by means of solutions of the diffusion equation subjected to the Born approximation and to Extrapolated Boundary Conditions [41,42].For the reconstruction, the volume is divided into (2.5 × 2.5 × 2.5 mm 3 ) cubic voxels.For the phantom case, a 104 × 80 × 45 mm 3 volume is considered, corresponding to the physical volume of the sample.For the in vivo measurements, a 100 × 100 × 30 mm 3 volume is considered, well representing the expected probed region by photons with this source-detector configuration.Due to the ill-posedness of the problem stated in Equation ( 1), a Tikhonov regularization is added to stabilize the solution with respect to noise.A detailed description of the software can be found in [21].

Simulation Approach
Simulations are performed using analytical temporal point-spread functions (TPSFs) that are obtained under the Diffusion Approximation of the radiative transfer equation for a semi-infinite medium [21], either homogeneous or embedding localized low absorbing perturbations (Born approximation) [41,42].In particular, we used a semi-infinite model and considered a homogeneous background with µ a = 0.01 mm −1 , µ s = 1 mm −1 and an internal refractive index of 1.5 (epoxy resin, i.e., the main constituent of the phantom used for the following experiments).
TPSFs are simulated up to about 2.1 ns, and then they are area-normalized to 10 6 counts, and finally Poisson noise is added to mimic realistic measurement conditions.Indeed, in simulations, the main source of noise is represented by the Poisson one, which is proportional to the overall number of counts.
The heterogeneous data are generated by placing a cubic absorbing perturbation in different positions.The simulated absorption value of the perturbation is derived using the approach of the Equivalence Black Volume [43].

Phantom Test
We subsequently tested the system on a solid mechanically switchable phantom [44] that was realized at the Physics Department of Politecnico di Milano in collaboration with the Physikalisch-Technische Bundesanstalt of Berlin.This phantom allows to rapidly change the optical properties from homogeneous to multiple degrees of heterogeneity.It is composed by a block of material (approximately a parallelepiped of 104 × 80 × 45 mm 3 ), the background, with a cylindrical hole from one lateral side to the opposite surface, and whose center is at 15 mm of depth and several cylinders, the rods, slightly lower in diameter than the hole, but longer than its height.These rods can be adjusted inside the background and interchanged between them.Since each rod embeds a different purely absorption inclusion, by modifying their position, a corresponding different perturbation is induced in the phantom itself and in changeable positions.The homogeneous background is in epoxy resin (optical properties at 690 nm: µ a = 0.01 mm −1 and µ s = 1 mm −1 ), while the inclusions are totally absorbing black polyvinyl chloride cylinders with different volumes.The suitability of the totally absorbing objects to mimic realistic optical inhomogeneities was rigorously demonstrated on Monte Carlo simulations for a wide range of experimental conditions [43].The constructed phantom revealed discrepancies up to 30% with respect to the theoretical expected perturbation, which was possibly ascribed to uncertainties in the fabrication process.Still, within this limitation in accurate quantitation, the phantom is effective in testing the sensitivity of a measurement system.
The probe is placed on the upper surface of the phantom, while the inclusion is translated at a depth of 15 mm, traversing the whole region beneath the detectors' circle.Measurements are taken at steps of 5 mm for a total of 20 positions, repeating the acquisition 10 times for 1 s at each position.A homogeneous measurement (i.e., with the inclusion far from the detector area) is always taken as reference for unperturbed state.

In Vivo Measurements
We monitored the absorption changes induced by the hemodynamic cortical changes on three healthy adult volunteers (male, age: 40, 45 and 50 years) during a motor task (right hand crossed finger tapping task, in which the thumb touches in sequence the index, the ring finger, the middle finger, and the pinkie).A simple protocol was considered consisting of five trials each lasting 60 s: 20 s rest, 20 s finger tapping at self-pace, and 20 s recovery to rest condition.A folding average of the five repetitions was then computed.All of the volunteers gave their written informed consent to participate, and the protocol was approved by the Ethical Committee for clinical investigations of Hospital Clinic de Barcelona and was conducted in compliance with the Declaration of Helsinki.A total of eight measurement points were recorded, by placing the neoprene holder in the left hemisphere around the C3 position over the motor area according to the 10-20 International System.The holder was held by Velcro stripes.After the positioning, we waited for about 15 min before starting the measurement, so that the detectors could properly thermalize in order to stabilize their breakdown voltage.In this way, background noise, quantum efficiency, and timing jitter are stable during measurements, thus avoiding possible artifacts in the tomographies.The average acquisition rate was about 800 kcounts/s, for subject 1, 1200 kcounts/s for subject 2, and 850 kcounts/s for subject 3.

Simulations
To test the capability of the system in identifying the presence of a perturbation in different conditions, we performed a tomographic reconstruction on simulations changing, one parameter at a time.The reference system and the optodes scheme is depicted in Figure 2.
finger tapping task, in which the thumb touches in sequence the index, the ring finger, the middle finger, and the pinkie).A simple protocol was considered consisting of five trials each lasting 60 s: 20 s rest, 20 s finger tapping at self-pace, and 20 s recovery to rest condition.A folding average of the five repetitions was then computed.All of the volunteers gave their written informed consent to participate, and the protocol was approved by the Ethical Committee for clinical investigations of Hospital Clinic de Barcelona and was conducted in compliance with the Declaration of Helsinki.A total of eight measurement points were recorded, by placing the neoprene holder in the left hemisphere around the C3 position over the motor area according to the 10-20 International System.The holder was held by Velcro stripes.After the positioning, we waited for about 15 min before starting the measurement, so that the detectors could properly thermalize in order to stabilize their breakdown voltage.In this way, background noise, quantum efficiency, and timing jitter are stable during measurements, thus avoiding possible artifacts in the tomographies.The average acquisition rate was about 800 kcounts/s, for subject 1, 1200 kcounts/s for subject 2, and 850 kcounts/s for subject 3.

Simulations
To test the capability of the system in identifying the presence of a perturbation in different conditions, we performed a tomographic reconstruction on simulations changing, one parameter at a time.The reference system and the optodes scheme is depicted in Figure 2.
We studied three different variables, that are (1) real Δ , (2) radial distance (along y axis) from the injection source, and (3) depth along the z axis.

Effect of the Perturbation
We investigated how the entity of a perturbation affects the tomographic reconstruction by changing the heterogeneity absorption coefficient Δ of 0.04, 0.017, 0.01, and 0.005 mm at the depth of 15 mm, which is a depth where a brain activation is mostly expected.The tomographic reconstructions for x-y images at specific z values are shown in Figure 3. Specifically, for the variations in Δ at 15 mm of depth, the system is able to identify the perturbation down to 0.017 mm , close to the expected x-y region.
Probably, lower Δ gives rise to signal variations that are too low to be efficiently detected by this system, which could explain the very pale and quite negligible darker shadows in the 2D representations for Δ = 0.01 mm .However, for all of the inclusions, depth is underestimated by about 5 mm, since the darker regions appear at shallower depths with respect to slices where a white We studied three different variables, that are (1) real ∆µ a , (2) radial distance (along y axis) from the injection source, and (3) depth along the z axis.

Effect of the Perturbation
We investigated how the entity of a perturbation affects the tomographic reconstruction by changing the heterogeneity absorption coefficient ∆µ a of 0.04, 0.017, 0.01, and 0.005 mm −1 at the depth of 15 mm, which is a depth where a brain activation is mostly expected.The tomographic reconstructions for x-y images at specific z values are shown in Figure 3. Specifically, for the variations in ∆µ a at 15 mm of depth, the system is able to identify the perturbation down to 0.017 mm −1 , close to the expected x-y region.
circle is depicted, highlighting slices where the inclusion is expected.Furthermore, artifacts always affect the superficial layers, showing the disposition of the optodes banana shapes.(columns) perturbations placed at z = 15 mm.White circles of radius = 10 mm are placed with the center on the actual position of the perturbation volume and are used to assess with more ease the accuracy in the localization of the retrieved perturbations.

Effect of Depth
We tested how deep the perturbation can be detected by the system.In this case, the chosen inclusion is the one with Δ = 0.017 mm , since it reasonably represents a localized brain activation.The depth was varied from 5 to 20 mm, increasing by 5 mm steps.In Figure 4, the tomographic reconstructions are reported.When we simulated the inclusion at increasing depths, we retrieved the perturbation quite well localized both in x and y.Disregarding the superficial artifacts, the most relevant Δ variations are localized correctly in z down to 10 mm, but reconstructed at a lower depth than desired for z = 15 mm.Probably, lower ∆µ a gives rise to signal variations that are too low to be efficiently detected by this system, which could explain the very pale and quite negligible darker shadows in the 2D representations for ∆µ a = 0.01 mm −1 .However, for all of the inclusions, depth is underestimated by about 5 mm, since the darker regions appear at shallower depths with respect to slices where a white circle is depicted, highlighting slices where the inclusion is expected.Furthermore, artifacts always affect the superficial layers, showing the disposition of the optodes banana shapes.

Effect of Depth
We tested how deep the perturbation can be detected by the system.In this case, the chosen inclusion is the one with ∆µ a = 0.017 mm −1 , since it reasonably represents a localized brain activation.The depth was varied from 5 to 20 mm, increasing by 5 mm steps.In Figure 4, the tomographic reconstructions are reported.When we simulated the inclusion at increasing depths, we retrieved the perturbation quite well localized both in x and y.Disregarding the superficial artifacts, the most relevant ∆µ a variations are localized correctly in z down to 10 mm, but reconstructed at a lower depth than desired for z = 15 mm.Nevertheless, this issue is quite common in tomographic reconstructions, since the sensitivity is lower at increasing depth [45,46].Unfortunately, for z = 20 mm, the perturbation seems invisible to the system.

Effect of Lateral Distance
In addition, we evaluated the tomographic capabilities of the system simulating different positions of the inclusion with Δ of 0.017 mm along the y axis, from 15 to 75 mm, with respect to the lateral side of the phantom (y = 0 mm) with steps of 10 mm and at 15 mm of depth (Figure 5).Regarding the distance, we can state that the inclusion is precisely detected until it remains within the superficial octagonal perimeter defined by the detectors (see Figure 2).Instead, as long as it is placed outside this boundary (i.e., a few mm), it is not properly detected anymore, as visible in the y = 75 mm position.Nevertheless, this issue is quite common in tomographic reconstructions, since the sensitivity is lower at increasing depth [45,46].Unfortunately, for z = 20 mm, the perturbation seems invisible to the system.

Effect of Lateral Distance
In addition, we evaluated the tomographic capabilities of the system simulating different positions of the inclusion with ∆µ a of 0.017 mm −1 along the y axis, from 15 to 75 mm, with respect to the lateral side of the phantom (y = 0 mm) with steps of 10 mm and at 15 mm of depth (Figure 5).Regarding the distance, we can state that the inclusion is precisely detected until it remains within the superficial octagonal perimeter defined by the detectors (see Figure 2).Instead, as long as it is placed outside this boundary (i.e., a few mm), it is not properly detected anymore, as visible in the y = 75 mm position.As a general overview about the results, we can affirm that the simulations are effective in showing perturbations of strength down to 0.017 mm −1 , with a satisfactory x-y localization and until the perturbation is under the superficial optodes coverage.In general, the perturbation depth is slightly underestimated (at most by 5 mm).These results can be considered quite encouraging for the proposed approach.In fact, the intersection of the light paths of several source-detector pairs are usually employed to obtain a precise localization, while here, despite that we only employed eight pairs, we can still attribute a specific position to the perturbation.Furthermore, the reduced number of optodes is probably the origin of the superficial artifacts that we have always observed on the reconstructed images at shallow depths.In the first reconstructed layers, the traces of the eight "banana shapes", which are the most probable regions of each photon path, are evident; they could probably result from the downsampling of the volume under the probe.In order to verify whether this effect was seen only in the presence of the inclusion, we tested the system without it and leaving the Poisson noise as the sole perturbation: the artifact is visible also in this case, but very similar superficial effects have been already reported in several works [1,45].As a general overview about the results, we can affirm that the simulations are effective in showing perturbations of strength down to 0.017 mm −1 , with a satisfactory x-y localization and until the perturbation is under the superficial optodes coverage.In general, the perturbation depth is slightly underestimated (at most by 5 mm).These results can be considered quite encouraging for the proposed approach.In fact, the intersection of the light paths of several source-detector pairs are usually employed to obtain a precise localization, while here, despite that we only employed eight pairs, we can still attribute a specific position to the perturbation.Furthermore, the reduced number of optodes is probably the origin of the superficial artifacts that we have always observed on the reconstructed images at shallow depths.In the first reconstructed layers, the traces of the eight "banana shapes", which are the most probable regions of each photon path, are evident; they could probably result from the downsampling of the volume under the probe.In order to verify whether this effect was seen only in the presence of the inclusion, we tested the system without it and leaving the Poisson noise as the sole perturbation: the artifact is visible also in this case, but very similar superficial effects have been already reported in several works [1,45].

Test on Phantoms
We subsequently tested the system on a solid switchable phantom that is capable of mimicking localized absorption perturbations in an otherwise homogeneous medium, as described in Section 2.3.
Figure 6 displays slices of three-dimensional (3D) reconstructions along horizontal x-y planes at different depths z (rows) for different choices of the effective absorption perturbation ∆µ a referred to a volume of 1 cm 3 (columns).We chose the y position of the inhomogeneity yielding the maximum contrast, which roughly corresponds to a location equidistant from each source-detector couple.The inhomogeneity can be reconstructed for all of the tested perturbations (down to ∆µ a = 0.005 mm −1 ) in the expected x-y location (white circle), but with a shallower depth, in particular, for the weaker ∆µ a .This issue is probably due to the strong decrease of the Jacobian coefficients with depth combined to the noise increase of signal coming from deep photons, as also observed in simulations.This combination tends to give more weight to superficial voxels in the reconstruction.

Test on Phantoms
We subsequently tested the system on a solid switchable phantom that is capable of mimicking localized absorption perturbations in an otherwise homogeneous medium, as described in Section 2.3.
Figure 6 displays slices of three-dimensional (3D) reconstructions along horizontal x-y planes at different depths z (rows) for different choices of the effective absorption perturbation ∆ referred to a volume of 1 cm 3 (columns).We chose the y position of the inhomogeneity yielding the maximum contrast, which roughly corresponds to a location equidistant from each source-detector couple.The inhomogeneity can be reconstructed for all of the tested perturbations (down to ∆ = 0.005 mm −1 ) in the expected x-y location (white circle), but with a shallower depth, in particular, for the weaker ∆ .This issue is probably due to the strong decrease of the Jacobian coefficients with depth combined to the noise increase of signal coming from deep photons, as also observed in simulations.This combination tends to give more weight to superficial voxels in the reconstruction.The upper slices confirm the presence of the superficial artifacts also in these phantom measurements, but this does not affect the retrieval of the real perturbation, since the darkest and widest shadows belong to deeper x-y planes.These results are quite consistent with findings on The upper slices confirm the presence of the superficial artifacts also in these phantom measurements, but this does not affect the retrieval of the real perturbation, since the darkest and widest shadows belong to deeper x-y planes.These results are quite consistent with findings on simulations, as reported in Figure 3.We can observe a good agreement for the first two inclusions, while the second two are slightly more visible than in the simulated ones (see Figure 3).This difference is a consequence of the less than ideal performance of a real measurement, whose additional noise and imprecision could make the inclusion appear bigger than it actually is, and the reduced number of optodes, which also determines less accurate identification of the object.Figure 7 (right) displays the reconstructed x-y slices, extracted at z = 16 mm for different ∆µ a values of the inhomogeneity (rows) and positions along the y axis (columns), as depicted in the scheme illustrated on the left panel.
simulations, as reported in Figure 3.We can observe a good agreement for the first two inclusions, while the second two are slightly more visible than in the simulated ones (see Figure 3).This difference is a consequence of the less than ideal performance of a real measurement, whose additional noise and imprecision could make the inclusion appear bigger than it actually is, and the reduced number of optodes, which also determines less accurate identification of the object.Figure 7 (right) displays the reconstructed x-y slices, extracted at z = 16 mm for different ∆ values of the inhomogeneity (rows) and positions along the y axis (columns), as depicted in the scheme illustrated on the left panel.
We observe an acceptable capability of the system to track the object location as far as it is within the detector's circle (with a radius of 30 mm, which stays between 10 < y < 65 mm), in agreement with the outcome of simulations (see Figure 5).Outside that region, the reconstructed perturbation is offplaced in the y direction.For the weakest perturbation, the detection gets worse, and it is basically lost at the edge of the detectors' circle.
For what concerns the repeatability of the measurements, we tested it for two different configurations and we acknowledged that the absorption reconstructions were always in the same positions, with an uncertainty of at maximum 1 voxel (2.5 × 2.5 × 2.5 mm 3 ), probably due to the intrinsic noise and/or regularizations during the reconstructions.

In Vivo Test
Concerning the in vivo measurements, for all of the subjects, Figure 8 reports Δ in the x-y plane at the fixed depth of z = 16.25 mm, for different times, each corresponding to integration over 5 s of measurement from the beginning to the end of the acquisition.For subject 1 and subject 2, we can clearly identify that the activation appears as a dark spot only after the beginning of the task, with about a 5 s delay, and then it disappears, progressively fading up to restoring the initial background condition in about 15-20 s after the end of the task.Unfortunately, the results for subject We observe an acceptable capability of the system to track the object location as far as it is within the detector's circle (with a radius of 30 mm, which stays between 10 < y < 65 mm), in agreement with the outcome of simulations (see Figure 5).Outside that region, the reconstructed perturbation is off-placed in the y direction.For the weakest perturbation, the detection gets worse, and it is basically lost at the edge of the detectors' circle.
For what concerns the repeatability of the measurements, we tested it for two different configurations and we acknowledged that the absorption reconstructions were always in the same positions, with an uncertainty of at maximum 1 voxel (2.5 × 2.5 × 2.5 mm 3 ), probably due to the intrinsic noise and/or regularizations during the reconstructions.

In Vivo Test
Concerning the in vivo measurements, for all of the subjects, Figure 8 reports ∆µ a in the x-y plane at the fixed depth of z = 16.25 mm, for different times, each corresponding to integration over 5 s of measurement from the beginning to the end of the acquisition.For subject 1 and subject 2, we can clearly identify that the activation appears as a dark spot only after the beginning of the task, with about a 5 s delay, and then it disappears, progressively fading up to restoring the initial background condition in about 15-20 s after the end of the task.Unfortunately, the results for subject 3 do not follow the expected time evolution: a darker variation in the absorption coefficient seems to be always present throughout the entire acquisition.The reasons are still under investigation.A 3D tomographic map of the changes in Δ during the activation is reported in Figure 9a for subject 1 and the time course of Δ for the three subjects is shown in Figure 9b.Three different regions of interest (ROIs) were identified to better estimate the spatial and temporal distribution of Δ : ROI_1, a 15 mm side cube centered over the activation area, as identified in Figure 9a, with the upper surface at a depth of 5 mm; ROI_2, similar to ROI_1 but far from the center of the activation; and, ROI_3, corresponding to the whole superficial layer down to a thickness of 5 mm.In Figure 9b, we notice that the time course of Δμa in ROI_3 is rather flat for all of the subjects.For subject 1 and subject 2, we notice a task related decrease of Δμa in ROI_1 during the task period, consistent with A 3D tomographic map of the changes in ∆µ a during the activation is reported in Figure 9a for subject 1 and the time course of ∆µ a for the three subjects is shown in Figure 9b.Three different regions of interest (ROIs) were identified to better estimate the spatial and temporal distribution of ∆µ a : ROI_1, a 15 mm side cube centered over the activation area, as identified in Figure 9a, with the upper surface at a depth of 5 mm; ROI_2, similar to ROI_1 but far from the center of the activation; and, ROI_3, corresponding to the whole superficial layer down to a thickness of 5 mm.In Figure 9b, we notice that the time course of ∆µ a in ROI_3 is rather flat for all of the subjects.For subject 1 and subject 2, we notice a task related decrease of ∆µ a in ROI_1 during the task period, consistent with the expected changes during brain cortex activation.Conversely, smaller and not task-related changes are observed in ROI_2.For subject 3, the results are very noisy in ROI_1 and ROI_2, therefore no task-related activation was detected in this case.

Conclusions
In conclusion, we have presented a compact eight-channel time-domain system for diffuse optical tomography, exploiting the latest advances in photonics technologies and the spatial sampling capability of the time-domain approach.In a first simple configuration, with eight SiPMs that are directly embedded on the probe in a circular arrangement centered on a single injection source at 689 nm, the system is capable to acquire full tomographic data at an acquisition rate of 1 Hz.A basic reconstruction procedure based on analytical solution for the Born approximation is applied to the diffusion equation [41,42], and a purely absorption change was adopted.
A first study on simulations with realistic noise added demonstrated that the proposed scheme can recover buried absorption perturbations with Δ ≥ 0.017 mm −1 , down to a depth z = 15 mm, and as far as the inhomogeneity is beneath the detectors' circle.
These findings were confirmed by experimental measurements on a solid switchable phantom.Tomographic images that are acquired at 1 Hz frame rate can detect the absorption inclusion (equivalent Δ ≥ 0.010 mm −1 over 1 cm 3 ) at a depth of 15 mm as far as the perturbation does not overcome the detectors' ring.
In vivo absorption tomographies on volunteers performing a standard motor-cortex activation

Conclusions
In conclusion, we have presented a compact eight-channel time-domain system for diffuse optical tomography, exploiting the latest advances in photonics technologies and the spatial sampling capability of the time-domain approach.In a first simple configuration, with eight SiPMs that are directly embedded on the probe in a circular arrangement centered on a single injection source at 689 nm, the system is capable to acquire full tomographic data at an acquisition rate of 1 Hz.A basic reconstruction procedure based on analytical solution for the Born approximation is applied to the diffusion equation [41,42], and a purely absorption change was adopted.
A first study on simulations with realistic noise added demonstrated that the proposed scheme can recover buried absorption perturbations with ∆µ a ≥ 0.017 mm −1 , down to a depth z = 15 mm, and as far as the inhomogeneity is beneath the detectors' circle.These findings were confirmed by experimental measurements on a solid switchable phantom.Tomographic images that are acquired at 1 Hz frame rate can detect the absorption inclusion (equivalent ∆µ a ≥ 0.010 mm −1 over 1 cm 3 ) at a depth of 15 mm as far as the perturbation does not overcome the detectors' ring.
In vivo absorption tomographies on volunteers performing a standard motor-cortex activation protocol (20 s finger tapping, 5 repetitions) showed a consistent brain activation response in two subjects out of three.In the positive cases, different task-related temporal evolutions of the absorption coefficient can be tracked in the cortex area as compared to the skull-scalp.
The system still suffers from some limitations.The first limitation is the single source, which implies a poor spatial coverage, leading to artifacts, particularly in the superficial layer, and weak spatial resolution.To cope with those problems, an increased number of source points can be easily inserted by adding an optical switch.However, the distribution of sources and detectors has to be carefully studied to have a sufficient equalization of recorded photon counts for each detector.Another limitation is represented by the single wavelength, since the system is not capable to disentangle deoxy-(HHb) from oxygenated-hemoglobin (HbO 2 ) variations.Also in this case the insertion of a switch coupled to lasers at different wavelengths can be a simple and effective solution to recover the concentration of both HHb and HbO 2 .However, the increase in the number of sources requires a time multiplexing of the DTOF curves, thus forcing it to use a shorter acquisition time for each injection point in order to keep the 1 Hz acquisition rate.This will result in a decrease of the overall signal and also of the signal to noise ratio.For this reason, a trade-off between the number of sources and the minimum signal that is acceptable has to be considered.In vivo results are not fully satisfactory, although not so different from the expected variability in these kinds of tests.Possibly, these results are consistent with the limiting sensitivity of ∆µ a ≥ 0.017 mm −1 at 15 mm depth observed in simulations.
The duplication of wavelengths and the increase in the number of injection sources are feasible, adopting e.g., fiber switches.Further, an increase in total count rate-and consequently in depth sensitivity-is possible with newest TDCs at 20 Mcounts/channel [47], though the technology is still under development for multi-channels devices.Substantial advancements can be anticipated in next few years, yet the actual system-within its limits of sensitivity and spatial resolution-is already an interesting tool for real-time tomography.

Figure 1 .
Figure 1.System architecture scheme (left) and photo (right).The acquired instrument response function of one channel (top left) shows a FWHM of 268 ps, while the tabled values highlights the compactness of the most innovative parts of the system.

Figure 1 .
Figure 1.System architecture scheme (left) and photo (right).The acquired instrument response function of one channel (top left) shows a FWHM of 268 ps, while the tabled values highlights the compactness of the most innovative parts of the system.

Figure 2 .
Figure 2. Reference system: in this representation of the upper surface of the medium, the zero coordinate of all the axes is the top left corner of the medium itself.Optodes scheme: the filled black dots represent the eight detectors, while the central one the source.

Figure 2 .
Figure 2. Reference system: in this representation of the upper surface of the medium, the zero coordinate of all the axes is the top left corner of the medium itself.Optodes scheme: the filled black dots represent the eight detectors, while the central one the source.

Figure 3 .
Figure 3. Tomographic reconstructions showing the effect of the perturbation: x-y images at specific z values (rows) of decreasing Δ(columns) perturbations placed at z = 15 mm.White circles of radius = 10 mm are placed with the center on the actual position of the perturbation volume and are used to assess with more ease the accuracy in the localization of the retrieved perturbations.

Figure 3 .
Figure 3. Tomographic reconstructions showing the effect of the perturbation: x-y images at specific z values (rows) of decreasing ∆µ a (columns) perturbations placed at z = 15 mm.White circles of 10 mm in radius are placed with the center on the actual position of the perturbation volume and are used to assess with more ease the accuracy in the localization of the retrieved perturbations.

Figure 4 .
Figure 4. Effect of inclusion depth on tomographic reconstructions: two dimensional planes taken at specific z (rows) show the perturbations arising from an inclusion of Δ = 0.017 mm placed at four increasing depths (columns), from 5 to 20 mm with steps of 5 mm.Two white circles of radius = 10 mm are placed in the layers just adjacent the depth at which the center of the inclusion should be retrieved.

Figure 4 .
Figure 4. Effect of inclusion depth on tomographic reconstructions: two dimensional planes taken at specific z (rows) show the perturbations arising from an inclusion of ∆µ a = 0.017 mm −1 placed at four increasing depths (columns), from 5 to 20 mm with steps of 5 mm.Two white circles of radius = 10 mm are placed in the layers just adjacent the depth at which the center of the inclusion should be retrieved.

Figure 5 .
Figure 5.Effect of lateral position of the perturbation on tomographic reconstructions: x-y reconstructed images at increasing z values (rows) of an inclusion of Δ = 0.017 mm placed at 15 mm of depth and at different positions along the y axis (columns).The two dashed horizontal lines in each slice correspond to the edges of the superficial coverage of the probe.The white circles indicate the first neighbor planes of the inclusion and trace its path along y axis.

Figure 5 .
Figure 5.Effect of lateral position of the perturbation on tomographic reconstructions: x-y reconstructed images at increasing z values (rows) of an inclusion of ∆µ a = 0.017 mm −1 placed at 15 mm of depth and at different positions along the y axis (columns).The two dashed horizontal lines in each slice correspond to the edges of the superficial coverage of the probe.The white circles indicate the first neighbor planes of the inclusion and trace its path along y axis.

Figure 6 .
Figure 6.Slices of three-dimensional (3D) reconstructions along horizontal x-y planes at different depths z (rows) for different choices of the effective absorption perturbation ∆ referred to a volume of 1 cm 3 (columns).

Figure 6 .
Figure 6.Slices of three-dimensional (3D) reconstructions along horizontal x-y planes at different depths z (rows) for different choices of the effective absorption perturbation ∆µ a referred to a volume of 1 cm 3 (columns).

Figure 7 .
Figure 7. (Left) Scheme of the measurement condition: the inclusion, represented by the black dot, is placed in different positions along y axis as shown by the bold dotted line, while the simple dotted lines qualitatively specify the superficial disposition of the probe.(Right) Tomographies showing horizontal x-y slices taken at z = 16 mm for different ∆ of the inhomogeneity (rows) and y positions (columns).The white circle identifies the expected "true" location, while the black dashed lines correspond to the boundary of the detectors' circle.

Figure 7 .
Figure 7. (Left) Scheme of the measurement condition: the inclusion, represented by the black dot, is placed in different positions along y axis as shown by the bold dotted line, while the simple dotted lines qualitatively specify the superficial disposition of the probe.(Right) Tomographies showing horizontal x-y slices taken at z = 16 mm for different ∆µ a of the inhomogeneity (rows) and y positions (columns).The white circle identifies the expected "true" location, while the black dashed lines correspond to the boundary of the detectors' circle.
follow the expected time evolution: a darker variation in the absorption coefficient seems to be always present throughout the entire acquisition.The reasons are still under investigation.

Figure 8 .
Figure 8.In vivo tomographies: x-y planes reconstructions at z = 16.25 mm at different times from t = 1 s to t = 60 s averaged in ranges of 5 s (rows) for three subjects (columns) during a finger tapping motor task.The grey rectangle marks the task.In this case the volume used is 100 × 100 × 30 mm 3 .

Figure 8 .
Figure 8.In vivo tomographies: x-y planes reconstructions at z = 16.25 mm at different times from t = 1 s to t = 60 s averaged in ranges of 5 s (rows) for three subjects (columns) during a finger tapping motor task.The grey rectangle marks the task.In this case the volume used is 100 × 100 × 30 mm 3 .

Figure 9 .
Figure 9. (a) 3D tomographic reconstruction of Δ for subject one, averaged only during the activation time course and indicating three selected regions: ROI_1 with solid lines and focused on the activation, ROI_2 with dashed lines and positioned far from the activation but at the same depth of ROI_1, and ROI_3 with dotted lines and considering the superficial 5 mm deep volume; (b) temporal time courses, averaged each 5 s, of the absorption properties of the three different regions for all three subjects.For each region, each Δ time value is calculated as the spatial average over the respective region.

Figure 9 .
Figure 9. (a) 3D tomographic reconstruction of ∆µ a for subject one, averaged only during the activation time course and indicating three selected regions: ROI_1 with solid lines and focused on the activation, ROI_2 with dashed lines and positioned far from the activation but at the same depth of ROI_1, and ROI_3 with dotted lines and considering the superficial 5 mm deep volume; (b) temporal time courses, averaged each 5 s, of the absorption properties of the three different regions for all three subjects.For each region, each ∆µ a time value is calculated as the spatial average over the respective region.

Table 1 .
Dimensions of the main building blocks composing the proposed system.