2.1.

ARPAE/US System

The AR-PAM/US system has been described earlier.⁸ Briefly, the system comprises a 1064 nm Nd: YAG laser (DPS-1064-Q,¹⁶, China) operating at 1 kHz with pulse energy of $\sim 8\mathrm{mJ}$, an US pulser-receiver, a function generator, a DAQ PC, and an imaging probe. The ultrasound transducer (US XDC) was custom made by Capistro Labs (San Clemente, California), using Polyvinylidene Fluorid (PVDF) material. The US transducer central frequency is 15 MHz with a 115% bandwidth. The US pulser receiver was a Panametric (currently Olympus, Tokyo, Japan) 5900PR, the function generator was a Tektronix (Beaverton, Oregon) AFG1022, and the DAQ was an Alazartech (Pointe-Claire, Quebec, Canada) ATS9462. To develop this system, the probe was modified to use a large 5 mm laser line mirror for better optical and mechanical stability. The function generator generates pairs of pulses with a short $16.3\mu \mathrm{s}$ delay: the first triggers the laser to generate a PA pulse, and the second triggers the US pulser-receiver to generate a US pulse. Both signals are amplified by the pulser-receiver before being received by the DAQ PC and digitized at 180 MHz.

Using a LabVIEW (National Instrument, Austin, Texas) code, the DAQ PC also synchronously drives a stepper motor that rotates the imaging head via a drive shaft embedded in the probe, shown in Fig. 1(c). Circular microsteps scan 1000 PA and US A-lines in a full rotation. The system is limited by the laser pulse repetition frequency at 1 kHz, giving a 1 kHz A-line scan rate and an $\sim 1\mathrm{Hz}$ B scan rate. The PA and US data obtained from scanning are processed in real time and displayed on the DAQ PC monitor for evaluation. The 6 dB US lateral resolution and axial resolution were characterized by $7\mu \mathrm{m}$ diameter carbon fiber targets to be 83 and $30\mu \mathrm{m}$, respectively, at the transducer focus. The PA resolution, characterized with the same target, measured $128\mu \mathrm{m}$ laterally and $68\mu \mathrm{m}$ axially at the US focal distance.

Fig. 1

AR-PAM/US system components (a) Probe body picture. (b) Probe body CAD. (c) Probe body with transparent parts. The stepper motor, at the extreme left, turns the drive shaft in small radial increments. The drive shaft and probe body house the US XDC wires (not shown) and the optical fiber, which extends along the center of the probe. (d) Cross-sectional view of the imaging head. The drive shaft is fixed to the imaging head with set screws. The optical fiber passes through a small opening in the waterproof housing holding the laser line mirror that reflects light into the short fiber fixed in the center hole of the US XDC.

2.2.

Spiral Phantom Preparation

A tissue-mimicking phantom with embedded spiral tubes was produced to simulate the imaging setting of a colorectal tissue scan with blood flow. Specifically, a two-pour approach and assemblies of three-dimensional (3D) printed phantom molds shown in Fig. 2(a) were used to first produce a tubular gelatin phantom and then to embed a spiral of tubing on the outside of the phantom. The phantom was made of 8%w gelatin solubilized in degassed deionized (DI) water to match the speed of sound in biological tissue’s while limiting optical absorbance. To retard degradation of the phantom, 1%v Germall Plus preservative was added.¹⁷ This mixture was degassed in a vacuum chamber to remove air bubbles formed during mixing. After the mold was configured as in Fig. 2(b), the gelatin mixture was added and allowed to cool to 10°C, forming a small hollow cylindrical phantom. Then, PTFE tubing with a $580\mu \mathrm{m}$ ID was spiraled onto the phantom [Fig. 2(c)]. The outer wall mold is mounted around the phantom, as shown in Fig. 2(d), then the gelatin mixture is poured into the mold and allowed to cool to give the final spiral flow phantom shown in Fig. 2(e).

Fig. 2

Drawings illustrating the preparation of a spiral flow phantom. (a) Individual molds. (b) Mold assembly for making a small cylindrical phantom. (c) Resulting small cylindrical phantom with spiral tubing on its exterior. (d) Mold with sleeve for encasing the spiral tubing to produce the final phantom. (e) Image of completed spiral flow phantom with embedded tubing.

2.3.

Doppler Processing

