2.1.

Vascular Tree Modeling

The tumor vascular tree is generated using a fractal geometry model.^24,34 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 $150\mu \mathrm{m}$ 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 $14\mu \mathrm{m}$ 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 $r^{\prime }$ and a point detector located at $r$.

The parameters of the chaotic vasculatures were acquired from literature^24,25 and are summarized here. The size of the daughter fragments is determined by the bifurcation index ($\beta$). The bifurcation index is used to correlate the left and right daughter segments through the equation below:

2.2.

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 $I$ and the absorption coefficient map of the tissue ${\mu }_a$. The product of these two parameters results in the heating function $H$. 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:^36,37

The heating function $H$ 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 $H(r_{\text{blood}},t)\gg H(r_{\text{tissue}},t)$ representing that the absorption of blood vessels is much larger than the absorption of the surrounding tissues, resulting in

2.3.

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 ($d$).

The diffused blood $(n)$ at location $x$ and time $t$ for a vessel boundary located at position $o$ (cylindrical symmetry was assumed) with initial concentration of $n_o$ was computed using the one-dimensional solution to the Fick’s law equation:

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 $30\mu \mathrm{m}$ 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)].

2.4.

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 $45\mu \mathrm{m}$ lateral resolution. The transducer was simulated for frequency ranges between 15 to 45 MHz (the $-6\mathrm{dB}$ 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 $\alpha$ was accounted for using the equation below:

2.5.

Experimental Validation

A vessel phantom was prepared using thin fishing line of $20\mu \mathrm{m}$ 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 $t=0$, 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.

2.6.

Photoacoustic Radiofrequency Spectral Analysis

The simulated PA signals were windowed using a moving Hann window of size $1.03\times 1.80\mathrm{mm}$ 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, $y$-intercept, and MBF. These parameters were averaged and compared for parameters derived from the vascular tree with and without extravasated blood.

3.1.

Modeling Vascular Bleeding from a Single Microvessel

A single microvessel $10\mu \mathrm{m}$ 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 $20\text{-}\mu \mathrm{m}$-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 $y$-intercept increases for this vessel geometry and ultrasound detection geometry/bandwidth.

3.2.

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 $20\mu \mathrm{m}$ is smaller than system axial resolution of $44\mu \mathrm{m}$. 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.

3.3.

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 $50\mu \mathrm{m}$, 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.

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. 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.

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 ($p<0.01$). The spectral slope is correlated to the size of the PA source.^11,16,41,42 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.^8,13 Figure 6(b) compares the $y$-intercept parameter, which is related to the concentration of the optical absorbers.^16,41,42 The $y$-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 ($p<0.05$). Figure 6(c) represents the MBF, which is a dependent parameter that combines the effect of the spectral slope and $y$-intercept. The results demonstrate a decrease in the MBF parameter 0.5 h after bleeding ($p<0.01$).

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.

Fig. 6

Spectral parameters acquired from vascular tree model simulation. The calculated parameters are the (a) slope, (b) $y$-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 ($p<0.05$).

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.^29,43,44 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.^45–47