HemoSYS: A Toolkit for Image-based Systems Biology of Tumor Hemodynamics

Abnormal tumor hemodynamics are a critical determinant of a tumor’s microenvironment (TME), and profoundly affect drug delivery, therapeutic efficacy and the emergence of drug and radio-resistance. Since multiple hemodynamic variables can simultaneously exhibit transient and spatiotemporally heterogeneous behavior, there is an exigent need for analysis tools that employ multiple variables to characterize the anomalous hemodynamics within the TME. To address this, we developed a new toolkit called HemoSYS for quantifying the hemodynamic landscape within angiogenic microenvironments. It employs multivariable time-series data such as in vivo tumor blood flow (BF), blood volume (BV) and intravascular oxygen saturation (Hbsat) acquired concurrently using a wide-field multicontrast optical imaging system. The HemoSYS toolkit consists of propagation, clustering, coupling, perturbation and Fourier analysis modules. We demonstrate the utility of each module for characterizing the in vivo hemodynamic landscape of an orthotropic breast cancer model. With HemoSYS, we successfully described: (i) the propagation dynamics of acute hypoxia; (ii) the initiation and dissolution of distinct hemodynamic niches; (iii) tumor blood flow regulation via local vasomotion; (iv) the hemodynamic response to a systemic perturbation with carbogen gas; and (v) frequency domain analysis of hemodynamic heterogeneity in the TME. HemoSYS (freely downloadable via the internet) enables vascular phenotyping from multicontrast in vivo optical imaging data. Its modular design also enables characterization of non-tumor hemodynamics (e.g. brain), other preclinical disease models (e.g. stroke), vascular-targeted therapeutics, and hemodynamic data from other imaging modalities (e.g. MRI).


Multicontrast optical imaging system
The imaging system employed three illumination sources: FL excitation of GFP was achieved via a 473 nm blue laser (100 mW, Cobolt AB, Sweden); a white light source (NI-150, Nikon Instruments Inc., NY) coupled with two bandpass filters (570±5 nm and 600±5 nm, Thorlabs Inc., NJ) provided illumination for imaging hemoglobin absorption based IOS, while a 632.8 nm He-Ne laser (0.8 mW, Thorlabs, NJ) was used for LS imaging.

Computation of hemodynamic variables (microvascular oxygen saturation, blood volume and blood flow)
For the computation of hemodynamic quantities based on IOS images, we first calculated a median intensity image from the 10 images acquired under 570 nm and 600 nm illumination, respectively, at each imaging time point. We repeated this for each image acquisition. We linearly resampled the median intensity image stacks for each wavelength with a fixed 30 second time resolution using a MATLAB script (Mathworks, MA). Next, using the ImageJ 1 MultiStackReg plugin 2 , we co-registered images in the time-series to the first image for each wavelength. Then, we used the 570 nm and 600 nm IOS images in a formulation of the modified Beer-Lambert law to compute the microvascular oxygen saturation (Hbsat) and total hemoglobin concentration (HbT) 3 . Briefly, concentrations of oxy-hemoglobin (HbO) and deoxy-hemoglobin (Hb) were calculated by solving the following system of equations: denote the absorption coefficients of Hb and HbO. These were calculated by using the tabulated molar extinction coefficients found in 4 , averaging the extinction coefficients over a ± 4 nm bandwidth around the two wavelengths (to account for finite filter bandwidth), and multiplying by a factor of 2.303 to convert extinction coefficients into absorption coefficients 4 . L(λ) is the differential path length at each wavelength λ. For simplicity, L(570 nm) was assumed to be unity.
A nominal value of 4 was used for L(600 nm).
Both the multiplication factor for incident light level calculation and L(600 nm) were estimated using a Monte Carlo simulation 5 of photon migration in a simple air-tissue model. Optical properties of the tissue were modeled using hemoglobin extinction coefficients as described above, an Hbsat level of 50% to indicate poor oxygenation status in tumor tissue, and an estimated anisotropy of 0.9 6 . A scattering coefficient of 10/mm was used. We first assumed the tissue to contain no blood vessels (i.e. a zero HbT level) to model tissue void of blood vessels. Then, using the Monte Carlo simulation, we computed the ratio between the amount of photons entering the tissue after specular reflection and the amount of photons diffusely reflected back from the tissue.
We obtained a value of 2.3 which was approximated to 2 to account for uncertainties in modeling parameters. We used this as the multiplication factor to convert the imaged intensity of light from a selected tissue area void of visible blood vessels into an estimate of the incident light level.
Next, we used an HbT level of 80 μM 6 to model tissue with blood vessels. Since we already denoted L(570 nm) as unity, L(600 nm) was estimated as the ratio of differential path lengths computed using the Monte Carlo simulation for 600 nm and 570 nm photons, respectively. The computed ratio of 3.6 was approximated to 4 to account of uncertainties in modeling parameters.
One could perform more detailed Monte Carlo simulations of photon migration in tissue to potentially arrive at more accurate estimates.
Next, we computed Hbsat and HbT according to: Hb sat = HbO HbT , HbT = HbO + Hb Assuming constant hematocrit, HbT served as a surrogate of blood volume (BV). We also computed images of relative blood flow (BF) from LS images by quantifying the speckle contrast (k) for each pixel. Briefly, speckle contrast is given by 7 : Here, µ and denote the mean and standard deviation of pixel intensities across the stack of images (40 images, acquired every 30s under 632.8 nm illumination) at each pixel location (x, y).
Under coherent laser illumination and orderly flow, speckle contrast (k) is approximately related to blood flow as 7 : These BF images were resampled to a 30 second time resolution, and co-registered to the IOS images.