Doppler flowmetry in the phantom and in tissue utilizes a correlation Doppler algorithm with extensive noise rejection features. The PA signal is first passed through an FIR filter with a 3 to 30 MHz passband, then the average intensity of the noisy regions is analyzed and used to apply amplitude-based thresholding to the signal to avoid noise correlation. The signals are interpolated by a factor of 3 to allow finer correlation shift calculation as:

For adjacent A-lines, a 25-point sliding window cross-correlation is calculated, and their correlation shift and correlation amplitude are recorded:

The correlation amplitude is used to filter out remaining noisy correlations. After filtering, the correlation shift can be converted to a velocity map based on Eqs. (1) and (2).¹⁸ Equation (1) calculates the Doppler velocity resolution, $V_{\mathrm{res}}$. Because the direction of the scanning and the direction of the flow are at an angle relative to one another, angle compensation is added to the velocity calculation by using Eq. (2).

Angle compensation is conventional in various Doppler velocimetry applications, where it is used to calculate the velocity component along the desired direction from the velocity measured along the A-line. Figure 3(b) shows the angle $\theta$ and the true versus measured velocities. We implemented angle compensation for in vivo images by calculating the angle between the imaging and flow direction, using the estimated time of arrival from Sec. 2.5 and the transducer surface offset of 5.57 mm, a process shown in Figs. 4(l)–4(n). The formula for the angle calculation is Eq. (3). Since the compensation factor approaches infinity as the angle between the flow direction and the scanning direction approaches 90 deg, the compensation angle is limited to 70 deg to prevent over-compensation at small angles:

Fig. 3

(a) Experimental flow phantom setup with linear tubing suspended in DI water (left) and with spiral tubing gelatin phantom fixed in DI water (right). Contrast agents suspended in DI water flow from the syringe pump through the tubing, which is imaged by the endoscopic imaging probe. (b) Illustration of the linear tubing being scanned. Upon laser exposure, optical contrast agents traveling in the tubing generate PA signals that are picked up by the US XDC. Angle $\theta$ for angle correction between the true velocity direction and measured velocity direction is illustrated in the zoomed-in figure.

Fig. 4

Time of arrival correction algorithm workflow and PA/US coregistration workflow. (a) The original US image is denoised at 10 dB above the noise floor to give (b). To compensate for tissue displacement, the rising edge in (c) is aligned with (d) the original PA image to yield (e) an aligned PA image, where the time of arrival is the same for all A-lines. Correlational Doppler processing yields (g). The Doppler data are subsequently backshifted (h) and coregistered with the original US image (j) to give (i) the Doppler and US coregistration image. This coregistered image is then converted to (k) polar display to better represent the scanning. (l) The rising edge time of arrival is changed to (m) polar display and used in angle calculation (n) for angle compensation.

To display the flowmetry, the velocity map is re-plotted with a color Doppler color scheme, in which black represents stationary regions, red to orange represents regions of increasing positive velocity, and blue to cyan represents regions of negative velocity with increasing speed. Finally, the color Doppler map is superimposed on the US image. For in vivo and spiral phantom images, the images are converted to a polar display for co-registered geometrical representation. For in vivo data, B scan images are uniformly divided into 20 regions, and the Doppler flow speeds in the regions are calculated for statistical analysis. One limitation of this correlational Doppler method is that it can detect flow only if absorbers move along the A-line direction. If they flow orthogonal to the A-line direction, they cannot be measured. In this research, the probe is almost always off-center, and the vasculature directions are highly variable, so the flow direction is usually not normal to the scan direction, which makes angle compensation particularly important in calculating flow velocities.

2.4.

Phantom Imaging Setup

Figure 3(a) shows the open loop flow setup for imaging the linear flow and spiral flow phantoms. The AR-PAM/US imaging system is mounted on a manual two-dimensional (2D) translation stage to allow fine adjustment of the scanning location. Because no appropriately sized microspheres absorbing at 1064 nm were available, we produced our own contrast agent by sonicating a mixture of 60%w $50\mu \mathrm{m}$ glass beads infused with 40%w hydrophobic black dye. A 2%w dilution of the contrast agent was found to be optimal for PA imaging, with only moderate photothermal damage to the flow tubing. The imaging setup and scanning geometry for the linear flow phantom experiment are shown in Fig. 3(a). The imaging head is held stationary for maximum consistency in imaging the $580\mu \mathrm{m}$ ID PTFE tubing, and both are submerged in degassed DI water for acoustical coupling. The axis of the imaging head is held at 60 deg to the tubing bore. To study the angle dependence of flow measurements, we wrote a LabVIEW program to rotate the imaging head, allowing fine adjustment of the scanning angle.

