Effects of time-gated detection in diffuse optical imaging at short source-detector separation

The adoption of a short source-detector distance, combined with a time-resolved acquisition, can be advantageous in diffuse optical imaging due to the stricter spatial localization of the probing photons, provided that the strong burst of early photons is suppressed using a time-gated detection scheme. We propose a model for predicting the effect of the time-gated measurement system using a time-variant operator built on the system response acquired at different gate delays. The discrete representation of the system operator, termed Spread Matrix, can be analyzed to identify the bottlenecks of the detection system with respect to the physical problem under study. Measurements performed on tissue phantoms, using a time-gated single-photon avalanche diode and an interfiber distance of 2 mm, demonstrate that inhomogeneities down to 3 cm can be detected only if the decay constant of the detector is lower than 100 ps, while the transient opening of the gate has a less critical impact.


Introduction
The study of light propagation into highly diffusive media like biological tissues (Diffuse Optics) is highly appealing due to the possibility to explore the medium non-invasively, deep beneath the surface and to recover information both on absorption (related to chemical composition) and on scattering (related to microstructure). Due to the relatively low absorption of most biological tissues in the 600-1100 nm wavelength range and to the inherent non-invasiveness of light at low power, new photonics tools have been devised for clinical diagnostics, such as optical mammography, functional imaging of tissue oxygenation and haemodynamics, diagnosis of osteoarticular diseases [1][2][3][4]. Yet, in the same spectral range most tissues are highly diffusive, causing a strong attenuation of the remitted light and blurring effects impairing spatial resolution. These effects result in a limited possibility of exploring biological structures deeply embedded in tissues.
A common way to increase sensitivity to deeper structures is to increase the distance ρ between the light injection (source) and collection (detector) points. Yet, what causes the deeper penetration is not the source-detector distance alone, but rather the photon arrival time, that is monotonically related to the mean visited depth [5][6][7]. Since for a larger ρ early photons are more severely attenuated by the effect of propagation, a larger mean photon transit time yields deeper penetration. Thus, the same depth can be reached even at short or null ρ provided that photons propagated for a sufficiently long time t are selected [8]. The physics of photon migration observed with null or small ρ is quite fascinating and yet largely unexplored. Ultimately, the depth is encoded only in the propagation time, while the stricter localization of photons detected at a small ρ leads to better spatial resolution, contrast and higher number of re-emitted photons [8].
Since at short ρ the number of photons re-emitted after a short time t (early photons) is huge, these concepts hold true if an efficient time-gating mechanism is applied, so to suppress such initial burst of photons that would otherwise wash out the few most meaningful photons arrived at a larger t. Such a scheme was indeed demonstrated using a Single-Photon Avalanche Diode (SPAD) operated with an ultrafast gate [9], achieving a record dynamic range of 8 decades [10,11], doubling the intrinsic dynamic range of the non gated detection chain. The chance to perform measurements at short ρ for in depth investigation of biological tissues has opened the way to new schemes, such as functional imaging of brain activation using a single optode [8], non-contact diffuse optical imaging [12,13], single-fiber interstitial spectroscopy [14], diffuse optical tomography using a compact probe [15]. Furthermore, the possibility to use any ρ for diffuse optical imaging makes it possible to devise dense and distributed source-detector coverage without the limitation of enforcing large source-detector distances [16,17].
The SPAD can be used for demonstrating these new concepts since it can be switched on and off in a short time due to the low driving voltage and small capacitance by applying a fast rising bias gate. Furthermore, the device is not damaged when exposed to high photon fluxes when in off state due to the fast electron-hole recombination rate, that prevents charge accumulation. However, the migration of charges beneath the depleted region and their subsequent diffusion cause a tail in temporal response with a characteristic time decay constant that causes small leakage of early photons even at much longer times (from ns to few µs [18]).
The ultrafast gated SPAD technology is at an early stage, with great room for further advancements, mainly towards massive parallelization via e.g. arrays of gated SPADs and more efficient gating capabilities. Still, modeling is urgently needed to understand the key limiting factors of the technology so to address the development phase, but also to better understand the true possibilities of this novel scheme. Yet, modeling a real time-gated diffuse reflectance measurement at small ρ is not trivial, for two main aspects. First, the measurement at a late gate can be severely contaminated by the burst of early photons depending on the specific properties of the non-ideal time-gated system and how it deals with out-of-gate photons. Second, the measurement system is inherently time-variant being dependent on the relative delay of photon arrival time with respect to the opening of the gate. Thus, although the linearity assumption still holds, the convolution theorem can not be applied and a diverse model must be accounted for.
In this paper we propose the modeling of a time-gated nonideal detection system applied to the physics of small ρ diffusive optical imaging. As far as possible, a general approach is adopted, valid for any time-gated system, yet with a practical application to time-gated SPAD detection. In the following, we will first recall the theoretical basis of time-variant system operators. Then we will discuss the matrix representation of the system operatorhere called Spread Matrix-and its physical meaning as temporal spread of collected photons caused by the non-ideal detection system. These concepts will be applied to the practical case of time-gated SPAD detection. Using phantom measurements both on homogeneous and inhomogeneous media, this approach will be first validated and then effectively used to identify the most crucial bottlenecks of the measurement system when applied to the detection of small optical inhomogeneities (e.g. brain activation areas, breast tumors) within a diffusive medium.

