Diffuse optical tomography to monitor the photocoagulation front during interstitial photothermal therapy: Numerical simulations and measurements in tissue-simulating phantoms

2.1 Simulations

To calculate the expected DOT optical signals for given geometries of source-detector positions and coagulation zone forward simulations were performed using NIRFAST [33], an open-source finite-element-method (FEM) software suite based on the diffusion approximation to photon radiation transport. The overall geometry is defined in Figure 1B, where a cylindrical coagulation zone, as generated by a linear diffusing fiber placed within the tumor target in the prostate to deliver the photothermal treatment [16, 32], expands in radius with its center fixed at a depth within the prostate. The rectum, prostate and all surrounding tissues of the coagulated lesion are considered as optically homogeneous with the same absorption and scattering coefficients. The computational mesh is large enough that the background can be considered essentially semi-infinite in extent, uniform and unchanging during the treatment. This is clearly an approximation to the real clinical situation but allows us to focus on the critical issue of the detectability and localization of the boundary of the coagulation zone without being concerned about the detailed structure of the rectum, the prostate gland and surrounding tissues. Tetrahedral FEM meshes of dimensions of 120 mm × 50 mm × 100 mm, a node-to-node distance of 2 mm and a total of 26,026 nodes were generated to simulate the homogeneous tissue background. A cylindrical anomaly representing the photocoagulated zone within the prostate, centered 10 mm from the ventral aspect of the mesh surface (representing the inner surface of the 2-mm thick rectal wall that is in contact with the DOT probe) and with a radius of 2 mm, was placed inside the homogeneous mesh to generate a mesh representing the lesion with 2 mm radius. This was repeated a further six times, increasing the radii from 3 to 8 mm in 1-mm increments, generating a total of 8 meshes. The smallest distance between the corresponding posterior coagulation boundary and the outer rectal wall (assuming a 2-mm thick rectum) then varies between 6 and 0 mm as the radius (denoted as “r”) increases from 2 to 8 mm. This distance is denoted as “d” throughout this paper (note that: d+r=8 mm). The optical properties assigned to the homogenous mesh as well as to the anomalies corresponded to the values used in the experimental phantoms for the different wavelengths employed (see Table 1) and were selected to be representative of typical soft tissues in this spectral range [29, 34].

Three different source-detector configurations were simulated with coronal and axial views schematically displayed in Figure 2. The first corresponded to the experimental set-up, as illustrated in Figures 2A and 3, with a line of 7 sources and a parallel line of 7 detectors with a 20 mm separation, which gives source-detector pair separations ranging from 20 to 63 mm. For this configuration all possible source-detector pairs were simulated, generating 42 measurements (not 49 measurements due to one bad detector channel as depicted in Figure 3A). Since this configuration was designed primarily for imaging the whole prostate, it provides a large depth of penetration and field of view and also has the advantage of being relatively insensitive to the exact position of the probe in the rectum. However, this may come at a cost of decreased sensitivity in monitoring the position of the coagulation front as it approaches the rectal wall.

Figure 2

Schematic representation of the simulated source-detector configurations relative to the coagulated lesion, showing both coronal and axial views for: (A) Experimental TRDOT probe [23] used in the phantom measurements; (B) Single line of source and detectors aligned in transverse direction; and (C) Single line of source and detectors aligned in longitudinal direction.

Figure 3

Experimental configuration used for phantom measurements.

(A) Schematic of the source-detector array and photograph of the experimental TRDOT probe [23]. Note that detector D3 was not functioning in the experimental system but is shown for completeness. (B) Illustration of the geometric relationship between the source-detector array, the coagulation zone, the prostate and the rectal wall in the solid phantom. (C) Photograph of the TRDOT probe in one of the solid phantoms showing the holes used to contain liquid phantom material that correspond to the coagulated region and (D) Configuration of the holes on the two solid phantoms.

The other two configurations – transverse and longitudinal (Figure 2B and C) – were therefore simulated also including smaller source-detector pairs (2 mm) than in the experimental configuration. In these configurations, the 12 sources and 12 detectors are placed along a single line, with sources occupying one half segment and detectors the other segment. This line is then placed along the transverse axis in the second configuration and along the longitudinal axis in the third configuration. Details of the source and detector locations for these both configurations are shown in Figure 4. While 12 sources and 12 detectors were available, only 12 source-detector pairs were simulated, with each pair being symmetrical to the middle point of the line (denoted as a black dot in Figure 2). Hence, the signal from source j is only measured at detector j. Using a 3.2 GHz processor, with 12 MB RAM on a 64-bit Windows 7 platform, these calculations required 5 s to perform.

