4: Time-Series MSOT Data Analysis

Time-Series photoacoustic image analysis is implemented for two common imaging methods in PATATO. Oxygen-Enhanced imaging data and Dynamic Contrast Enhanced imaging data can be processed. Full details of the methods involved can be seen here: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5562224/. To briefly summarise, mice with subcutaneous tumours were imaged with MSOT and two separate scans were taken:

• Oxygen-enhanced (OE) scan: A single slice was imaged continuously. To start with, the breathing gas was 100% air and then changed to 100% oxygen after 5 minutes.

• Dynamic contrast-enhanced (DCE) scan: A single slice was imaged continuously. After five minutes of baseline imaging, a bolus of indocyanine green (ICG), a contrast agent, was injected intravenously.

Image reconstruction (backprojection) and linear spectral unmixing was applied. For the OE scan, spectral unmixing was applied to estimate the relative concentration of deoxyhaemoglobin and oxyhaemoglobin. For the DCE scan, spectral unmixing was applied to estimate the relative concentration of deoxyhaemoglobin, oxyhaemoglobin and ICG.

import patato as pat
import matplotlib.pyplot as plt
import numpy as np

Firstly, we load the example data provided with patato. The pre-calculated reconstructed images are extracted, along with the pre-drawn ROIs. Example data can be found in our example data repository here.

# Change this to be a OE dataset instead of a DCE dataset.
pa_so2 = pat.PAData.from_hdf5("data/so2-timeseries-data.hdf5")
pa_dce = pat.PAData.from_hdf5("data/icg-timeseries-data.hdf5") # We'll get the ROIs from here.
pa_so2.set_default_recon(("Reference Backprojection", "0"))
pa_dce.set_default_recon(("Reference Backprojection", "0"))
pa_dce.default_unmixing_type = "ICG"
pa_dce.external_roi_interface = pa_so2

# Get the oe reconstructions, so2 and delta icgs:
rec = pa_so2.get_scan_reconstructions()
so2 = pa_so2.get_scan_so2()
icg = pa_dce.get_scan_unmixed()[:, 2:] # ICG is the 3rd dataset in the unmixed data.

times = pa_so2.get_timestamps()[:, 0]
times -= times[0]
times /= 60

icg_times = pa_dce.get_timestamps()[:, 0]
icg_times -= icg_times[0]
icg_times /= 60

roi_tumour_right = pa_so2.get_rois()["tumour_right", "0"]
roi_tumour_left = pa_so2.get_rois()["tumour_left", "0"]
roi_reference = pa_so2.get_rois()["reference_", "0"]

Firstly, we show an example of the photoacoustic \(sO_2^{MSOT}\) estimate when breathing air and when breathing oxygen. A clear increase in the tumour oxygenation is observed on oxygen vs. on air.

fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(10, 8))

# Show initial slice
rec[0].imshow(ax=ax1)
im = so2[0].imshow(roi_mask=[roi_tumour_right, roi_tumour_left, roi_reference], clim=(0.2,0.8), ax=ax1)
ax1.set_title("Air")

# Show final slice
rec[-1].imshow(ax=ax2)
im = so2[-1].imshow(roi_mask=[roi_tumour_right, roi_tumour_left, roi_reference], clim=(0.2,0.8), ax=ax2)
ax2.set_title("Oxygen")

plt.colorbar(im, ax=[ax1, ax2], label="Unmixed $sO_2^{MSOT}$")

# Show initial slice
rec[0].imshow(ax=ax3)
im1 = icg[0].imshow(roi_mask=[roi_tumour_right, roi_tumour_left, roi_reference], ax=ax3, cmap="Greens")
ax3.set_title("Baseline")

# Show final slice
rec[30].imshow(ax=ax4)
im2 = icg[75].imshow(roi_mask=[roi_tumour_right, roi_tumour_left, roi_reference], ax=ax4, cmap="Greens")
im1.set_clim(im2.get_clim())
ax4.set_title("After ICG Injection")
plt.colorbar(im2, ax=[ax3, ax4], label="Unmixed $ICG$ [a.u.]")

plt.show()

We can also plot a trace of the blood oxygenation or ICG concentration over time in the three different regions.