Theory
Let us consider a time-resolved diffuse optical measurement. A short laser pulse is injected at the surface of a highly scattering medium and photons are re-emitted at a given point r with a temporal distribution N(t). The measurement system can be described by a linear operator Ŝ which transforms the real photon distribution N(t) into the measured one =m t N t S ( ) ( ). The system is time-invariant if the effect of Ŝ does not depend on the arrival time t of photons and can be described by the instrumental transfer function S(t) that is the response of the system to a delta-like photon distribution. Then, the convolution theorem holds true and the measurements can be derived as the convolution of the system response with the photon time-distribution: ). This is the case-for instance-of a Time-Correlated Single-Photon Counting (TCSPC) instrument [19] operated under standard conditions (i.e. with no gating). Conversely, the system can be time-variant, meaning that the system behavior is not the same for any arrival time t. This always happens for a time-resolved system operated with a temporal gate. In the ideal case, the system is completely blind (m = 0) if the photons arrive while the gate is closed (off), while it is active for photons impinging when the gate is open (on). In the real case, the opening of the gate is not instantaneous, with possible oscillations and the suppression factor is not infinite, leading to changing responses of the detector under different operating conditions along the gate. In the most general case the response of the system S(t, τ) under a delta-like temporal illumination depends on the temporal position τ of the opening gate.
The convolution theorem can not be applied since the system is time-variant. The measurements can be derived using a pseudo-convolution, where the system term is timevarying [20]: The index τ on the system operator ̂τ S signifies that the actual system action and thus the measurement m depends on the temporal position τ of the opening gate. Equation (1) can be interpreted as the superposition of the effects of all photon packets N(t′)dt′. Each of them contributes at t with a weight given by the response function S(t − t′, τ − t′) produced by a burst of photons at a time-delay t − t′. The response function must correspond to the actual delay τ − t′ of the opening gate seen by photons impinging at t′. In the following, t′ will be addressed as the real time, i.e. the actual time when the photon escapes from the medium surface, while t is the measured time, i.e. the time-channel when the photon is classified by the detection electronics. Thus, the measurement system, described by the operator ̂τ S , transforms the photon distribution from the real time to the measured time.
Moving now from the continuous to the discrete case, both the measurement m(t) and the photon distribution N(t) can be expressed as vectors m = {m i } and N = {N i } where the index i spans all the M temporal channels of the acquisition system, each of them with width Δt. Thus, equation (1) can be recast into matrix form: where the matrix S τ can be derived using the system responses obtained for all possible time gate delays and has the following form: where S i, j = S (iΔt, jΔt) represents the response function evaluated at time iΔt acquired for a gate delay set at time jΔt. We have neglected integration of the functions over the time interval Δt, assuming it is sufficiently small (typically <50 ps). It is worth noting that, to build S τ , the system response must be recorded both for positive and negative delays of the time-gate, since actual photons can arrive both before and after the opening of the gate. Normally, a large set of delays is acquired [10], so that the S τ matrix can be constructed for any delay τ of interest.
To better understand the physical meaning of the S τ matrix, it is useful to look at the measurement at a given time channel i: Thus, the i-th row of the S τ matrix contains the weights by which the elements of N contributes to produce the i-th measurement channel. If the system were ideal and with no time-gating, than only elements m i, i are non-null and the S τ matrix is diagonal. Thus, the S τ matrix represents the spread of photons among measurement channels due to the effect of the non-ideal system. For this reason, we will refer in the following to S τ as Spread Matrix.

