A light-fluence-independent method for the quantitative analysis of dynamic contrast-enhanced multispectral optoacoustic tomography (DCE MSOT).

MultiSpectral Optoacoustic Tomography (MSOT) is an emerging imaging technology that allows for data acquisition at high spatial and temporal resolution. These imaging characteristics are advantageous for Dynamic Contrast Enhanced (DCE) imaging that can assess the combination of vascular flow and permeability. However, the quantitative analysis of DCE MSOT data has not been possible due to complications caused by wavelength-dependent light attenuation and variability in light fluence at different anatomical locations. In this work we present a new method for the quantitative analysis of DCE MSOT data that is not biased by light fluence. We have named this method the two-compartment linear standard model (2C-LSM) for DCE MSOT.

Introduction

A variety of biomedical imaging modalities, including magnetic resonance imaging (MRI), computed tomography (CT), and positron emission tomography (PET), have been used to measure the Dynamic Contrast Enhancement (DCE) after intravenous injection of a non-specific contrast agent [[1], [2], [3]]. These DCE imaging methods can be utilized to evaluate vascular characteristics in patients with cancer and in mouse tumor models [4]. These vascular characteristics can be used to assess tumor angiogenesis and changes in tumor vasculature caused by anti-angiogenic and anti-vascular therapies [5]. DCE imaging methods can also gauge the vascular patency of some tumors to predict potential for intravenous drug delivery to the tumor. Therefore, DCE imaging methods have substantial value for pre-clinical and clinical diagnoses [6,7].

DCE imaging results can be analyzed with moderate-to-advanced pharmacokinetics models to quantitatively estimate K, the transport rate of a contrast agent from plasma to the extracellular extravascular space within a tumor [8]. These analysis methods can also quantitatively estimate k, the rate of transport of an agent from extracellular space to the plasma. However, accurate estimates of K and k require rapid image acquisitions to provide fast temporal resolution of the DCE results [9]. Unfortunately, routine MRI, CT, or PET protocols often have low temporal resolution, or sacrifice spatial resolution to gain temporal resolution [10]. Advanced MRI or CT protocols that can obtain fast temporal resolution while retaining fine spatial resolution often require expensive instrumentation and substantial expertise. For these reasons, DCE imaging methods have typically only been used for qualitative assessments of tumors rather than providing quantitative measurements of vascular characteristics.

In this work, we sought to introduce a new model for quantitative DCE MultiSpectral Optoacoustic Tomography (MSOT) that can provide DCE imaging with high temporal resolution and spatial resolution. This relatively new imaging technique directs pulses of light into tissue, which are absorbed by chromophores in the tissue. The absorption of light leads to transient heating and subsequent thermoelastic expansion, giving rise to ultrasound waves that can be detected by an ultrasound transducer [11]. The signal creation is rapid, allowing for DCE MSOT studies at high temporal resolution. Due to the detection of relatively low-scattered ultrasound waves, MSOT provides multi-spectral optical contrast at ∼2-3 cm depth and high spatial resolution, overcoming the traditional optical imaging limitation of reduced resolution in vivo [12,13]. Although MSOT has a lower depth of view relative to MRI and CT, the ∼2-3 cm depth of view is still sufficient for imaging the whole mouse torso and for imaging solid tumors in some patients who have cancer.

Although qualitative and semi-quantitative DCE MSOT studies have shown utility for specific applications [14,15], the development of a quantitative DCE MSOT analysis method has not yet been demonstrated. So far, the quantitative analysis of optoacoustic images has been aimed at correcting for optical light fluence, and calculating the absolute concentration of the contrast agent from the acoustic pressure maps formed by a reconstruction algorithm [16]. Furthermore, a true assessment of response to therapy with DCE MSOT requires a quantitative representation of the observed kinetics [17].

In this report, inspired by methods used in quantitative DCE PET and DCE MRI [[18], [19], [20]], we deviate from the current paradigm and present a new model for the quantitative analysis of DCE MSOT that does not require the estimation of the absolute concentration of a contrast agent in a certain pixel as an input, and can overcome the problem of light fluence [21]. We have called this model the two-compartment linear standard model (2C-LSM) for DCE MSOT.

Theory

Derivation of the two-compartment model for DCE MSOT

In this work we have adapted the notation and symbolic conventions previously described for DCE MRI and DCE PET studies assuming a two-compartment model (Fig. 1A) [22,23]. A differential equation describes the temporal evolution of a single contrast agent in a tissue of interest (TOI),where C is the concentration (nM) of the contrast agent in the TOI at time t, C is the concentration of the contrast agent in plasma at time t, K is the rate (min−1) at which the contrast agent leaves the plasma and enters the extravascular extracellular space of the TOI, and k is the rate (min−1) at which the agent returns from the extravascular extracellular space to the plasma [22]. Eq. (1) assumes that the concentration of the contrast agent is known. However, the optoacoustic signal S(λ,r,μ) does not directly represent the concentration of the contrast agent, and it is instead a linear combination of the product of light fluence and absorption scaled by the Grüneisen parameter [24],

