Photoacoustic temperature imaging based on multi-wavelength excitation.

Building further upon the high spatial resolution offered by ultrasonic imaging and the high optical contrast yielded by laser excitation of photoacoustic imaging, and exploiting the temperature dependence of photoacoustic signal amplitudes, this paper addresses the question whether the rich information given by multispectral optoacoustic tomography (MSOT) allows to obtain 3D temperature images. Numerical simulations and experimental results are reported on agarose phantoms containing gold nanoparticles and the effects of shadowing, reconstruction flaws, etc. on the accuracy are determined.

Introduction

In the framework of cancer treatment, hyperthermia therapy can improve the therapeutic outcome of radiotherapy and chemotherapy by synergistic combination of heating and radiation based attack of malignant cells and antineoplastic drug cytotoxicity. Meanwhile, through modulation of immune activity, hyperthermia can elicit a cellular immune response [[1], [2], [3]]. If the temperature increase is pushed high enough, as in thermal ablation approaches, then the heating can be used to destroy tumor cells. In order to tune the heating dosage and optimize the outcome of the treatment while minimizing side effects, non-invasive, real-time monitoring of temperature with high contrast and resolution during therapy is of essential importance.

Until now, many techniques used for anatomical imaging have been investigated for their potential extension to real-time temperature imaging with high resolution and sensitivity. Thermocouples are commonly used in temperature measurement in view of their low cost and ease of use. However, direct contact with the tissue is necessary, making this approach quite invasive [4]. Infrared thermography allows real-time and high-resolution temperature measurement, but it only provides superficial detection [5]. Ultrasound (US) imaging can penetrate tissue up to several centimeters depth because long wavelength acoustic waves are only slightly affected by scattering. However, US image contrast is determined by the strength of gradients in the mechanical properties of tissues [6], which only weakly depend on temperature. Attempts to exploit the temperature dependence of the speed of sound have been explored, but it is quite challenging to extract the very small related effects from US images with high spatial resolution [7]. MRI offers the advantage of high sensitivity and spatial resolution, but it is rather time-consuming, which puts some constraints on its use for real-time monitoring [8]. Moreover, MRI remains a costly method, which limits its use to a limited number of large scale clinical facilities.

Photoacoustic imaging (PAI), also called optoacoustic imaging, was initially considered to overcome the diffusion limit encountered in pure optical imaging. This emerging method combines the strong contrast of optical imaging with the high spatial resolution and deep penetration of US imaging [[9], [10], [11], [12], [13], [14], [15], [16], [17]]. It has been demonstrated that PAI is able to provide functional imaging [[18], [19], [20], [21], [22], [23], [24], [25]] through excitation of hemoglobin and DNA&RNA. Moreover, with the help of exogenous contrast agents, PAI can even target reporter genes and biomarkers to achieve molecular imaging [[26], [27], [28], [29], [30], [31]].

PAI has also been used for temperature measurements [[32], [33], [34], [35], [36], [37], [38], [39]]. By exploiting the roughly linear temperature dependence of the Grüneisen parameter, which expresses the thermomechanical conversion efficiency and thus the photoacoustic (PA) signal magnitude, PA amplitude changes are associated with temperature-induced Grüneisen parameter variations, and thus with variations of the tissue temperature. Although recently quite promising results have been obtained using this approach [[40], [41], [42], [43]], it also has some limitations. Beyond 50 °C, due to protein denaturation and coagulation, the linear relation between the Grüneisen parameter and temperature does not hold anymore and tissue optical properties are altered [43]. Also, the Grüneisen parameter is tissue dependent, so that a calibration on one tissue cannot be applied to a different one. In principle, tissues of different patients can also exhibit different Grüneisen parameters, in spite of being of the same type. In addition, the local laser light fluence varies with the depth and position probed, leading to spatial heterogeneities in energy deposition. As a consequence, the absolute magnitude of the PA signal cannot be used to extract temperature information, and acquisition of a reference PA image at a known temperature is necessary. In summary, further steps are needed to establish an efficient and cost-effective absolute temperature monitoring for deep tissue.

In this article, we investigate to what extent the unique possibilities MSOT can be exploited to achieve temperature monitoring, without the above-mentioned limitations. We exploit the difference in temperature dependence of the optical absorption coefficient at different optical wavelengths, rather than making use of the temperature dependence of the Grüneisen coefficient. In this way, provided that sufficient knowledge on this dependence is available, differences between PAI images obtained at different optical excitation wavelengths can be translated into temperature images. In case of uniform illumination, such relative images are independent of the absolute magnitude of the PA signal, and thus independent of the concentration of chromophore and laser light fluence.

In the following, we verify to which extent MSOT can be applied to temperature imaging without the need of a reference image, and how accurately the resulting temperature images compare to the ones obtained using single wavelength PAI.

A numerical analysis is performed in which finite difference (FD) simulations are used to generate artificial ultrasonic signals in a tomographic transducer array. FD is also used to reconstruct PA source images from transducer signals. All simulations are done for different optical wavelengths, assuming temperature dependences of the optical absorption coefficient at different wavelengths comparable to that observed in experimental data.

