Differential diffuse optical tomography

Abstract: We formulate a perturbative solution for the heterogeneous diffusion equation which demonstrates how to use differential changes in diffuse light transmission to construct images of tissue absorption changes following contrast agent administration. The analysis exposes approximations leading to an intuitive and simplified inverse algorithm, shows explicitly why transmission geometries are less susceptible to error than the remission geometries, and why differential measurements are less susceptible to surface artifacts. These ideas about differential diffuse optical tomography are not only applicable to tumor detection and characterization using contrast agents, but also to functional activation studies with or without contrast agents and multi-wavelength measurements. ©1999 Optical Society of America OCIS codes: (110.6960) Tomography; (170.6960) Tomography; (170.3010) Image reconstruction techniques; (170.5280) Photon migration; (170.3660) Light propagation in tissues; (170.5270) Photon density waves; (999.9999) Contrast agents


Introduction
The use of contrast agents in Diffuse Optical Tomography [1] (DOT) for disease diagnostics and for probing tissue functionality follows established clinical imaging modalities such as Magnetic Resonance Imaging (MRI) [2,3,4], Ultrasound [4,5,6] and X-ray Computed Tomography [4,7,8] (CT).Contrast agent administration provides for accurate difference images of the same heterogenous tissue volume under nearly identical experimental conditions.This approach yields superior diagnostic information.In MRI of breast cancer for example, it is the relative signal enhancement compared to the baseline image, and the architectural features of this enhancement that aid in characterizing lesions as benign or malignant [9,10].While architectural features are unlikely to improve breast cancer diagnosis using DOT in current implementations, due to the low resolution achieved [11,12], the relative signal enhancement after contrast agent administration could significantly increase sensitivity and specificity.
The theory presented in this paper is motivated by clinical DOT measurements of the human breast [13].In these measurements optical data are obtained before and after intravenous administration of the optical contrast agent Indocyanine Green (ICG).ICG is an absorber and a fluorophore in the Near Infrared (NIR); here we focus on its absorption function.In principle DOT images of the breast before and after ICG administration may be reconstructed and subtracted.However, in practice a more robust approach derives images of the differential changes due to the extrinsic perturbation.In the latter case experimental measurements use the exact same geometry within a few minutes of one another, thereby minimizing positional and movement errors and instrumental drift.Furthermore the use of differential measurements eliminates systematic errors associated with the different medium required to calibrate operational parameters of the instrument or to provide a baseline measurement for independent reconstructions of the pre-and post ICG images.In addition the effect of surface absorbers such as hair or skin color variation is also minimized.
The main analytical difficulties of the differential approach arise because the media are inhomogeneous.Thus the total diffuse light field in the ICG perturbation problem does not separate into a homogeneous background field and a scattered field in a straightforward way.Furthermore the average background properties, particularly the absorption, change as a result of contrast agent administration.
In this paper we derive a perturbation solution that can be used to yield differential images of induced absorption perturbations.This is the primary effect of many molecular optical contrast media that are administered in low concentration.The solution is also applicable for imaging the differential variations in tissue scattering.In obtaining these results we will expose the full consequences of inhomogeneous tissue volumes, and identify the conditions wherein these effects are small.We demonstrate the superiority of the method compared to the typical perturbation approach using simulated data of contrast-enhanced breast.The method is also suitable for imaging differential changes as arise in measurements of functional activity such as brain activation or muscle exercise (the only requirement is that baseline measurements of the same tissue volume are taken during the procedure), or in measurements at multiple wavelengths.In the following analysis the term pre-ICG breast is used to indicate the "baseline" breast before the contrast agent injection and the term post-ICG breast marks the breast following contrast agent administration.

