Exploiting Nonlinear Photoacoustic Signal Generation in Gold Nanospheres for Selective Detection in Serial 3 D PA Tomography

The photoacoustic (PA) signal amplitude measured in gold nanosphere suspensions has been shown to increase nonlinearly with the incident excitation fluence. In this work, this effect is exploited to recover the spatial distribution of gold nanoparticles in tomographic 3D photoacoustic (PA) images against the background contrast provided by absorbers that exhibit a linear relationship between the PA signal amplitude and the fluence. Serial tomographic PA images of a tissue phantom containing gold nanospheres and a tissue-mimicking absorber were acquired. By assessing the linearity of the PA intensity voxel by voxel, the spatial distribution of the gold nanosphere suspension was recovered. The method is shown to enable the robust detection of gold nanoparticles.


Introduction
Photoacoustic (PA) tomography is a noninvasive, non-ionizing imaging technique that combines optical specificity and sensitivity with high ultrasound resolution and penetration depth [1].Biomedical PA imaging relies on the absorption of short optical pulses by endogenous tissue chromophores such as hemoglobin as well as exogenous contrast agents in certain applications [2].Due to thermoelastic expansion, PA waves are generated inside the tissue and detected at the tissue surface.3D images of the initial pressure are then reconstructed from the time-resolved PA signals.
The use of gold nanoparticles as exogenous contrast agents offers several advantages to applications in biomedical PA imaging.Gold nanoparticles are characterized by a high absorption cross-section due to localized surface plasmon resonance, which in turn provides strong PA contrast [3], including in vivo [4][5][6][7][8].Furthermore, gold nanoparticles have a high biocompatibility, and studies have shown that gold nanoparticles in certain applications are nontoxic to human and animal cells [9].Finally, gold nanoparticles can be functionalized by conjugation to other molecules, tumor ligands, or monoclonal antibodies, thus enabling detection and visualization of early-stage tumors and in vivo tracking of targeted particles [10].However, the unambiguous detection of gold nanoparticles against the overwhelming background absorption by endogenous tissue chromophores remains a challenge.
Recent theoretical models describing the PA signal generated in gold nanospheres suggest that its amplitude is nonlinear with respect to fluence.This is due to the high temperatures that are reached in a nanoscaled volume in the immediate vicinity of the gold nanosphere, which cause a non-negligible change in the thermal expansion coefficient (and hence Grüneisen parameter) of water [11].Early studies investigating this effect had only been able to observe the nonlinearity at a temperature of 4 • C (where the Grüneisen parameter for water is zero) for gold nanospheres of 40 nm in diameter [12].Recently, Pang et al. [13] presented experimental evidence that the predicted nonlinearity of the PA signal amplitude with respect to fluence can be observed at room temperature in suspensions of large gold nanospheres with a diameter of up to 150 nm.This nonlinearity leads to a greater PA signal amplitude in gold nanospheres compared to that measured in a homogeneous dye solution of identical absorption coefficient and linear fluence dependence of the PA signal amplitude [1].
This study provides experimental evidence that the aforementioned nonlinear fluence dependence of the PA signal amplitude can be exploited to detect and visualize the spatial distribution of gold nanospheres against the background contrast that is provided by absorbers that exhibit a linear PA signal with fluence, such as typical endogenous tissue chromophores.A phantom containing a gold nanosphere suspension and an ink-based absorber was imaged using serial 3D PA tomography at increasing excitation fluences.Two image analysis methodologies were developed and evaluated to obtain a robust method for recovering the distribution of absorbers that exhibit a nonlinear PA signal with respect to fluence.