MSOT experiments in the wavelength range 670 nm–1255 nm are carried out on an agarose phantom containing gold nanoparticles. PA signals from a quasi-circular array of 256 transducers are used to reconstruct PA images at different optical wavelengths, and therefrom temperature images. The actual accuracy of the method is assessed, the reasons for its current limitations are sought, and the possibilities for future improvements are suggested.

Temperature and wavelength dependence of the PA signal

When matter is illuminated by laser light, a part of the energy is converted to heat by optical absorption and followed by relaxation via non-radiative pathways. The sudden thermal expansion resulting from pulsed light illumination generates acoustic waves: this is referred to as the PA effect. Typically, the excited initial PA pressure signal, p [Pa] can be expressed aswhere η [-] is the fraction of absorbed light that is converted to heat (the rest is re-emitted as fluorescent/phosphorescent light), [m−1] is the optical absorption at wavelength [nm] and temperature T [°C], and [J. m-2] is the local laser fluence at position (x,y,z). Γ(T)=β/(кρCv) is the dimensionless Grüneisen parameter, which, for water has been shown to vary with temperature T [°C], with a variation of about 2.7%/°C with respect to 36 °C [43], i.e., body temperature:with β [oC−1] the volume thermal expansion coefficient, к [Pa−1] the isothermal compressibility, ρ [kg/m3] the density and Cv [J. kg−1.oC−1] the specific heat capacity at constant volume.

The proportionality of the PA signal with the Grüneisen parameter has been exploited to extract temperature information from the PA signal magnitude, by extrapolating the validity of Eq. (2) to body tissue. However, as stated above, protein denaturation and coagulation above 50 °C affect Eq. (2), absorption and scattering coefficients, thus limiting the accuracy of the temperature monitoring techniques beyond this point.

When multiple signals are acquired at multiple laser excitation wavelengths (e.g. two wavelengths, λ1 and λ2), assuming that the quantum yield and the thermoelastic conversion efficiency are equal for both acquisitions, the ratio between the signals is then determined by the ratio of the respective optical absorption coefficients and the local light fluence:

If the ratio of the local light fluence at two wavelengths is independent of temperature, and if an exogeneous chromophore is sufficiently encapsulated so as not to be chemically influenced by its biological environment, while still sensing the temperature of the surrounding tissue, then the temperature of the PA ratio depends (quasi) only on the temperature dependence of the local absorption characteristics of the chromophore, allowing for temperature measurement with a simple calibration without the aforementioned drawbacks. In case of PA imaging, if the above conditions are fulfilled, then Eq. (3) can be used for every location (x,y,z), so that the temperature image can be calculated from the ratio of the PA images obtained at the different wavelengths. Moreover, the more images at different wavelengths are available, the more effective the wavelength and temperature dependence of PA images can be exploited, leading in turn to a more accurate and robust reconstruction of temperature images. In the following, we assess the feasibility of the proposed temperature imaging method. Effects of wavelength dependent shadowing of the fluence in deep regions due to temperature and wavelength dependent optical absorption in shallow regions (so that the fluence ratio in Eq. (3) on deep regions becomes affected by the shadowing by more shallow regions) are addressed. For the sake of clarity, we first present a simulation-based analysis, which allows to illustrate the essentials of the method and of shadowing related artefacts and ways to deal with them in a straightforward and systematic way.

Description and optical characterization of model sample

The goal of this work is to verify by means of numerical simulations and an experimental case study, to which extent the feasibility of the proposed approach is confirmed. We also assess to which extent deviations from the assumptions enabling Eq. (3) affect the temperature reconstruction accuracy on ideal, simulated data, and on experimental data. For the sake of clarity and simplicity, a cylindrical agarose gel phantom (1 wt% Agar from Sigma-Aldrich® in ultrapure water), with a diameter of 24 mm, encapsulated in a 4 mm thick polydimethylsiloxane (PDMS) cylinder. The phantom size is similar to the adult mouse main body cross-section. PDMS was chosen because of its good acoustic impedance match with water (ZPDMS = 1.08 MRayl, Zwater = 1.57 MRayl), which reduces acoustic reflections at the interfaces with the enclosed hydrogel and the surrounding water bath, as well as reducing the alterations the water would have on the gel. Gold nanorods (GNRs) were used as exogeneous chromophores, in order to increase the optical absorption and thus PA excitation efficiency in a short wavelength range, and to enrich the spectral variability of the temperature dependence of the optical absorption coefficient, and thus maximize the temperature sensitivity of the PA ratio in Eq. (3). GNRs were synthesized following a seed-mediated protocol. All syntheses were carried out using aqueous solutions. A gold seed solution was prepared by vigorous mixing of 4.7 mL of 0.1 M CTAB solution and 25 μL of 50 mM HAuCl4 with 300 μL of 10 mM NaBH4 at 30 °C. The seed solution was kept at 25 °C for 2 h prior to further use. The GNR growth solution was prepared by mixing 45 mL of 0.2 M CTAB with 9 mL of 5 mM HAuCl4 solution and 112.5 μL of 0.1 M AgNO3. After 10 min of stirring, 5.25 mL of 10 mM ascorbic acid was added to the solution. This was followed by an additional 10 min stirring before addition of 21 μL of the seed solution under vigorous stirring. The resulting solution was aged for 12 h at 25 °C, before purification of the nanorods by three rounds of centrifugation at 8000 rpm and re-dispersion in 5 mM CTAB.