Time-gated SPAD
A couple of SPADs were tested. Both the detectors were developed at Politecnico di Milano ( [21]). They differ for the production run, active area diameter and for the device cross section, in fact depletion region thickness and dopant diffusions are not the same. Differences in structure give differences in performances: timing response of the two detectors is different, in particular for the decay time constant of the response to a pulsed laser [22]. The first detector (SPAD1) has a decay time of 85 ps, an active area of 100 µm and a breakdown voltage of 22 V, while the second detector (SPAD2) has a decay time of 240 ps, an active area of 200 µm and a breakdown voltage of 36 V.

Diffuse spectroscopy system
The system set-up is shown in figure 1. The system is based on a supercontinuum fiber laser (SC450, Fianium Ltd., UK) as laser source. It generates white-light picosecond pulses at a repetition rate of 40 MHz. Wavelength selection is performed by means of an Acousto-Optic Tunable Filter (Neos Technologies, FL, USA), not shown in figure 1. Laser pulses at 750 nm, with a spectral width of about 5 nm, were selected. Light is delivered and collected from the sample by means of 1 mm core step-index fibers. The collected light is focused on the detector by means of an aspheric duplet lens system with a demagnifying factor of 1.5. A small fraction (5%) of the source laser light is split-off and delivered to a photodiode (PD) in order to generate a synchronization signal, to monitor slow laser power fluctuations occurring during data acquisition and compensate them during data analysis. A constant fraction discriminator (CFD) is employed in order to reduce time jitter in the evaluation of the laser pulse time position after the photodiode. A stack of 3 metallic variable attenuators (VA) with an achievable maximum attenuation of 120 dB is used in order to keep the total count rate within single-photon counting statistics to avoid pile-up distortion [19]. This attenuation stage is placed where laser light is coupled in the injecting fiber. A fully calibration of the stack of attenuators was performed in order to normalize all the acquisitions to the same equivalent input energy during the data analysis. A second VA is placed in front of the PD in order to adjust the signal level.
The set-up is completed by a TCSPC board (SPC130, Becker and Hickl, Germany) for data acquisition and by an ultra-fast pulse generator for the time gating of the detector. Both TCSPC board and time gating electronics were synchronized by means of a signal derived directly from the PD monitoring the laser source. The detector can be turned on in less than 200 ps. The gate window (temporal window in which the detector is in the on-state) can be adjusted in time with a precision of tens of picoseconds thanks to an homemade passive delayer based on RF relays and transmission lines of different lengths. The delayer is controlled via USB by a home made software written in C language. The gate window width is adjustable from less than 1 ns up to 10 ns in 10 ps steps. When a photon is detected the avalanche signal from SPAD is sensed by means of a homemade high speed comparator based on fast ECL (Emitter-Coupled Logic) electronics and sent to the CFD input of the TCSPC board for the arrival time estimation. In the meanwhile the high speed comparator provides a reset pulse to SPAD in order to restore the operative condition after the photon detection.

