Photoacoustic simulations of microvascular bleeding: spectral analysis and its application for monitoring vascular-targeted treatments.

Solid tumors are typically supplied nutrients by a network of irregular blood vessels. By targeting these vascular networks, it might be possible to hinder cancer growth and metastasis. Vascular disrupting agents induce intertumoral hemorrhaging, making photoacoustic (PA) imaging well positioned to detect bleeding due to its sensitivity to hemoglobin and its various states. We introduce a fractal-based numerical model of intertumoral hemorrhaging to simulate the PA signals from disrupted tumor blood vessels. The fractal model uses bifurcated cylinders to represent vascular trees. To mimic bleeding from blood vessels, hemoglobin diffusion from microvessels was simulated. In the simulations, the PA signals were detected by a linear array transducer (30 MHz center frequency) of four different vascular trees. The power spectrum of each beamformed PA signal was computed and fitted to a straight line within the −6-dB bandwidth of the receiving transducer. The spectral slope and midband fit (MBF) based on the fit decreased by 0.11  dB  /  MHz and 2.12 dB, respectively, 1 h post bleeding, while the y-intercept increased by 1.21 dB. The results suggest that spectral PA analysis can be used to measure changes in the concentration and spatial distribution of hemoglobin in tissue without the need to resolve individual vessels. The simulations support the feasibility of using PA imaging and spectral analysis in cancer treatment monitoring by detecting microvessel disruption.

Introduction

A class of cancer therapeutic drugs targeting blood vessels, known as vascular disrupting agents (VDAs), induces intertumoral hemorrhage. The drugs are activated at the location of the tumor to disrupt endothelial and/or intravascular cells (e.g., red blood cells), which results in hemoglobin leakage into the surrounding tissues. The outcome is vessel disruption followed by tumor cell death and extensive hemorrhagic necrosis within hours of treatment delivery.,

Photoacoustic (PA) imaging has been proposed to monitor tumor response to various types of cancer therapies. PA imaging can be used to acquire images with high contrast and greater penetration depth, compared with other optical and ultrasound modalities. The most common type of analysis in PA imaging is the time-domain amplitude analysis of the PA signals. However, microvessels deep in tissue cannot be resolved using acoustic-resolution PA imaging. Ultrasound frequency analysis of the PA signals has been shown to detect changes smaller than the system resolution, and this analysis may permit assessment of microvessel bleeding. The ability to detect and quantify the occurrence of bleeding can have an important role for PA modeling in cancer therapy.

The potential of spectral analysis of the PA signals has not been fully investigated, especially for cancer treatment monitoring. Utilizing spectral analysis can improve the sensitivity of PA imaging to surpass the acoustic-resolution approaches. Changes in the frequency content of the PA signals are associated with the size, concentration, and spatial distribution of dominant optical absorbing chromophores. This is analogous to quantitative ultrasound analysis, which is associated with the dominant ultrasound scattering source.,

PA spectral analysis provides system-independent parameters that are correlated to the structure of optically absorbing chromophores., To extract these spectral parameters, the time-domain radiofrequency (RF) PA signals are transformed into the frequency domain. System-dependent parameters are removed by normalizing the PA signal to a reference power spectrum. The slope, -intercept, and midband fit (MBF) parameters are extracted from the line of best fit of the normalized power spectra and are correlated to tissue microstructures. The advantages of PA spectral analysis have been used in assessing tumors, liver conditions, and osteoporotic patients.,,

For cancer monitoring, a decrease in the PA spectral slope at multiple wavelengths (45% at 750 nm and 73% at 850 nm) has been reported 2 h after introducing a VDA (HaT-DOX). The changes in the spectral slope have been linked to hemorrhaging of blood vessels. However, the impacts of bleeding due to intertumoral hemorrhaging on the generated PA signals have never been simulated.

Simulations were undertaken to model the generation of PA signals from vascular trees. This would allow us to examine the feasibility of detecting the effects of VDAs on the PA signals from tissue. The tumor microvessel geometry was simulated using a vascular tree model. The model uses cylinders that bifurcate into two daughter vessels, including a bleeding model as a result of vascular collapse, which has been applied in modeling the extravasation of blood outside the vessels. The detection was simulated for widely available linear array transducers. In this work, we perform simulations that can provide insights into PA spectral analysis parameters to assess the effectiveness of cancer therapy. This could potentially improve the therapeutic outcome and reduce the tumor reoccurrence rate noninvasively and without the need of contrast agents, to provide timely feedback on the effectiveness of the drug/therapeutic approach.