Theory
Photon transport in highly scattering low absorbing media such as tissue can be approximated as a diffusional process [14].For heterogeneous media, the diffusion equation in the frequency domain is [15]: Here U( r r ) is the photon fluence [W⋅cm -2 ], ω is the source modulation frequency, c is the speed of light in the medium [cm⋅s -1 ] , ) ( a r r µ′ the medium absorption coefficient [cm -1 ] and ) (r S r the source term [W⋅cm -3 ].In this work we consider solutions of the heterogeneous diffusion equation in the frequency domain employing the perturbation approach [16,17].
The first order perturbation expansion divides the absorption ( ) (

′
, is the field that would have been detected from the same medium if these heterogeneities were not present.These terms are incorporated into the diffusion equation, whose formal solution can be expressed as an integral equation using the appropriate Green's function for the geometry implemented.For simplicity we present this analysis for an infinite medium.This theory however can be easily extended to other simple geometries such as semi-infinite or slab, using weights derived with the method of image sources [18] and the appropriate extrapolated boundary condition [19,20,21].
The first order perturbative solution of the heterogeneous diffusion equation yields represents the absorption (scattering) weight of the voxel at position r r , due to a source at s r r and a detector at d r r .In the infinite medium these weights are: Following the administration of a contrast agent the background optical properties change.The new, post-ICG, total field can be written in a similar form: .The first order perturbative solution of the heterogeneous diffusion equation yields where W a (W s ) represents the absorption (scattering) weight of voxels at position Combining Eq.2 with Eq.7 we obtain the relative scattered field, i.e We will show that ) , , (  During the administration of an absorption contrast agent the scattering properties of tissue are not expected to change.Therefore we assume 0 0 0 . Substitution of Eq.3 and Eq.8 into Eq.12 yields: , , ( Here we have also assumed . This is a very good approximation when the average absorption change due to the contrast agent is small or in the transmission geometry (see Appendix and discussion for absorption variations below).
r δµ be the total absorption perturbation due to the ICG injection that includes both position-independent and position-dependent contributions.Then The quantity represents the position-dependent absorption heterogeneities induced by the contrast agent.The relative scattered field is computed by substitution of Eq. 15 into Eq.13.It depends on contrast agent induced absorption heterogeneities and on pre-ICG tissue absorption heterogeneities.

∫ ∫
The second integral in Eq.16, describes the influence of the pre-existing (intrinsic) absorption heterogeneity of the breast on the relative scattered field.Since the intrinsic heterogeneity is weighted by the difference the influence of this term can be quite small.Using the analytical forms of the weights (Eq.4 and Eq. 9) we can write out Eq.17 explicitly, i.e.

(
) where The term in Eq.18 is approximately unity and S≈0 when the average absorption increase due to the ICG injection is very small (i.e. ). Usually however . For example the recommended ICG dosage for humans (0.25mg/kg) introduces an average µ a increase within the interval [0.005-0.015]cm -1 depending on breast vascularization [13].Fig. 1a and b show the amplitude and phase of the term  when α increases).However, the probability for photons to pass through these "distant" perturbations decreases exponentially via the weight a W ′ ′ in the integrand of Eq.18.Hence accumulated contributions of the heterogeneities at large α are small.Figure 2  by taking S=0. Figure 2a depicts the ratio of the amplitude detected with no approximation to the amplitude detected assuming S=0.Similarly Fig. 2b depicts the phase shift between the phase detected with no approximation and the phase detected assuming S=0.The error is plotted for a single perturbation at different positions α for the geometry depicted in Fig. 1c.The values assumed in Eq. 16 were The simulation of Fig. 2 explicitly shows that the errors introduced because of the approximation S=0 are very small for physiologically relevant optical properties (i.e.relatively small ) ( a r r µ δ ′ ) provide the most probable photon paths.The same behavior is exhibited for the scattering weights as shown in the Appendix.Eq.16 thus becomes Our conclusions do not change when image sources are invoked to satisfy more complex boundary conditions such as semi-infinite or slab geometries.In these cases will appear in all the terms corresponding to image sources.Note however that the assumption that S≈0 is best suited for slab geometry where for the most probable photon paths.This condition is not always true for reflectance geometry.
Notice that the differential measurements are insensitive to surface artifacts such as small skin absorbers and hair under a certain source or detector.The term in Eq.18 could be used to approximate surface heterogeneities by taking r r to be close to medium surface, near to the corresponding source or detector.The influence of such terms is virtually zero since in such a geometry 0 ) ( ≅ r R r and subsequently S=0.For image reconstruction, Eq. 20 is discretized into a sum of volume elements [22] (voxels) and the scattered field is obtained at each modulation frequency ω for every sourcedetector pair ω.For n voxels and m=o×p×q measurements, where o is the number of sources, p is the number of detectors and q is the number of frequencies employed, the discretization yields a set of coupled linear equations Inverting the weights' matrix determines the spatial map of absorption due to contrast agent injection.

3.Results and Discussion
Although equations 20 and 21 resemble the result of typical perturbation analyses, there are fundamental differences and constraints that must be considered when using them.First the parameter imaged is the synthetic perturbation term . This term expresses the change in the incident field due to the average absorption coefficient increase of the post-ICG breast.Its use in Eq.12 leads to significant reconstruction improvements [23].Note that the term can be analytically calculated for simple geometries such as infinite, semi-infinite or slab or calculated numerically for more complicated geometries if we know the average optical properties of the pre-and post-ICG breast.
Using simulated data we have examined the performance of Eq.20 to image the contrast-enhanced breast as compared to: 1) The typical Rytov approximation which assumes This comparison captures a critical feature of our work.The classical perturbation approach that does not consider the average absorption increase due to the extrinsic contrast.
2) Eq.17 including the term S, namely This comparison investigates the effect of the approximation S=0 assumed in Eq.20.
3) A difference image produced by subtracting the post-ICG image from the pre-ICG image.Both pre-and post-images were produced using Eq.3, assuming a homogeneous medium as baseline.The optical properties of the homogeneous medium were 0 a µ ′ =0.03cm -1 and 0 In order to perform the comparisons we obtained two MRI coronal slices of a human breast: one before and one after contrast enhancement.Figure 3a depicts the T 1 -weighted MR image.This image depicts structure.White regions correspond primarily to adipose (fatty) tissue while dark regions correspond to parenchymal (glandular) tissue.Figure 3b depicts the signal enhancement of the same T1-weighted image due to injection of the MRI contrast agent Gd-DTPA.The Gd-DTPA enhancement is superimposed in color.An infiltrating ductal carcinoma (shown in yellow) demonstrated the highest signal enhancement.Gd-DTPA and ICG have similar distribution patterns.Here we assume that the Gd-DTPA distribution reflects the ICG distribution.We converted the MR images to optical property maps, separating four structures based on the image intensity information (by applying appropriate thresholds).The cancer is assumed having two states: pre-and post-ICG contrast.The structures selected and the corresponding absolute optical properties are shown in TABLE I.The optical properties are chosen to simulate breast optical properties as obtained from our breast clinical measurements [13].Scattering has been assumed constant for all structures for simplicity.The resulting absorption maps are shown in figure 4. The medium surrounding the breast was arbitrarily simulated as a highly absorbing diffuse medium (µ a =0.30cm -1 µ s '= 8cm -1 ).The average absorption of the pre-and post-ICG breast were found to be   For reconstruction purposes, the region of interest (indicated in Fig. 5 as a green double line) was segmented into 35x25 voxels.The inversion was performed using the Algebraic Reconstruction Technique (ART) [16], with relaxation parameter λ=0.1.Convergence was assumed when an additional 100 iterations did not change the result more than 0.1%.The simulated image and the reconstructed   The superior performance of Eq.20 compared to the typical perturbative formulation (Eq.22) can be evaluated by examining Fig. 6a and b.Although both methods resolve the cancerous lesion with comparable positional and size accuracy, the typical formulation (Fig. 6b) yields several strong artifacts close to the boundaries.These artifacts illustrate that in the presence of distributed absorption, the perturbative method converges preferentially to localized regions of high absorption.This is often true when inverting underdetermined systems.Eq.20 on the other hand removes the "average absorption increase" from the measurement vector.Therefore Fig. 6a images weaker perturbations introduced by the ICG injection, relative to the average absorption increase.Since by construction the perturbation method works especially well for weak perturbations [16], it is expected that the use of Eq. 20 will more accurately image the heterogeneous medium.The same behavior is expected for a Born-type [16] perturbative formulation.We note that Fig. 6b 6a and 6c.Therefore it is reasonable that the reconstructed value for cancer in Fig. 6b is higher than the value reconstructed in Fig. 6a and 6c.The difference in reconstructed values equals approximately the average absorption increase in the post-ICG breast ( 0 0 a a µ µ ′ ′ ′ − =0.0116 cm -1 ). Figure 6c has been produced after correcting the measurement vector with S instead of setting it to zero as in Fig. 6a.Only minor differences exist between the two images as had been predicted in Figure 2. In this simulation the pre-ICG cancer had a contrast of 2:1 to the average pre-ICG background value.This contribution has most likely resulted in the minor differences observed between the two images, especially in the structures close to the boundaries.
Figure 6d is the result of subtracting an image of the post-ICG breast from an image of the pre-ICG breast.Here the ) (r ICG a r δµ is imaged.The magnitude of the cancer is slightly overestimated and its size is significantly overestimated.Similarly to Fig. 6b, strong artifacts appear close to the boundary.A distributed absorption is also reconstructed which does not correspond to the ICG distribution and is also an artifact.Compared to the other approaches the subtraction yields the most artifacts. In these simulations the average optical properties were known by simple integration over the optical property map.In our clinical implementation [13] the average optical properties of the pre-ICG breast are calculated by fitting the experimental timedomain data obtained to the appropriate solution of the diffusion equation for the geometry used.Furthermore an algorithm has been developed [24] that calculates the difference ) with an accuracy of the order of 10 -3 cm -1 .
To conclude we have described a formulation of perturbation theory, which is particularly well suited for image reconstructions of differences in the absorption properties of tissues as a result of optical contrast agent administration.Importantly, these results enable the extraction of differential contrast agent absorption even within media that are heterogeneous in the absence of the contrast agent.Our primary result is an intuitive equation, which is valid over a large range of conditions.We have shown explicitly what these corrections are and how these corrections can be included in more careful analyses.The results should be applicable for a broad range of other DOT applications wherein baseline and "stimulated" measurements are available, particularly functional imaging in brain and muscle.for the geometry of Fig1a and for the same absorption variation of the post ICG-breast.
Figure A1 shows that the amplitude and phase of R remains very small for the most probable photon paths, (i.e. when the distance α is small).For higher α, R increases but similarly to the arguments that reduced Eq.16 to Eq.20, the contribution of terms for higher α is weighted less.Therefore for the small absorption changes considered here and for small variations of the diffusion coefficient the scattering terms can be ignored when reconstructing the absorption perturbation only due to contrast agent injection.When introducing boundaries, the assumption of small α compared to the source-detector distance for the most probable photon paths works better for transmittance geometry.

