In diffuse optical tomography (DOT), light is employed to interrogate soft tissue organs to recover their optical absorption and scattering coefficient distributions. These optical properties carry useful diagnostic information to detect pathological conditions such as cancer.1, 2, 3, 4 Since light diffuses in a highly scattering medium like tissue, spatial resolution in the reconstructed images is poor. To get around the difficulty of insufficient localization of optical properties, a dual-model approach has been adopted, bringing in another interrogating radiation, the ultrasound. Since ultrasound suffers comparatively low scattering in tissue, it can be focused to narrow regions to tag photons traversing through the focal region. These ultrasound-tagged photons are used to recover a qualitative image of the insonified region. This branch of dual-mode imaging is known in literature as ultrasound-assisted optical tomography (UAOT).5, 6, 7
To move from a qualitative image to a quantitative recovery of optical properties, models bringing out the interaction of ultrasound with the tissue and of light with insonified regions are required. These models are discussed in a number of past publications.6, 7, 8 The main points brought out by these investigations are that the ultrasound brings about 1. a modulation in the refractive index of the insonified tissue material, and 2. a periodic oscillation of the particles of the tissue. The oscillation and the refractive index modulation impart phase modulations to the interacting light. Sakadzic and Wang7 have considered the light interaction model for the general case when the tissue has a scattering anisotropy. They have derived expressions for the amplitude autocorrelation of the light exiting the object, which show forth a periodic modulation, owing to the periodic phase modulation picked up by the light during its interaction with the insonified region. The measurement carried out on the boundary is the intensity autocorrelation , which is related to . It has been shown that the modulation depth in is related to the local absorption coefficient of the ultrasound-modulated region, and this modulation depth can be used to get a qualitative map of this distribution. (In principle, it can be used to recover quantitative information on absorption coefficients as well.) It has been pointed out recently9 that the UAOT readout [i.e., the modulation depth in ] is affected not only by the optical absorption coefficient in the region but also by its viscoelastic properties.
Since the measurement is affected by both mechanical and optical properties, the question arises as to how one can separate out the contributions from amplitude of vibration (which depends on mechanical property) and the optical absorption coefficient. In this context it was shown in a recent study that the direction of displacement of the vibrating particles in the focal region of an ultrasound transducer is almost along its axis.10 For moderately scattering objects [typically the (poly) vinyl alcohol (PVA) phantom used in our experiments, for which the reduced scattering coefficient varies from ], the transport mean path is large enough that light reaching the insonified region does not become fully diffuse. In such cases, if the interrogating light is launched perpendicular to the ultrasound transducer axis, the contribution of displacement to phase modulation is small, which means that the effect of variation of the local storage modulus on the modulation depth is minimal (but not zero). The possibility of quantifying and eliminating this residual contribution from the overall modulation depth is also mentioned in that work.10 The major fall-out of this observation is that by using transverse interrogation geometry, it is possible to ensure that the contribution to modulation is more or less owing to variation. This paved the way for a possible quantitative recovery of from the measured modulation depth in the case of moderately scattering objects. When the reduced scattering coefficient increases beyond , which is typical for human tissue, this ability to discriminate the properties based on the direction of illumination is lost. As shown in the works of Sakadzic and Wang,11, 12 the modulation depth is affected not by the component of the vector displacement but by the mean square of its magnitude over the transport mean free path of light. This is similar to an expression for derived in the context of diffusing wave spectroscopy (DWS) for light propagated through a multiple scattering system driven by temperature fluctuation, wherein the field autocorrelation is affected by the mean-square displacement and not the vector displacement. Therefore, in practical cases of interest where scattering is high, we need to resort to other means to effect a separation of these properties from measurement. In the rest of this work, we deal with relatively low-scattering objects, of which the PVA phantom of our experiment is a typical representative.
For moderately scattering objects such as the ones mentioned before, since there is a distinct possibility of recovery from transverse interrogation, it would be tempting to see if one can complement this by recovering the storage modulus [uncorrupted by the local variation of ] from the measured axial interrogation modulation depth. Since the optical property is isotropic, it influences the axial as well as transverse modulation depths equally. Therefore it is essential, for quantification of the displacement contribution, to separate from the axial modulation depth the contribution from . This is possible, for we already have this contribution from the transverse interrogation modulation depth, which, after proper adjustment for the larger interaction length in the axial case, can be used to arrive at the contribution of . However, if this correction is not made, the recovered displacement component and the storage modulus will be erroneous. For example, we show in this work that a mechanical inhomogeneity with a low optical absorption coefficient can give the same axial modulation depth as an optical inhomogeneity with a low storage modulus. In the present work we have not yet demonstrated the separation of optical and mechanical properties from the measured axial and transverse modulation depths. However, through simulations and experiments, we show the major differences in the information carried by the two interrogations. For example, we show that for a storage modulus inhomogeneity, the contrast in an image constructed with the transverse modulation depth is very poor compared to that from an axial interrogation. For inhomogeneity, the transverse illumination contrast is as large as the one from the axial illumination. We have also shown the experimentally measured modulation depth variations from axial as well as transverse measurements from objects with mechanical and optical inhomogeneities. We hope these experiments provide data that delineate the route to be followed in future experiments aimed at separate recovery of these properties through ultrasound assisted optical elastography (UAOE). A brief summary of the work is as follows. In Sec. 2 we arrive at analytical expressions for the expected contributions to phase modulation from displacement and refractive index modulations. These are used to derive expressions for . Alternatively, Monte Carlo simulation is employed to take photons through an insonified tissue-mimicking object. The accumulated average phase is used to compute . The modulation depth in [and ] is computed when the storage modulus and in the insonified region are varied. Section 3 describes the experiments done on tissue-mimicking phantoms, made of PVA gel, to verify the theoretical simulations of Sec. 2. The experimentally obtained modulation depths are compared with the corresponding simulation results. Here we also discuss the case of an object with a much larger scattering coefficient , where our method of separation of the properties based on varying the illumination direction fails. Section 4 discusses the results and also puts forth the concluding remarks.
Pressure at the Focal Region of an Ultrasound Transducer
Our objective here is to compute the force exerted on the tissue because of the ultrasound insonification in the focal region. Force can be computed from the average intensity distribution of the ultrasound as , where is the acoustic absorption coefficient of the material of the object and is the acoustic velocity . The time average acoustic intensity can be obtained from the pressure distribution in the ultrasound focal region using the relation . Here, is the density of tissue. In this section we calculate the pressure distribution in the focal region of the transducer by solving the Westervelt equation, which models the propagation of sound in the medium.13 The Westervelt equation is given byis the velocity of sound in the medium, is the sound diffusivity, and is a constant related to the material nonlinearity. The parameter is usually expressed as , where is a factor related to nonlinearities in media. is taken as 6.2 for breast tissue.14 Equation 1 is solved numerically using the scheme given by Kamakura, Ishiwata, and Matsuda.13
In the simulations here and elsewhere in this work, we assume the object properties and dimensions mimic those used in the experiments described in Sec. 3. The characteristics of the ultrasound (US) transducer are also assumed to be identical to those of the transducer employed in the experiments. The acoustic properties of the object are: 1. for the background material the average sound , sound absorption at , density , and acoustic ; and 2. for the inclusion the average sound , sound absorption at density, , and acoustic .
The details of the US transducer are as follows. It is a continuous-wave-type transducer operating at a central frequency of . It has a focal length of , and aperture radius of , giving a full aperture angle of . We used a custom-made transducer (manufactured by Roop Telsonic Ultrasonix Limited, Mumbai, India) with a hole at the center for allowing axial illumination of its focal region. The coupling coefficient for the transducer is 0.622. The electrical input to the transducer is adjusted such that the pressure at the transducer surface is .
With these parameters for the transducer and the phantom material, we have solved the Westervelet equation for pressure in the focal region of the transducer. The pressure distribution is found to be ellipsoid in shape with an eccentricity of . The typical value of the pressure in the center of the ellipsoidal region is found to be approximately . The cross sectional plots (both radial and axial) of the pressure through the center of the ellipsoidal region are shown in Figs. 1a and 1b . The pressure distribution is converted to average intensity, which is used to compute the force exerted by the US transducer in its focal region.
Computation of Displacement in the Focal Region of the Ultrasound Transducer
The force distribution computed in the previous section is used here to arrive at the displacements suffered by the tissue particles. For this, the equilibrium equations of elasticity are set up and solved. The governing equation for tissue displacement loaded by a sinusoidal force along the transducer axis (taken as axis here) is given by10and are the Lame parameters, and the commas in the subscripts denote partial differentiation with respect to the indices that followed. Summation is implied with respect to the repeated indices.
A discretized version of Eq. 2 is solved using Ansys (Ansys®, Incorporated, Canonsburg, Pennsylvania). We considered a object, which is meshed by 29,791 nodes to create 27,000 equal volume cubic elements. The location of the focal point of the US transducer is selected to coincide with the center of the object. A typical solution of Eq. 2 for , with force from the transducer used in the experiments and representative values of elastic modulus , Poisson’s ratio , and for the phantom, is given by , , and . We point out that the displacement is almost along the US transducer axis. With the values of , , and suitably modified, the algorithm is used to compute the displacement in the inclusions as well.
Interrogation of the Insonified Object With light
We employ a coherent light beam to interrogate the focal region of the US transducer. The force distribution in the focal region results in periodic movement of the tissue particles as well as a periodic modulation of the refractive index . These two together result in a periodic phase modulation of light passing through the region, which results in a periodic modulation in , the amplitude autocorrelation of the exiting light. This is observed in the measured intensity autocorrelation . The force exerted in the focal region of the transducer makes a typical (say the ’th) optical path length into , where is the component of the time-dependent displacement of the ’th scatterer in the direction of , and is the time-dependent refractive index modulation. The path length modulation introduced by the US force is approximately equal to , which gives a phase modulation , the contributions from refractive index modulation and displacement, respectively. Denoting , the field autocorrelation of the detected light at the boundary corresponding to a particular path can be derived.7 This autocorrelation is given by
Hereand are derived in Ref. 7, which are made use of to compute is the probability density function for path of the photon. The periodic variations of and contribute to a periodic modulation on , whose depth is affected by both and . This modulation can be thought of as a carrier, which is amplitude modulated by the optical absorption coefficient of the insonified region. Considering that does not vary much within the object, the depth of modulation is affected primarily by the variation in and . If can be evaluated [for geometrically well-defined objects like the slab used in the experiments, analytical expressions exist for ], one can derive an expression for by substituting this in Eq. 5. From , can be obtained using the Siegert relation15 . Here is a constant related to the collection optics.
The field autocorrelation function can also be computed using Monte-Carlo (MC) simulation to launch and track a large number of photons through the insonified tissue-mimicking object. The procedure for doing this is explained in Ref. 7, and is adapted successfully in Ref. 10, to quantify the contributions of and in . We make use of both the analytical expression for and its MC simulated equivalent to compute the modulation depth of in the case of light transmitted through a typical tissue-mimicking object. Light is to intercept the insonified region, and in the simulations both and storage modulus of this region are varied, and the power spectral density of at the insonified frequency is computed. Our objective is to study variation in with and when the illumination is either along or perpendicular to the US transducer axis.
In keeping with the properties of the phantom used in the experiment, whose elastic properties cannot be controlled independent of its scattering properties, we varied storage modulus from , which automatically ties its reduced scattering coefficient to vary from . With this range of moderate scattering values, we have computed modulation depths for both axial and transverse scans through an object with storage modulus inhomogeneity, and another with an absorption coefficient inhomogeneity. The background properties of both the object are and . The first object had a storage modulus inhomogeneity of and the second an absorption coefficient inhomogeneity of . We have observed that the contrast in modulation depth computed from the axial illumination was larger compared to the transverse illumination for both the objects. There are two possible reasons for this. 1. The amplitude of vibration of the tissue particles has a large contribution to the phase modulation and the consequent modulation. Since the direction of vibration is almost axial, it has a large component in the mean direction of light transport, assuming that the light propagation has not yet fully become diffuse in the insonified region. Therefore a change in will result in a large change in the phase modulation and the consequent modulation.2. Since the insonified region is an ellipsoid with a large aspect ratio favoring the axial direction, interaction length is larger for light launched along the axis; this alone should explain the larger contrast from axial illumination for the object with absorption inhomogeneity. When the background reduced scattering coefficient is increased to , the difference between the computed modulation depths for axial and transverse illuminations became much less pronounced for the same contrast in . Interestingly, when the interaction lengths along these directions also were made equal, the difference in modulation depths became almost negligible, indicating that light propagation has already become diffuse.
This is the summary of the simulations done. These results are also verified through experiments. In addition, we have conducted experiments on composite phantoms with both absorption and storage modulus inhomogeneities. The results, which bring forth the need to quantify and separate the contributions to modulation from and en route to a quantitative recovery of these properties, are discussed in the sections on experiments and results.
Fabrication of the Phantom
The phantoms used in the experiments are slabs of dimensions with inclusions of size . Both the background material and the inclusions are made of PVA gel, whose mechanical, optical, and acoustic properties can be tailored to mimic breast tissue in normal conditions as well as in disease. The way the PVA gel is formed and its properties tailored are described in Ref. 16. Altogether we have employed four phantoms in the experiments. All of them have a background storage modulus and absorption coefficient of and , respectively. The first object has an inclusion of with the same as that of the background. The second has an absorption inhomogeneity of with in the inhomogeneity the same as that of the background. The third object has two inclusions, one a storage modulus inhomogeneity of value and the other an absorption inhomogeneity of value . The fourth object also has two inclusions, the first one a storage modulus inhomogeneity of (with absorption in it at the background level) and the other an inhomogeneity in both and with values of and , respectively.
The elastic property of the phantom is tailored by varying the number of freezing-thawing cycles the PVA gel is subjected to, which also alters the physical cross-linking of the polymer, which in turn controls the scattering coefficient of the gel. When the number of freezing thawing cycles was two, was found to be and . When this number increased to six, and became and , respectively. The absorption coefficient in the inhomogeneities is tailored by adding appropriate quantities of India ink while preparing the gel. The anisotropy factor of the PVA gel is found to be around 0.89, irrespective of the number of freezing-thawing cycles the gel underwent.
When the mechanical properties are tailored by adjusting the number of freezing-thawing cycles, the acoustic properties also varied. The acoustic properties of the PVA phantom after two freezing-thawing cycles (with ) are: 1. average acoustic , 2. acoustic , 3. , and 4. sound absorption at . The corresponding values after six freezing-thawing cycles (with ) are: 1. average acoustic , 2. acoustic , 3. , and 4. sound absorption at .
Intensity Autocorrelation Measurement
The setup in Fig. 2 is used to measure the intensity autocorrelation.10 PVA samples in the form of rectangular slabs are fabricated with rectangular inclusions, which had either a higher storage modulus or optical absorption coefficient or both when compared to the background. The typical values measured for the fabricated phantoms are: is for the background, and for inclusions it is , which were also used in earlier simulations. To match the objects used in the simulations, we have fabricated objects: 1. with either a or a inhomogeneity; 2. with two inclusions, one of and another of ; and 3. one which had inclusions with higher and .
The object, held in position in a water-filled glass tank (water provides an appropriate coupling medium for delivering the US beam into the object) (Fig. 2), is insonified by radiation from a focusing US transducer. The tissue particles in the focal region are set in vibratory motion by the force of the US beam, which, as pointed out earlier, is along the transducer axis. The focal region is interrogated by a laser beam sent along the axis of the US transducer, through its central hole, with a view to pick up a component of the amplitude of vibration [Fig. 2a]. The transmitted light is detected by a single-mode fiber, which ensures a high signal-to-noise ratio in the detected intensity by picking up a single speckle for the photon counting detector. The photon counting system is a single unit with a photomultiplier tube (PMT) and a pulse amplifier discriminator (PAD) (Hamamatsu H7360-03). The output from the photon counting unit is given to an intensity correlator (Flex 021d, see http://www.correlator.com), which helps us get the intensity autocorrelation of the exiting photons from the phantom. The correlator output is Fourier transformed, and the power spectral amplitude at is measured, which is proportional to the modulation depth in .
The phantom, which was mounted on an translation stage, was moved so that the midpoint of the US focal region scanned a transverse plane of the slab. Thus a point-by-point scanning was done for different positions of the slab. For each position of the slab as well as its depth of modulation are measured The latter is smaller for the inhomogeneous inclusion (for which and are larger) and larger for the background. For transverse geometry, the photon packets are injected in a direction perpendicular to the US axis, as shown in Fig. 2b. The experiments described are repeated with the transverse illumination as well.
Modulation depths in the measured with focal regions of the US transducer scanning across cross sections of these objects are obtained and plotted. These are shown in Figs. 3 and 4 . These results are further discussed in the following section.
Results and Discussions
Figure 3a shows the variation of the measured power spectrum amplitude as the focal region of the US transducer scans across the first object. Here the simulation results are compared with the experimental values, which show a close match. Because the object is only moderately scattering, the experimentally measured has reasonably high signal-to-noise ratio, which partly explains the good match between the simulation and experimental results. It is seen that the contrast from the axial data is much larger compared to that from the transverse data. This is only to be expected, since the inhomogeneity is in the storage modulus, which alters only the axial component of displacement. The small contrast discernable in the transverse plot is because of the light diffusion and the nonzero component of the displacement in the direction of some of the diffuse photon paths.
Figure 3b is the experimentally measured modulation depth variation for the second object, which has an absorption inhomogeneity. (There is no comparison here with simulation results.) The axial scan still shows forth a larger contrast in comparison to the transverse scan. This is owing to the larger interaction length of the photons along the axis of the ellipsoidal focal region. We see from this result that the axial scan provides large contrast in measurements for an object with only the absorption coefficient variation.
The experimental result shown in Fig. 4a is for the third object, which has two inclusions. We see that the axial measurements do show a contrast between an absorption inhomogeneity of in an background, and a storage modulus inhomogeneity of in a background absorption coefficient, which we are not in a position to ascribe to either the or the change. In fact, an increase in in a low absorption coefficient background produces a lowering of similar to that from an increase in in a low background storage modulus, so that a inhomogeneity in this case is almost indistinguishable from a inhomogeneity. It is interesting to note that the inclusion with higher absorption coefficient (in a low storage modulus background) appears with a negative contrast. Assuming, as is done in UAOT experiments, that contrast is only due to variation, negative contrast implies the showing forth of a higher corresponding to a larger than for a lower . However, the transverse scan can distinguish the absorption inhomogeneity with the right contrast. This clearly shows that the contrast in the axial scan is predominantly due to displacement caused by the US force (which has resulted in the negative contrast in this experiment), whereas in the transverse scan, it is due to the absorption coefficient.
The last result shown in Fig. 4b is for the fourth object. The object has an inclusion with a higher storage modulus and absorption coefficient, which represents the change in a cancerous inclusion in a practical scenario, where the stiffness as well as optical absorption increases with disease. The second inclusion, which is used for the sake of comparison, has only a storage modulus inhomogeneity. We see once again that the axial scan fails to bring out a large enough contrast between the two inclusions. The larger contrast in the first inclusion (compared to the background) is primarily owing to the storage modulus change. In our phantom experiment, from the difference in contrast in the axial scan at the two inclusions, one can possibly quantify the contribution of the absorption coefficient. Once the contribution of is accounted for, the information from the axial scan can be used to recover the storage modulus. In a similar way, the transverse scan can be used to recover the absorption coefficient. For detection of a cancerous lesion, which is distinguished both by its increased optical absorption and storage modulus compared to the background, the axial scan provides a much larger contrast distinguishing the lesion from the healthy background tissue.
Finally, we have simulated and its modulation depth for an object with using both axial and transverse scans. In the simulations, we have considered an object with a inhomogeneity (background of and of inhomogenous inclusion of ). The results in Fig. 5 show that the contrast, which is evident in earlier simulations and experiments (for example, Fig. 3) using lower scattering coefficient objects, has become much smaller. This is attributed to the spreading of light because of diffusion and the consequent randomization of propagation. This is evident from comparing curves, A and C of Fig. 5. It is also seen that the difference is not negligible, especially as seen from the modulation depths from the background region. This could be because of the difference in interaction lengths of the axial and transverse light with the focal region of the US transducer, which, as mentioned earlier, has an ellipsoidal shape. Also in the background, the stiffness is smaller, which means that the displacements of the scattering centers are larger, introducing a pronounced difference in phase modulation through the interaction length difference between the two scans. Curve B in Fig. 5, which is from a transverse scan with the modification that the US focal region has been constrained to be spherical, is almost indistinguishable from curve A of the axial scan. This last set of simulation results makes it clear that for a highly scattering medium like tissue, discrimination based on illumination direction to separate contributions from optical and mechanical properties is ineffective.
We demonstrate through experiments the capability of UAOE to image objects with mechanical as well as optical property inclusions. For objects that only moderately scatter (i.e., with less than ), we show that the axial scan is predominantly sensitive to storage modulus variations, whereas the transverse scan is to optical absorption coefficients. The experiments reported here are meant to be the preliminary investigations en route to finally employing UAOE to map the mechanical and optical properties of such objects quantitatively. It is clearly shown that the axial and transverse data gathered from UAOT experiments require compensation to recover either the mechanical or optical property without one influencing the other.
We also show through simulation that the present method fails when the object becomes highly scattering, say becomes greater than . If a computation-friendly Monte Carlo simulation (such as the perturbation Monte Carlo) can help us calculate the derivative of the measurement with respect to the absorption coefficient and storage modulus in the focal region of the US transducer, then the axial data itself can be used to reconstruct the optical and mechanical properties separately. This is quite feasible, which will render an axial UAOT scan capable of reconstructing both the mechanical and optical properties.