To simulate the in vivo imaging conditions, we imaged the spiral flow phantom, as shown in Fig. 3(a). The AR-PAM imaging probe is mounted on a two-axis translational stage, with its imaging head submerged in degassed DI water. The processing for the spiral flow phantom is special because the imaging direction is perpendicular to the spiral flow direction, thus constant angle correction was applied. Due to tissue motion and the handheld probe consideration, our application of Doppler velocimetry employs limited averaging of 100 A-lines to the spiral flow phantom and in vivo data.

2.5.

Time of Arrival Correction for In Vivo and Spiral Flow Phantom Data Processing

Time of arrival correction, shown in Fig. 4, compensates for shifts in the signal due to shifts in the physical location of the rectal tissue, allowing correlational Doppler to more accurately estimate shifts in the optical absorber locations. Assuming the PA and US A-lines overlap, we developed an algorithm for aligning the individual PA A-lines. Because the physical boundary between the coupling DI water and the phantom/tissue is well defined in the US image, rising edge-based US time of arrival detection was applied to the US A-lines, yielding Fig. 4(c), in which the red line represents the time of arrival. To eliminate noise points in the rising edge detection, the threshold was manually set to 10 dB above noise floor, yielding Fig. 4(b). To remove outliers, the calculated time of arrival locations were passed through an outlier rejection algorithm. Based on their times of arrival, PA A-lines are shifted in time for correlational Doppler calculations and shifted back to their original location for coregistration with the US image, as shown in Figs. 4(e)–4(i). The PA Doppler image co-registered with the US image is then warped for display in polar coordinates, as shown in Fig. 4(k).

2.6.

In Vivo Imaging of Post Treatment Colorectal

For in vivo imaging, patients with rectal adenocarcinoma who had undergone radiation and chemotherapy were imaged pre-operatively before rectal tissue resection. Patients were first put under anesthesia per the standard of care, a colonoscopic exam was performed to localize the cancer, and our AR-PAM/US probe was then used to scan the patient’s rectum. The findings from five patients are included in this paper: three patients (P1, P2, P3) with residual adenocarcinoma, one with complete response (P4) of no remaining cancerous tissue, and one patient (P5) with a stricture below the tumor so only lesion-free tissues were imaged. This study was approved by the Institutional Review Board and was Health Insurance Portability and Accountability Act (HIPAA) compliant (ClinicalTrials.gov identifier NCT04339374), and all patients provided written informed consent.

3.1.

Linear Phantom Flow Experiment

A set of flow speeds, ranging from 0.63 to 6.3 mm/s in increments of 0.63 mm/s were measured with the linear flow imaging setup to observe the flow profile. Twenty US and PA coregistered B Scans were acquired for each flow speed setting, with the results shown in Fig. 5(a). The measurements demonstrate good linearity across the flow speeds in the measurement interval, and even though a slight discrepancy between the theoretical and measured values can be observed, the differences are within 15%.

Fig. 5

(a) Flow speed measurements from 0.01 to $0.1\mathrm{ml}/\min$. A fitted line with $R^2$ value of 0.965 and theoretical velocity is shown. (b) Experimental angle compensation results, where blue is the calculated velocity, red is the compensated velocity, and yellow is the actual flow velocity.

To study the effect of angular changes on the flow speed measurement, we varied the angle between the imaging head and flow direction from 70 deg to 50 deg, in 10 deg increments, while maintaining the flow speed at $6.3\mathrm{mm}/\mathrm{s}$. As Fig. 5(b) shows, the angle compensation procedure can compensate for large angles between the flow direction and imaging direction, but at angles $<70\mathrm{deg}$, over-compensation occurs because the compensation factor increases exponentially. Consequently, we limited the angle compensation to 70 deg.

3.2.

Spiral Doppler Flow Experiment

Figure 6(a) shows a 3D PA scan of the spiral flow phantom. This image allowed us to identify the turn of the tubing spiral that overlapped the most with the 2D imaging plane and consequently provided the most flow speed data. It also confirmed that the tube had no obstructions or kinks to cause flow speed irregularities. The flow speed in the spiral phantom was calculated from AR-PAM/US system data. Because the PA signals from the inner surface of the tubing are strong, as shown in Fig. 6(b), only a weak signal can be observed inside the tubing. The Doppler image coregistered with the US image, Fig. 6(c), appears quite similar, with measurable flow near the inner tubing surface. At the same time, in Fig. 6(c), we observe flow inside the tubing highlighted by Doppler processing, in contrast to the low signal amplitude in Fig. 6(b). These sets of images emphasize Doppler processing’s strength in providing perfusion information with little dependence on signal amplitude, which is valuable for in vivo data. Due to the difficulty of aligning the PA data with the low US signal, PA Doppler could be calculated only for regions with a strong US signal, but this should not be an issue for in vivo data since the water–tissue interface generates a strong US reflection.