Figure 4

Relative positions of sources and detectors used in the simulations of the transverse and longitudinal configurations. Measurements are acquired at source-detector pairs symmetric to the center of the configuration, i.e., light from S6 is measured only at D6. The view between S8 to D8 is expanded to provide greater details of the relative source-detector positions.

2.2 Tissue-simulating phantoms

Tissue-simulating phantoms were designed and fabricated, as shown in Figure 3B–D. The main solid volume was made from silicone impregnated with India ink as an optical absorber and titanium dioxide powder to provide optical scattering. This corresponds to an optically homogenous semi-infinite volume representing the rectum, the prostate and surrounding tissues (Figure 3B), which for this purpose were assumed to be the same. Each block contains a 30 mm×30 mm square hole in the center running the length of the block, into which the DOT probe (which has a flat working face housing the light source-detector array) could be inserted (Figure 3C). This inner surface then represents the luminal (or inner) surface of the rectal wall. Seven cylindrical holes of circular cross section and of different diameter, from 2 to 8 mm in 1 mm increments, centered at 10 mm from the closest flat face of the central square hole and running parallel to it, were made to represent the growing photocoagulation zones (Figure 3D). Hence, one of the phantoms contained four such optical anomalies, while the other had three of them. The separation between these holes is sufficient to minimize any influence on the detected light field when any one of them is in use. It is of note that, as currently practiced clinically, PTT uses cylindrical diffusing fibers that are typically 1–3 cm long and so generates near-cylindrical coagulation volumes except at the ends. It can also be assumed that the DOT measurements are sampling the cylindrical region and that the treatment fibers are parallel to the rectal wall, which is a reasonable approximation.

Intralipid^® (Baxter, Canada) at three different dilutions of 20% stock solution and containing India ink as an absorber were used to fill the cylindrical holes, generating three sets of optical scattering and absorption contrast, as per Table 1. The absorption spectrum, μ_a (λ), of India ink was measured in a spectrophotometer and added to the intrinsic absorption of the silicone and Intralipid^® [35]. The reduced scattering spectrum, μ′s(λ), of the solid phantom was measured using an established frequency-domain diffuse reflectance technique [36, 37], while the reduced scattering coefficients of various dilutions of the 20% Intralipid^® stock solution at 690 nm were measured using an in-house CW-DOT system that will be the basis for the first clinical prototype and will be described in a future publication. The values at each of the three wavelengths were then calculated using the following equation, which is based on Mie scattering theory [38]:

in which, for a phantom with a given dilution of a 20% Intralipid^® stock solution, the reduced scattering value at wavelength λ₂ can be calculated if the value at wavelength λ₁ is known.

2.3 DOT experiments

The phantom experiments carried out for cross-validation with the simulations used a CW-TRDOT system that has been previously reported by Piao and colleagues [24] for in vivo canine prostate cancer imaging studies. The source and detector fibers are mounted on a near-flatface substrate attached to a cylindrical 2 cm diameter ultrasound probe. It is of note that one of the detector channels (D3) was malfunctioning during these experiments, so that the corresponding simulation data have been excluded from the analysis. The probe was inserted into the central hole (representing the rectum), which was made square to ensure good contact with the sources and detectors on the probe face.

The measurements were performed over three separate sessions. At the beginning of every session, a set of measurements was first made with the sources and detectors facing the homogenous side, i.e., away from any of the anomalies representing the coagulation zone. These three sets of data served to calibrate all the sets of heterogeneous measurements made in the same session, as described below. After these homogeneous measurements were completed, the probe was rotated to face each one of the optical anomalies in turn. Measurements were then taken sequentially for each of three Intralipid^®/ink composition filling the anomaly. In the first session, a set of homogenous measurements and a set of measurements of hole #7 (r=8 mm) were taken. In the second session, a set of homogenous measurements and then three sets of measurements with holes #5 (r=6 mm), #3 (r=4 mm) and #6 (r=7 mm) were taken. In the final session, a set of homogeneous measurements and three sets of measurements with holes #4 (r=5 mm), #2 (r=3 mm) and #1 (r=2 mm) were taken (see Figure 3D for locations of the holes). For each probe orientation, measurements were made sequentially with each of the three solutions of liquid phantoms (contrast #1 – contrast #3) without moving the probe. For every combination of anomaly size and optical contrast, each of the seven sources (each multiplexed with three wavelengths, 705, 785, and 850 nm) was activated in turn and the signals from each of the six functioning detectors was recorded by the charge-coupled device (CCD) at all three wavelengths simultaneously. Thus, for three sets of contrast ratios and seven anomaly sizes 21 sets of heterogeneous measurements were acquired, in addition to three sets of homogeneous measurements. Measurements were also made at 940 nm, however these data were inconsistent for reasons that are not well understood so that they are not included in the analysis. The data were subsequently analyzed off-line.