def get_trace(time_series_data, region, filter_bad=False):
    mask, _ = region.to_mask_slice(time_series_data)
    timeseries = time_series_data.raw_data.T[mask.T].T[:, 0, :]
    # Unmixed SO2 can be out of the interval [0, 1] when there is low SNR.
    if filter_bad:
        timeseries[(timeseries>1) | (timeseries<0)] = np.nan
    return np.nanmean(timeseries, axis=-1)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4))

ax2.plot(icg_times, get_trace(icg, roi_reference), label="Reference")
ax2.plot(icg_times, get_trace(icg, roi_tumour_right), label="Right Tumour")
ax2.plot(icg_times, get_trace(icg, roi_tumour_left), label="Left Tumour")
ax2.set_xlabel("Time [minutes]")
ax2.set_ylabel("ICG [a.u.]")
ax2.set_title("Dynamic Contrast Enhanced (ICG)")

ax1.plot(times, get_trace(so2, roi_reference, True), label="Reference")
ax1.plot(times, get_trace(so2, roi_tumour_right, True), label="Right Tumour")
ax1.plot(times, get_trace(so2, roi_tumour_left, True), label="Left Tumour")
ax1.set_xlabel("Time [minutes]")
ax1.set_ylabel("$sO_2^{MSOT}$")
ax1.set_title("Oxygen Enhanced ($sO_2^{MSOT}$)")
fig.legend(*ax1.get_legend_handles_labels(), loc="upper center", ncols=3, bbox_to_anchor=(0.5, -0.01))
plt.show()

The built in GasChallengeAnalyser and DCEAnalyser are then applied to extract the change in the ICG content and \(sO_2^{MSOT}\) from baseline to peak.

gca = pat.GasChallengeAnalyser(display_output=False)
dso2, _, (baseline_so2, baseline_sigma_so2) = gca.run(so2, pa_so2)

dce = pat.DCEAnalyser(display_output=False)
dicg, _, [baseline_icg, baseline_sigma_icg] = dce.run(pa_dce.get_scan_unmixed(), pa_dce)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4))

# Draw delta sO2
rec[0].imshow(ax=ax1)
im = dso2.imshow(roi_mask=[roi_tumour_right, roi_tumour_left, roi_reference], ax=ax1, clim=(-0.04,0.1))
ax1.set_title(r"$\Delta sO_2^{MSOT}$")
plt.colorbar(im, ax=ax1, label=r"$\Delta sO_2^{MSOT}$")

# Draw delta ICG
rec[0].imshow(ax=ax2)
im = dicg.imshow(roi_mask=[roi_tumour_right,roi_tumour_left,roi_reference], ax=ax2)
ax2.set_title("$\Delta ICG$")
plt.colorbar(im, ax=ax2, label=r"$\Delta ICG$ [a.u.]")
plt.show()

By plotting \(\Delta sO_2^{MSOT}\) against \(\Delta ICG\), we can see that the oxygen enhanced experiment is a good metric of perfusion.

def get_pixel_values(delta_values, region, filter_bad=False):
    mask, _ = region.to_mask_slice(delta_values)
    delta = delta_values.raw_data.T[mask.T].T

    if filter_bad:
        delta[(delta>0.1) | (delta<-0.1)] = np.nan
    return delta

from scipy.stats import pearsonr

plt.scatter(np.log(get_pixel_values(dicg, roi_reference)), get_pixel_values(dso2, roi_reference, True), alpha=0.2, label="Reference")
plt.scatter(np.log(get_pixel_values(dicg, roi_tumour_right)), get_pixel_values(dso2, roi_tumour_right, True), alpha=0.2, label="Right Tumour")
plt.scatter(np.log(get_pixel_values(dicg, roi_tumour_left)), get_pixel_values(dso2, roi_tumour_left, True), alpha=0.2, label="Left Tumour")
plt.xlabel(r"$\log(\Delta ICG)$")
plt.ylabel(r"$\Delta sO_2^{MSOT}$")
plt.title(r"$\Delta sO_2^{MSOT}$ correlates with $\log (\Delta ICG$)")
plt.legend()
plt.show()

summary_measurements = pa_so2.summary_measurements(metrics=["so2", "dso2", "thb", "icg"])
# This is a pandas data frame containing measurements averaged over regions of interest.

No icg channel found
Skipping metric dso2
Skipping metric dso2
Skipping metric dso2