Four hundred milligrams of agarose together with 40 mL of ultrapure water were added to a sealed vial and heated in a 90 °C water bath until complete agarose dissolution. After this, 200 μL of concentrated GNR solution were mixed with 40 mL of the hot agarose solution and this suspension was used to completely fill the PDMS cylinder. The cylinder was then left to cool in a 20 °C water bath for 30 min to allow complete gelification.

According to Eq. (1), at a given temperature and light fluence, taking into account that the quantum yield of gel and GNRs is 0% so that η = 1, the PA signal at a given wavelength is simply proportional to the optical absorption coefficient. In view of this, anticipating the MSOT experiment, we have measured the optical absorption of the involved materials, the GNR - agarose sample material using a simple optical transmission setup. Measurements were performed at different temperatures by combining a white light source, a temperature controller and an Ocean Optics USB4000® spectrometer as shown in Fig. 1. In the experiment, the light from an Ocean Optics HL2000® lamp exiting a broadband emission fiber was collimated by a lens and sent through a 1 × 1 cm2 Sevenlight® quartz cuvette that contained the sample material. The cuvette was placed in a copper holder that contained a platinum resistor and a resistive heater for monitoring and controlling the temperature with 0.1 °C resolution, by making use of a Red Lion TCU® temperature control unit, which was operated by a homemade National Instruments Labview® program with resolution of 0.1 °C. The program synchronized the acquisition of the spectrum of the transmitted light Itransmitted(λ,T) with the temperature control. Additional empty cuvette measurement Iempty(λ,T = 22 °C) was performed to correct for the absorption and reflection of the empty cell, and to determine the incident light power. The optical absorption coefficient μa [m−1] of the sample was then determined as:with d = 10 mm, the optical path length of the light transmitted through the sample.

Fig. 1

Schematic representation of the optical transmission setup used to determine the temperature dependence of the optical absorption spectrum of the sample.

Fig. 2(a) shows the temperature-dependent absorption spectra of the GNRs – agarose phantom calculated based on Eq. (4) from 25 to 40 °C with a step of 1 °C. The peak at ∼700 nm arises from the longitudinal surface plasmon resonance (LSPR) of the GNRs while the second peak results from water absorption. The LSPR peak is very stable and shows almost no variation with temperature, while the absorption of water monotonically increases with increasing temperature. The ratio R between the absorption coefficients 965 nm (∼water peak) and 700 nm (∼Au peak) increases linearly with temperature (Fig. 2(b)). The line fitted through R(T) has a slope of 0.005 and intercept of 0.43. A sensitivity s = 0.82%/°C can be determined with respect to ratio at 36 °C.

Fig. 2

(a) Optical absorption spectra of the GNRs phantom sample measured at different temperatures (varied between 25 and 40 °C in steps of 1 °C). A larger variation is observed around 965 nm for water molecules than for the longitudinal surface plasmon resonance peak of gold nanorods at 700 nm. (b) Ratio between absorption coefficients at 965 nm and 700 nm as a function of temperature.

If all conditions validating Eq. (3) are satisfied, this implies that the temperature of this material can be mapped by taking the ratio of the PA signal values at the two wavelengths. The expected temperature reconstruction uncertainty for an image reconstruction accuracy of say Δp = 1% is then expected to be  Δp/s = 1.2 °C. This implies that, among other aspects discussed further in this paper, the future clinical implementation of temperature tomography is mainly determined by the magnitude of the difference in temperature sensitivity at different experimentally accessible wavelengths.

The magnitudes of the optical absorption coefficients, between 0.1 cm−1 and 0.4 cm−1 are such that sufficiently large PA signals can be expected, while the signal to noise ratio (SNR) is not hampered by reduction of the intensity of the excitation light reaching central regions, due to shadowing by more eccentric regions.

Numerical simulation based analysis of the feasibility of PA thermotomography

Combining Eq. (3) with the linear temperature dependence of the ratio of the PA signals at λ1 = 700 nm and λ2 = 965 nm, in ideal circumstances, temperature images could be obtained as

In Eq. (5), the human body temperature 36 °C was chosen ad hoc as reference temperature, and s = 0.82%/°C, is the temperature coefficient evaluated at 36 °C. In the following, by performing numerical simulations, we systematically verify to what extent real life, non-ideal conditions affect the validity of Eq. (5), and the accuracy of the resulting temperature map. The investigated influential factors are:

Image reconstruction errors

Non-equal and temperature dependent light fluence between the two PA signals

The numerical study was performed by considering a 2D gel disk with 20 mm diameter, with the optical characteristics of Fig. 2, surrounded by water, and based on a FD approach and time reversal (TR) image reconstruction.

Considering that PA amplitudes generated in vivo are relatively small and typically in the order of 10 kPa, wave behavior within the soft tissue can be described by three basic linear equations of motion, continuity and state,with [m/s], ρ [kg/m3], and p [Pa] respectively the photoacoustically induced change of the velocity vector, the density and the pressure. ρ is the ambient density, and c is the speed of sound accordingly. The speed of sound and density of sample are approximated by the ones of water. For the FD emulation, these equations were discretized (index i in the x-direction and j in the y-direction) as the following update equation for the discretized velocity and pressure field:

The photoacoustic source was modeled via the starting condition:

The size of the calculation domain was 70 × 70 mm2, with 1061 × 1061 pixels and the sample disk centered in the middle. The boundaries were modeled as perfectly matching layers (PMLs) of 7 mm thickness. The time step was taken to be Δt = 20 ns, thus respecting the stability condition Δt≤0.5min[Δx/c,Δy/c].

In the actual experiment, a MSOT small animal imaging system (inVision 256-TF; iThera Medical GmbH) was used [44]. Briefly, a tunable optical parametric oscillator (OPO) pumped by a Nd:YAG laser provided excitation pulses with a duration of 9 ns at wavelengths comprised between 660 nm and 1300 nm at a repetition rate of 10 Hz with a wavelength tuning speed of 10 ms and a peak pulse energy of 90 mJ at 720 nm. Ten arms of a fiber bundle provided even illumination of a ring-shaped light strip of approx. 8 mm width. For ultrasound detection, a ring detector, consisting of 256 toroidally focused ultrasound transducers with a center frequency of 5 MHz (60% bandwidth, filtered between 500 kHz to 7 MHz), organized in a concave array of 270-degree angular coverage and a radius of curvature of 4 cm, were used. The spatial resolution of images reconstructed by the inVision 256-TF is 150 μm. The sample was placed in the center of the transducer array, immersed in water, and water temperature was varied between, in the setup, 25 °C and 37 °C in steps of 1 °C.

In the forward numerical simulations, the laser illumination was assumed to originate from around the sample, in concordance with the experiment, and was set to have a uniform angular distribution. The local light fluence was modeled according to Lambert-Beer law, by setting the initial PA pressure distribution (t = 0) to decay exponentially within the sample diameter, as depicted in Fig. 3(a). For the sake of simplicity, optical absorption by the water surrounding the sample and light scattering were neglected. In accordance with the real life 256-element imaging ring array, throughout the FD-calculation, the pressure signals at the 256 transducer positions, pkl with k = 1…N the time index and l = 1…256 the transducer index, were recorded for later use in reconstruction calculations. As mentioned above, the absorption coefficients used in the FD calculations were taken from the temperature-dependent optical absorption measurement at 700 nm and 965 nm (Fig. 2).

Fig. 3

(a) Illustration of the assumed angle independent and radially exponential decaying distribution of PA-source strength at 25 °C and 700 nm (μa = 0.4295 cm−1) used for forward FD-calculations. (b) TR-FD-reconstructed source strength distribution extracted from 256 virtual transducer signals at the indicated locations (white circles) corresponding to Fig. 3(a). (c) and (d) Maps of the slope and intercept of the temperature dependence of the ratio of the TR-FD-reconstructed maps at 965 nm and 700 nm. Without shadowing effect, the slope and intercept values are 0.005 and 0.43 respectively. (e) and (f) Evolution of the spatial average of the reconstructed temperature map with the set temperature used in the forward calculations, with (e) implementation of local slope and local intercept and (f) mean slope over the sample region and a calibration based on a reference measurement at 36 °C.

A TR algorithm [45,46] was used to reconstruct images simulated in different scenarios. In short, the 256 signals recorded in the forward FD calculation were time reversed and implemented synchronously as time dependent pressure conditions at the respective transducer locations. Given the number of time steps used in the forward calculation, N, the TR-FD-calculated pressure pattern after N steps was taken as the PA source reconstruction. Each reconstruction (one per wavelength and temperature) was done on the basis of one laser shot, within the time between 2 laser shots (100 ms).

Fig. 3(b) shows the resulting TR-FD reconstructed image, based on the signals of 256 virtual transducers, recorded during the forward FD calculation with the PA source strength distribution p(x,y,t = 0) of Fig. 3(a). The virtual transducers are indicated as white point circles in Fig. 3(a) and (b). In view of determining the feasibility of MSOT-based temperature mapping, we have performed forward simulations at 16 temperatures (spaced by 1 °C) between 25 °C and 40 °C, and at 965 nm and 700 nm wavelength. For each forward calculation, virtual transducer signals were recorded and used for a TR-FD reconstruction. Next, the ratio of reconstructed PA source strength maps at 965 nm and 700 nm was then taken for each temperature. For each position, linear regression was done, yielding the intercept and slope of the best linear fit of R(T) the ratio versus temperature function. The fitted slope (Fig. 3(c)) and intercept values (Fig. 3(d)) (corresponding to the ratio at 0 °C) are position dependent. The R(T) slope values are distributed in the range of 0.0045-0.005, while the R(T) intercepts are spread from 0.43 to 0.55 as compared to the fitted slope 0.005 and intercept 0.43 shown in Fig. 2(b). In this case, the temperature coefficient is between 0.63-0.82%/°C at 36 °C, as compared to 0.82%/°C, the coefficient obtained through the fitted line in Fig. 2(b).

Looking at Eq.(5), the discrepancies between the TR-extracted temperature sensitivity (slope and intercept) values and the original ones are most likely due to a spatially non-uniform temperature dependence of the ratio between the light fluence values at the two wavelengths. This can be caused by deeper regions in the sample being partially shadowed by shallower regions, with an associated position dependent shadowing attenuation that is determined by Lambert-Beer law, and by the temperature dependences of the optical absorption coefficient at the two wavelengths.