Simulation Approach

The vascular tree structure of tumor vessels was simulated using a fractal model with cylindrical bases. The extravasation of blood outside the vessels was simulated using Fick’s law of diffusion, as previously done to model bleeding. The PA signals of an omnidirectional point detector were generated using the solution to the PA wave equation with a custom simulator of the Green’s function approach accounting for the directivity and detector bandwidth. Two-dimensional PA images were generated using delay and sum of the PA signals. Spectral analysis was performed by extracting the spectral slope, -intercept, and MBF.

Vascular Tree Modeling

The tumor vascular tree is generated using a fractal geometry model., In this paper, vascular tree geometrical parameters were selected to model chaotic vessels similar to vessels found in tumors [Fig. 1(a)]. The simulations were performed in three dimensions, starting with the third branch of in diameter and 1 mm in length simulated at the origin point (0, 0, 0) mm, leading up to the 12th branch with an average diameter of and length of 0.35 mm. These values were chosen to mimic measured tumor vasculatures of mammary carcinoma.

Fig. 1

(a) An illustration of the fractal tree geometrical parameters used for generating the vascular tree. The parent segment bifurcates into two daughter segments with different diameters and lengths. (b) A schematic representation of the detected PA signal from an absorber located at and a point detector located at .

The parameters of the chaotic vasculatures were acquired from literature, and are summarized here. The size of the daughter fragments is determined by the bifurcation index (). The bifurcation index is used to correlate the left and right daughter segments through the equation below: where the value range is . The left daughter segment () was assumed to be larger than right daughter segment () to represent the diameter asymmetry. Solving for and is achieved using the hemodynamic energy minimizing constrains. The was fixed to a value of 0.95 for this simulation. The distance () was used to correlate the length () of the daughter and parent segments and distance was set to 0.9 for this simulation. Since the simulations are done in 3D, there were two angles for each daughter segment. The first angle is called the branching orientation () and ranges from 0 deg to 360 deg, while the other angle is called the branching angle () and ranges from 25 deg to 140 deg. These parameters were used to control the branching angle of the daughter segments with respect to the parent segment.

Photoacoustic Signal Generation

When short laser pulses irradiate optically absorbing objects, the absorbers emit PA waves that encode their geometric information. The PA signal can be computed for a defined laser profile and the absorption coefficient map of the tissue . The product of these two parameters results in the heating function . The heating function was used to compute the shape of the generated PA signal. The forward solution of the PA wave equation was computed based on the free-space Green’s function:, where is the thermal coefficient of the volume expansion, is the specific heat capacity of the tissue at constant pressure, is the speed of ultrasound in water set at , is the location of the point detector, and is the location of the absorber, as demonstrated in Fig. 1(b). The above equation is used to simulate the PA signal generated from the point detector with a small aperture as a response of the heat deposition function . The equation assumes constant density and speed of sound between the vessels and the background tissue. The equation applies spherical integration with a radius determined by the acoustic time of flight. The constant only contributes to the amplitude of the PA signal and was set to 1 for these simulations.

The heating function was computed for individual vessels, and generated pressure signals were then superimposed for the entire tree with the appropriate geometrical time delays. Each cylindrical vessel has a defined radius, start position, and end position specified by the fractal tree geometry. The model assumes uniform illumination of the laser light with representing that the absorption of blood vessels is much larger than the absorption of the surrounding tissues, resulting in

Modeling the Diffusion of Blood outside Microvessels

To model VDAs, the diffusion of extravasated blood outside the vessels was simulated using Fick’s law of diffusion. Fick’s law relates the amount of bleeding that we anticipate to the duration of time elapsed since vascular disruption. It states that the rate of change of the concentration at a point in space is proportional to the second derivative of concentration in space with constant diffusion coefficient ().

The diffused blood at location and time for a vessel boundary located at position (cylindrical symmetry was assumed) with initial concentration of was computed using the one-dimensional solution to the Fick’s law equation: where is the diffusion of hemoglobin with a value of and erfc is the complementary error function. To test the applicability of this approach, Eq. (4) was computed for a single microvessel in diameter and 2 mm in length. The diffused blood extends to the region in which hemoglobin is present. This results in an increase in apparent vessel size of full width at half maximum of 20, 34, 50, for time intervals of , 0.1, 0.5, and 1.0 h. The time intervals that were chosen as preliminary studies from our group suggest that vascular hemorrhaging occurs early after the administration of treatments that target blood vessels.,

