Fluorescence time-resolved imaging system embedded in an ultrasound prostate probe.

Ultrasound imaging (US) of the prostate has a low specificity to distinguish tumors from the surrounding tissues. This limitation leads to systematic biopsies. Fluorescent diffuse optical imaging may represent an innovative approach to guide biopsies to tumors marked with high specificity contrast agents and therefore enable an early detection of prostate cancer. This article describes a time-resolved optical system embedded in a transrectal US probe, as well as the fluorescence reconstruction method and its performance. Optical measurements were performed using a pulsed laser, optical fibers and a time-resolved detection system. A novel fast reconstruction method was derived and used to locate a 45 µL ICG fluorescent inclusion at a concentration of 10 µM, in a liquid prostate phantom. Very high location accuracy (0.15 cm) was achieved after reconstruction, for different positions of the inclusion, in the three directions of space. The repeatability, tested with ten sequential measurements, was of the same order of magnitude. Influence of the input parameters (optical properties and lifetime) is presented. These results confirm the feasibility of using optical imaging for prostate guided biopsies.

1. Introduction

Diffuse optical tomography is an emerging and promising technique for biological tissue screening [1]. The main advantage of this optical method compared to others in molecular imaging is that it remains non invasive and non-ionizing.

These assets are particularly beneficial for prostate imaging applications. Indeed the current diagnosis is based on prostate specific antigen (PSA) dosage and digital rectal examination [2,3]. If one of these tests is abnormal, an ultrasound (US) guided biopsy may be performed. Even if it is considered as the gold standard for now, this technique presents several drawbacks. In particular, the lack of sensitivity and specificity of sonography alone impacts significantly the localization of tumor targets. Biopsies are thus performed randomly yielding a global false negative rate of 30-40% [4,5], and implying repeated biopsies. In addition to patient discomfort and potential complications [6], these biopsies have an even higher risk of false negative diagnosis [7,8].

To improve the specificity on prostate imaging by guiding biopsies to tumors, it has been proposed to combine near infrared diffuse optical tomography with ultrasound [9]. The assumption is that tumors are more vascularised than the surrounding tissue. This technique was successfully demonstrated in vivo on canine prostate to discriminate tumor nodules and normal tissue [10]. The authors used a sagittal fiber array coupled to a continuous wave source.

Here, we propose to use fluorescence optical tomography to further increase the contrast between tumors and healthy tissue. We also chose to perform time-resolved acquisitions to accommodate the prostate reflectance geometry. Fluorescent time-resolved experiments have been performed in reflectance mode on other organs such as the brain [11], or on tissue-like phantoms [12]. Gao et al. [13] presented a full tomographic device in transmittance mode on a small animal sized phantom.

Optical measurements on prostate are challenging because of the high absorption coefficient (µ = 0.3-0.4 cm−1 at 786 nm [14]) and the small contact surface between the probe and the tissue. Our group had previously presented a prostate dedicated device combining fluorescence tomography and ultrasound technologies [15]. Based on this device, we developed a novel faster detection system based on time-correlated single photon counting (TCSPC). We also improved our reconstruction algorithm and characterized the device resolution.

In Section 2, our time-resolved experimental setup is described. We also present the forward model we employed within the diffusion approximation to obtain the Green’s functions. The derivation of our original method is explained subsequently and preprocessing of the experimental data is detailed. In Section 3, we show the ability of the system to localize a single inclusion in a prostate mimicking phantom.

2. Materials and methods

2.1. Experimental setup

The main parts of the experimental system (Fig. 1 ) are described in detail in reference [15]. Here a similar laser was operated but the detection part has been changed to speed up data acquisition [16]. The gated CDD camera has been replaced by four photomultiplicators (PMTs) coupled to a TCSPC system.

Fig. 1

Illustration of the experimental setup.

The laser source was a pulsed Ti:sapphire laser (Chameleon, Coherent, Santa Clara, USA) tuned at 775 nm, at 80 MHz repetition rate. The light was filtered by a band pass excitation filter (FB770, Thorlabs) then was injected to an input excitation fiber mounted on a displacement plate in order to rout sequentially the light through a six fiber bundle and the six sources on the probe. It allowed the medium to be probed by light from six different locations. The light reemitted by the probed medium was then collected by four plastic detection fibers. Each detection fiber was coupled to its own photomultiplier (PMC-100-20, Becker&Hickl, Germany). A set of emission filters (Notch NF01-785U-25, Razoredge LP02-785RS-25, Semrock and RG850, Itos) was placed in front of each PM head to select fluorescent light. A four-channel TCSPC module (SPC-134, Becker&Hickl, Germany) analyzed the signals from the four PMTs. The synchronization signal was derived from the main beam after the laser exit and was measured with a photodiode. An adequate cable length ensured the correct synchronization delay regarding our system.