Fig. 1

The two-compartment linear standard model (2C-LSM).

A. The model depicts a tissue of interest (tumor, C) and the arterial input function (plasma, C). B. Curves for C and C were simulated with an injection speed of 10 s, K = 0.25 s−1 and k = 0.64 s−1. The concentration is not directly detected by MSOT, but the optoacoustic signal is proportional to concentration.

A more rigorous version of Eq. (2) could be presented but all additional terms would eventually be canceled during our derivation [11]. In Eq. (2), Γ is the Grüneisen parameter, Φ(r,λ) is the light fluence (J/cm2), μ(λ,r) is the absorption coefficient (cm−1), while εj(λ) and Cj are the molar extinction coefficient (L mol−1cm−1) and concentration (mol L−1) respectively for each absorber j for a set of N total absorbers. The relative position and wavelength are given by r and λ repectively. Eq. (2) can be rewritten as Eq. (3) where the optoacoustic signal has been separated into its background and contrast agent components, assuming that pharmacokinetics studies with MSOT are acquired over several minutes at a fixed excitation wavelength and using a single contrast agent:

Assuming that the background optoacoustic signal and fluence remain constant over time in a given pixel for the duration of the experiment, Eq. (3) can be rearranged to obtain Eq. (4), which relates the change in optoacoustic signal as a function of time to the concentration of the contrast agent:

The approach described above can be used for C or C. Thus, replacing C and C in Eq. (1) by their corresponding equations derived from Eqs. (4b) yields (5), the two compartment model for MSOT:

This equation can be simplified further by multiplying both sides by ΓΦ(λ,r and removing the explicit dependency on λ, r, and t,where r and r indicate the location of the TOI and plasma, respectively. This simplification is possible because Γ remains fairly constant in the near infrared wavelength range of 680–980 nm [25]. ΔS and ΔS are the change in the optoacoustic signal in the TOI and plasma, respectively, and K is the true K scaled by Φ(r In the next section, we will describe how to make this model independent of Φ(r

Solutions of the 2C-LSM for DCE MSOT

Eq. (6) can be solved by obtaining its Laplace transform and solving for ΔS assuming that the change in optoacoustic signal at time zero (ΔS) is equal to zero, where s is the independent variable in the frequency domain:

The inverse Laplace transform is then applied,where ⊗ denotes a convolution.

The effect of fluence on DCE MSOT analysis

When Φ(r  = 1, then K reduces to K, resulting in estimated k and K that are unbiased and independent of Φ(r and Φ(r even though the analysis is performed using the change in optoacoustic signals instead of concentrations as in Eq. (1). We used Eq. (8b) when we performed simulations to evaluate the version of the two compartment model that does not remove the effect of fluence. When Φ(r the estimated k is unbiased, while the estimated K is scaled by the relative fluences (Φ(r) resulting in the biased estimator K. Unfortunately, correcting this biased estimator is not possible without modeling or measuring the relative fluence, which are tasks that are still under investigation [24]. Thus, a method that removes the effect of fluence on DCE MSOT would obviate the need for correcting this estimator.

A ratiometric approach to estimate fluence-independent PK constants for DCE MSOT

Our approach to remove the effect of fluence on DCE MSOT scales ΔS and ΔS by their corresponding maximums at time t (time to reach the maximum C or C signals):

Next, Eqs. (9) and (10) are solved for ΔS and ΔS and substituted into Eq. (8b). This substitution eliminates Φ(λ,r and Φ(λ,r from the model, and K is now scaled by a fluence-independent factor C:

As described in Eqs. (6a) and (6b), K=Φ(r and thus Eq. (11b) can be reduced,which contains the true K scaled by the ratio of the maximum concentration at C and CTOI. It is impossible to estimate C based on the experimental data. Thus, Eq. (11c) can be reduced,to obtain the two-compartment non-linear standard model (2C-NSM) for DCE MSOT, where K is K scaled by the ratio of the maximum concentrations of C and CTOI, in other words K = K.

Pharmacokinetics parameters

The pharmacokinetics parameters K and k can be calculated by non-linear least square curve fitting (NLSQ) of Eq. (12) to the normalized changes in the optoacoustic signal at the TOI and plasma. However, the analysis of dynamic imaging data by NLSQ has a slow processing time due to its iterative nature, and is highly dependent on the initial estimated parameters [18,20]. A linear version of Eq. (12) is expected to overcome these limitations [20], and can be obtained by applying the Laplace transform to Eq. (12) and multiplying both sides by (s + k:

Next, Eq. (13) is divided by s and the inverse Laplace transform is applied to the resulting equation,wherewhich represents the two-compartment linear standard model (2C-LSM) for DCE MSOT [26​]. Eq. (14b) denotes the linear algebra notation of Eq. (14a), where the rows of columns one and two in matrix X are the cumulative sums of the normalized changes in the optoacoustic signal at the tissue of interest (Δσ) and plasma (Δσ) respectively. Eq. (14) can be used to estimate the elements of b (K and k) using standard linear least-square curve fitting (LLSQ). The goodness-of-fit for all fitting was estimated using the coefficient of determination r [27​].

Materials and methods

Simulations

Computer simulations were used to test the effect of fluence, temporal resolution, and noise on the accuracy and precision of estimating Kapp and kep with the 2C-LSM for DCE MSOT. All simulations and data analyses were performed with MATLAB® 2015b (Mathworks Inc., Natick, MA) using a Windows 7 PC with a 3.5 GHz INTEL i5 processor, 16 GB of memory at 1600 MHz. All dynamic concentration curves were simulated using the 2C-NSM. The pharmacokinetics constants Kapp and kep were estimated using the 2C-LSM. The MATLAB function lsqnonneg was used for linear nonnegative least-squares curve fitting of the 2C-LSM to the simulated data. No bounds or initial estimates were needed for this linear fitting method.

Simulations of DCE MSOT data

The accuracy and precision of PK modeling is highly dependent on the pharmacokinetics profile of a contrast agent’s concentration in plasma (also known as an arterial input function or AIF). Therefore, we simulated AIF curves with three different injection speeds based onand the parameters listed in Table 1, which cover a practical range of experimentally validated AIFs [28,29]. Eq. (15) includes time constants that describe the uptake (B,E) and wash-out (C,F) of the agent in tumor tissue.

Table 1

Parameters of the simulated Arterial Input Function (AIF).

AIF name	Injection speed (sec)	A (nM sec−B)	B	C (min−1)	D (nM)	E (min−1)	F (min−1)
AIF-1	10	30.0	1.0	4.00	0.65	5.0	0.04
AIF-2	40	120.0	3.0	4.34	0.80	1.0	0.07
AIF-3	80	55.0	5.5	3.98	0.80	1.0	0.07

Each simulated AIF was then used to calculate the change in optoacoustic signal as a function of time for the TOI (ΔS):

Finally, all AIFs were also converted to changes in optoacoustic signal as a function of time:

Effect of injection speed and temporal resolution on the expected value of Kapp

The expected Kapp was simulated (without noise) as a function of temporal resolution (Δt) and injection speed (10, 40, and 80 s) using the three AIFs described in the previous section [8]. All curves were initially calculated from t = 0 to t = 5 min at Δt = 1.0 s, and subsequently undersampled from their initial Δt = 1 to Δt = 30 in steps of 1.0 s. Kapp was then estimated for each injection speed and temporal resolution, and compared to Kapp at Δt = 1.0 s.

Effect of temporal resolution, fluence and pharmacokinetics constants on the systematic errors in the estimated Kapp and kep

The Δt and Φ(r were explored as sources of systematic error when estimating K and k using simulated data that did not include noise. First, a series of S and S curves were calculated from t = 0 to t = 30 min at Δt = 1.0 s using Eqs. ((15), (16), (17)) at an injection speed of 40 s; with K = 0.1 to 0.5 min−1 in increments of 0.10 min−1; with k = 0.25–1.25 min−1 in increments of 0.25 min−1; and with Φ(r = 10.0 to 100.0 mJcm-2 in increments of 10.0 mJ cm-2. K and k where then estimated for each combination of Φ(r and K by fitting the simulated data to the 2C-LSM as described in Section 3.1.1. Next, S at Φ(r = 10 mJ cm-2, and S at Φ(r = 20 mJ cm-2 were undersampled from their initial Δt = 1.0 s to Δt = 15 s in steps of 1.0 s, and the fitting routine to the 2C-LSM was performed for each Δt. All simulations were repeated assuming an injection speed of 10 s.

Effect of the signal-to-noise ratio (SNR)

The AIF with the best performance determined by systematic error analysis was used to study the combined effect of SNR and Δt on the accuracy and precision of the 2C-LSM to estimate K and k. First, all simulated ΔS and ΔS curves were recalculated with Δt = 1, 5, 10, and 15 s, Φ(r = 50 mJ cm−2, and Φ(r = 25 mJ cm−2. Second, the MATLAB function awgn was used to measure the power of each curve and to add white Gaussian noise to produce signals with SNRs from 5 to 30 dB. The SNR was defined as SNR = 10*log10(P/P), where P and P are the power for the signal and noise respectively. Third, each simulation was repeated 1000 times for each condition of SNR and Δt with different white Gaussian noise.

The predictive performance of the 2C-LSM for the estimation of K and k was evaluated in terms of the root mean squared error (RMSE),where P is the jth estimated value of the ith parameter, P is the true value of the ith parameter, and N is the number of samples/repetitions (N = 1000).

Materials

Genhance 680 was purchased as a powder (Perkin-Elmer, Inc., Waltham, MA, USA), and diluted to a final concentration of 0.91 mM in saline buffer. All mice were injected intravenously through a tail vein catheter with 91 nmol or 4.5 μmol/kg mouse body weight, using 100 μL of the stock solution, followed by a flush of the catheter with 100 μL of saline. The total injection was performed over ∼30 s.

Mouse model

Three severe combined immune deficiency female mice weighing approximately 20 g bearing MDA-MB-231 mammary carcinoma xenograft flank tumors were housed and maintained under pathogen-free conditions in accordance with the guidelines of the American Association for Laboratory Animal Care. All experiments met the current regulations and standards of the U.S. Department of Agriculture, the U.S. Department of Health and Human Services, the National Institutes of Health, the Institutional Animal Care and Use Committee of the University of Louisville, and the University of Louisville Animal Care Facility. Each mouse was injected subcutaneously in the mammary fat pad with 1 × 106 MDA-MB-231 cells. MSOT studies were performed when tumors reached a diameter of ∼7 mm.

MultiSpectral optoacoustic tomography (MSOT)

Optoacoustic measurements were acquired using a commercially available MSOT inVision 256-TF animal scanner (iThera Medical GmbH; Munich, Germany). This instrument uses nanosecond near-infrared laser pulses at wavelengths between 680–980 nm and a 10 Hz repetition rate to illuminate an axial cross-section of a mouse. Optoacoustic signals were detected using a toroidally focused, ring shaped piezocomposite detector with 256 detection elements covering 270° around the sample and a center frequency of 5 MHz. Together with sensitive electronics that simultaneously sample all detector elements, the setup enables cross-sectional tomographic imaging with an in-plane resolution of 150 μm and a 0.1 s temporal resolution at the 10 Hz laser repetition rate.

To prepare for MSOT scans, a mouse was anesthetized with 1.5% isoflurane in oxygen gas. Using 100% oxygen gas increased the ratio of oxyhemoglobin to deoxyhemoglobin, which reduced the absorbance of 680 nm light by hemoglobin, thereby increasing contrast-to-noise during the DCE MSOT study. A 27 G catheter was inserted into the tail vein, and the mouse was allowed to acclimate in the MSOT scanner for 10 min. An anatomical localizer scan was acquired with 31 contiguous image slices, each with 1 mm thickness, collected around the expected location of the tumor using 800 nm excitation and 5 averages, which took ∼4 min to acquire. This scan was performed while waiting for the mouse to acclimatize to the temperature of the system, so that this ∼4 min localizer scan did not extend the scan session. After the mouse acclimatized to the temperature of the system, multispectral baseline images were acquired with three contiguous image slices, each with 1 mm thickness, centered on the tumor, and were collected after excitation at 680, 700, 715, 730, 760, 800, 850 and 900 nm, using 10 averages at each wavelength. DCE MSOT studies were performed by interleaving the acquisition after excitation at 680 nm (maximum for Genhance 680) and 900 nm (control frequency) for the same three contiguous slices, using 10 averages for each wavelength, resulting in a temporal resolution of 13.1 s/image. A wavelength of 900 nm was selected for the control frequency, because the extinction coefficient of Genhance 680 at 900 nm is negligible [30]. All DCE MSOT studies were initiated by collecting the same 3 slices at each wavelength 10 times before injection of Genhance 680. Data acquisition was continued for 30 min. A moving average was used to improve the contrast-to-noise of the pharmacokinetics profiles of the temporal changes in MSOT signal.

Image reconstruction and analysis

Images were reconstructed using standard filtered back projection along with impulse response correction with a pixel size of 150 μm, a field of view of 25 mm and a set of bandpass filters with a passband from 50 kHz to 6.5 MHz [31]. Laser energy fluctuations were automatically corrected prior to reconstruction. The 2C-LSM for DCE MSOT (Eq. (14b)) was implemented in MATLAB to estimate K and k in each voxel. DCE MSOT analysis without removing the effect of fluence (Eq. (8b)) was also performed with Matlab 2015b. Hemoglobin and deoxyhemoglobin maps were estimated using ViewMSOT (iThera Medical, GmbH, Munich, Germany) with the multispectral data acquired before the injection of Genhance 680.

Temporal correlation analysis

The 2C-LSM for DCE MSOT assumes that the contrast agent accumulation in one region or voxel did not reduce the optoacoustic signal of other voxels. Inspired by resting state functional MRI [[32], [33], [34], [35]], we used the temporal correlation of voxels and regions after injection of a contrast agent to determine if this assumption was fulfilled. A seed region was manually defined, followed by the calculation of the Pearson correlation coefficient (ρ) between the temporal changes of the average signal of the seed region and each voxel in the dataset,where ΔS(t) is the spatial-average of the change in optoacoustic signal in the seed region (normalized to the first time point), ΔS(t) is the normalized change in optoacoustic signal of a voxel, Cov is the time-dependent covariance, and σ and σ are the standard deviations of the normalized changes in optoacoustic signal in the seed region and voxel respectively. The MATLAB function corr was used to estimate ρ and the p-value under the null hypothesis that no correlation is present (ρ = 0). The p-value was estimated using a two-sided t-distribution test.

The tumor and AIF regions of interest were used as seed regions to identify voxels with a negative correlation coefficient that are not spatially contiguous, which would indicate that the change in optoacoustic signal at the tumor or AIF resulted in an effective reduction in the optoacoustic signal at the pixel being analyzed. Conversely, a positive correlation coefficient was an indication of similar dynamic behavior that was attributed to the pharmacokinetic uptake of the agent in the voxel and region.

Results

The 2C-NSM was used to simulate DCE MSOT data under various conditions. The 2C-LSM was used to estimate K and k by nonnegative least square curve fitting of Eq. (14b) to the simulated data. Noise-free data was simulated to study the systematic error in the estimated K caused by changes in injection speed and temporal resolution (Fig. 2A,B). These simulations were performed because K depends on C and C, and t may be dependent on injection speed and temporal resolution, and therefore K may be dependent on the injection speed and temporal resolution of DCE MSOT. However, our simulations demonstrated that K values estimated by 2C-LSM at injection speeds of 10, 40, and 80 s are equal to the true K if data are acquired every 7 s or faster (Fig. 2C). A systematic underestimation of K of up to 20% is observed for injection speeds of 10 and 40 and seconds when data are acquired every 15 s or slower. Notably, an injection speed of 80 s did not result in a systematic underestimation of K This result agreed with prior DCE MRI studies that show improved accuracy when using a “slow bolus” injection [18].

Fig. 2

K is independent of temporal resolution.

A. C was simulated at injection speeds of 10, 40, and 80 s, and a temporal resolution of 1.0 s. B. A similar simulation was performed for C. C. K was estimated by undersampling C and C from the initial Δt = 1.0 to 30.0 s in steps of 1.0 s. The ratio of the true K at each temporal resolution was compared to K estimated at Δt = 1.0 s.

A second simulation study was conducted to test the systematic error in the estimated K and k as function of fluence and pharmacokinetics constants. The estimated K and k are highly dependent on fluence when the optoacoustic signal is used for the pharmacokinetics analysis using Eq. (8b) that does not remove the effects of fluence (Fig. 3B,D). However, when the optoacoustic signal is analyzed with the 2C-LSM using Eq. (14b), the estimated K and k closely match the values of K and k used to generate the simulated data, regardless of the absolute values of K and k (Fig. 3A,C).

Fig. 3

The estimated K and k are independent of fluence when the 2C-LSM is used.

Optoacoustic signals were simulated with fluence at the AIF fixed to 20 mJcm−2, while fluence at the TOI was varied from 5 to 20 mJ cm−2. The dashed horizontal lines show the expected values for K (1.44, 1.86, 2.16, 2.42, 2.68 min-1) and k (0.25, 0.50, 0.75, 1.00, and 1.25 min-1) used to simulate the data. Panels A and C show the estimated K and k when the data was analyzed using the 2C-LSM (Eq. (14b)), while panels B and D show the estimated K and k when the data was analyzed without removing the effect of fluence (Eq. (8b)). The squares and circles are colored-coded, with Ktrans/Kapp/kep values for blue = 0.1/1.45/0.25 min-1; red  = 0.2/1.86/0.50 min-1; yellow = 0.3/2.16/0.75 min-1; purple = 0.4/2.42/1.00 min-1; green = 0.5/2.68/1.25 min-1. The symbols denote the estimated K and k values. Thus, if symbols fall on a dashed line, then the estimated values match the expected value.

The simulations in Fig. 3 were obtained at a temporal resolution of 1.0 s, which may not be attainable experimentally. Thus, a third simulation was performed to test the effect of temporal resolution in combination with fluence. The estimated K using the 2C-LSM was independent of fluence and was affected only by the injection speed (as expected; Fig. 4A,E). If the DCE MSOT data were analyzed without removing the effect of fluence, then the estimated K did not match the simulated K, and was underestimated by as much as 100% (Fig. 4B,F). The estimated k matched the expected k at high temporal resolution (≤5 s) regardless of the injection speed and was independent of fluence (Fig. 4C,G). k also showed a smaller systemic error than K when the data was analyzed without removing the effect of fluence (Fig. 4D,H).

Fig. 4

The estimated K and k are robust with respect to temporal resolution when the 2C-LSM is used.

Optoacoustic signals were simulated with Φ(rp) = 20 mJ cm−2, Φ(rTOI) = 10 mJ cm−2, and a temporal resolution of 1.0 sond/image. The data were then undersampled from their initial temporal resolution of 1.0 s to 15 s in steps of 1.0 s. The dashed horizontal lines show the expected values for K and k. Panels A, C, E, and G show the estimated K and k when the data were analyzed using the 2C-LSM that accounts for fluence (Eq. (14b)), while panels B, D, F, and H show the estimated K and k when the data were analyzed without consideration for fluence (Eq. (8b)). Panels A–D show results with 10-second injection speed, and panels E-H show results with 40 s injection speed. The symbols are colored-coded, with Ktrans/Kapp(10 s injection)/ Kapp(40 s injection)/kep values for blue = 0.1/1.45/1.00/0.25 min-1; red = 0.2/1.86/1.35/0.50 min-1; yellow = 0.3/2.16/1.60/0.75 min-1; purple = 0.4/2.42/1.83/1.00 min-1; green = 0.5/2.68/2.07/1.25 min-1. The symbols denote the estimated K and k values. Thus, if squares or circles fall on a dashed line, then the estimated values match the expected value.

To investigate the precision and accuracy of the model, the 2C-LSM was used to evaluate simulated data that varied in SNR (5, 10, 15, 20, 25, and 30 dB) and temporal resolution (Δt = 5, 10, 15, 20, and 40 s), assuming an injection speed of 40 s. The RMSE was determined for the estimated K and k values of 1000 simulations that included white Gaussian noise relative to the K and k used to simulate the data (Fig. 5). The RMSE in the estimated K and k decreased as the SNR increased. The effect of temporal resolution was also reduced as the SNR increased, with the exception of Δt = 40, which systematically underperformed compared to Δt = 5, 10, 15, and 20 s. These results indicated that high temporal resolution is critical to achieve reliable estimates at low SNR. For example, changing the temporal resolution from 5 to 20 s increased the RMSE of K by 50%, while the same change caused a 100% increase in the RMSE for k (Fig. 5A).

Fig. 5

The effect of the signal-to-noise ratio (SNR) and temporal resolution on the estimated K and k.

A. The root-mean-squared error (RMSE) of the predicted K as a function of SNR and temporal resolution of 5 (circle), 10 (square), 15 (diamond), 20 (triangle), and 40 s (cross). B. The RMSE of the predicted k at the same SNR and temporal resolution of panel A. C. Representative simulated data at SNR = 5 dB. D. Representative simulated data at SNR = 15 dB.

Experimental results

Based on our simulation results, we acquired experimental DCE MSOT data with three mice bearing xenograft flank tumors of MDA-MB-231 mammary carcinoma. Genhance 680 was injected at a speed of ∼30 s, and MSOT images were acquired every 13.1 s (Fig. 6). All DCE MSOT data was acquired interleaving excitation at 680 nm, the wavelength of maximum absorbance for Genhance 680 in vitro, and a control wavelength of 900 nm where Genhance 680 does not absorb light (Fig. 6A). Most voxels in the tumor were not enhanced 13.1 s after injection (Fig. 6C). At 4.9 min after injection, multiple voxels showed a 50% increase in the optoacoustic signal relative to the signal observed before injection, while a low enhancing region was observed in the center of the tumor. For comparison, the optoacoustic signal in the tumor remained constant after excitation at 900 nm (Fig. 6B). These results demonstrated that the observed changes in the optoacoustic signal were due to the accumulation of Genhance 680.

Fig. 6

DCE MSOT of Genhance 680 in a mouse model of breast cancer.

A. The absorbance spetrum of Genhance 680 shows strong absorbance at 680 nm and negligible absorbance at 900 nm. B. The change in average signals (relative to pre-injection signal) at 680 nm and 900 nm for the ROI were plotted for all time points. C. The tumor is shown in a red region-of-interest (ROI). D–F. Parametric maps of the signal in the ROI at 680 nm relative to the signal before injection. G-I. Parametric maps of the signal at in the ROI at 900 nm relative to the signal before injection (scale bar = 2 mm).

The quantitative analysis of DCE MSOT data using the 2C-LSM requires measuring the change in the optoacoustic signal in the AIF after injection of a contrast agent (Fig. 7A,B). The main vessels in the mouse were first localized by visualizing oxygenated and deoxygenated hemoglobin maps of all voxels. The average optoacoustic signal in the AIF followed a similar behavior to AIFs measured by DCE MRI and DCE CT [10,19], increasing by ∼65% within the first minute and decreasing to only a 15% enhancement relative to baseline in less than 5 min. This result demonstrated that the AIF can be measured using DCE MSOT. This behavior was observed only with excitation at 680 nm, but not with excitation at 900 nm. The average change in the AIF after excitation at 900 nm was 0.3% with a standard deviation of 3%. Thus, excitation at a wavelength for which the contrast agent does not absorb can be used as an internal control for delivery of the agent into the tissue.

Fig. 7

The temporal correlation of DCE MSOT signals can evaluate if fluence is affected by the contrast agent injection.

The 680 nm signal was used for this analysis. A. Each voxel within an ROI of the entire mouse (red) was compared to the average signal in an AIF seed region (white). The temporal Pearson correlation coefficients (r between the signals of each pixel and AIF seed region were estimated and displayed as a parametric map. Most pixels within the AIF seed region were highly correlated with the average AIF seed region signal (yellow pixels). This analysis also identified other regions that are highly correlated (r ≥ 0.80) with the average AIF (purple pixels). No pixels with r < 0.0 were found (scale bar  = 4 mm). B. The temporal changes in the average signal for the AIF seed region and the average signal for voxels with r ≥ 0.90 demonstrate this correlation. Signals are shown relative to pre-injection signal. C. Using the same analysis method with the tumor as the seed region (blue) yielded regions of positive correlation (yellow) and negative correlation (purple). Pixels with | r | ≥ 0.90 are shown in color. A negative correlation indicated that the accumulation of Genhance 680 in the tumor region reduced the optoacoustic signal in another region, and thus the fluence was not constant in these purple regions with respect to the tumor (scale bar  = 4 mm). D. The temporal changes in the average signal for the tumor seed region and the average signal for voxels with r ≤ −0.90 demonstrate this negative correlation.

Before applying the 2C-LSM for DCE MSOT, we verified our assumption that the background optoacoustic signal remained constant and that the amount of light delivered to a particular voxel remained the same before and after injection of a contrast agent. Our strategy is based on estimating the correlation coefficient of a seed region and every other voxel in the DCE MSOT image (Fig. 7). A negative correlation coefficient for two regions is an indication that one region is limiting the delivery of light into the other region.

This analysis was performed using the AIF as the seed region (Fig. 7A). No significant negative correlation was observed with any other voxel in the DCE MSOT image. Other major vessels showed a correlation greater than 0.90 with the AIF, which is an indication of similar dynamic behavior (Fig. 7B). Most importantly, all of the voxels in the tumor had a p-value greater than 0.15, and the median r value for all voxels was 0.49. This important result indicated that the assumptions made during the derivation of the 2C-LSM are fulfilled.

A similar analysis was performed using the average optoacoustic signal in the tumor as the seed region (Fig. 7C). Some contiguous voxels in the center of the mouse showed a negative correlation with the tumor, which suggested that the dynamic changes in the tumor potentially “masked” the signals in these interior tissue regions, causing a change in fluence in these interior regions (Fig. 7C,D). Thus, these regions with negative correlations relative to tumor could not be used in a quantitative DCE MSOT analysis.

We analyzed all voxels in the tumor region using the 2C-LSM (Eq. (14b)), and repeated the analysis without removing the effect of fluence (Eq. (8b)). The median r2 value for voxels in the tumor area was greater than 0.7 when removing the effect fluence (Fig. 8F), and the median r value for the same voxels was less than 0.55 without removing the effect of fluence (Fig. 8C). This result indicated that the normalized change in optoacoustic signal for DCE MSOT was well-described by 2C-LSM, while analyzing the optoacoustic signal without removing the effect of fluence was inferior. Furthermore, the parametric maps of K and k while removing the effect of fluence (Fig. 8D,E) were uniform compared to the parametric maps without removing the effect of fluence (Fig. 8A,B).

Fig. 8

Quantitative DCE MSOT analysis while removing the effect of fluence.

The same ROIs for the tumor and AIF shown in Fig. 8 were selected for quantitative analysis. A–C. Parametric maps of K, k, and r2 (goodness of fit) were determined without avoiding the effect of fluence. D–F. Parametric maps of the same parameters were determined using Eq. (14) that avoids the effect of fluence. Thresholding was applied to each parametric map which was overlaid on a grayscale optoacoustic image. The red dashed boundary shows the tumor ROI.

The estimated K and k values using a model that does not remove the effect of fluence were outside the range observed for tumors using DCE MRI [28], while the estimated K and k using the 2C-LSM lie within the range observed for tumors. For example, the mean k estimated by fitting data without removing the effect of fluence was 2.3 min−1, indicating that the mean transit time (MTT = 1/kep) of Genhance 680 in the tumor was 26.1 s and the contrast agent quickly returned to the blood, which was inconsistent with the observed mean signal in the tumor after injection of Genhance 680. The mean kep and MTT in the tumor region after removing the effect of fluence were 0.07 min−1 and 14.3 min, respectively. These estimates were consistent with the observed change in optoacoustic signal. Thus, consideration for fluence was needed to achieve physically meaningful estimates for these pharmacokinetics constants.

Discussion

We have presented the theoretical derivation and application of a new model, the 2C-LSM, that performs quantitative analysis of DCE MSOT studies. The pharmacokinetics analysis of DCE data required the fitting of the estimated concentration of a contrast agent as a function of time to a two-compartment model while removing the effects of fluence. The derivation of pharmacokinetics models for MRI, CT, and PET was achieved decades ago, and we have shown the equivalent model for quantitative DCE MSOT in this work.

Our work provides a theoretical foundation to interpret the analysis of dynamic optoacoustic imaging presented in the literature. First, we demonstrated theoretically and experimentally that using the absolute optoacoustic signal after injection of a contrast agent leads to unreliable estimates of the pharmacokinetics constants K and k of the two-compartment model used in this work. Second, Eqs. (6b) and (8b) indicate that a simple scaling of the OA signal relative to the signal before injection makes it possible to estimate k without knowing the fluence. Third, it is not possible to estimate the absolute K without knowing the true fluence, so instead we proposed a method known as the 2C-LSM for DCE MSOT to estimate a new PK constant, K If the change in the optoacoustic signal relative to the signal before injection is scaled, then the estimated K and k are not affected by differences in light fluence, but only by the SNR and temporal resolution of the dynamic data. In addition, our simulations provide a foundation for additional simulations that may aid in selecting experimental acquisition parameters for specific DCE MSOT studies performed in the future.

The performance characteristics of the optoacoustic imaging instrument used in this study were critical for accomplishing accurate DCE MSOT analyses. Rapid image generation with MSOT enables high temporal resolutions as fast as the laser repetition rate, performed at 10 Hz per wavelength in this study. We obtained sufficient SNR for our analyses with only 10 signal averages for each wavelength, or 2.0 s/image for 2 images with absorbance at 680 nm and 900 nm in this study, which resulted in rapid image acquisition. This detection sensitivity also facilitated the interrogation of a cross section of the entire mouse torso, which simplified the identification of the AIF within the interior of the mouse. Instrumentation with a more shallow depth of view may have more difficulty in obtaining the AIF signal within the mouse torso. Finally, the 270° tomographic detection around the mouse may aid in avoiding the “masking” of deeper voxels by uptake of agent in a small region of shallow voxels. Instrumentation with a linear array ultrasound detector that is placed nearest the shallow voxels will always suffer from this “masking” effect when detecting deeper voxels below the shallow voxels.

Our theoretical derivation included the assumption that fluence only changed in pixels where the contrast agent accumulated, and otherwise the endogenous fluence was constant, which is the primary potential limitation of this method. Other DCE imaging methods rely on the same assumption of static endogenous properties, such as a static T2 relaxation time and x-ray absorbance in endogenous tissues during DCE MRI and DCE CT, respectively. To test our assumption with DCE MSOT, we introduced a new method inspired by resting state fMRI to evaluate correlations in MSOT signal between pixels that accumulated the agent (i.e., the tumor tissue and AIF) and all other pixels. Our results showed that the signals of some voxels near the center of the field of view are negatively correlated with the signals in voxels near the surface of the mouse due to accumulation of a contrast agent and the resulting change in light fluence at depth. An unintended consequence of this analysis was the discovery of voxels with high positive correlation with the AIF that were not contiguous to the selected AIF, allowing the identification of several additional vessels. The DCE profile of almost all other image voxels were not correlated with voxels in the tumor region. Although this lack of correlation suggests that accumulation of agent in the tumor did not “mask” other voxels, the possibility still exists that an increase in MSOT signal in a deep pixel due to agent accumulation was exactly canceled by the “masking” due to agent accumulation in a more shallow pixel. Therefore, future studies are warranted to further explore this assumption of static fluence. In particular, future studies should investigate deeper tissues, such as orthotopic tumors or internal organs including the liver, where changes in fluence caused by injection and accumulation of an agent are expected to have a much greater impact on photoacoustic imaging signals.

Changes in mouse physiology may cause changes in fluence. To address this potential limitation, we ensured that the temperature of the mouse had reached the temperature of the water bath in the MSOT instrument. We also ensured that the mouse physiology was stable for 10 min before performing our DCE MSOT study. Our slow-bolus injection over 30 s, followed by a slow saline flush, has been shown to avoid large fluctuations in vessel volume and composition [18]. Our linear analysis methods also reduce the influence of first-pass kinetics due to a fast bolus injection, relative to non-linear analysis methods that are weighted towards rapidly changing regions of pharmacokinetics curves. We limited the DCE MSOT study to 30 min, which avoided longer-term physiological changes that are often encountered after on hour of anesthesia. Yet future studies should evaluate the stability of fluence to further confirm this assumption.

Additional experimental studies can be performed test the practical utility of the 2 C-LSM for DCE MSOT. For example, k has been used as a good predictor of response to neoadjuvant chemotherapy in breast cancer using DCE MRI. A similar study can be conducted to assess the availability of DCE MSOT to predict or detect response to chemotherapy. Additionally, the functional pharmacokinetics constants derived from DCE MSOT can be combined with multispectral imaging of hemoglobin to aid the selection of chemotherapy targeted to hypoxic regions of tumors. As another potential application, DCE MSOT studies may evaluate vascular perfusion in tissues other than tumors, including normal organs such as kidney and liver, or tissues undergoing wound healing. DCE MSOT can be compared with other imaging methods such as DCE MRI, which may further evaluate the merits of each modality.

Conclusion

In this work, we presented the derivation and experimental verification of a two compartment linear standard model (2C-LSM) for the quantitative analysis of dynamic contrast enhanced multispectral optoacoustic tomography (DCE MSOT). The 2C-LSM allows for the estimation of the apparent vascular permeability (K) and washout rate (k) in a tumor mouse model using linear regression, even when light fluence is variable. Thus, our 2C-LSM offers the possibility of combining molecular information from MSOT with functional information about malignancies in a single imaging session. This approach may lead to new opportunities to noninvasively assess tumor vascular characteristics and the response to anti-vascular therapies.

Conflict of interest

The author declare that there are no conflicts of interest.