The diffusion of blood outside the vessels was then applied to the entire vascular tree. The diffusion of blood was simulated using Fick’s law for vessels smaller than in diameter after 0.5 and 1 h of vessel disruption to represent newly generated tumor vessels. The generated PA signals from a vascular tree simulated with extravasated blood were compared with the same vascular tree without extravasated blood. The comparison involved generating PA B-mode images and performing spectral analysis on the resulting beamformed PA signals [Eq. (2)].

Photoacoustic Signal Modeling and Image Beamforming

The simulated PA signals were generated by considering several features of a high-frequency commercial ultrasound transducer used in our laboratory (VevoLAZR LZ550, Fujifilm-Visual Sonics, Toronto, Canada). It has the capability to image tumors at a depth of more than 1 cm with lateral resolution. The transducer was simulated for frequency ranges between 15 to 45 MHz (the bandwidth) using a Butterworth bandpass filter of the third order, which was applied to the generated PA signals. The PA signals were computed for 256 elements (14.08 mm in length). The distance of the transducer to the original parent vessel segment was set to 11 mm (the approximate focus of the laser beam).

The signals were beamformed using a conventional delay-and-sum method for every 64 elements. The directivity of the transducer was accounted for using the equation below: where is the angle to the vector line, is the diameter of a single element in the transducer, and is the ultrasound wavelength at the central frequency of the transducer ( for the speed of sound of ). Finally, apodization of the transducer was applied using a Hamming filter. The generated B-mode images of the vasculature tree with and without extravasated blood were then visually compared.

Experimental Validation

A vessel phantom was prepared using thin fishing line of in diameter embedded in 10% porcine gelatin skin (Sigma-Aldrich Inc., Canada). The fishing line was removed and Sudan black dye (Sigma-Aldrich Inc., Canada) dissolved in methanol was injected into the vessel. The phantom was then imaged using the VevoLAZR system (FujiFilm-VisualSonics Inc., Canada) coupled to the LZ-550 transducer at 700 nm. Signals of the vessel phantom were acquired at , 1, 5, and 10 min after setup. A reference signal was acquired from a 20-nm-thick gold film. The measured RF signals were extracted from the VevoLAZR system and beamformed using the conventional delay-and-sum method. The beamformed RF signals were windowed and transformed into the frequency domain. The calculated power spectra were normalized to the gold film reference signals.

Photoacoustic Radiofrequency Spectral Analysis

The simulated PA signals were windowed using a moving Hann window of size with 50% overlap. The windows were applied for elements 32 to 224 at the location of the vascular tree. The power spectra of the windowed signals were computed and averaged. The averaged power spectrum was fitted to a straight line for the simulated frequencies. The parameters calculated from the fitted line were the spectral slope, -intercept, and MBF. These parameters were averaged and compared for parameters derived from the vascular tree with and without extravasated blood.

Results and Discussion

Modeling Vascular Bleeding from a Single Microvessel

A single microvessel in diameter and 2 mm in length was modeled. The PA signal at a point detector 11 mm away from the microvessel source was computed using Eq. (2). Hemoglobin diffusion into the interstitium was modeled using Eq. (4) to account for microvessel damage due to the VDA. Spatial and frequency analyses of generated PA signals from the single microvessel are presented in Fig. 2. Figure 2(a) is a schematic representation of the geometry of a simulated microvessel. Figure 2(b) shows images of the cross section of two simulated microvessels before vessel disruption and 1 h after. This demonstrates the changes in the distribution and content of hemoglobin after bleeding. For bleeding, there is a gradient change in the hemoglobin content once it enters the tissues. In contrast, before bleeding there is a sharp change in the hemoglobin content outside the vessel. This will have an impact in the generated PA signals through increasing the effective absorber size and thereby affecting the spatial distribution of the optical absorbers.

Fig. 2

(a) A schematic representation of the simulated vessel geometry of -thick vessel. (b) A schematic representation of the cross section of the vessel in (a) before bleeding and after 1 h of bleeding simulated using Fick’s law. (c) Line profile of the blood distribution, generated PA signal in (d) spatial and (e) frequency domains at different time points after vessel disruption.