2.4 Calibration

The calibration of these measurements is critical, since there are several factors that need to be taken into account, including the:

• differences in the output power of each source at each wavelength and the variations of these over the 2-day course (three sessions) of the experiments;

• corresponding varying sensitivities of the detector channels; and

• variations in repositioning of the probe between sets of experiments.

The simulations for homogeneous meshes at each wavelength were used to normalize the experimental data acquired within each measurement session, with the calibration factors calculated according to:

and

I_inhom (i)=M_inhom (i)·C_i,(3)

where I_simu (i) is the simulated intensity (particle fluence in mm^-2) for the homogenous mesh, M_hom (i) is the corresponding experimental measurement, M_{in hom} (i) is the measurement on the phantom containing the anomaly and I_{in hom} (i) is the corresponding corrected measured intensity for source-detector pair i. Thus, a set of 42 calibration factors, C_i, were calculated for each wavelength and then applied to all experimental measurements in a given measurement session.

As a measure of the agreement between the simulations and measurements, the following metric was calculated:

where N=42 is the number of source-detector pairs or measurements within a single set.

3.1 Comparison of experimental and simulated data

Figure 5 shows representative plots of measurements and the corresponding simulated data for the dependence of the optical signal on the detector number when specific sources were activated, and how this varies as the photocoagulation front approaches the rectal wall. The specific cases shown in Figure 5 are for contrast #3 at 850 nm when either source S5 or source S7 was activated. As expected in both cases, the detected signal falls off monotonically with increasing source-detector distance, with the gradient being highest for smallest separations. It can be seen that the measurements and simulations show good agreement in how the detected signal changes, both as the scattering anomaly grows (and the posterior edge approaches the rectal wall) and as the particular source-detector pairs are varied.

Figure 5

Comparison of experimental data (left) and corresponding simulations (right) for the signal from each detector for contrast #3 at 850 nm when either source S5 (A, B) or source S7 (C, D) is activated. The lines are simply to guide the eye. The background values of optical properties are for the homogenous phantom case, i.e., no scattering or absorption anomaly. Note that detector D3 was not functioning, so that these data are omitted. The legend indicates the distance, d, of each anomaly to the outer rectal wall.

Figure 6 shows summary plots for 3 representative anomaly sizes of the experimental vs. simulated signals for all the data (42 in a set) at 850 nm for contrasts #1 (Figure 6A) and #3 (Figure 6B). Overall, the two sets of data are in good agreement, but there is evidence of a systematic deviation at larger anomaly sizes.

Figure 6

Measured vs. simulated signals at 850 nm for 3 representative sized anomalies for contrast #1 (A) and contrast #3 (B). The legend indicates the anomaly radius, r, while the line of equality is also plotted.

This systematic difference is evident also in the goodness-of-fit values presented in Tables 2 and 3, which show that the agreement between the simulated and experimental data is also best for the lowest scattering contrast. Several possible reasons may explain this trend. First, poor contact between the probe and rigid phantom when performing measurements at several positions may cause an increased measured signal. This potential effect was only discovered during analysis after all data were collected. There may also be a discrepancy between the measured and true values of the optical properties of the anomaly, the effect of which would become more evident as the size of the anomaly increases: the measured signals are higher than predicted by the simulation, which would suggest that either the optical attenuation of the experimental anomaly is too low or that the limited dynamic range of the instrument cannot measure the expected attenuation. Also, since the optical properties of the dilutions of 20% Intralipid^® solutions were measured at only one wavelength (690 nm) and then calculated for all respective wavelengths, any error in the use of Eq. (1) would propagate. To determine the change in optical properties required to match the measured and simulated data, simulations were repeated using a range of optical properties. Full agreement was achieved when the background absorption was 2× greater than the absorption of the anomaly, significantly higher than the measured contrast (∼1:1) and beyond any measurement error in the optical properties of the silicone phantom [36, 39]. The limited mesh resolution (2 mm node-to-node distance) may also contribute to the discrepancy; the large volume required to model the experimental set-up limited the resolution. Increasing the resolution for the same volume would increase the number of nodes, but this was limited in our case by the available computer memory. Alternative probe configurations with smaller source-detector separations and meshing algorithms that increase the node density near the surface would reduce this effect. Regardless, the trends in the simulations are more consistent with the expected physical behavior of the signal.