Phantoms
Calibrated phantom mimicking optical parameters of homogeneous diffusive media were obtained using a water solution of Intralipid and Indian Ink (Higgins) [23]. The liquid phantom was poured in a black PVC tank (120 × 140 × 70 mm) with holes in the middle for hosting the launching and collecting fibers at an interfiber distance ρ = 2 mm. The inhomogeneous phantom [24] was obtained placing in a liquid phantom a black PVC cylinder (volume = 100 mm 3 ), equivalent [25] to a 1 cm 3 inhomogeneity with an absorption change Δμ a = 0.2 cm −1 . The PVC cylinder is placed along the fiber axis, equidistant from the two fibers by means of a rigid metallic wire (0.5 mm music wire) painted white and moving it at increasing depths z beneath the fiber tips plane. The absorption (μ a ) and the reduced scattering coefficient (μ′ s ) at 750 nm were 0.05 and 10 cm −1 , for the homogeneous phantom and 0.1 and 10 cm −1 , for the inhomogeneous phantom.

Effects of time gating
An example of the system response function obtained for different delays of the opening gate is shown in figure 2. The injection and collection fibers were placed in contact, while the enabling gate was shifted in time via a programmable delayer, so to progressively cut out the laser pulse, that is assumed to be of negligible width with respect to the timejitter of the electronics. The count rate was always kept fixed by rotating the variable attenuator. For the first delays the laser pulse impinges onto the detector while the gate is fully open, while upon increasing the delay, the laser pulse gets out of the open gate and only the tail of the SPAD response is detected. The injected power was increased from ≈10 nW for the first delay up to ≈1 mW for the last delay. The gate is far from being ideal since the SPAD response is not just sliced, but severely altered by the transient sub-optimal biasing of the gate and the ringings of the gating pulse. Just for these measurements, the data shown were obtained with a commercial fast pulse generator (AWN-W3, Avtech Electrosystems, USA) used in a previous study [9], where the non ideal and time-variant effects of the gate are markedly shown. With the dedicated fast pulser developed specifically for the SPAD in use this effect is less pronounced, but still present. It's worth nothing that the effect of the gate can not be simply described multiplying the SPAD response under non-gated operation (free running) for the gate pulse profile, since the SPAD behavior is dramatically altered for intermediate values of the bias voltage with a large increase of temporal broadening.

Constructing the spread matrix
The assembling of the Spread Matrix is described with reference to figure 3. The instrument is prepared so as to detect the system response without the presence of the scattering medium, that is for an ideal delta-like response. For fibre-based systems this is obtained by directly facing the injection and collection fibres with the insertion of a Teflon foil in-between to excite all propagation modes in the collecting fibre, thus simulating the angular distribution of a diffusive medium. A wide set of responses are acquired spanning all possible arrival times (t′) of the photon with respect to the opening of the gate (τ). In practice, this is accomplished more easily by moving, instead, the gate temporal position τ for a fixed photon arrival time t. These experimental curves are properly shifted to reproduce the response to a variable photon arrival time and stored in the Spread Matrix. Light injection power is always adjusted so as to keep within optimal single-photon statistics, thus data stored in the Spread  Matrix data are first multiplied by the source attenuation. The finesse of the delay steps must match exactly the required resolution on the measured time axis. Normally a step-size of 25 ps will do. Since the discretization of the delayer control and of the TCSPC acquisition board are different, a resampling of the curves using an interpolation is needed. In principle, the Spread Matrix should be constructed for the specific value of the gate delay τ involved in the measurement. Nevertheless, it is sufficient to acquire a general matrix by spanning τ and then extract a submatrix out of the general one by choosing properly the first row and column depending on the desired value of τ required by the specific application.

Using the spread matrix in a forward solver
To show the aptness of the proposed approach to predict real measurements, time-gated diffuse reflectance curves were acquired on the homogeneous phantom with an interfiber distance of 2 mm, using SPAD1 (left pane) or SPAD2 (right pane). Figure 4 shows the experimental data (blue dashed line), compared to the theoretical prediction (red solid line) obtained implementing equation (2) with the results of a Monte Carlo simulation of the photon distribution N(t) for the same homogeneous medium. Generally, the model succeeds in describing the experimental data and in particular the deformation of the time-resolved curve with the marked increase of the peak. There are some discrepancies between data and model that could be ascribed to (i) the jitter of the raising edge of the gate, (ii) the not easy calibration of the variable attenuators over 8 decades, (iii) some inaccuracies in the temporal calibration of the passive delayer. Nonetheless, the shape of the gated measurements are well tracked.