Figure 2 shows the setup of the fibers at the outer end of the US probe and its orientation. The fibers are housed between the plastic shell of the probe and the transducer core.

Fig. 2

(a) Optical fibers configuration at the end of the probe (b) Position of the probe and axis representation (the basis is clockwise).

We designed a liquid phantom mimicking the prostate optical properties. It was made of Intralipid at a concentration of 1% as the scattering medium and black India ink (Dalbe, France) as the absorber. The relative concentrations were adjusted to reach levels of = 0.3 cm−1, = 12 cm−1 at 775 nm which are typical prostate optical parameters [14]. The optical properties were controlled by a standard TCSPC techniques [17,18]. A cylindrical plastic tank (diameter 11 cm; height 8 cm) was filled with the phantom solution. The probe was immersed 1 cm below the surface to simulate the pressure of the probe on the rectum wall inside the prostate.

We simulated a tumor targeted by a biomarker by using a 45 µL volume of Indocyanine Green (ICG) encapsulated in nanodroplets [19] (peak excitation and emission wavelength at 796 and 820 nm respectively) placed in a glass tube (inner diameter: 3 mm, outer diameter: 4.4 mm and length: 5 mm). The glass tube was positioned perpendicularly to the sources-detectors line, along the x axis, its outer extremity hold on the x<0 side. The lifetime of the fluorophore (τ = 0.55ns) was also measured by standard time-resolved techniques.

2.2. Forward model

Diffusion equation

Light propagation in the diffusive medium is described by the diffusion equation [20]:

where c is the light velocity, is the diffusion coefficient, is the absorption coefficient, is the Dirac distribution and is the photon density resulting from the Dirac source excitation (Green’s function).

For media with arbitrary boundaries or arbitrary optical properties, numerical calculation must be performed to obtain G. However, solving directly Eq. (1) implies huge memory and computation requirements since one must solve N times a steady state diffusion equation (N being the number of time steps of the calculations). In order to reduce memory and computation time, we used the moments of G. Here, the n order moment of a function f(t) is defined as:

As shown in [21], moments of G are easily reachable since they follow the equation:

where for n = 0, and = 0 otherwise, and = 0 for n<0. Equation (3) shows that is the solution of steady state diffusion equation with a source term including . After meshing our volume, the Green’s function list,…,can therefore be numerically determined by inverting n + 1 times an equation LX = S where operator is an operator whose expression can be found by using the finite element method (FEM) or the finite volume method (FVM). Partial current boundary conditions were applied to the probe boundary and the medium’s edges. Finally, we obtain a table with the Green functions values on each point of the mesh.

Excitation and emission measurements

Two types of time-domain measurements were performed to optically probe a medium of interest: excitation measurements, i.e. measurements at the wavelength of the laser source, and fluorescence or emission measurements, i.e. measurements at a wavelength above the laser source obtained by using a filter. Furthermore, multiple sources (indexed by s) and detectors (indexed by d) were used to provide a set of measurements allowing tomographic reconstruction. In addition, optical properties were assumed to be the same at the excitation and emission wavelength.

Excitation measurements (named ) are described by three physical processes: the time shape of the source (the laser source and propagation inside the source fiber), the propagation inside the diffusive medium (named as the model-predicted data) and the detection time response (detector fiber and detector). Some noise could be added to this model but we will neglect it here, considering only high signal to noise measurements. As these processes are independent, emission measurements are expressed as below,

where IRF is the instrument response function (IRF, convolution of the time response of the source and detection parts) and * stands for the time convolution operator. would be the excitation response of an ideal instrument, that is to say for a Dirac instrument response function.

Emission measurements may be written as the convolution product of five processes: the source emission, the propagation from the source to a fluorescence site, the fluorescence exponential decay, the propagation at the fluorescence wavelength from the fluorescence site to the detection point and the detection. As these processes are independent, emission measurements follow this equation:

where is the fluorescence time-response at location , is the Green’s function for a Dirac located at the detector position. IRF is the same function as in Eq. (4) as we assumed that detection responses are equivalent at the excitation and emission wavelength. would be the emission response of an ideal instrument. Here, we suppose that only one kind of fluorophore with lifetime τ is present, therefore F is modeled as where is related to the concentration of fluorophore at position and to the cross section of the fluorophore. We are interested in reconstructing the map from the optical measurements.

2.3. Reconstruction method