Fig. 3(e) shows the reconstructed temperature, spatially averaged over the 20 mm diameter disk region. The temperature was obtained by applying the R(T) slope and intercept shown in Fig. 3(c) and (d). In spite of the location, temperature and wavelength dependence of the light fluence, the correspondence with the true temperature is very good, with maximum standard deviations smaller than 0.02 °C. For the considered angle independent geometry, the temperature and wavelength dependent shadowing effect can be estimated by the factor exp(-μa(T,λ1)(Rc-r))/exp(-μa(T,λ2)(Rc-r)), with Rc = 10 mm, the cylinder radius and r, the radial coordinate at the considered location. Correcting for this factor, the temperature reconstruction error can be pushed further down to below 0.002 °C. However, it should be mentioned that for real life conditions with non-uniform distributions of endogenous thermochromic entities with strong optical absorption in the wavelength range of interest, and for non-uniform distribution of the illumination, this kind of correction is not possible due to a lack of symmetry and a priori knowledge about those distributions.

As an alternative to the above approach, temperature maps can also be extracted by using the spatial average of the R(T) slope over the 20 mm diameter disk region in Fig. 3(c) and the image ratio R(Tref) at a known reference temperature Tref. Using Tref = 36 °C, the body temperature, would be a logical choice for temperature mapping to monitor hyperthermia treatment, before the start of the treatment of a patient. In this way, the intercept error is largely cancelled.

Implementing this procedure forces the spatially averaged temperature reconstruction error to zero at 36 °C. The error for other temperatures is found to be increasing with around 0.023 °C per degree distance from 36 °C (Fig. 3(f)).

Measurement based analysis of the feasibility of photoacoustic thermotomography

PA images were acquired at 120 wavelengths between 660 nm and 1255 nm and at 13 temperatures between 25 °C and 37 °C. For each wavelength, one image was acquired at a rate of 10 fps without averaging after the setting temperature was reached. The laser was then scanned to the next wavelength. Since the tuning of laser wavelength to any desired one in the range can be completed within less than 100 ms and does not affect 10 Hz laser trigger rate, multispectral images at 10 fps can be guaranteed. In practice, with promising nanoparticles and ratio obtained at 2 wavelengths at 200 ms, temperature reconstruction can be acquired within one second. In order to remove possible influences of laser intensity fluctuations, every image was normalized by the laser intensity, which was measured synchronously with the ultrasound signal acquisition.

Fig. 4(a) shows a PA source strength map on the Au nanoparticle dispersed PDMS contained agarose gel sample obtained by pulse illumination at 670 nm, recording of 256 transducer signals, and ViewMSOT v3.8 software (iThera Medical) image reconstruction. The temperature of the water surrounding the sample was maintained at 30 °C. The elliptic lines in the image are other non-uniformities originate mainly from the absence of transducers along a sector of 90 ° non-uniform illumination and from shadowing for central sample regions by optical attenuation of light trespassing more eccentric regions. The dashed red circle highlights roughly the boundaries of the sample region and is here defined as region of interest (ROI), over which some of the quantities are spatially averaged in the analysis that followed. The SNR in this case is about 4.5. This value was determined as the ratio between the average signal in the ROI and the standard deviation in circular region with radius from 16.8 mm to 24 mm. Given a concentration of GNR in the sample of 38 pM, this gives a noise equivalent sensitivity of around 9 pM.

Fig. 4

(a) Experimentally obtained PA signal image at 670 nm and 30 °C. The dashed circle highlights roughly the boundary of the phantom sample which is taken as region of interest (ROI). (b) Comparison of the absorption spectrum (blue) and PA spectrum (red) between 660 nm–1000 nm. Both spectra were normalized for the sake of easier comparison. The PA spectrum was obtained by calculating the mean value inside of the dashed circle, ROI. (c) Comparison between PA spectra of the sample, PDMS and water. The dip at 920 nm in the PA spectrum in the sample is clearly caused by shadowing of the surrounding PDMS container, which has a spectral peak at that wavelength.

In Fig. 4(b), the normalized spatial average over the ROI of the PA amplitude of the reconstructed images obtained between 660 nm and 1000 nm at 30 °C is compared with the respective normalized optical absorption coefficients obtained from optical transmission measurements at the same temperature. According to Eq. (3), in the absence of shadowing effects, the PA spectrum should be isomorphic with the spectrum of the optical absorption coefficient of the sample. This is clearly not the case. From the PA signal spectra in the water surrounding the phantom, PDMS and phantom sample (Fig. 4(c)), which were determined by analyzing the spatially averaged spectrum in the respective regions (268 mm2 annular region for PDMS, 20 mm2 rectangular region close to the bottom laser fiber for water, and 452 mm2 ROI region for sample) at 30 °C, it is obvious that the dip in the PA spectrum in the sample region between 870 nm to 930 nm coincides with the peak in the PA spectrum in the PDMS region, and thus is caused by shadowing due to strong absorption by PDMS. Although other temperature dependent shadowing effects, e.g. of inner regions by outer regions within the same material, are less obvious, they can have some disrupting effects on the extraction of temperature values from PA signal ratios.