Understanding the spread matrix
So far, the Spread Matrix was used to simulate the deformation produced on a given photon distribution. In the following, we will exploit the Spread Matrix as a powerful tool to understand the key features of a time-gated system. With reference to equation (1), we recall that t′ represents the real time, i.e. the actual arrival time of the photon at the detector. Conversely, t represents the measured time, i.e. the temporal channel where a given photon is classified by the TCSPC electronics. Thus, with reference to equation (3), a certain row of the Spread Matrix contains the distribution of real time t' that corresponds to a given measured time t, that is it tells when the photons recorded at a given temporal channel have actually reached the detector. This concept is better elucidated in figure 5 that displays a few rows of the Spread Matrix of the two systems described above. Each curve corresponds to a different measured time of the detection electronics, taken at intervals of 100 ps, starting from the measurement point at the early onset of the gate (noisy curve on the left) up to the measurement points taken when the gate is stably open (collapsed curves on the right). The x-axis is the time-shift, that is the lead or lag caused by the detection system. Thus, an ideal time-invariant system with an instantaneous response would appear in figure 5 as a delta-distribution centered at t shift = 0.
When the gate is stably opened, the system gets time-invariant and reaches the optimized response, thus all curves collapse to the same shape that is very similar to the time-reversal of the System response of the SPAD operated in free running, i.e. without gate, but for a much larger dynamic range. Conversely, upon approaching the gate edge, photons get more and more dispersed. The Full Width at Half Maximum (FWHM) increases from 50 ps up to 150 ps (SPAD1) and from 75 ps up to 300 ps (SPAD2). Furthermore, there is a negative shift in the peak position, meaning that the measured signal originates from earlier real times (negative time shift) due to the slower response of the detector. Moving at even earlier measurement times the curves drop to zero, implying that the system does not classify any photons before the gate is opened due to the low bias voltage preventing the avalanche multiplication process. Yet, this must not be confused with the suppression factor of early photons. Indeed, no photons are classified at an early measurement time, but even for a relatively long measurement time there is a non-negligible abundance of early photons in the case of SPAD2 [18].
In summary, the key features that can be seen from the Spread Matrix signature of the system under study are (i) a long decay slope (around 80 ps for SPAD1 and 250 ps for SPAD2) causing a delayed contamination of early real time even at a late measurement time and (ii) a local spreading of photons around the peak position across the rising edge of the gate (still <300 ps in the worst case), (iii) an increase in the background noise at high laser power inpinging onto the detector due to the memory effect [18].