We derived a method based on the first two “moments” of the measured fluorescent signal (intensity and mean time), with an explicit account of the excitation measurement, which is also called the Born approximation. The appliance of the excitation moments allowed us to get rid of the instrument response (sources and detectors specific characteristics: sources intensity, sources delays, detectors efficiencies and detectors delays, which would require cumbersome calibration procedures if we wanted to evaluate them).

Let us consider the signal: where represents the deconvolution operator. We define the intensity and the mean time t(f), defined as . From now, we will omit the dependence with the time t to simplify the writing. Using the moment’s rules, the contribution of moments of the IRF disappears and we get:

This indicates that these two quantities are the same for the model-predicted data and the experimental data, and as a matter of fact, are independent of the IRF. In practical terms, these quantities are 1) the ratio of the fluorescence intensity on the excitation intensity, 2) the difference of the fluorescence and excitation mean time. Our main hypothesis was that we were searching for a unique fluorescent inclusion. In the case of a fluorescent dot located at position , the models of intensity and mean time of S, respectively I ant t have simple expressions:

where . In particular, Eq. (9) shows that the mean time of S is independent of the fluorescence concentration.

To find the position of the fluorophore, each position m of the fluorophore was tested. Here is a summary of the different steps in our reconstruction method:

Step 1: Forward model calculation: numerical calculation of the Green functions ,,and for each source-detector couple, and for every voxel m.

Step 2: For every voxel, calculation of the optimal (Eqs. (10) and (11)) by comparing the intensities.

Step 3: For every voxel, calculation and optimization of a global criterion taking into account both intensities and mean times, with calculated at step 2:

where the operator var designs a modeled variance of the signal, taking into account the photon noise and laser intensity random variations.

Step 4: A threshold was applied on the criterion to select a localization map.

Step 5: The localization map was scaled to minimize the criterion in Eq. (11) to quantify the fluorescence.

2.4. Signal preprocessing

Several steps were required before the measurements could be handled for reconstruction. We call a measurement a set of three measures: i) an excitation measurement, ii) a fluorescence measurement, iii) a fluorescence measurement without the fluorescent inclusion. The third measurement corresponds to a parasite signal featuring leaks from the excitation and autofluorescence of the medium and system. This measurement was subtracted to the raw fluorescence data [22]. This was possible as we were working on phantom experiments.

Step 1: For each signal, a mean background value was subtracted. It was calculated on the first part of the temporal curve before the peak.

Step 2: For each raw fluorescent signal, the corresponding parasite signal was subtracted.

Step 3: For intensity calculation, the fluorescent intensity was normalized by the excitation intensity.

Step 4: For mean time calculation, a mask was applied to the data in order to remove the effects of the background. The time interval was chosen between 1% of the peak value on both edges of the peak.

2.5. Evaluation of the method

The reconstruction accuracy was studied by varying the location of the inclusion in the three directions under the imaging probe. In the first experiment, the inclusion was placed at different depths under the probe (z change). The depth of the inclusion was varied from 6.2 mm under the probe up to 18.2 mm (center of the inclusion), with a 3 mm step. In the second set of experiments, the inclusion was moved perpendicularly to the lines of sources and detectors (x change) within a 16 mm range, with a 4 mm step, for a depth of the inclusion of 9.2 mm under the probe. In the third set of experiments (same depth), the inclusion was moved parallel to the lines of sources and detectors (y change) within a 16 mm range, with a 4 mm step. Accuracy was reported through the difference of the reconstructed values and the true value that is called the residual. The median (δ), the standard deviation (δ) and the maximum absolute deviation (δ) of the residual were computed to evaluate the reconstruction.

Then we checked the precision of the method by performing ten sequential sets of measurements. Finally, the input optical parameters were manually changed in order to study the influence of false parameters on the results. The absorption coefficient, the diffusion coefficient and the fluorescence lifetime were changed on a 20% range ( ± 10% and ± 20%) and the reconstruction was processed for different z locations. We chose these values as ± 10% corresponds usually to the error made when measuring the optical properties. Thus the ± 20% range increases that range and also may account for heterogeneity inside the prostate.

3. Results

3.1. Varying the inclusion location

In the first experiment where the depth of the inclusion varied, a very high accuracy was obtained on the z location with = 0.06 cm. The results can be represented visually by plotting the calculated values and the true values, in the (y,z) plane as in Fig. 3 (a) . Figure 3 (b) shows the linear regression for these values, with a slope close to 1. The model we applied assumed a point-like inclusion, which is not the case in practice. With this model, we estimated a depth which is close to the center of the inclusion. andwere equal to 0.24 and 0.16 cm respectively. This may be explained by the experimental errors when we centered the inclusion under the probe.