Generation of tumor masks
The tumor extent was manually identified using FL images. Corresponding BV images were then used to identify an in-focus FoV for each tumor. These FoVs were used to generate unique masks for each tumor xenograft (Supplementary Fig. 4) which were used for subsequent analyses.

Temporal and spatial filtering of multicontrast optical signals
We applied a mean temporal filter (1 minute step size, 5 minute long kernel) to reduce noise levels for each of the three hemodynamic variables BV, Hbsat and BF. This also ensured that all HemoSYS modules operated on hemodynamic time-series with 1 minute temporal resolution. Moreover, the 5 min filter kernel reduced the number of time points available for analysis. For example, a 60 min imaging session yielded 56 time points at 1 min temporal resolution, while a 30 minute session yielded 26 time points at the same temporal resolution.
In addition, all hemodynamic images were subjected to a mean spatial filter with a discrete spatial kernel of 5050 pixels ensuring that all HemoSYS analyses were performed in 5050 pixel sub-regions. This helped minimize the computational burden and provided resilience against moderate motion that cab occur in a soft tissue vascular bed such as that of the mammary fat pad. (The propagation analysis module of HemoSYS was the only exception, and used hemodynamic images that underwent mean filtering with a continuous spatial kernel of 5050 pixels). illumination source was manually switched ON/OFF as it was used only once at the beginning of the experiment. The IOS and LS illumination sources were autonomously gated. The IOS illumination source was fitted with a filter wheel containing two bandpass filters (570 nm ± 5 nm and 600 nm ± 5 nm) and a light stop slot. The LS illumination source was fitted with a custommade light stop switch using a microcontroller actuated solenoid. These were controlled by the master control program running on a personal computer (PC). For clarity, these two controls are not indicated in the schematic. Light from these illumination sources was shone onto the animal imaging window. A lens set was used to focus images on a CCD image sensor. Additionally, a 496 nm longpass filter prevented blue excitation light from entering the imaging light path. The master control program also controlled image acquisition and storage. Dashed arrows indicate light paths, and the hollow bi-directional arrows indicate digital control and data transfer. Supplementary Fig. 3: Schematic of the multicontrast pre-processing pipeline.  where analyses were restricted to selected discrete spatial locations of the TME. The green colored raw shows a single instance where analysis was performed on data that was continuous in the spatial and temporal domains. However, this analysis was limited to oxygenation dynamics.