3.2 Sensitivity to the approach of the photocoagulation front

Figure 7 shows the same measured and simulated data depicted in Figure 6 for contrast #3 but now as relative changes in signal. In Figure 7A, the relative change in signal is plotted for anomalies of 2, 3, 5 and 8 mm radius for all source-detector separations. For the r=2, 3 and 5-mm anomalies, the measured and simulated trends are similar, although the magnitude of the changes does not match (see above). The impact of poor probe-phantom contact can be seen in the measured data for the r=8-mm anomaly; here the measured data show a decreasing relative change as the source-detector separation increases while the simulations show attenuation generally increasing with larger separation. The same data set is plotted in Figure 7B but as the anomaly increases in size and hence as it approaches the rectal wall. Measurements at source-detector separations of 20 and 22 mm are shown, while the remaining lines show simulated results. In all cases, the signal decreases by approximately 10% per mm. Despite the discrepancy between the measured and simulated signals in Figure 7A, the large change in signal as the anomaly progresses towards the rectal wall indicates the likely utility of the transrectal measurements for tracking the evolution of the photocoagulation front.

Figure 7

Relative change in signal, as the photocoagulation anomaly approaches the rectal wall, for 850 nm and contrast #3.

Measured data (exp) are represented by solid lines and black symbols while simulations (sim) are represented by dashed lines and white symbols. (A) Changes for 2, 3, 5 and 8 mm diameter anomalies as a function of source-detector separation. Identification of each set of data is provided on the right of the figure. For the smaller anomalies, both measured and simulated data show greater attenuation with increasing separation. For the 8-mm anomaly, poor contact between the probe and phantom surface generates an artifact of increasing signal change. (B) Attenuation for each source-detector separation as the anomaly approaches the rectal wall.

The geometry of the above phantom measurements and corresponding simulations is that the line of sources and line of detectors are parallel to each other and to the longitudinal axis of the prostate. This configuration is designed to probe the whole prostate in which the shortest source-detector separation, 20 mm, is larger than the radius of any expected anomaly created by photothermal therapy. Two alternative configurations (Figure 2B and C) were considered to examine their sensitivity to the expanding coagulation front. In these configurations, the source-detector separation ranges from 2 to 60 mm (Figure 4), with all sources and detectors lying on the same line and with every active source-detector pair symmetric to the same midpoint. As would be expected for these geometries, the signal is strongest when the source and detector lie directly opposite to one another, i.e., when their separation is smallest, and falls rapidly as the separation increases (as shown for the experimental configuration in Figure 5). For any given source-detector separation, the signal also falls monotonically as the radius of the anomaly increases, but not in a regular step-wise fashion: for small anomalies (r≤3 mm) there is little change in the signal, even for the closest source-detector pair. In this case the fraction of detected photons that pass through the anomaly is small, due in part to its greater depth within the prostate. As the scattering anomaly grows and approaches the probe (i.e., towards the posterior edge of the prostate and the rectal wall), the signal depend on size (and depth) of the anomaly increases.

Figure 8A and C shows, for the transverse and longitudinal configurations respectively, the fractional change in the detected signal with increasing source-detector separation, as the scattering anomaly grows larger and the photocoagulation front gets closer to the DOT probe. Figure 8B and D shows the same data but with signal changes plotted as a function of anomaly size for different source-detector separations. These plots essentially reflect how much of the diffuse light field is intersected by the scattering anomaly. Where this relative intersection is largest, the signal experiences greater attenuation. At source-detector distances in the range of 6–10 mm, the detected signal initially increases as the lesion grows (Figure 8B and D) and then sharply decreases when the distance between the lesion and rectal wall reaches d=2 mm. At these distances, the primary diffuse field does not interact with the small highly-scattering lesion; however the lesion may also act as a partial diffuse reflector and increase the signal. As the lesion grows into the diffuse field for these distances, the signal is attenuated as it traverses the lesion. The increasing signal seen at small source-detector separation is likely due to the increase in the backscatter of light from the anomaly being greater than the loss of light within it. Small anomalies (which are also most distant from the probe) are essentially undetected at small source-detector separations, since they intersect a minimal portion of the diffuse light field. Hence, it would require very high signal-to-noise ratio to detect the initial growth of the scattering anomaly if this lies deep within the prostate. Conversely, as the scattering anomaly becomes large and its posterior boundary approaches the probe, the signal becomes extremely sensitive to its growth, especially for the smaller source-detector separations, since the light field is significantly perturbed. From these plots, there are no striking differences between the transverse and longitudinal source-detector configurations. The changes in relative signal are larger (10–30%) for the longitudinal configuration, since more of the anomaly falls within the diffuse field for most source-detector pairs. The large change in signal response at different source-detector positions in these configurations as the lesion grows (as seen in Figure 8B and D) raises the possibility to recover the location of the posterior edge of the lesion boundary without performing full tomographic reconstruction of the lesion.