Fig. 3

Results of the reconstruction (diamonds) for different true depths of the inclusion under the probe (square). The real locations of the inclusion are represented with squares, from 14 to 26 mm with a 3 mm step. (a) in the (y,z) plane (b) z calculated as a function of z true.

In the second set of experiments where the inclusion was moved perpendicularly to the lines of sources and detectors, the reconstruction gave a very good accuracy, with standard deviations of 0.15 cm and 0.07 cm for x and y respectively (Fig. 4 ). However, the x estimation was not as good for positions 0.4 and 0.8 cm as for positions −0.4 and −0.8 cm. Even if these positions were supposed to be symmetric, the difference was due to the glass tube holding the fluorophore. For positions 0.4 and 0.8 cm, the empty glass tube was centered under the probe and thus guided the light toward the negative y. This distorted part of the reconstructed location and explains that the regression on Fig. 4 (b) as well as the standard deviation on z are not as good as in the previous experiment.

Fig. 4

Results of the reconstruction (diamonds) for different lateral positions of the inclusion (constant depth). The inclusion was moved perpendicularly to the lines of sources and detectors. The real locations of the inclusion are represented with squares, from −0.8 to 0.8 mm with a 4 mm step. (a) in the (x,y) plane (b) x calculated as a function of x true.

In the third set of experiments (inclusion moved parallel to the lines of sources and detectors), the y position is well estimated with = 0.11 cm (Fig. 5 (a) ). Instead of a zero centered position, the x result is slightly shifted on negative values with = 0.15 cm. This is due again to the geometry of the glass tube we used. This was proved with the tube placed on the other side (x>0) and the shift followed on the positive values (results not shown here). The calculated depth has a good accuracy as well ( = 0.07 cm) similar to the one that was obtained for the first experiment.

Fig. 5

Results of the reconstruction (diamonds) for different lateral positions of the inclusion (constant depth). The inclusion was moved parallel to the lines of sources and detectors. The real locations of the inclusion are represented with squares, from −0.8 to 0.8 mm with a 4 mm step. (a) in the (x,y) plane (b) y calculated as a function of y true.

A summary of all the results is given in Table 1 .

Table 1

The median, the standard deviation and the maximum absolute deviation from the residual of the reconstructed values for x, y and z, for the three displacements

	z variation (cm)				x variation (cm)				y variation (cm)
	δ_med	δ_rms	δ_max		δ_med	δ_rms	δ_max		δ_med	δ_rms	δ_max
x	0.15	0.17	0.24		0.15	0.15	0.22		0.06	0.09	0.15
y	0.16	0.12	0.16		0.01	0.07	0.13		0.07	0.11	0.17
z	0.02	0.06	0.11		0.11	0.17	0.26		0.12	0.07	0.23

3.2. Repeatability

Using the ten sequential sets of measurements, the repeatability standard deviations were found to be 0.10, 0.11 and 0.06 cm for x, y and z, respectively. These precision values are of the same order of magnitude as the previous accuracy values.

3.3. Influence of background optical properties determination

Varying the absorption coefficient in a ± 20% range around the measured value has very little influence on the calculated values. The results can be seen on Fig. 6 (a) for the first experiment with different depths. The maximum standard deviation was = 0.13 cm for = 0.24 cm−1. Even if it is remaining small ( = 0.18 cm), the impact was maximal for the highest depth z = 2.82 cm for which we can clearly observe the z dispersion versus the absorption coefficient. By checking the speed’s formula in an infinite medium () for example, we observe that it increases with the absorption coefficient. This is in agreement with the plot where the calculated depth increased with the absorption coefficient. A larger path is making up for a higher speed.

Fig. 6

Calculated depth (z) as a function of theoretical depth (a) for different absorption coefficients around = 0.3 cm−1 ± 10% and ± 20% (b) for different diffusion coefficients around = 12 cm−1 ± 10% and ± 20%.

In a similar way, the impact on the calculated values was even smaller when the diffusion coefficient changes (Fig. 6 (b)). The maximum standard deviation was = 0.13 cm for = 14.4 cm−1 (highest simulated value). For the highest depth (z = 2.82 cm), we notice that the calculated depth decreased when the diffusion coefficient increases, which is in agreement with the above discussion.