Applying the spread matrix to diffuse reflectance imaging
The Spread Matrix can be used to understand the bottlenecks of a systems with respect to a specific application and to derive useful indications for further development. In the following we present a powerful methodology to understand the effect of the time-gated system on the problem of imaging a localized inhomogeneity within a diffusive medium, as in the case of functional imaging of brain activity. Figure 6 shows the simulated temporal composition (real time) of the signal detected at different measured times for ρ = 2 mm using SPAD1 (left column) and SPAD2 (right column). These plots were derived using Monte Carlo simulations for a homogeneous medium with μ a = 0.1 cm −1 and μ′ s = 10 cm −1 . According to equation (4), each curve corresponds to a given measured time and represents the product of the Spread Matrix row with the real photon distribution N, plotted against the real time. Both a logarithmic (top row) and a linear (bottom row) scale are represented, to highlight the suppression of early photons and the effective temporal composition of the measured signal, respectively.
A first clear indication to be drawn is that the SPAD1 System succeeds in extracting information from late photons, while the SPAD2 System-with a longer decay constant-carries information only from early photons and thus is less effective for in-depth measurements. Secondly, the temporal spreading in the time-variant region, caused by the transient opening of the gate, is reflected in a broadening of the photon temporal distribution. Nonetheless, this is a local effect, with a well confined spread of information around a given photon arrival time, not exceeding 200 ps (FWHM) as for figure 6(a). This broadening is still lower than the width of the integrating gate that is often applied at analysis time to increase the SNR of the data, thus it can be neglected. We underline here that the largest signal derives from the time-variant region (the main peak in figure 6(a)), as expected since the diffuse reflectance is almost exponentially decreasing with time. Thus, the possibility to effectively use also the most abundant time-variant fraction has important implications for in vivo applications.
An experimental validation of these inferences is presented in figure 7, that shows the contrast produced by a local optical perturbation set beneath the fibers at a variable depth z within a homogeneous medium, as described in section 3. The contrast is defined as: where m i 0 and m i (z) are the measured signal for the homogeneous and inhomogeneous medium, respectively. The measured signal is integrated over a time-window starting at time t w with width Δt w .
Two contrast plots as a function of z are displayed for various time windows with different starting time t w and the same width Δt w = 0.5 ns, using SPAD1 (left pane) and SPAD2 (right pane). No effective contrast can be detected using the SPAD2 because of the longer decay constant, in agreement The mild effect caused by the local temporal spreading during the gate opening permits to use also the time-variant section of the data. This is confirmed by experimental data    figure 8, showing the contrast (top row) and counts (bottom row) using SPAD1 at different values of the time windows t w (columns) for different delays τ (legend). Basically, as far as there are enough counts, the choice of the delay does not affect the contrast, that is dependent solely on the time window. As a consequence, the delay that yields the maximum counts for the desired t w is to be preferred because of best signal-to-noise ratio. Results reported in figure 7(a) were obtained applying this criterion, which corresponds to select the first time-variant portion of the temporal curve. A consequent implication is that we can devise a practical strategy for the choice of the optimal delay when using a gated system in combination with a TCSPC acquisition to reach maximum depth. The delay is increased up to the point when-given the maximum usable injected power-the count rate is at the upper limits for either single-photon counting statistics or board electronics. Increasing any further the delay will result in loss of useful photons without increasing late counts. Conversely, there is no need to decrease the delay so as to shift the time-variant portion of the curve towards early times, at the expenses of reducing the injected light, since the timevariant portion is just as effective in retrieving best contrast as the time-invariant one. More care must be taken if the exact time-profile of reflectance is to be recorder for e.g. model-fitting. In that case, the delay must be reduced at the expenses of signal-to-noise ratio in order to retain the portion of interest in the photon time-distribution within the time-invariant region.

Conclusions
In conclusion, we have proposed the modeling of time-gated detection in time-resolved diffuse optical imaging. Due to the time-variant nature of the detection system, the convolution theorem is not applicable and the analytical expression of the system operator must include the information on the system response for multiple gate delays. Further, the discrete expression of the system operator, termed Spread Matrix, is presented that transforms the real time photon distribution into the measured time distribution.
This approach was applied to the study of time-gated SPAD detection for diffuse reflectance measurements at a short source-detector distance. The analysis of the Spread Matrix revealed two separate effects, that are the spread of temporal information due to the transient opening of the gate and the delayed detection of photons due to the diffusion decay constant of SPAD. While the former has little impact on diffuse optical imaging since leads to a localized spread of photons, the latter can completely vanish the information from deep structures if the detector decay constant is not fast enough to damp the contamination of early photons. These concepts were confirmed by phantom measurements on heterogeneous media showing that diffuse optical imaging in reflectance geometry with an interfiber distance of 2 mm can be achieved only if the SPAD decay constant is <100 ps.
The proposed modeling and analysis can be applied also to other set-ups and problems. A first area is the use of fastgated Intensified CCD camera for diffuse optical imaging [26,27] and tomography [16,28]. A further application is optical biopsy using a single interstitial fiber for both injection and collection [14]. Also in this case early photon suppression is key to provide insight well beyond the fiber tip. Another interesting field is time-resolved fluorescence spectroscopy for molecular imaging [29][30][31][32][33], for fluorescence lifetime imaging [34,35] or Forster Resonance Energy Transfer imaging [36].  Gated detection is growing interest also in Raman spectroscopy through diffusive media [37,38], where gated acquisition is a plus to gain information beneath the surface. The proposed method can be extended to all these fields providing better insight on the system design and on the physical information that can be gained by it.