The nature of the PA spectrum is quite different from the dual peak (670 nm and 965 nm) optical transmission spectrum that was used as a basis for the straightforward temperature mapping based on simulated data. We believe that this can be attributed to the difference in photoacoustic energy conversion between the gel and the nanoparticles. In order to assess the reliability of using features of the PA spectrum to reconstruct a temperature map with reduced sensitivity to detailed image noise, additional preprocessing of the data and a different choice of wavelengths was needed. First each of the PA images ROI was divided into 300 regions of equal surface, obtained by cross-sectioning 30 angular sectors (index θ) of 12° angular width and 10 rings (index r) of 1.5 mm2 surface. For each region, the spatial average of the PA spectrum was determined. In order to suppress stochastic fluctuations further and to assess trends of interest we rearranged the 4D data matrix P( θ,r,λ,T) into a matrix of spectra P(X,λ), with X(r,θ,T) an index that is representative for the position and temperature at which the spectrum was obtained. We then performed principal component analysis (PCA) to extract the mean spectrum and the eigenspectra [47]. We know from Fig. 4 that the PA spectrum of the GNR-agarose sample was affected by the absorption of water and PDMS, but the shadowing of the phantom sample region by the surrounding water and PDMS only has a minor influence on the temperature dependent character analysis by PCA.

For each spectrum, the principal component that accounted for the largest common variability, 84%, across temperatures and positions, was determined. Fig. 5(a) shows PCA determined mean PA spectrum for all regions and temperatures, and illustrates how the determined PA spectrum (mean PA spectrum + first eigen PA spectrum) for the sector region with θ=[0,12°] and r=[0,3.8 mm], and for 25 °C and 36 °C represents the essential features of the original spectrum while removing the noise. The difference in magnitude of the PCA spectra at 25 °C and 36 °C can be mainly attributed to the temperature dependence of the Grüneisen parameter. Their less pronounced shape difference results from the difference in temperature dependence of the absorption at different wavelengths. For every temperature and wavelength, also the spatially averaged value for the 300 subregions in the ROI was determined. Fig. 5(b) illustrates that the PCA-processed PA spectra are monotonically increasing with increasing temperature for every wavelength. The overall increase is about (2.2 ± 0.5)%/°C (evaluated at 36 °C), which is of the order of the temperature coefficient of the Grüneisen parameter of water (2.7% at 36 °C, cfr Eq. (2)). Fig. 5(c) shows that the temperature dependence of the PCA-processed PA amplitude for wavelengths in the range 670 nm – 870 nm is roughly linear. In view of finding a spectral parameter that is linked to the ratio rather than to the average amplitude of the PA signal spectrum, we have also determined the temperature dependence of the ratio between the principle component of the PA spectra at the two extremes of the reliable part of the spectrum, 670 nm and 840 nm (many of the PA images above 840 nm were distorted because of too high absorption or too much shadowing), as shown in Fig. 5(d). Therefore, unlike the simulation based analysis in section 4, which made use of the ratio between images at 965 nm and 700 nm wavelength, here we have made use of the R(T) image ratio between 840 nm and 670 nm. Fig. 5(e) shows the statistics of that ratio at 25 and 36 °C. The spatial averages, in the ROI, inside the red dashed circle as in Fig. 4(a) and excluding the outlier tails of the distributions, are different by about 6.3% (per 36-25 = 11 °C, with an average of R = 0.61 ± 0.01 at 25 °C and R = 0.65 ± 0.01 at 36 °C). This validates the possibility to use the ratio to determine absolute temperature. Note that the temperature sensitivity of R of 0.52% per degree with respect to the ratio at 36 °C is of the same order but not equal to the one that was used in the simulation study (0.77% per degree for wavelengths 965 nm and 700 nm). This is because (i) the simulation made use of different wavelengths, and, (ii) in the simulation it was assumed that the PA signal magnitude was proportional to the optical absorption coefficient as derived from optical transmission data.

Fig. 5

(a) Comparison between the mean PCA spectrum, the PA spectra at 25 and 36 °C, spatially averaged over one region (θ=[0,12°] and r=[0,3.8 mm]), and their largest principal component. (b) PCA-processed spectra in this region at different temperatures. (c) Corresponding amplitude variation with temperature at selected wavelengths. (d) PA ratio between 840 and 670 nm variation with temperature. (e) Statistical distribution of the PA ratio after PCA at 36 °C and 25 °C.

