Applying dynamic contrast enhanced MSOT imaging to intratumoral pharmacokinetic modeling

Examining the dynamics of an agent in the tumor microenvironment can offer critical insights to the influx rate and accumulation of the agent. Intratumoral kinetic characterization in the in vivo setting can further elicudate distribution patterns and tumor microenvironment. Dynamic contrast-enhanced Multispectral Optoacoustic Tomographic imaging (DCE-MSOT) acquires serial MSOT images with the administration of an exogenous contrast agent over time. We tracked the dynamics of a tumor-targeted contrast agent, HypoxiSense 680 (HS680), in breast xenograft mouse models using MSOT. Arterial input function (AIF) approach with MSOT imaging allowed for tracking HS680 dynamics within the mouse. The optoacoustic signal for HS680 was quantified using the ROI function in the ViewMSOT software. A two-compartment pharmacokinetics (PK) model constructed in MATLAB to fit rate parameters. The contrast influx (kin) and outflux (kout) rate constants predicted are kin = 1.96 × 10−2 s-1 and kout = 9.5 × 10-3 s-1 (R = 0.9945).


Introduction
Multispectral optoacoustic tomography (MSOT) is a novel in vivo imaging modality with high spatial and temporal resolution at clinically significant detection depths [1][2][3]. Near-infrared (NIR) light pulses are absorbed by chromophores in the tissue, leading to heat production and thermoelastic expansion. This results in acoustic signal production which can be detected by an ultrasound transducer. MSOT imaging takes advantage of the optoacoustic (OA) effect with ultrasound detection of optical excitation [4]. NIR light with wavelength between 600-900 nm do not attenuate in most tissues, resulting in high resolution images at improved depth than optical imaging [1]. The use of OA overcomes the main limitations of optical and ultrasound imaging including limited tissue penetration (< 8 mm), and lack of potential for molecular information using a wide array of contrast agents [5,6].
Because each chromophore absorbs light and emits a distinct optoacoustic signal, it is possible to detect multiple chromophores simultaneously in high-resolution images as well as to obtain biological information, such as oxygenation state and tissue inflammation [2,[7][8][9][10]. This ability to simultaneously track multiple chromophores based on absorption spectra allows for one multispectral scanning session to yield multiple images depicting the biodistribution of each chromophore [11][12][13]. Thus, endogenous and exogenous agents can be used in combination to facilitate the identification of tumors. One hallmark of tumors is a region of tissue hypoxia from abnormal microvessel formation and decreased diffusion, leading to lack of vascular flow [14]. HypoxiSense 680 (HS680) targets and binds to carbonic anhydrase IX, a marker for tumor hypoxia [15]. HS680 absorbs at 670 nm and emits at 685 nm ( Fig. 1), which can be tracked using MSOT imaging along with oxygenated hemoglobin (HbO 2 ).
In addition to detecting distinct spectral signatures, optoacoustic imaging is an intrinsically fast imaging technique [16]. High temporal resolution allows the advantage of visualizing dynamic, biological events, which is difficult in ultrasound imaging due to motion artifacts [17][18][19][20]. One application is a pharmacokinetics (PK) study. Pharmacokinetic modeling maps out the kinetics of an administered drug or contrast agent in the plasma of an organism, thus playing a predictable role in the dose-response relationship in toxicity assessment [21][22][23].
Here we used a 2-compartmental PK model to simplify the body of an organism into one or two compartments, a commonly used approach [24,25]. in vivo data are collected and fitted against a PK model with estimated rate parameters.
Combining the functionality of contrast agents with a fast tracking imaging modality yields a powerful tool that can elucidate contrast kinetics and agent biodistribution in vivo [18]. DCE-MSOT has not been used previously to investigate PK with a contrast agent targeting tumor hypoxia [18,26]. This study utilizes HS680-MSOT contrast-imaging to study intratumoral pharmacokinetics in the in vivo setting.

Materials
The powder form of HS680 (Perkin-Elmer Inc., Waltham, MA, USA) was diluted to a concentration of 6.1 mM in saline solution. The absorbance and fluorescent emission of HS680 were measured using a FP-8300 spectrofluorometer (Jasco, Oklahoma City, OK, USA) and the OA spectra was measured using an MSOT inVision 256-TF (iThera Medical GmbH, Munich, Germany) within a tissue mimicking phantom as in previous studies ( Fig. 1) [8,27]. A tail vein catheter was set up to inject each mouse with 100 μL of HS680 and a concentration of 4.5 μmol/kg mouse and followed by a catheter flush with 100 μL of saline solution. Each injection took ∼30 s.

Mouse breast tumor model
Three female, athymic mice, weighing 20 g each, were injected with 1 × 10 6 MDA-MB-231 cells in the mammary fat pad. After injection, mice were housed under pathogen-free conditions following the guidelines of the American Association for Laboratory Animal Care. When tumors were 7 mm in diameter (45 days post injection), PK and MSOT imaging studies were performed in accordance with the current standards of the University of Louisville Animal Care Facility, the Institutional Animal Care and Use Committee of the University of Louisville, the National Institutes of Health, the U.S. Department of Agriculture, and the U.S. Department of Health and Human Services.

Multispectral optoacoustic tomography (MSOT) imaging
A commercially available MSOT inVision 256-TF (iThera Medical GmbH, Munich, Germany) was used for optoacoustic measurements. The MSOT emitted near-infrared laser pulses at 10 Hz and wavelengths between 680 to 900 nm A 256 element, toroidal, ring-shaped transducer received optoacoustic signals. This 270°transducer setup allowed for cross-sectional tomographic imaging with an in-plane resolution of 75 μm which has been further confirmed in our previous manuscript [28], and a center frequency of 5 MHz. This arrangement provides clear visualization of the spatial distribution of injected agents in the animal, avoiding the "masking" of deeper pixels by contrast uptake in a shallow region. In addition, implantation of the orthotopic breast tumor in the mammary fat pad avoided high signal noise level from internal organs.
A mouse was anesthetized with 1.5% isoflurane and 0.8 l medical air with 0.1 l O 2 gas in preparation for imaging. A 27 gauge catheter was inserted in preparation for tail vein injection. To limit physiological changes in the animal and manual introduction of vessel volume fluctuation, we ensured that mouse physiology was stable for 10 min while the mouse was located within the MSOT unit before data collection and performed bolus injection over 30 s followed by a slow saline flush [29]. MSOT signal data were collected 5 min prior to injection of HS680, providing a baseline reading in the tumor region that represents the hemoglobin level. Multispectral images were acquired with three contiguous image slices centered on the tumor for 20 min following injection of HS680 using three wavelengths 680, 760, and 900 nm and 5 pulses per wavelength.
Advanced Molecular Imager (AMI-1000X) was used to detect the biodistribution of HS680 in the mouse model 24 h after injection.

Image reconstruction and spectral unmixing
A standard back projection algorithm from the ViewMSOT 3.6 software suite (iThera Medical) was used to reconstruct the images [30]. The linear spectral unmixing method was used to detect and separate unique signals for contrast agents (i.e. hemoglobin and HS680) [31,32] by comparing the optoacoustic signals measured at each pixel with the known absorption spectra of the target contrast agent(s), allowing for separation and visualization of the distribution of each agent within the animal. In the current study, this method was utilized to distinguish between signal differences in HS680 and hemoglobin in the tumor and blood vessels.

Pharmacokinetic modeling
The PK model for this study is a 2-compartmental PK algorithm custom-designed using MATLAB. The two compartments represent the plasma and the breast tumor microenvironment (Fig. 2). The mathematical approach is similar to the 2-compartmental models previously utilized in previous DCE-MRI and DCE-PET studies [33,34]. A differential equation is used to describe the kinetics of HS680 within the tumor microenvironment (Eq. (1)). At a given time point, C tumor (t) is the concentration (nM) of HS680 in the tumor microenvironment; C plasma (t) is the concentration of contrast in the plasma. The rate constant k in (s −1 ) describes the rate at which HS680 leaves the plasma and is taken up into the tumor space and k out describes the rate at which the agent returns from tumor space into the plasma. The differential equation assumes a known working concentration, which is 0.91 mM of HS680.
This in vivo experiment measures optoacoustic signal of HS680. The absorbance and fluorescent emission of HS680 was determined using a spectrofluorometer. (B) The optoacoustic spectrum was determined by inserting HS680 into a tissue mimicking phantom and measured using MSOT. The optoacoustic spectrum was utilized to identify HS680 in vivo as detected using MSOT.
MSOT imaging provides continuous OA signal acquisition at selected near-infrared wavelengths using HS680, acquiring OA signal for background and contrast shown in Eq. (3). Then, OA signal for HS680 is isolated and substituted with the HS680 concentration expression in Eq. (2). Isolation of C HS680 derives an expression for HS680 concentration in relation to observed difference in OA signal (ΔS) in Eq. (4), which can be substituted into Eq. (1) (See Eq. (5)). Simplification of the resulting equation, along with a major assumption on light fluence, yields an expression for HS680 kinetics in terms of OA signal rather than HS680 concentration (Eq. (6)). Light fluence between the aorta and tumor "region of interest" can be considered to be the same because both are approximately the same distance from the surface of the mouse, 2.51 mm and 2.35 mm (Fig. 3).

HS HS
680 680 HS680 was administered in the form of a 100 μL single bolus, IV injection via the tail vein catheter (100% bioavailability) of the female athymic mouse. Contrast signal in the plasma was measured in the aorta and quantified with the "region of interest" function in ViewMSOT software, and referred to as the arterial input function (AIF). Given the quantification of AIF, these known values were used to solve contrast concentration within the tumor microenvironment by substituting AIF for ΔS plasma in Eq. (6). A MATLAB function was written to solve the differential equation (Eq. (6)) for the contrast tumor time-concentration profile (C tumor ) given AIF OA signal intensity values and estimated rate constant parameters. The experimental MSOT data for HS680 within the tumor was used to estimate compartment model rate constant parameters, which was performed with the nonlinear least squares (lsqnonlin) algorithm in MATLAB. Lastly, the algorithm calculates the correlation value (R value) for the PK model against the experimental C tumor data.

Isolating target absorber in images
Biodistribution of HS680 within the breast cancer mouse model using MSOT imaging was visualized via multispectral unmixing. This process separated the unique HS680 spectral signature peaked at 680 nm from the background signals, identifying HS680 signal and differentiating signal intensity within the animal. Cross-sectional images of the mouse were taken during the imaging process at the location of the breast tumor. The first row showed MSOT images of grayscale background. The second row showed the intensity of oxygenated hemoglobin (HbO 2 ) biodistribution within the animal. The third row showed the intensity of HS680 biodistribution. The third row showed an image overlay of HS680 signal gradient on top of a singlewavelength background (Fig. 4). The background optoacoustic image was produced by a single-wavelength signal at 900 nm. HbO 2 reading provided the optimal method to define the anatomical structures of the animal. The optoacoustic image overlay process visualizes HS680 biodistribution and signal intensity with respect to live mouse anatomy.

Pharmacokinetic imaging for intratumoral HypoxiSense 680 dynamics
MSOT imaging was selected for continuous imaging of HS680 dynamic within the tumor microenvironment of a living mouse. For the purpose of the intratumoral PK study, cross-sectional MSOT images were captured in the region of the tumor over time after intravenous injection of HS680. Spectral unmixing of the images allowed for isolation of contrast signals. Overlaying these signals on top of the singlewavelength grayscale background provided visualization of HS680 kinetics within the tumor over time. Signal increased in the tumor and the tumor environment as time progressed. Accumulation of HS680 was confirmed with near infrared fluorescent imaging using Advanced Molecular Imager (AMI-1000X) 24 h post-injection ( Supplementary  Fig. 1).
"Region of interest" (ROI) tool in the ViewMSOT software was used to define regions. The three regions identified (Fig. 3) represented the tumor, aorta, and a region outside of the mouse. Defining ROIs allowed for calculation of the mean contrast signal intensity within these regions using spectral unmixing. The ROIs shown have areas 3.01 mm 2 for the tumor ROI and 0.82 mm 2 for the AIF ROI. "Region of interest" of the AIF and tumor were drawn with the intention to keep the area of the corresponding ROIs constant in other mice. The empirical data points calculated from the tumor ROI represented the experimental concentration-time profile for HS680 within the tumor of the mouse model. A 2-compartmental PK algorithm was written using MATLAB in which Fig. 2. The pharmacokinetic model schematics. The basic structure of the PK model consists of two compartments: plasma compartment and tumor compartment. The two parameters, k in and k out , represent the rates at which a contrast agent fluxes in and out of the tumor compartment, respectively. Based on the pictorial schematic, the two-compartmental PK model assumes homogeneous distribution and constant volume in each compartment. Parameters are derived from curve-fitting the PK model curve to empirical MSOT data.  the compartments were used to account for blood plasma and tumor microenvironment. The model was fitted to the empirical data points by optimizing the relevant rate constants (Fig. 5).

Discussion
The non-invasive nature and capability to perform live imaging with MSOT allow advanced novel tracking capabilities of intratumoral contrast dynamic and biological processes within a native environment. Decreasing variability between animals combined with the individual AIF model provides more potential patient-specific applications in the future. MSOT image acquisition is in the range of 0.1 s while MRI image acquisition is in the range of 10 s to 100 s [36,37]. Significantly faster image acquisition makes MSOT a better imaging option for a PK study by providing more data points in a given period of time (Fig. 5). In addition, the lower level of photon scattering from acoustic waves contributes to high spatial resolution at clinically relevant depth compared to optical imaging, computed tomography, and MRI [1].
MSOT images of the same slice of tumor at different time points after HS680 injection (0, 5, 10, 15, 20, 25 min.). It demonstrates the unique advantage of collecting information in regard to multiple chromophores acquired during a single scan [12]. The overlay provides spatial confirmation of the location of the tumor as the HS680 signal is located where HbO 2 signal is depleted. In addition, visualization of the increased and sustained HS680 signal at the tumor location confirms the quantified signal data obtained for the PK study. This study with HS680 provides previously unknown quantification of intratumoral kinetics with hypoxia-targeting agents, which could be compared with in vivo kinetics of hypoxic-targeting therapeutics and enable further understanding of dosing and time-concentration profiles of therapeutic agents.
Another unique aspect of this investigation is the utilization of HS680 to target carbonic anhydrase IX, which regulates hydrogen ion flux in cell membrane and is a known tumor hypoxia marker [15]. This investigation identifies the HS680 kinetics and accumulation within the tumor in comparison to the plasma (Fig. 3). Selection of a hypoxia targeting contrast agent differentiates this study from previous PK studies performed with nonspecific blood-pool contrast agents, which are more suitable for applications such as vascular imaging [38] and lymph node mapping [39].
Selection of HS680 as the contrast agent contributes to the time delay observed (Fig. 5). The increased vascular permeability based on the enhanced permeability and retention (EPR) effect in most solid tumors indicates that diffusion is the main method of contrast transport into the tumor microenvironment [40]. Thus, size is a main determinant of kinetics within the tumor. HS680 consists of a modified acetazolamide conjugated to IR-680 fluorochrome, having a molecular weight of 1500 g/mol. In comparison, ICG has a molecular weight of 774.96 g/mol while gadolinium-based MRI contrast agents have molecular weights ranging from 36 to 480 g/mol [41]. Thus, HS680 diffusion into the tumor microenvironment is delayed compared with these commonly used smaller contrast agents. This possible explanation is affirmed by previous PK studies using ICG and MSOT imaging, but revealed no time delay between tumor and AIF signals [18,26]. This 2-compartmental PK model requires only two rate parameters, contrast influx into (k in ) and efflux out (k out ) of tumor microenvironment (Fig. 2). The model in this investigation was constructed using a similar AIF approach as the Flow-Limited Model for DCE-MRI [33], which calculated HS680 concentration within the tumor microenvironment. Avoiding possible error derived from overfitting due to unnecessary parameters is the main reason for selecting a 2-compartmental PK model over a more detailed and complicated PBPK model. A PBPK model consists of multiple compartments that correspond to different tissues in the body, which are not necessary in this specific kinetic modeling approach using an AIF [42].
The AIF approach is utilized in quantitative DCE-MRI analyses. It uses an input measurement to describe arterial concentration of a contrast agent. An image-derived AIF approach must satisfy two criteria: visibility in target field of view and sufficiently high temporal resolution [43]. Alternative approaches have produced similar results using a population AIF in preclinical setting [44] and clinical setting [45,46] or using the reference region (RR) approach [47][48][49], which replaces input AIF with a differential equation that quantifies contrast agent accumulation in a well-characterized tissue (i.e. muscle) over time. The individual AIF approach is selected for this study with MSOT because it generates the most accurate data for each individual animal and it satisfies both crucial criteria for using individual AIF approach. First, the aorta was selected because the vessel was present in the same field of view as the tumor (Fig. 3). More importantly, MSOT has the temporal resolution required to capture the rapid dynamic of the contrast agent for PK curve characteristics. Both the population AIF or RR approach are alternative approaches in DCE-MRI because of limitations to MRI acquisition time. MSOT avoids this temporal resolution limitation and can utilize an individual AIF approach. The modeling results of the AIF approach are shown in Fig. 5. The deviation between PK model curve and experimental data is likely due to the high noise level of the AIF (Fig. 5). The AIF model quantified signal from the ROI representing the mouse aorta, which is prone to changes in vessel dimension, resulting in high level of noise. Breathing could have contributed to the cyclic noise patterns through flux of oxygenated and deoxygenated hemoglobin. Because the AIF is used directly to calculate the tumor concentration of HS680 (Eq. (1)), noise in the AIF is amplified in the PK model curve. Noise does not significantly affect the initial spike because of the signal spike from the bolus injection. A motion compensation algorithm can be applied to offset the rhythmic fluctuation of noise [50]. Another possibility is to attempt the reference region (RR) model rather than the individual AIF model. However, the RR method yielded significantly more error than the AIF method in a similar type of study using DCE-MRI [51].
The characterization of intratumoral kinetics using MSOT is not only helpful in tracking a contrast or drug in a mouse model but is also meaningful for clinical applications. Pharmacokinetic models have been developed to predict anticancer drug dosing and tumor growth [52], aide in development and optimization of anticancer therapy [53], and have been used to characterize antineoplastic drug interactions in human prostate cancer models [54]. We expect increasing clinical applications for the combination of target specific contrast agent and MSOT imaging such as for drug delivery. Numerous studies have investigated the effectiveness of various antineoplastic drug combinations in treating malignancy [55][56][57]. Drug regimens can benefit from the knowledge of peak intratumoral kinetics of the anticancer drugs, thus allowing for optimization of the synergistic effects of drug combinations. Another critical potential area of clinical application is for intraoperative detection of tumor for surgical resection [1]. Accurate prediction of peak intratumoral kinetics of a contrast agent during surgical resection would allow for improved surgical efficacy, minimize unnecessary surgical procedures, and increase the chances of successful tumor resection.

Conflict of interest
All contributors are associated with Wake Forest Baptist Comprehensive Cancer Center and Wake Forest University School of Biomedical Engineering.
Ted Xiao completed his Bachelor in Science at the California Institute of Technology and is a MD/PhD candidate at Wake Forest School of Medicine. He is currently a graduate student in the Virginia Tech-Wake Forest University School of Biomedical Engineering and Sciences. He is a member of the McNally Lab and focuses on nanoparticle construction and MSOT imaging in his ongoing PhD thesis. He is interested in effective cancer treatments and detection methods through application of nanotechnology and imaging techniques.   He is the author of over 200 peer-reviewed publications and has been awarded three patents on medical image analysis.
Umit Topaloglu, Ph.D., is the Associate Director for Informatics at the Wake Forest University Comprehensive Cancer Center and an Associate Professor of Cancer Biology and Biostatistics. He develops and implements a Biomedical Informatics Vision that supports and advances the clinical and translational research needs around the scalable, interoperable, findable, sharable data frameworks. He has worked on pathological and radiological image management and secure access including leading statewide systems for trauma and telemedicine network implementation. Additionally, he has focused automated or semi-automated identification of features utilizing variety of machine learning and deep learning methods in spatio-temporal domain, which has the promise for detecting agents to accomplish the proposed tasks. McNally's laboratory develops actively targeted contrast agents, theranostic nanoparticles, and novel treatments with the overall goal to improve patient outcomes. She uses multispectral optoacoustic tomography to track contrast agents and nanovehicles in both the pre-clinical and clinical setting.