The effect of a vessel’s bleeding on the spatial and frequency contents of generated PA signals is presented in Figs. 2(c)–2(e). Figure 2(c) shows the distribution of hemoglobin before (0 h) and at different time intervals after bleeding (0.1, 0.5, and 1.0 h). Upon modeling vascular disruption, the previously localized hemoglobin now is progressively distributed outside the vessels, with the hemorrhaging diameter increasing as a function of time post bleeding. Figure 2(d) shows the simulated PA signals in the spatial domain as time progressed after the vessel damage; the amplitude of the PA signals decreased while the length of the PA signals increased. The decrease in the amplitude of generated PA signals is due to introducing a gradient change in the hemoglobin content. The temporal length of the PA signals increased with time after bleeding due to the increase in the spatial distribution of the hemoglobin occupying a larger area within the interstitium.

The changes of the PA signals in the spatial domain also translate into changes in the power spectra of the signals, as shown in Fig. 2(e). Within the transducer bandwidth simulated in this experiment (15 to 45 MHz), the power spectra exhibit significant changes as a result of hemoglobin diffusion into the interstitium. Specifically, the power spectra contain features (minima and maxima) that change with time as the hemorrhaging area increases. Such features are characteristic of increasing absorber sizes. As the bleeding progresses, the slope of the power spectra decreases, and the -intercept increases for this vessel geometry and ultrasound detection geometry/bandwidth.

Vessel Phantom

The B-mode PA images of the prepared vessel phantom with Sudan black dye injected into the vessel at two time points are shown in Fig. 3(a). From the B-mode images, it is difficult to detect the diffusion of the injected Sudan black dye from the vessel to the gelatin matrix. This is expected as a vessel size of is smaller than system axial resolution of . Figure 3(b) displays the power spectra of the RF signals of the vessel phantom at different time points. The power spectra reveal changes in the frequency component as the dye diffuses into the gelatin matrix with the passage of time. The changes in the frequency components include a decrease in the amplitude of the higher frequency components and an increase in the amplitude of the lower frequency components. These changes are similar to the simulated results of the single vessel shown in Fig. 2(e).

Fig. 3

(a) A vessel phantom imaged using the VevoLAZR injected with Sudan black dye at time equal 0 and 5 min and (b) the average power spectra of the imaged vessel phantom normalized to a thin gold reference slide.

Modeling Photoacoustic Signals of Vasculature before and after Bleeding

After the single vessel is simulated, PA signals from the entire vascular tree were modeled before and after vascular bleeding. An example of simulated tumor vasculature is shown in Fig. 4(a). The tumor vasculatures are chaotic and dense as is the case for experimentally measured tumor vasculatures. The location of the transducer is denoted by the blue line for 256 elements. The beamformed images for the selected bandwidth of 15 to 45 MHz are presented in Figs. 4(b) and 4(c) for vascular trees without bleeding and 0.5 h after bleeding, respectively. Although the ultrasound attenuation was not modeled in these simulations, at this frequency range, the ultrasound attenuation is higher compared with clinical ultrasound (1 to 5 MHz). The higher frequencies will result in more negative values in the reported spectral slope measurements due to frequency-dependent attenuation. However, the effect of the ultrasound attenuation will be similar for with and without microvessel bleeding. The results in Fig. 4 demonstrate the challenge of detecting the bleeding from nonresolvable microvessels in acoustic resolution PA imaging. Since the size of the vessels is smaller than the system resolution of , it is difficult to extract image features that would suggest changes in the vascular structure. Visual inspection of the B-mode images shows little difference in the images generated from the bleeding and nonbleeding vessels.

Fig. 4

(a) A schematic representation of the simulated vascular tree. The blue line represents the location of the linear array transducer considering of 256 simulated point source detectors. The beamformed image of the vascular tree (b) without microvessel bleeding and (c) 0.5 h after microvessel bleeding.

The parameters of the vascular tree model are used to approximate the vascular geometry of real tissue vessels. The vascular tree models have several assumptions. Vessels are perfectly cylindrical segments of constant radius that bifurcate into two daughter segments. Hence, there are no trifurcations and no “dead-end” junctions where only one child segment is formed. Actual vessels have a varying radius and follow a curved trajectory. Another assumption is that the tissues have only arterioles of constant oxygenation with uniformly distributed branching angles. Tissue mechanical properties are assumed acoustically homogeneous without any dispersion. Finally, it is assumed that the tissues have optical properties with uniform illumination from a single laser pulse. Regarding the light illumination, the loss of photons with depth due to absorption and scattering affects the frequency content of the power spectra to result in a more negative spectral slope.