Tissue Phantom
The tissue phantom consisted of polymer tubes (inner diameter 460 µm, outer diameter 500 µm) (Paradigm Optics, Vancouver, BC, USA) immersed in a scattering lipid suspension (diluted milk, µ s ' ≈ 1 mm −1 ) as shown in Figure 1a,b.The tubes were parallel to each another and positioned 1 mm apart at a depth of approximately 2-3 mm below the illuminated surface of the lipid suspension (Figure 1a).The tubes were filled with either a tissue-mimicking absorber or a colloidal suspension of monodispersed gold nanospheres.India ink was chosen as the tissue-mimicking absorber which is known to exhibit a linear PA signal with respect to fluence.The gold nanospheres colloids (BBI Solutions, Crumlin, UK) were of 150 nm diameter and suspended in water.Experiments were carried out at the undiluted colloidal concentration, which was given by the manufacturer as 1.66 × 10 9 particles/mL, i.e., an optical density of OD1 at 532 nm.India ink was diluted with water to an optical density similar to that of the gold nanosphere colloid.After filling the polymer tubes with the gold nanosphere suspension, the particles were allowed to settle to a form a thin layer at the bottom of the tubes as illustrated in Figure 1b.

PA Imaging Setup
A schematic of the experimental setup for acquiring serial photoacoustic (PA) images is shown in Figure 1a.An all-optical 3D PA scanner based on a Fabry-Pérot polymer film ultrasound sensor [14,15] was used to acquire 3D image data sets in forward mode.A wavelength tunable OPO laser system (Spitlight 1000 OPO, Innolas, Munich, Germany) provided PA excitation pulses of 8 ns duration at a wavelength of 532 nm and a repetition rate of 30 Hz.The output of the OPO coupled into a 1.5 mm multimode fiber in order to homogenize the beam.The fiber output was collimated and directed onto the tissue phantom, which was placed on top of the Fabry-Pérot ultrasound sensor.The beam diameter was approximately 9 mm at the phantom surface.
A tunable diode CW laser (Yenista Tunics T100S, Lannion, France) with a wavelength range of 1500 to 1630 nm was employed to interrogate the Fabry-Pérot sensor.The acoustically modulated optical power reflected by the sensor was detected using a photodiode, the output of which was recorded using a digitizer card (PCI-5124, National Instruments, USA).By raster-scanning the interrogation laser across the sensor surface, the spatial and temporal distribution of the PA pressure field was mapped in 2D from which 3D PA images of the initial pressure distribution were reconstructed using a fast Fourier transform (FFT)-based image reconstruction algorithm [16].

Serial PA Tomography
A total of eight PA image data sets of the tissue phantom were acquired at fluences ranging from 7.7 mJ/cm 2 to 34.6 mJ/cm 2 in steps of approximately 4 mJ/cm 2 .Each image data set was acquired in forward mode over a scan area of 6 mm x 15 mm with a step size of 100 μm.To increase the signal-to-noise-ratio, PA signals were averaged over three acquisitions at every scan position.The speed of sound was estimated for every tomographic scan by employing an autofocus function [17], which yielded an average sound speed of 1501 m/s.

Recovering the Spatial Distribution of Gold Nanospheres from 3D PA Images
A method based on determining the degree of nonlinearity in the PA image intensity as a function of fluence was developed to recover the gold nanosphere distribution on a per-voxel-basis.To minimize computation time, voxels corresponding to the absorber-filled tubes were segmented from the background by applying a linear fluence-dependent image intensity threshold to each of the datasets acquired from the serial PA tomographic imaging procedure.The threshold level was chosen ad hoc to remove the image background without affecting the image intensity of the tube volumes.The segmented tube subvolumes represent voxels where the intensity of the majority of the serial images was found to be greater than the threshold.This included data points where low absorption and/or low fluence produced image intensities below the threshold.Voxels with an intensity below the threshold irrespective of the fluence were set to zero.The segmented subvolumes were subsequently analyzed voxel by voxel.The image intensity as a function of fluence was extracted to which two functions were fitted.
Method I involved using a linear regression analysis according to Equation 1 for the voxel-by-voxel fitting analysis after the image segmentation.

Serial PA Tomography
A total of eight PA image data sets of the tissue phantom were acquired at fluences ranging from 7.7 mJ/cm 2 to 34.6 mJ/cm 2 in steps of approximately 4 mJ/cm 2 .Each image data set was acquired in forward mode over a scan area of 6 mm × 15 mm with a step size of 100 µm.To increase the signal-to-noise-ratio, PA signals were averaged over three acquisitions at every scan position.The speed of sound was estimated for every tomographic scan by employing an autofocus function [17], which yielded an average sound speed of 1501 m/s.

Recovering the Spatial Distribution of Gold Nanospheres from 3D PA Images
A method based on determining the degree of nonlinearity in the PA image intensity as a function of fluence was developed to recover the gold nanosphere distribution on a per-voxel-basis.To minimize computation time, voxels corresponding to the absorber-filled tubes were segmented from the background by applying a linear fluence-dependent image intensity threshold to each of the datasets acquired from the serial PA tomographic imaging procedure.The threshold level was chosen ad hoc to remove the image background without affecting the image intensity of the tube volumes.The segmented tube subvolumes represent voxels where the intensity of the majority of the serial images was found to be greater than the threshold.This included data points where low absorption and/or low fluence produced image intensities below the threshold.Voxels with an intensity below the threshold irrespective of the fluence were set to zero.The segmented subvolumes were subsequently analyzed voxel by voxel.The image intensity as a function of fluence was extracted to which two functions were fitted.Method I involved using a linear regression analysis according to Equation (1) for the voxel-by-voxel fitting analysis after the image segmentation.
where p is the PA image intensity, Φ is fluence, and a 1 is the linear regression coefficient.The coefficient of determination to the linear regression (R 2 lin ) was used to distinguish the voxels representing the gold nanospheres from those representing the background absorption provided by the India ink.
Method II involved fitting a second-order polynomial according to Equation (2) to the fluence-dependence of PA image intensity: where a 1 and a 2 are the first and second-order fit coefficients, respectively.a 2 was used as a measure of nonlinearity to distinguish voxels containing gold nanospheres from those belonging to the linear absorber (India ink) or the background.

Reconstructed PA Signal Intensity Images
Figure 1c shows a representative reconstructed 3D PA image acquired at an excitation fluence of 17.8 mJ/cm 2 .The image region shown will be used for all subsequent analyses.The images of the phantom show the location of the tubes filled with gold nanospheres and India ink.It should be noted that the gold nanospheres sedimented at the bottom of the tubes as seen in Figure 1c.By contrast, no sedimentation was observed in the tubes filled with India ink and a more homogeneous initial pressure distribution can be seen across the tube lumen.However, irrespective of the phantom geometry, a PA image data set acquired at a single excitation fluence does not provide enough information for differentiating the two absorbers.

Nonlinear Fluence Dependence of the Signal Amplitude Measured in Gold Nanospheres
Figure 2a shows a set of xy-plane maximum intensity projections (MIPs) of the 3D images from serial measurements at eight different excitation fluences.The PA signal from all absorbing regions increases with fluence which is clearly visible in the images.It can also be seen in Figure 2a that the increase in the PA image intensity with fluence is more pronounced in the gold nanosphere-filled tubes compared to in the India ink.To provide a quantitative comparison of the dependence of PA image intensity on fluence, the mean image intensities of two regions of interest (ROI) as indicated in Figure 2a, corresponding either to a subvolume within a tube filled with gold nanospheres or India ink, respectively, were determined.The mean PA signal intensity of each ROI is plotted as a function of fluence in Figure 2b.The error bars represent the standard deviation of all voxels within the respective ROI. Figure 2b shows that the image intensity in the tubes filled with India ink increases linearly with fluence, while nonlinear fluence dependence can be observed for the gold nanospheres.This was quantitatively confirmed by fitting a linear and a second-order polynomial function to each plot from which the coefficient of determination (R 2 ) was obtained.For the linear regression (Equation ( 1)), the plot corresponding to India ink resulted in R 2 lin = 1.00 while that of the tubes filled with gold nanospheres yielded R 2 lin = 0.59.The PA image intensity measured from the tubes filled with gold nanospheres has a lower coefficient of determination to a linear regression, confirming the nonlinear fluence dependence of the PA signal (with a continuously increasing rate of change).However, fitting a second-order polynomial function to the data corresponding to the gold nanospheres shows better agreement than the linear regression as evidenced by a coefficient of determination, R 2 2p , of 0.99.The fitted functions are shown as dashed lines in Figure 2b.Additional non-linear functions were examined to fit the nonlinear PA signal in the gold nanospheres, including a third-order polynomial (R 2 3p = 0.98) and an exponential function (R 2 exp = 0.99).However, because neither of these alternatives provided a significantly better agreement with the data, and because the nonlinearity is expected to follow a quadratic function [13], the second-order polynomial was chosen for subsequent analysis.

Recovering the Gold Nanosphere Distribution from Serial PA images
Figure 3 shows a schematic of the steps involved in the image analysis process of Method I as described in Section 2.4.The 3D image of resulting from applying Method I to the serial PA imaging data sets is shown in xy-and xz-MIPs in Figure 3a.The voxels representing the India ink are approximately unity, as expected because of the linear relationship between the image intensity and fluence.By contrast, the voxels representing the gold nanospheres exhibit a lower value due to the nonlinear fluence dependence.Figure 4b presents a box plot comparing the values of the image regions that correspond to the two absorbers.Note that the outliers primarily result from voxels close to the tube boundaries.A two-sample t-test reveals a significant difference in (p < 0.005) between the two types of absorbers.Linear regression of the relationship between fluence and PA image intensity can be applied to distinguish gold nanospheres from background absorbers.The results of Method II are shown in Figure 4 a as xy-and xz-MIPs of , the coefficient to the second order term in Equation 2. The voxels corresponding to the regions containing gold nanospheres show a higher coefficient than those corresponding to India ink, with the latter displaying a near zero coefficient, which is expected as the PA signal from the India ink is linear  lin values of the image regions that correspond to the two absorbers.Note that the outliers primarily result from voxels close to the tube boundaries.A two-sample t-test reveals a significant difference in R 2 lin (p < 0.005) between the two types of absorbers.Linear regression of the relationship between fluence and PA image intensity can be applied to distinguish gold nanospheres from background absorbers.

Recovering the Gold Nanosphere Distribution from Serial PA images
Figure 3 shows a schematic of the steps involved in the image analysis process of Method I as described in Section 2.4.The 3D image of resulting from applying Method I to the serial PA imaging data sets is shown in xy-and xz-MIPs in Figure 3a.The voxels representing the India ink are approximately unity, as expected because of the linear relationship between the image intensity and fluence.By contrast, the voxels representing the gold nanospheres exhibit a lower value due to the nonlinear fluence dependence.Figure 4b presents a box plot comparing the values of the image regions that correspond to the two absorbers.Note that the outliers primarily result from voxels close to the tube boundaries.A two-sample t-test reveals a significant difference in (p < 0.005) between the two types of absorbers.Linear regression of the relationship between fluence and PA image intensity can be applied to distinguish gold nanospheres from background absorbers.The results of Method II are shown in Figure 4 a as xy-and xz-MIPs of , the coefficient to the second order term in Equation 2. The voxels corresponding to the regions containing gold nanospheres show a higher coefficient than those corresponding to India ink, with the latter displaying a near zero coefficient, which is expected as the PA signal from the India ink is linear The results of Method II are shown in Figure 4a as xy-and xz-MIPs of a 2 , the coefficient to the second order term in Equation ( 2).The voxels corresponding to the regions containing gold nanospheres show a higher a 2 coefficient than those corresponding to India ink, with the latter displaying a near zero a 2 coefficient, which is expected as the PA signal from the India ink is linear with fluence.In Figure 4a, the tubes containing the India ink cannot be seen clearly because the second order term a 2 within the background was set to zero and is therefore similar to the near-zero a 2 coefficients corresponding to India ink.A boxplot of a 2 for the image regions corresponding to the two absorbers is shown in Figure 4b.The results of a two-sample t-test reveal a significant difference in a 2 (p < 0.005) between regions with gold nanospheres and India ink.Therefore, Method II is capable of differentiating between gold nanospheres and absorption sources that exhibit a linear PA signal with respect to fluence in 3D PA images.
second order term within the background was set to zero and is therefore similar to the near-zero coefficients corresponding to India ink.A boxplot of for the image regions corresponding to the two absorbers is shown in Figure 4b.The results of a two-sample t-test reveal a significant difference in (p < 0.005) between regions with gold nanospheres and India ink.Therefore, Method II is capable of differentiating between gold nanospheres and absorption sources that exhibit a linear PA signal with respect to fluence in 3D PA images.

Definition of Contrast-To-Noise Ratio for Distinguishing the Gold Nanoparticle Signal from the Ink Signal
In the previous section, the parameters R 2 lin and a 2 were shown to provide simple and robust metrics for recovering the spatial distribution of gold nanospheres.Both R 2 lin and a 2 exhibit a significant difference (p < 0.005) when comparing the values from the gold nanospheres and India ink.To analyze which of these methods provides better contrast in the identification of gold nanospheres in a measured 3D PA image, a contrast-to-noise ratio (CNR) metric was defined as where x is the mean of the chosen parameter (either parameter R 2 lin or a 2 ), σ x is the standard deviation of the chosen parameter, and the subscripts "Gold" and "Ink" indicate from the image regions corresponding to gold nanospheres and Indian ink, respectively.The CNR for Method I (linear regression) was found to be 2.6 whereas it was 1.6 for Method II (polynomial fit).These values are also included in Figures 4 and 5.This indicates that Method I provides a better CNR when the analysis is carried out on eight PA data sets.

Definition of Contrast-To-Noise Ratio for Distinguishing the Gold Nanoparticle Signal from the Ink Signal
In the previous section, the parameters and were shown to provide simple and robust metrics for recovering the spatial distribution of gold nanospheres.Both and exhibit a

Robustness of the Method against Background Noise
Method I and Method II were also applied to the detection of gold nanospheres without using a threshold to segment the tubes from the background (removing the step shown in Figure 3b or simply setting the threshold intensity to zero).This was done to determine the robustness of the methods against background noise, which may adversely affect the results.3D maps of R 2 lin and a 2 obtained using Method I and Method II without image segmentation are shown in the xy-and xz-MIPs in Figure 6a,b, respectively.Figure 6a indicates that the image noise confounds the detection of the gold nanospheres using Method I (linear regression) in the xy-MIP of R 2 lin .While the typically low R 2 lin values for gold nanospheres can be seen, they are surrounded by substantial noise in the determined parameter.This is confirmed in Figure 6c, which shows a cross-sectional xy-image of the 3D image at the tube depth (as indicated by the dotted red line in Figure 6a,b).The typically low R 2 lin values for gold nanospheres are recognizable but could not be separated unambiguously from the background.By contrast, the results of Method II (polynomial fit) suggest that it enables the robust detection of gold nanospheres even in the presence of noise as shown in in Figure 6b,d, respectively.The tubes filled with India ink merge with the background because both regions show a negligible second-order a 2 coefficient.
To compare the contrast from Method I and Method II against the background noise in a quantitative manner, a contrast-to-background-noise ratio (CNR bg ) metric was defined as where the subscript "bg" indicates the background region between the tubes.The denominator of Equation ( 4) will be referred to as background noise.The CNR bg of Method I (linear regression) was found to be 3.0 whereas for Method II (polynomial fit) it was 21.6 (see also Figure 6).This indicates that Method II may be a generally applicable approach to detecting absorbers with a nonlinear fluence-dependent PA signal amplitude.Method II is robust against background noise and yields a higher CNR bg compared to Method I.

Minimum Number of Measurements and Maximum Permissable Exposure
The methods were compared in terms of CNR and CNRbg as a function of the number of measurements and maximum fluence.Of particular interest were the image data sets acquired at fluences below the maximum permissible exposure (MPE), which is 20 mJ/cm 2 for the laser wavelength and pulse duration used in the current experiments reported here [1].The number of images included in the analysis, , was varied, starting from the lowest fluence (7.7 mJ/cm 2 ) up to a specific maximum fluence.Figure 7a shows the resulting CNR as a function of for Methods I and II.The corresponding maximum fluence is also indicated.For Method I, CNR appears to plateau until 5 (i.e. a maximum fluence of 23.4 mJ/cm 2 ) and decreases for 5.This would suggest that acquiring additional images at fluences higher than the MPE reduces the CNR obtained using Method I.By contrast, the CNR of Method II does not show a strong dependence on N. The biggest discrepancy between Methods I and II occurs at N < 6.For all values of , the CNR from Method I is higher than from Method II.The highest CNR of 6.4 was found using Method I on a subset of images acquired at five fluences up to 23.4 mJ/cm 2 .Method II yields a maximum CNR of 2.5 for a subset of four images acquired at fluences up to 17.8 mJ/cm 2 .3D images of and obtained using Methods I and II, respectively, from images acquired at fluences below the MPE are shown in Figure 8a,b and illustrate that sufficient contrast for gold nanosphere detection is obtained without exceeding the MPE.The CNRbg obtained using Method II is greater than that of Method I for all .Method II is therefore a more robust and generally applicable method irrespective of N compared to Method I for recovering the distribution of gold nanospheres from tomographic PA images.

Minimum Number of Measurements and Maximum Permissable Exposure
The methods were compared in terms of CNR and CNR bg as a function of the number of measurements and maximum fluence.Of particular interest were the image data sets acquired at fluences below the maximum permissible exposure (MPE), which is 20 mJ/cm 2 for the laser wavelength and pulse duration used in the current experiments reported here [1].The number of images included in the analysis, N, was varied, starting from the lowest fluence (7.7 mJ/cm 2 ) up to a specific maximum fluence.Figure 7a shows the resulting CNR as a function of N for Methods I and II.The corresponding maximum fluence is also indicated.For Method I, CNR appears to plateau until N = 5 (i.e., a maximum fluence of 23.4 mJ/cm 2 ) and decreases for N > 5.This would suggest that acquiring additional images at fluences higher than the MPE reduces the CNR obtained using Method I.By contrast, the CNR of Method II does not show a strong dependence on N. The biggest discrepancy between Methods I and II occurs at N < 6.For all values of N, the CNR from Method I is higher than from Method II.The highest CNR of 6.4 was found using Method I on a subset of images acquired at five fluences up to 23.4 mJ/cm 2 .Method II yields a maximum CNR of 2.5 for a subset of four images acquired at fluences up to 17.8 mJ/cm 2 .3D images of R 2 lin and a 2 obtained using Methods I and II, respectively, from images acquired at fluences below the MPE are shown in Figure 8a,b and illustrate that sufficient contrast for gold nanosphere detection is obtained without exceeding the MPE.The CNR bg obtained using Method II is greater than that of Method I for all N. Method II is therefore a more robust and generally applicable method irrespective of N compared to Method I for recovering the distribution of gold nanospheres from tomographic PA images.

Nonlinearity of the Fluence-Dependent PA Image Intensity in Gold Nanosphere Suspensions
The nonlinearity of the fluence-dependent PA image intensity measured in the tubes filled with gold nanospheres (Figure 2) is in agreement with the findings of Pang et al [6].Suspensions of nanospheres with a diameter of 100 nm and 150 nm were shown to exhibit a nonlinear quadratic PA signal amplitude with respect to fluence.A similar behavior was also seen for even larger diameters of 200 nm by Egerev and Oraevsky [18].By contrast, gold nanospheres of less than 100 nm in diameter were not found to exhibit nonlinearity, except at temperatures near 4 °C where the Grüneisen parameter of water is zero [12].In each of the aforementioned studies, PA measurements were made at a single point.
Theoretical models have succeeded in describing the appearance of a nonlinear fluence-dependence of the PA signal amplitude through the temperature dependence of the thermal expansion coefficient.The initial PA pressure, p, is proportional to the fluence, i.e., Φ Γ • • Φ, where Γ is the Grüneisen parameter, T is the temperature, and is the absorption coefficient [1].
It is also well accepted that the PA signal generation in nanoparticle suspensions occurs primarily from thermal expansion in the surrounding fluid and not in the nanosphere itself [19].For gold nanospheres, theory suggests a temperature-induced change in  within a nanoscaled volume surrounding the nanosphere is the primary cause for the nonlinear pressure generation [11].If 

Nonlinearity of the Fluence-Dependent PA Image Intensity in Gold Nanosphere Suspensions
The nonlinearity of the fluence-dependent PA image intensity measured in the tubes filled with gold nanospheres (Figure 2) is in agreement with the findings of Pang et al [6].Suspensions of nanospheres with a diameter of 100 nm and 150 nm were shown to exhibit a nonlinear quadratic PA signal amplitude with respect to fluence.A similar behavior was also seen for even larger diameters of 200 nm by Egerev and Oraevsky [18].By contrast, gold nanospheres of less than 100 nm in diameter were not found to exhibit nonlinearity, except at temperatures near 4 °C where the Grüneisen parameter of water is zero [12].In each of the aforementioned studies, PA measurements were made at a single point.
Theoretical models have succeeded in describing the appearance of a nonlinear fluence-dependence of the PA signal amplitude through the temperature dependence of the thermal expansion coefficient.The initial PA pressure, p, is proportional to the fluence, i.e., Φ Γ • • Φ, where Γ is the Grüneisen parameter, T is the temperature, and is the absorption coefficient [1].
It is also well accepted that the PA signal generation in nanoparticle suspensions occurs primarily from thermal expansion in the surrounding fluid and not in the nanosphere itself [19].For gold nanospheres, theory suggests a temperature-induced change in  within a nanoscaled volume surrounding the nanosphere is the primary cause for the nonlinear pressure generation [11].If   The nonlinearity of the fluence-dependent PA image intensity measured in the tubes filled with gold nanospheres (Figure 2) is in agreement with the findings of Pang et al. [6].Suspensions of nanospheres with a diameter of 100 nm and 150 nm were shown to exhibit a nonlinear quadratic PA signal amplitude with respect to fluence.A similar behavior was also seen for even larger diameters of 200 nm by Egerev and Oraevsky [18].By contrast, gold nanospheres of less than 100 nm in diameter were not found to exhibit nonlinearity, except at temperatures near 4 • C where the Grüneisen parameter of water is zero [12].In each of the aforementioned studies, PA measurements were made at a single point.
Theoretical models have succeeded in describing the appearance of a nonlinear fluence-dependence of the PA signal amplitude through the temperature dependence of the thermal expansion coefficient.The initial PA pressure, p, is proportional to the fluence, i.e., p(Φ) = Γ(T) • µ • Φ, where Γ is the Grüneisen parameter, T is the temperature, and µ is the absorption coefficient [1].It is also well accepted that the PA signal generation in nanoparticle suspensions occurs primarily from thermal expansion in the surrounding fluid and not in the nanosphere itself [19].For gold nanospheres, theory suggests a temperature-induced change in Γ within a nanoscaled volume surrounding the nanosphere is the primary cause for the nonlinear pressure generation [11].If Γ increases approximately linearly with temperature, and thus also linearly with fluence, then the PA signal would exhibit a nonlinear change with approximately a quadratic dependence.If the temperature change is negligible after photoacoustic excitation, as in the case for most tissue chromophores, Γ remains constant and the PA signal is linear with respect to fluence.Γ is the Grüneisen parameter and µ is the absorption coefficient [1].It is also well accepted that in nanoparticle suspensions, the majority of the PA signal generation occurs from thermal expansion in the surrounding fluid and not in the nanosphere itself [19].For gold nanospheres, the theory suggests a change in Γ within a nanoscaled volume surrounding the nanosphere is the primary cause for the nonlinear pressure generation [11].If Γ increases approximately linearly with temperature, and thus also linearly with fluence, then the PA signal would exhibit a nonlinear change with approximately a quadratic dependence.If the temperature change is negligible after photoacoustic excitation, as in the case for most tissue chromophores, Γ remains constant and the PA signal is linear with respect to fluence.

Sedimentation of the Gold Nanosphere Colloid
The gold nanospheres were found to sediment at the bottom of the tube, which differs from the experiments reported in Pang et al. [13] on nanospheres suspended in water.To avoid any temporal changes in the local concentration of the gold nanospheres, imaging acquisitions were started at least 30 min after the preparation of the phantom to allow the particles to settle.This was confirmed using an approximation of the sedimentation rate using Stoke's law [20], which gives a sedimentation time of 30 min for 150 nm diameter particles in water and the diameter of the tubes.
The possibility that close proximity or aggregation of gold nanospheres may have affected the fluence-dependence of the PA signal amplitude was investigated by making independent measurements in agar gel suspensions of 150 nm gold nanospheres using a 1D backward-mode PA sensor (based on the one used by Pang et al. [13]).The agar gel prevented particle movement, and thus sedimentation, and the results confirmed the nonlinear fluence-dependence reported in the literature [12,13,18].It was therefore assumed that negligible aggregation occurred in the settled nanospheres used in this study.Agar gel was not used in the phantom due to the practical implementation of filling thin tubes without introducing bubbles.Also, in real biomedical applications, stable suspensions may not be found since gold nanoparticles can accumulate in close proximity at tumors or other specific cell sites [21,22].It is therefore advantageous to validate detection methods, such as those reported in this paper, on phantoms that do not provide gold nanospheres in a stable colloidal suspension.

Potential for Biomedical PA Imaging Applications
This work demonstrates that the nonlinear fluence-dependence of the PA signal measured in gold nanospheres can be used to recover their spatial distribution from 3D images acquired below the MPE.Previous studies have also shown that the nonlinear PA behavior of gold nanospheres is seen with serial PA excitation to fluence values well under the MPE (8 mJ/cm 2 , Pang et al. [13]; 6 mJ/cm 2 , Simandoux et al. [12]), indicating that an even lower fluences may be sufficient.This suggests that the current methods may be applicable to recover the spatial distribution of gold nanospheres in vivo even if they are located at greater depths below the surface compared to those in the phantom.
According to the literature, gold nanospheres of 100 nm in diameter and larger have been found to exhibit a nonlinear fluence-dependent PA signal amplitude [13,18].However, nanospheres of up to 80 nm in diameter have been shown to have a linear dependence [13].While small nanoparticles show a more widespread distribution in various organ systems after intravenous injection [23,24], the proposed methods may not be used to detect nanospheres below 100 nm, unless a nonlinear PA signal with respect to fluence manifests itself through other mechanisms such as aggregation or bubble formation [25,26].Nanospheres of up to 200 nm in diameter have nevertheless been detected in blood, liver, lung, and spleen.Minute concentrations of large gold nanoparticles have also been detected in the brain, stomach, and pancreas [24].Therefore, although larger size spheres, such as the 150 nm diameter nanospheres used in the current study, may have a limited biodistribution, they may be applied in in vivo biomedical imaging studies targeting specific organs.
Gold nanospheres exhibit an absorption maximum in the green wavelength region, which coincides with strong absorption in blood.For biomedical PA tomography, excitation wavelengths in the visible and near infrared part of the electromagnetic spectrum are preferred since they offer penetration depths of up to several centimeters [1].Exogenous contrast agents absorbing in this optical window, such as gold nanorods, are therefore of interest [10].To date, the linearity of the PA signal with respect to fluence has not yet been thoroughly investigated in gold nanorod suspension.In principle, the methods presented in this paper can also be used for the detection of gold nanorods provided if certain sizes or shapes of gold nanoparticles are found that exhibit a nonlinear fluence dependence of the PA signal amplitude.

Conclusions
This study demonstrates an experimental imaging methodology for recovering the spatial distribution of gold nanospheres from tomographic 3D images.The methodology exploits the nonlinear fluence-dependent PA signal amplitude observed in nanospheres with a diameter above 100 nm.Strong contrast between gold nanospheres and the India ink was observed using methods based on fitting a linear function and a second-order polynomial function to the PA image intensity as a function of fluence voxel-by-voxel.The results also show that PA imaging using fluences up to the MPE is sufficient for the robust recovery of the spatial distribution of gold nanospheres.In cases where regions containing the absorbers of interest can be segmented from the rest of the image, the method based on linear regression produced the best contrast-to-noise ratio.If the entire image is used, the method based on a second-order polynomial fit was shown to be superior.The methodology may be applied for selective detection of gold nanospheres, or any other contrast agent exhibiting a nonlinear fluence-dependence of the PA signal, in biomedical PA imaging in vivo against the contrast produced by endogenous tissue chromophores.As a result, gold nanoparticles may be better imaged in applications in in vivo molecular PA imaging of superficial tumors and their vascular supply, in monitoring of drug delivery, and in evaluation of treatment efficiency.

Figure 1 .
Figure 1.(a) Schematic of the experimental setup.The tissue phantom was placed on the Fabry-Pérot sensor and illuminated from the top by the output of the excitation laser.(b) Schematic of the tissue phantom, which consisted of parallel polymer tubes filled with either gold nanospheres (diameter 150 nm) or diluted India ink.The gold nanospheres were allowed to settle before the PA imaging was performed.(c) 3D photoacoustic (PA) image of the phantom acquired at a fluence of 17.8 mJ/cm 2 .

Figure 1 .
Figure 1.(a) Schematic of the experimental setup.The tissue phantom was placed on the Fabry-Pérot sensor and illuminated from the top by the output of the excitation laser.(b) Schematic of the tissue phantom, which consisted of parallel polymer tubes filled with either gold nanospheres (diameter 150 nm) or diluted India ink.The gold nanospheres were allowed to settle before the PA imaging was performed.(c) 3D photoacoustic (PA) image of the phantom acquired at a fluence of 17.8 mJ/cm 2 .

Figure 2 .
Figure 2. (a) xy-plane maximum intensity projections (MIPs) of the reconstructed PA images of the 3D tissue phantom acquired at fluences ranging from 7.7 mJ/cm 2 to 34.6 mJ/cm 2 .(b) Mean PA image intensity in two regions of interest (ROIs) corresponding to gold nanospheres or India ink as a function of fluence.The error bars indicate the standard deviation.The dashes lines show the fitted linear and polynomial functions.

Figure 3 .
Figure 3. Flow chart of Method I.After image reconstruction (a), a linear, fluence-dependent threshold is applied to remove the image background (b).(c) The fluence-dependent image intensity is extracted from the segmented voxels, which correspond to the location of absorber-filled tubes.A linear function is fitted to the image intensity voxel by voxel, and the coefficient of determination, , is plotted in the final image (d).All images shown are MIPs.

Figure 2 .
Figure 2. (a) xy-plane maximum intensity projections (MIPs) of the reconstructed PA images of the 3D tissue phantom acquired at fluences ranging from 7.7 mJ/cm 2 to 34.6 mJ/cm 2 .(b) Mean PA image intensity in two regions of interest (ROIs) corresponding to gold nanospheres or India ink as a function of fluence.The error bars indicate the standard deviation.The dashes lines show the fitted linear and polynomial functions.

3. 3 .
Figure 3 shows a schematic of the steps involved in the image analysis process of Method I as described in Section 2.4.The 3D image of R 2 lin resulting from applying Method I to the serial PA imaging data sets is shown in xy-and xz-MIPs in Figure 3a.The voxels representing the India ink are approximately unity, as expected because of the linear relationship between the image intensity and fluence.By contrast, the voxels representing the gold nanospheres exhibit a lower R 2 lin value due to the nonlinear fluence dependence.Figure 4b presents a box plot comparing the R 2lin values of the image regions that correspond to the two absorbers.Note that the outliers primarily result from voxels close to the tube boundaries.A two-sample t-test reveals a significant difference in R 2 lin (p < 0.005) between the two types of absorbers.Linear regression of the relationship between fluence and PA image intensity can be applied to distinguish gold nanospheres from background absorbers.

13 Figure 2 .
Figure 2. (a) xy-plane maximum intensity projections (MIPs) of the reconstructed PA images of the 3D tissue phantom acquired at fluences ranging from 7.7 mJ/cm 2 to 34.6 mJ/cm 2 .(b) Mean PA image intensity in two regions of interest (ROIs) corresponding to gold nanospheres or India ink as a function of fluence.The error bars indicate the standard deviation.The dashes lines show the fitted linear and polynomial functions.

Figure 3 .
Figure 3. Flow chart of Method I.After image reconstruction (a), a linear, fluence-dependent threshold is applied to remove the image background (b).(c) The fluence-dependent image intensity is extracted from the segmented voxels, which correspond to the location of absorber-filled tubes.A linear function is fitted to the image intensity voxel by voxel, and the coefficient of determination, , is plotted in the final image (d).All images shown are MIPs.

Figure 3 .
Figure 3. Flow chart of Method I.After image reconstruction (a), a linear, fluence-dependent threshold is applied to remove the image background (b).(c) The fluence-dependent image intensity is extracted from the segmented voxels, which correspond to the location of absorber-filled tubes.A linear function is fitted to the image intensity voxel by voxel, and the coefficient of determination, R 2 lin , is plotted in the final image (d).All images shown are MIPs.

Figure 4 .
Figure 4. Results of Method I (linear regression).(a) xy-and xz-MIPs of the coefficient of determination of a linear regression, (Equation 1) applied to the PA image intensity as a function of fluence for each voxel of a 3D PA image.The contrast-to-noise ratio (CNR) computed with is also reported.(b) Boxplot of for the tubes filled with gold nanospheres compared with those filled with India ink.

Figure 5 .
Figure 5. Results of Method II (polynomial fit).(a) xy-and xz-MIPs of the coefficient from a second-order polynomial fit according to Equation 2, applied to the PA image intensity as a function of fluence for each voxel in the PA image.The contrast-to-noise ratio computed with is also reported.(b) Boxplot analysis of corresponding to the tubes filled with gold nanospheres and those filled with India ink.

Figure 4 .
Figure 4. Results of Method I (linear regression).(a) xy-and xz-MIPs of the coefficient of determination of a linear regression, R 2 lin (Equation (1)) applied to the PA image intensity as a function of fluence for each voxel of a 3D PA image.The contrast-to-noise ratio (CNR) computed with R 2 lin is also reported.(b) Boxplot of R 2 lin for the tubes filled with gold nanospheres compared with those filled with India ink.

Figure 4 .
Figure 4. Results of Method I (linear regression).(a) xy-and xz-MIPs of the coefficient of determination of a linear regression, (Equation 1) applied to the PA image intensity as a function of fluence for each voxel of a 3D PA image.The contrast-to-noise ratio (CNR) computed with is also reported.(b) Boxplot of for the tubes filled with gold nanospheres compared with those filled with India ink.

Figure 5 .
Figure 5. Results of Method II (polynomial fit).(a) xy-and xz-MIPs of the coefficient from a second-order polynomial fit according to Equation 2, applied to the PA image intensity as a function of fluence for each voxel in the PA image.The contrast-to-noise ratio computed with is also reported.(b) Boxplot analysis of corresponding to the tubes filled with gold nanospheres and those filled with India ink.

Figure 5 .
Figure 5. Results of Method II (polynomial fit).(a) xy-and xz-MIPs of the a 2 coefficient from a second-order polynomial fit according to Equation (2), applied to the PA image intensity as a function of fluence for each voxel in the PA image.The contrast-to-noise ratio computed with a 2 is also reported.(b) Boxplot analysis of a 2 corresponding to the tubes filled with gold nanospheres and those filled with India ink.

Figure 6 .
Figure 6.Results of Method I and II when omitting the image segmentation step.The 3D images are displayed as xy-and xz-MIPs showing (a) from Method I (linear regression) and (b) from Method II (polynomial fit).Contrast-to-background-noise ratios (CNRbg) are also reported.A dotted red line shows the z-plane of the xy cross-sectional images shown in (c) for from Method I and (d) from Method II.

Figure 6 .
Figure 6.Results of Method I and II when omitting the image segmentation step.The 3D images are displayed as xy-and xz-MIPs showing (a) R 2 lin from Method I (linear regression) and (b) a 2 from Method II (polynomial fit).Contrast-to-background-noise ratios (CNR bg ) are also reported.A dotted red line shows the z-plane of the xy cross-sectional images shown in (c) for R 2 lin from Method I and (d) a 2 from Method II.

Figure 7 .
Figure 7. (a) Contrast-to-noise ratio (CNR) and (b) contrast-to-background-noise ratio (CNRbg) of Methods I and II as a function of the number of images N included in the analysis.The maximum fluence is shown in brackets.CNR and CNRbg are defined by Equation 3 and Equation 4.

Figure 8 .
Figure 8. 3D image of (a) (Method I, linear regression) and of (b) (Method II, polynomial fit) using N = 4 images, all acquired with fluence values below the MPE.

Figure 7 .
Figure 7. (a) Contrast-to-noise ratio (CNR) and (b) contrast-to-background-noise ratio (CNR bg ) of Methods I and II as a function of the number of images N included in the analysis.The maximum fluence is shown in brackets.CNR and CNR bg are defined by Equations (3) and (4).

Figure 7 .
Figure 7. (a) Contrast-to-noise ratio (CNR) and (b) contrast-to-background-noise ratio (CNRbg) of Methods I and II as a function of the number of images N included in the analysis.The maximum fluence is shown in brackets.CNR and CNRbg are defined by Equation 3 and Equation 4.

Figure 8 .
Figure 8. 3D image of (a) (Method I, linear regression) and of (b) (Method II, polynomial fit) using N = 4 images, all acquired with fluence values below the MPE.

Figure 8 .
Figure 8. 3D image of (a) R 2 lin (Method I, linear regression) and of (b) a 2 (Method II, polynomial fit) using N = 4 images, all acquired with fluence values below the MPE.

1 .
Nonlinearity of the Fluence-Dependent PA Image Intensity in Gold Nanosphere Suspensions