Fig. 6(a) (c) and (b) (d) shows respectively the statistical and spatial distribution of the slope and intercept of linear fits of the temperature dependence of the ratio, R(T). In the northeast and southwest regions, outliers appear due to the image reconstruction artefacts as discussed above. In order to discard these regions from the evaluation, from hereon, we limit the analysis to the 75% most frequently occurring points of reconstructed data in ROI and consider the remaining 25% as outlier zones. The green curves shown in Fig. 6(a) and (b) represents the best fitting Student t-distribution, which we have chosen because of its suitability to characterize the ratio of random quantities, where denominator variations due to outliers may be amplified, leading to longer tails. The black lines demarcate the interval of the 75% most frequently occurring values (from hereon referred to as the 75% region). The white cross hatched regions in Fig. 6(c) and (d) depict the discarded outlier regions. The same procedure was applied to the temperature coefficient, determined from the R(T) slope and the image ratio at 36 °C. The used slope was averaged over the 75% region as described in Fig. 6(c). Fig. 6(e) compares the statistical variation of the obtained temperature coefficients within the 75% region, in which the temperature coefficient is found to vary between 0.51%/°C and 0.55%/°C while the coefficient in whole ROI is changing between 0.51%/°C and 0.66%/°C (for comparison: in the simulation study, the coefficient is between 0.63%/°C and 0.82%/°C). Fig. 6(f) shows the spatial distribution of the temperature coefficient. White hatching is used to depict the discarded 25% outlier regions. In this case, the spatial average of the temperature coefficient is around 0.52%/°C, with a standard deviation is 0.01%/°C. The standard deviation on the temperature coefficient in the ROI including the outlier zones is 0.02%/°C, with a spatial average of 0.53%/°C. While in the simulated case, standard deviation of 0.05%/°C and average of 0.77%/°C are obtained. Larger average value is extracted because of more temperature sensitive wavelength selection at 965 nm in the simulation.

Fig. 6

(a), (b), (c), (d) Statistical (a, b) and spatial (c, d) distribution of the R(T) slope (a, c) and intercept (b, d) within the gel cylinder. The green curves in (a) and (b) are t-distribution fits to the 75% most frequent values indicated by the black dashed lines. (e), (f) Statistical (e) and spatial (f) distribution of the temperature coefficient, based on the spatially averaged value of the R(T) slope and the local image ratio at 36 °C. The slope is spatially averaged over the 75% region. The white hatching in (c), (d) and (f) indicates the 25% outlier region.

In order to make the analysis more tangible, we have analyzed temperature reconstructions using the local slope and local intercept in Fig. 6(c) and (d). Fig. 7(a) and (b) shows temperature images for the water bath at 26 °C and 36 °C. In the reconstructed region without outliers, the reconstructed temperature is deviating between 0.5–1 °C from the bath temperature. Fig. 7(c) and (d) compare the spatially averaged reconstructed temperature values with the bath temperature and depicts the standard deviation within the determined 75% better reconstructed points. In most of the temperature cases, the match is better than 0.5 °C.

Fig. 7

Pseudo color temperature maps obtained at 26 °C (a) and 36 °C (b), within the gel cylinder. Plots of (c) spatially average reconstructed temperature versus set water bath temperature and (d) reconstructed temperature deviation distributions. The figures were obtained from the local slope and intercept of R(T), as shown in Fig. 6(c) and (d) respectively. The error bars in (c) and deviation in (d) are based on the 75% most frequently occurring values, excluding the outlier regions that are indicated by white hatching in (a) and (b).

In clinical implementations it would not be feasible to perform a spatially resolved calibration as the one used to perform the image reconstructions in Fig. 7(a–d). One can imagine that in practice a quick calibration to obtain the slope of the ratio versus temperature curve, averaged over the best reconstructed part of the region of interest is performed in a limited temperature range that is not harmful for the tissue. Alternatively, a calibration over a large temperature range could be done on a specimen from a tissue bank, which has similar characteristics as the tissue of interest. In order to test this approach, we have determined the average slope of the above data in the 75% region. Next, based on a reference image at 36 °C and taking the same approach as the one used to obtain the simulated result in Fig. 3(f), the ratio between PA signal maps at 670 nm and 840 nm was converted into a temperature map. Fig. 8(a) shows the obtained temperature map for a bath temperature of 26 °C. In the 75% region, the temperature prediction is quite uniform and acceptable. The root mean square reconstruction error is about 1.2 °C at 26 °C (Fig. 8(b–c)) and tends to 0 °C as the temperature approaches the reference temperature 36 °C.

Fig. 8

Pseudo color temperature map obtained at 26 °C (a). Plots of (b) spatially averaged reconstructed temperatures versus set water bath temperature and (c) reconstructed temperature deviation distributions. The curves were obtained from the slope average over its 75% most frequent values in the ROI and with a reference temperature of 36 °C. For Figs. 7(d), 8 (c) and 9 (b), a bin size of 0.5 °C was used to calculate the distributions, and the black line at the base of the plots indicates a 0 °C deviation.

Fig. 9

Plots of (b) reconstructed average temperatures versus bath temperatures and (c) reconstructed temperature error distributions based on the PA signal at 670 nm as in a traditional Grüneisen parameter-based thermography approach. The analysis was done in the region of 75% most frequently occurring values.

For the technique to be suitable for thermotherapy monitoring, an accuracy of 1 °C seems feasible: tissue coagulation typically occurs at 50–55 °C [32,43], and irreversible damage to tissue is expected to occur when it is heated to 50 °C for 6–12 min or longer [48,49] Nevertheless, for more refined applications, such as detection of locally increased or decreased metabolism due to tumor growth or inflammation in organs, 3D temperature tomography with an accuracy of 0.1 °C would be needed. Moreover, since the experimental results presented were obtained on an ideal sample in terms of uniformity of material composition, optical density and temperature, it is important to understand the cause of the imperfections in the above temperature reconstructions. In the simulation section we already discussed the effect of wavelength and temperature dependent shadowing of incoming light in deep regions as a consequence of temperature and wavelength dependent absorption by shallow regions. This effect is difficult to correct for, but it can be minimized by choosing a thermochromic material with a low absorption to avoiding shadowing effects and with a large difference in temperature coefficient between two wavelengths in a wavelength range for which the tissue is as transparent as possible. Also imaging artefacts, e.g. caused by a not perfectly symmetric transducer geometry, should be avoided. A translucent ring of scattering material could be placed around the sample in order to provide a more uniform distribution of incoming laser light.