We model bleeding effects as a function of time elapsed after vessels have been disrupted. Spatial and frequency domain signals of simulated PA signals generated from the vascular tree model are presented in Fig. 5. The time-domain PA signals simulated before bleeding and 0.5 h after bleeding are shown in Fig. 5(a) for bandlimited signals. Again, it is difficult to distinguish between the two cases by visual inspection of the PA signals. Figure 5(b) presents the average power spectra and the line best fitted before bleeding and 0.5 h after bleeding. The line fitted to the power spectra suggests a change between the two groups and that these changes are more pronounced at higher frequencies.

Fig. 5

(a) Representative PA signals acquired from simulating vascular tree without microvessel bleeding and 0.5 h after microvessel bleeding. (b) The average power spectra of four different vascular trees and the average line best fitted with its standard deviation.

A comparison of the selected spectral parameters between the three time points is presented in Fig. 6. These three groups represent the analysis of PA signals before bleeding, 0.5 h after bleeding, and 1 h after bleeding. Th e spectral slope is presented in Fig. 6(a). The slope decreased as early as 0.5 h of bleeding (). The spectral slope is correlated to the size of the PA source.,,, For our work, this corresponds to the hemorrhaging vessels. A decrease in the spectral slope at the measured frequencies suggests an increase in the size of the PA source. This is due to the diffusion of hemoglobin into the interstitium that generates an overall larger PA source. The decrease in spectral slope following vessel destruction has been observed in vivo., Figure 6(b) compares the -intercept parameter, which is related to the concentration of the optical absorbers.,, The -intercept parameter increased as bleeding progressed. However, these changes are not as significant as the changes to the spectral slope analysis. Changes were detected after 1 h of bleeding (). Figure 6(c) represents the MBF, which is a dependent parameter that combines the effect of the spectral slope and -intercept. The results demonstrate a decrease in the MBF parameter 0.5 h after bleeding ().

Fig. 6

Spectral parameters acquired from vascular tree model simulation. The calculated parameters are the (a) slope, (b) -intercept, and (c) MBF for the three groups of simulated vascular trees to represent an intact tumor vasculature, 0.5 h after bleeding and 1 h after bleeding. * Statistical significance compared with the without bleeding group ().

Spectral analysis of simulated PA signals could be used to differentiate between vasculatures before and 0.5 h after bleeding through analyzing the changes in the spectral slope parameter. The decrease in the spectral slope correlates to an increase in the effective vessel diameter after bleeding due to the diffusion of blood outside of the vasculature. The MBF could also be used to detect microvessel bleeding; however, it is a dependent variable on the spectral slope.

In this study, we modeled hemorrhaging of tumor blood vessels post VDA to assess the potential of PA RF analysis for monitoring this process. Other changes due to vascular bleeding could potentially be used as a marker to monitor cancer therapy. In future work, this could be confirmed by further simulation accounting for the distribution of multiple chromophores. For example, as red blood cells escape the vessels, the concentration of deoxyhemoglobin and methemoglobin rises due to changes in the surrounding environment.,, These changes can lead to alteration in PA signals that depend on the optical wavelength. We confined these simulations to early time points, so we would not have to consider the new chromophores that are expected to be created due to the oxidation of hemoglobin (such as methemoglobin and hemichrome). The chromophore distribution as a result of bleeding could also be modeled using a vascular tree, and its effect on generated PA signals can be analyzed at different wavelengths. Finally, additional validation of these simulations can be performed by comparing our results with in-vivo results obtained from a tumor mouse model treated with different types of VDAs. Our group is currently conducting detailed studies to examine in-vivo changes of PA spectral parameters measured during cancer treatment. The combination of ultrasound backscatter to monitor changes in cellular structure and PA to monitor changes in vascular structure provides a unique approach to monitor structural changes in the most important elements of tumors, which are cancer cells and the tumor vasculature.

Conclusion

This study demonstrates how PA spectral analysis can be used in principle to detect structural changes due to bleeding at the micron scale. According to the simulations performed, changes due to microvessel bleeding can be detected by quantifying changes in the spectral slope and MBF of the power spectra in ultrasonic resolution PA. This demonstrates a potential for the use of PA spectral analysis for cancer treatment monitoring and for treatments that target the vasculature. The early clinical feedback that could be potentially obtained by detection and quantification of microvessel bleeding may result in a significant improvement in treatment monitoring and, eventually, outcome. It may also result in a reduction in the reoccurrence rate of tumors through feedback provided on the efficacy of the treatment and through using bleeding as a biomarker of tumor response.