For five lifetimes between 0.44 and 0.66 ns (τ = 0.55 ns ± 10% and ± 20%), we calculated the corresponding depths for the first experiment. These changes had a few impact on the (x,y) location, indeed the standard deviation remains constant, around 0.15 cm and 0.11 cm for x and y respectively. As expected, the changes were more important for the z location. The results are shown on Fig. 7 . The calculated depth decreases as the input lifetime increases. As a matter of fact, for a same mean time, the distance has to decrease in order to compensate for a larger lifetime. was minimal at 0.06 cm for τ = 0.55 ns (measured value) and was equal to 0.22 and 0.28 cm for τ = 0.44 ns ( = 0.31 cm) and 0.66 ns ( = 0.5 cm) respectively.

Fig. 7

Calculated depth (z) as a function of theoretical depth for different lifetime around τ = 0.55 ns ± 10% and ± 20%.

4. Discussion

In this paper, we have presented a new device and method for fluorescence tomography of the prostate. The optical fibers were integrated in a US probe. The acquisition time was reduced to a few minutes, which fits the constraints of the clinical environment. We successfully applied our reconstruction method to an optical phantom containing a 3 × 3 × 5 mm fluorescent inclusion. The reconstruction method was based on the search for a unique inclusion and managed a normalized Born approach. A very good accuracy was obtained both in depth and lateral displacement of the inclusion. We put some errors down to the glass tube that hold the inclusion. Further studies on very small tumors would require another type of inclusion (smaller glass tube or fluorophore embedded in a gel ball or a sponge for instance). However, we chose this size (0.045 cc) as it is in agreement and even smaller than the mean size of non palpable cancers, around 0.1 cc [2]. Besides, we restrained ourselves to a region up to 1.82 cm under the probe, assuming that the probe is squeezed one cm into the prostate. This is in accordance with the fact that most of the clinically important prostate cancers (those that are potentially morbid or lethal) are found in the peripheral zone whereas the clinically unimportant are found in the transition zone [2].

Different hypothesis may be discussed. Firstly, optical properties and fluorophore lifetime were assumed to be known. Indeed, prior information on the prostate optical properties may help the optical reconstruction. Future work will investigate two different ways to obtain the medium’s optical properties. The first one would be based on information from measurements at the excitation wavelength. The second one could use information for the US measurement. According to the results, even a 20% error on the optical properties measurement implies a small error on the inclusion localization ( = 0.18 cm). Concerning the fluorescent lifetime, we underlined that this parameter is more predominant in the depth resolution. With an average speed of light around 4 cm/ns, a 0.1 ns error may induce a 0.2 cm error in the inclusion localization. Our choice was to use the value we had measured previously in a dedicated measurement. This was justified as the fluorophore did not change its environment in our glass tube inclusion. Nevertheless, this will not be the case anymore for an in vivo experiment and it is well-known that this parameter is very sensitive to local tissue environment [23,24]. Different methods may be explored to overcome this a priori knowledge. Keeping the basis of our reconstruction method, a differential method as in reference [25] could be adapted. Other groups have also chosen different more elaborate [13,26] or simple [27] reconstruction methods in which the lifetime is an additional unknown. Another limitation for in vivo measurements is the subtraction of the parasite signal.

For an actual patient, it can be envisioned to measure this signal before injection of the fluorescent marker but one can imagine both patient discomfort and potential calibration problems. Another possibility would be to find a spot on the tissue without markers. However, we are currently working on a feasible method to establish a model for the parasite signal and apply a temporal filter to get rid of this signal.

Furthermore, we presented an algorithm dedicated to the search of a unique inclusion as a first step in our studies. Prostate cancer tends to be multifocal in 50-76% in all cases of radical prostatectomy [28], moreover, they tend to be bilateral for small-volume prostate cancers [29]. The average number of prostate tumors is between 2 and 3 [28,30]. These facts justify the search for a unique tumor in each of the bilateral parts of the prostate. We are currently working on algorithms for detection of two tumors to ensure an even better reliability. Besides, we did not address the quantification problem for now as we assumed that this was not crucial for biopsies guiding.

In the results presented here, background fluorescence was set to zero. As fluorescent markers are not developed yet for the prostate, very few information is available on the expected levels of specificity. Ongoing work is studying the effects of non specific fluorescence on reconstruction and new methods to overcome this signal.

5. Conclusion

A novel device for time-resolved imaging of the prostate was presented. Our original method was derived and evaluated. High accuracy was shown on localization of a 45 µL fluorescent inclusion of ICG in an optical phantom mimicking the prostate. We obtained accuracy better than 0.15 cm in the three directions of space and a precision around 0.1 cm. This is compatible with the size of early stage tumors. The ultimate goal of this device is to work simultaneously with the US imaging. After phantom validation for systems coregistrations, further investigations will be then needed on canine animal model for example to demonstrate the benefits provided by an optical diagnosis.