′.
) coefficients of the pre-ICG breast into spatially varying ( Throughout this paper a single ′ denotes pre-ICG tissue volumes.In the Rytov approximation[16] the total photon density wave measured at position d r r due to a source at position s r r is written as the product of two components, i.e.
component scattered from the post-ICG heterogeneities (i.e.
primarily to perturbations created by the contrast agent injection.
from the average optical properties of the pre-and post-ICG breast (see discussion section).
as a function of the post-ICG breast absorption coefficient for a source detector separation =6cm, using the geometry of Figure1c.The background µ a =0.05 cm -1 and the background µ s '=10cm -1 .

Figure 1 .
Figure 1.(a) Amplitude of the term ) ( ) ( r R k k i e

Figure 2 .
Figure 2. (a) Amplitude deviation and (b) Phase shift introduced in the field measured in Eq. 16 when S is assumed zero.The calculation is done as a function of the distance α., for the geometry shown in Fig. 1c, assuming 200MHz, background a0 actual pre-ICG and post-ICG measurements, and the multiplicative term

)
the arguments that led on the elimination of S from Eq.16 cannot be applied to this term since |

Figure 3 .
Figure 3. (a) T1-weighted MR coronal slice of a human breast, (b) Gd-DTPA distribution (in color) of the same coronal slice.A ductal carcinoma appears in yellow.

′
=0.0589 cm -1 so that average absorption increase due to the ICG is

Figure 4 . 5 .
Figure 4. Absorption maps used for the simulation.(a) pre-ICG breast (b) post-ICG breast.The maps of figure 3 served as an input to a finite-differences implementation of the frequency-domain diffusion equation.The simulation assumed 7 sources and 21 detectors as

Figure 5 .
Figure 5. Sources and detector arrangement used for the simulation.The region reconstructed is outlined with a green double line.
cases examined are shown in Figure6.

Figure 6 .
Figure 6.Reconstruction results imaging the region indicated with the green double line on Fig.5 assuming an infinite slab geometry.a)The reconstructed image using Eq.20.b)The reconstructed image using Eq.22.c)The reconstructed image using Eq.23.d) The result of subtracting an image of the post-ICG breast (reconstructed relative to a homogeneous baseline medium) from an image of the pre-ICG breast (reconstructed relatively to the same baseline medium).The optical properties of the homogeneous medium were 0 a µ ′ =0.03cm -1 and s µ′ = 8cm -1 .

Figure A1 .
Figure A1.(a) Amplitude and (b) Phase of R as a function of the average absorption coefficient of the post-ICG breast, plotted for the geometry of Fig 1c.

TABLE I
Absolute optical properties of the different structures used for the simulations.