For the sake of completeness, we have also performed “traditional” photoacoustic temperature mapping, based on a single wavelength, and thus exploiting the temperature dependence of the product of the Grüneisen parameter and the optical absorption coefficient at that wavelength. The average slope over the 75% region and a reference image at a temperature of 36 °C were used to predict temperature images in the range 25–37 °C. Fig. 9 shows that the average resulting reconstruction accuracy at 26 °C is limited to 0.8 °C. Discarding the 25% of outlier regions, with which the standard deviation on the reconstructed temperature is about 3.6 °C at 26 °C, the reconstruction uncertainty is about 1.2 °C in the 75% region of the spectrally resolved method depicted in Fig. 8. In spite of this, we are convinced that spectrally resolved information can be more robust under many real-life circumstances, mainly in cases of non-uniform tissue composition and thus variability in Grüneisen parameter properties in the region of interest and for hyperthermia therapy going beyond 50 °C. In the latter cases, the temperature dependence of the PA amplitudes is affected, while the temperature dependence of the ratio of PA amplitudes at different optical wavelengths remains reliable. Moreover, by designing thermochromic entities with non-trivial temperature induced spectral changes [50,51], neural network recognition of the full spectrum can be used to give additional robustness to the reconstructed temperature images [[52], [53], [54]].

Conclusion

A novel temperature imaging method was proposed, based on multi-wavelength laser excitation and PA amplitude ratio calculation. The possibilities and limitations of the method are assessed, both on the basis of a simulation and of experimental data, by using the ratio between the PA signal magnitudes at two optical wavelengths.

Despite limitations introduced by ultrasonic imaging imperfections and temperature and wavelength dependent shadowing effects due to the setup, and although the experimentally accessible temperature control range was limited to the range 25 °C–37 °C, the experimental results show a promising potential for this MSOT method in the context of pre-clinical applications and clinical monitoring during hyperthermia therapy, complementary to conventional Grüneisen-parameter-based photoacoustic temperature tomography. For the experimental results on a phantom sample containing GNRs, PCA was applied for noise reduction and a temperature reconstruction was obtained based on the photoacoustic signal ratio between 840 nm and 670 nm. The temperature reconstruction error at 26 °C was found to be of the order of 1.2 °C and 0.3 °C, respectively with and without using of a reference ratio map at 36 °C, which was 0.22 °C and 0.01 °C in the simulation respectively. These experimental values were determined after excluding the 25% most deviating values that occur in regions containing image reconstruction artefacts.

The above results were obtained with a limited amount of time averaging. Relaxing the acquisition rate thus gives room for further improvement of the accuracy. Further research is needed to find thermomarkers with a highly temperature-sensitive shape of the optical absorption spectrum. By making use of response surfaces and neural network recognition, multiple wavelengths can be used to increase the robustness of the information on temperature. Thermomarkers should also exhibit high optical contrast with other absorbers so that dependence of the temperature reconstruction parameters on the optical properties of the tissue can be minimized.

The main intrinsic limitation of the current approach is that the light fluence in deep regions (further from the light source) is influenced by temperature dependent optical absorption in more shallow regions (closer to the light source). This issue can be compared to the problem of shadowing in ultrasound echography imaging, where the echoes coming from deeper layers are weakened due to reflection of energy at more shallow layers. One can envisage to correct the reconstruction of deeper layers by using the PA magnitude and temperature at the shallow layers that affect the fluence, but further research is needed to validate the feasibility.

One feature of the method that limits its applicability in clinical applications is that it relies on a calibration of the temperature dependence of the PA spectrum. The question is whether calibration on a particular sample tissue would remain valid for other tissues or similar tissue of other patients.

In order to cope with this, region of interest targeting nanoparticles with high optical absorption with respect to the tissue of interest are needed. In that way, the PA amplitude is dominantly determined by the known nanoparticles, and less influenced by the particular tissue. If moreover the nanoparticles have a strong and linear temperature sensitivity, then a calibration of their temperature dependence in a limited temperature range (e.g. within a harmless temperature increase of 1 °) can be extrapolated to higher temperatures. The uncertainty on the extrapolation can be used to set a sufficient safety margin on the monitored temperature in case of hyperthermia treatment monitoring.

For 3D photoacoustic temperature tomography, encapsulation of the thermomarkers will also be needed, so that their optical absorption remains insensitive to their biochemical environment. We note that also without encapsulation and with the tissue itself as optical absorber, multispectral photoacoustic imaging can still be used for early detection of tissue damage, via the accompanying changes in features of the optical absorption spectrum.

Conflict of interest

The authors declare that there are no conflicts of interest.