Figure 8

Percent change in the detected signal relative to the signal with no anomaly (i.e., before the start of treatment) plotted for different scattering anomaly sizes (2 to 8 mm) as a function of the source-detector separation (left) and the distance between the posterior edge of the anomaly and the nearest surface of the rectal wall (right) for the source-detector separations listed in the legend. (A, B) Transverse source-detector geometry, and (C, D) longitudinal geometry. Optical properties are for contrast #3 at 850 nm (see Table 1).

For the largest anomaly investigated here (r=8 mm, corresponding to the photocoagulation front touching the inner surface of the rectal wall, i.e., d=0 mm and 2 mm from the face of the DOT probe), the slope of the signal vs. anomaly radius is ∼8% per mm for contrast #3 at 850 nm, similar to that observed with the experimental configuration described above (∼10% per mm). Thus, it should be feasible to meet the clinical goal of determining the position of the posterior border of the coagulation zone within ±∼1 mm as the coagulation approaches the posterior capsule of the prostate, especially if the scattering contrast due to photocoagulation is high.

The changes in optical signal with increasing anomaly size were generally invariant to the background and anomaly absorption properties, as seen in Figure 9, where the background absorption of the background tissue and anomaly is μ_a=0.03 mm^-1, representing an increase by a factor of five from the data depicted in Figure 8. As with the low absorption data, at the longest source-detector separations the rate of attenuation is larger than with the lower absorption values. At shorter source-detector separations, scattering effects dominate; for 2 and 4-mm separations, the signal increases as the photocoagulation anomaly nears the rectal wall, while for separations of 6–10 mm, the signal first increases then decreases.

Figure 9

Percent change in the detected signal relative to the signal with no anomaly (i.e., before the start of treatment) plotted for different scattering anomaly sizes (2 to 8 mm) as a function of the source-detector separation (left) and the distance between the posterior edge of the anomaly and the nearest surface of the rectal wall (right) for the source-detector separations listed in the legend. (A, B) transverse source-detector geometry, and (C, D) longitudinal geometry. Scattering properties are identical to those used in Figure 8, while absorption was increased by a factor of five to μ_a=0.03 mm^-1.

3.3 Optimal clinical configuration

While the three alternative DOT configurations showed sensitivity to a growing coagulation lesion, each has advantages and disadvantages in the clinical setting. The experimentally tested configuration was originally designed for whole prostate imaging [23–25]. With its large field of view, lesion monitoring is relatively insensitive to probe placement. However, its sensitivity to the lesion growth is realized only as relative decrease in intensity. The rate of attenuation is similar for all source-detector pairs (see Figure 5), and hence this configuration loses sensitivity as the lesion nears the rectal wall.

In contrast, the two simulated configurations (transverse vs. longitudinal) using much smaller source-detector separations were especially sensitive to the lesion approaching the rectal wall. Conceptually, the longitudinal geometry samples along the length of the roughly cylindrical coagulation lesion and so may be more sensitive to misalignment of the probe with respect to the shape and location of the coagulation zone. The transverse geometry samples across the coagulation boundary and so may be more sensitive to the placement of the probe in the superior-inferior axis. It would also sample only one transverse plane and could potentially miss the lesion growing past the outer rectal wall in another transverse plane. In practice, the probe could be placed at the transverse plane where the distance between the treatment fiber and the rectal wall is the smallest. Assembly of a probe using the transverse geometry would be easier than the longitudinal configuration, since optical fibers could run parallel to the direction of the probe; in the longitudinal configuration, fibers would need to curve into a single longitudinal line. The ideal probe may be one that combines elements of both geometries by probing, for example, three transverse planes separated by several millimeters.