Fig. 6

Spiral flow phantom imaging. (a) 3D AR-PAM endoscopic image of the spiral phantom. Contrast dye deposited on the inner radius of the spiral tubing wall generated a higher PA signal than dye on the outer radius. (b) PA coregistered with the US image of the scan plane that maximally intersects the tubing bore (c) Doppler image coregistered with US image. (d) Variation of Doppler-calculated flow speed with actual speed. 100 A-line ROI averaging is applied to the data, and a fitted curve with $R^2$ value of 0.913 is shown.

Statistically, spiral flow signals exhibit linear trends similar to those linear flow, but with larger deviations, which is expected since only 100 A-line ROI averaging was applied to the data as compared to averaging across 1000 A-lines in Fig. 5(a), and the signal-to- noise ratio is poorer because a portion of the tube does not exhibit a strong PA signal. Still, the actual speed and estimated speed are approximately linearly correlated with an $R^2$ value of 0.913, demonstrating the utility of flowmetry in spiral imaging scenarios.

3.3.

In Vivo Processing Results

Doppler processing was applied to in vivo images to test velocity estimation. Six rectal images were processed for five different patients: three patients with residual cancer (P1-3), one treatment responder (P4) with no residual cancer, and one patient (P5) with normal rectal tissue.

Figures 7(a)–7(c) show a polar scan, co-registered US and AR-PAM of a tumor bed, and a flow map of the same region, respectively, demonstrating that correlational Doppler can highlight low flow in a tumor bed (the darker regions in the red box) as distinct from the nearby higher flow regions. The residual tumor in Figs. 7(a)–7(c) displays perfusion on the surface, with no perceivable flow inside the tumor, possibly due to necrotic tissue, but higher velocity and densely packed vessels are observed adjacent to the residual tumor. In the treatment responder [Figs. 7(e)–7(g)] and normal tissue [Figs. 7(i)–7(k)], a more uniform blood flow profile across the rectal cross section is observed. For a more quantitative analysis, we averaged the Doppler flow velocity of the cancerous patients’ images in cancerous and normal regions, the treatment responder’s images in treatment scar and normal regions, and the normal patient’s images in normal regions, with a sample size of 6 for each patient and type of tissue. A Student’s $t$-test was used to evaluate the difference of P1 to P4’s cancerous or responding tumor bed, versus normal tissues. As shown in the box plot Fig. 7(m), there are clear statistically significant lower velocities in the cancer regions compared to normal regions in all three cancer patients. The responding tumor bed regions with scarring exhibited higher flow velocities than normal tissue regions, though not statistically significant. And even though other confounding factors including inter-patient variability and tissue heterogeneity influence flow velocities in patient rectums, the flow velocities of normal regions and responding tissue region are generally higher across all five patients.

Fig. 7

Examples of PA images coregistered with US, Doppler coregistered with US, and flow speed across a rectal cross-section. (a) Polar and (b) linear US coregistered with PA scan of a patient with cancerous tissue bordered in the red rectangle. (c) US coregistered with a Doppler scan of the same patient (d) Doppler flow speed bar chart across rectal tissue of the same patient, with cancerous tissue bordered in red. (e) Polar and (f) linear US coregistered with a PA scan of a patient with treatment scarring, bordered in green (g) US coregistered with Doppler scan of the same patient with treatment scar, labeled in green (h) Doppler flow speed bar chart across rectal tissue of patient with treatment scar, labeled in green (i) polar and (j) linear US coregistered with PA scan of patient with normal rectal tissue (k) US coregistered with Doppler scan of the patient with normal tissue (l) Doppler flow speed bar chart across rectal tissue of the same patient. (m) Box plot of average Doppler velocity across patients, analyzed with a $t$-test, the cancer regions exhibit statistically significant lower velocity than normal regions in all three cancer patients (P1–P3), and the responding tumor bed region exhibits slightly higher flow velocities than normal regions, but these are not statistically significant.