Volumetric tumor delineation and assessment of its early response to radiotherapy with optical coherence tomography

: Texture analyses of optical coherence tomography (OCT) images have shown initial promise for differentiation of normal and tumor tissues. This work develops a fully automatic volumetric tumor delineation technique employing quantitative OCT image speckle analysis based on Gamma distribution fits. We test its performance in-vivo using immunodeficient mice with dorsal skin window chambers and subcutaneously grown tumor models. Tumor boundaries detection is confirmed using epi-fluorescence microscopy, combined photoacoustic-ultrasound imaging, and histology. Pilot animal study of tumor response to radiotherapy demonstrates high accuracy, objective nature, novelty of the proposed method in the volumetric separation of tumor and normal tissues, and the sensitivity of the fitting parameters to radiation-induced tissue changes. Overall, the developed methodology enables hitherto impossible longitudinal studies for detecting subtle tissue alterations stemming from therapeutic insult.


Introduction
There is emerging evidence that detailed analysis of diagnostic and therapy-monitoring image characteristics in oncology may provide valuable insights into cancer detection and prognosis of clinical outcome. Hundreds of first-, second-and higher-order quantitative features, like tumor shape or image texture descriptors can be extracted from computed tomography (CT), magnetic resonance (MR) and positron emission tomography (PET) images and analyzed by the radiomics approach [1][2][3][4][5][6][7][8]. These features in radiomics are compared to one another to find any similarities and to correlate with the clinical, molecular or genetic data. In early cancer stages, however, when tumors are smaller or close in size to MR, CT or PET-CT imaging resolution (millimeter scale), diagnostic tests do not detect skin and other cancers very well. In cases of highly metastatic cancers (e.g., pigmented melanoma), it is crucial to diagnose and initiate the treatment before the cancer has spread, therefore microscopic non-invasive imaging modalities like optical coherence tomography (OCT) are gaining increasing preclinical and clinical use [9,10].
OCT was introduced three decades ago as an effective label-free high-resolution 3D imaging technique with a moderate penetration depth of 1-3 millimeters [11]. OCT allows in situ, safe, noninvasive real-time investigation of micromorphology, pathology, and often functionality without tissue removal. Its image analysis methods are continuing to be developed for early cancer detection, tumor delineation, and noninvasive optical biopsy of dangerous aggressive cancers [12]. Many of these methods rely on OCT speckle, which can cause image noise artefacts but also contains useful information about underlying tissue structures.
Image speckle originates from the interference of backscattered light from sub-resolution scatterers in the tissue being imaged [13]. Although often referred to as a "speckle noise", it is not random but pseudo-deterministic. That is, for a given spatial distribution of scattering tissue structures, the speckle pattern carries information of the underlying scatterers' average properties, such as size, shape, and distribution, without individually resolving them [14]. There is emerging evidence that image analysis of OCT speckle patterns can extract additional and valuable tissue information, potentially enabling detection and delineation of various tissue types [15].
Numerous first order (mean, variance, kurtosis, energy, etc.) and second order (correlation, inertia, entropy, inverse difference, etc.) texture features have been proposed in the literature utilizing various approaches (principal component analysis [16], generalized Gamma (G) distribution [17], K-distribution [18], morphological analysis [19], among others) to extract useful information from OCT images. Most of these approaches to date have been limited to oneand two-dimensional analyses, working with axial and cross-sectional scans and avoiding the more challenging volumetric analysis problem. One of the most important desired outcomes of these explorations is a robust method to differentiate between different types of tissues (e.g., tumor vs normal) for early disease detection, diagnosis and treatment response monitoring.
Among the variety of texture analysis methodologies, fitting the speckle pattern statistics to the generalized G-distribution has previously demonstrated excellent sensitivity to changes in the tissue scattering properties using high-frequency ultrasound [20,21] and OCT [17,22]. These approaches examine the signal intensity distributions within a given speckle pattern region of interest (ROI) in the image, fitting it using a non-linear regression least-squares method to the G-distribution function [23]: where k and θ are referred to as the shape and scale parameters, respectively, of the G-distribution, and I is the signal intensity. Within a given speckle pattern, kθ represents the mean of the distribution, essentially the "average signal brightness" [24]. In ultrasound research, parameter k was found to be related to the effective tissue scatterer number density [21]. θ was considered to be proportional to backscattered signal intensity I and related to cellular changes (e.g., normal size nuclei of healthy cells vs increased size nuclei of tumor cells) that affect the scattering properties of tissue. In our recent OCT microsphere phantom and Monte-Carlo simulation study [25], fitting parameter θ decreased with the increasing scattering strength (cross section) of the individual scatterer while parameter k showed essentially concentration-independent behavior, in contrast to ultrasound-based findings in ref. [21]. We've also reported that a spatial speckle analysis algorithm based on G-distribution fit signal histograms shows initial promise for differentiation of normal and tumor tissues [17]. Farhat et al. found a strong correlation of coefficients of this fit with the survival rate of acute myeloid leukemia (AML) cells treated in-vitro with the chemotherapeutic agent cisplatin that caused apoptosis [26]. In Pires et al. [22,27], we have quantitatively demonstrated the high sensitivity of this method to changes in OCT signal scattering in melanomas following tissue optical clearing agent application [22] and photodynamic therapy [27]. Building upon our previous OCT research and that of others [17,18,22,25,26,27], and supplemented by high-frequency ultrasound [20,21] and magnetic resonance imaging [28,29] studies that attempt to utilize the G-distribution fitting, here we develop an automated and user-independent objective methodology for volumetric delineation of tumor and normal tissues otherwise not distinguishable in structural OCT images. We validate this new approach in two important pre-clinical xenograft models, using epi-fluorescence microscopy, combined photoacoustic microscopy / high-frequency ultrasound imaging and histology. Further, we demonstrate the utility of automatic 3D tumor segmentation and Gamma-fit parameter extraction in a longitudinal treatment monitoring assessment pilot study of early radiobiological response of tumors to a 15Gy single dose of ionizing radiation therapy (RT). All animal procedures were performed in accordance with appropriate standards under protocols approved by the University Health Network Institutional Animal Care and Use Committee in Toronto, Canada (animal use protocols #3256 and #4401). Two important and representative tumor models were studied, towards demonstrating the generalizability of the method: dorsal skin window chamber grown pancreatic adenocarcinoma (to explore the developed method performance in a well-established controlled environment) and cutaneous pigmented melanoma (to test the methodology in bare skin imaging conditions).
Human BxPC-3 pancreatic cancer cells labeled with Discosoma sp. Red Orange Mushroom (DsRed) derived red fluorescent proteins were purchased from AntiCancer Inc. (San Diego, CA, USA) and cultured in Roswell Park Memorial Institute 1640 medium supplemented with 2 mM L-glutamine, 10% fetal bovine serum and 1% Penicillin-Streptomycin at 5% CO 2 and 37°C. Tumors were generated by injection of 3×10 5 cells prepared in 10 µL of 1:1 PBS:Matrigel (BD Biosciences, ON, Canada) solution into the dorsal skin of seven-to eight-week-old NRG mice using 30Ga needles. The dorsal skin window chamber (DSWC) surgery was performed several weeks post-injection of BxPC-3 cells (after tumors reached 3-5 mm diameter as shown in Fig. 1(b),(c)). Titanium-framed window chambers (APJ Trading Co., Inc., Ventura, CA, USA) were surgically implanted into the dorsal skinfold of anesthetized (mixture of 80 mg/kg of ketamine and 5 mg/kg of xylazine) mice following the procedure described previously [30]. Animals were kept in micro-isolator cages with access to food and water ad libitum. Imaging was performed after a recovery period of three to five days post-DSWC installation. Tumors were irradiated after full vascularization was confirmed with the OCT microangiography.

OCT imaging
Tumor-bearing mice were imaged with the swept-source OCT system based on a quadrature interferometer to suppress the complex conjugate artifact [31][32][33]. Its light source (HS2000-HL, Santec, Japan), with 20 kHz rotating-polygon-based tunable filter and wavelength centered at 1320 nm, had a sweep range of 110 nm and an average output power of 10 mW. The axial and lateral resolutions in air were 8 and 15 µm, respectively.
For OCT imaging, animals were anesthetized with either isoflurane (5% for knock-out, 2% for maintenance) or ketamine/xylazine mixture and placed under the imaging probe on a heating pad (melanoma-bearing mice) or a custom 3D printed animal restrainer with built-in heater (animals with window chamber). Imaging ( Fig. 1(d),(e)) was performed in conventional OCT speckle variance (svOCT) mode [34] to get microstructural and derive microvascular images of tissue (field of view 6 mm x 6 mm, with ∼1.5 mm in depth), with eight microstructural images being taken from the same location with 25 ms inter-frame interval. Speckle variance contrast of 10% above noise floor was considered as blood microvasculature. For visual representation, microvascular images were depth-encoded in RGB color space [35] with 256 color gradations (green = top layers below the tissue surface and black = deepest tissues). Matlab (Mathworks, MA, USA) software was used to process and visualize the data.

Radiation therapy
Prior to treatment, two window chamber bearing animals were anesthetized with isoflurane delivered through a mask. Tumors were irradiated using a commercial small animal X-ray micro-irradiator system XRad225Cx (Precision X-Ray Inc., North Branford, CT, USA) through the custom-built 12×6mm 2 semicircular collimator ( Fig. 1(f)). The innovative half-chamber irradiation geometry was chosen to enable direct side-by-side comparison of full-dose tissue response versus adjacent un-irradiated control in the same tumor. The X-ray system delivered single focal radiation beams (225 kVp, 13 mA, added filtration of 0.32 mm Cu) directly to BxPC-3 tumors, with a dose rate of 2.63 Gy/min to the total dose of 15Gy. Treated area and dose levels were confirmed with calibrated Gafchromic EBT-2 film (ISP Inc., Wayne, NJ, USA) consisting of a radiosensitive monomer that polymerized and changed color with absorbed dose (Fig. 1 g,h).

Ultrasound imaging, photoacoustic tomography, fluorescence microscopy and histology
Cutaneous melanoma tumors were imaged with an ultrasound Vevo 2100 system (FUJIFILM VisualSonics Inc., Toronto, Canada) operating at 32-56 MHz. The system incorporated a tunable pulsed laser (680-970 nm, 6 ns pulse length, 20 MHz pulse repetition rate) to enable photoacoustic imaging at 700 nm that was combined with the ultrasound images (photoacoustic signal confirmed the presence of melanin throughout the pigmented tumor volume). Images were fused together and analyzed using ImageJ software (NIH, Bethesda, MD, USA). DsRed-labeled pancreatic adenocarcinoma cell fluorescence images (535 nm excitation, 580 nm emission) were obtained using an epi-fluorescence microscope (Leica MZ FLIII, Leica Microsystems, Richmond Hill, ON, Canada) with consistent exposure settings to map out the general shape (two-dimensional projection) of the tumors. For histological staining, animals were euthanized by isoflurane followed by cervical dislocation. Tumors were resected, fixed in 10% formalin and processed for hematoxylin and eosin (H&E) staining to view cellular morphology. Slides were scanned by Aperio Scanner and viewed using Aperio ImageScope software (Leica Biosystems, Concord, ON, Canada).

OCT image processing and texture analysis for tumor delineation
For tumor delineation, images were first prepared for analysis. Figure 2 shows cross-sectional OCT images (B-scans) of BxPC-3 tumor grown in the dorsal skin window chamber (left column of panels) and B16-F10 tumor grown in the flank of a nude mouse (right column).
In the case of the window chamber, images ( Fig. 2(a)) were prepared for further processing by correcting for the OCT probe tilt and removing the window chamber glass ( Fig. 2(b)) and compensating for exponential depth-decay of the OCT signal ( Fig. 2(c)). In our earlier proof-of-principle study [17], only one depth of an OCT cross-sectional image (B-scan) was chosen for statistical analysis to minimize complicating effects due to optical signal attenuation with depth. Attenuation possessed considerable challenge to two-and three-dimensional tumor delineation method development and required accurate compensation to keep the relative intensities comparable between hundreds of B-scans in a volumetric tissue image. For this, phantom-based measurements of attenuation coefficients specific to our particular OCT system were performed following the procedure described in detail in [36]. For cutaneous melanoma case (pre-processed OCT image in Fig. 2(d)), the skin surface was defined with a sliding Gaussian filter (kernel 5 × 5 pixels) as labeled with the red curve in Fig. 2(e), with the signal depth-decay compensation performed from the surface of the tissue into deeper layers (Fig. 2(f)).
At the next processing step, 3D images were divided into voxel-by-voxel sliding volumes of interest (VOI) as demonstrated on the structural OCT image in Fig. 3(a). The sliding VOI was 24 × 24 × 3 voxels (lateral x × lateral y × depth z, corresponding to a physical size of 180 × 180 × 20 µm 3 ), chosen large enough to allow for robust histogram statistics but sufficiently small to retain tissue spatial heterogeneity information. The OCT pixel intensity distribution in the VOI was then represented as a histogram, with the variable bin size (ranging from 20 to 100 intensity units) to find the optimal least squares histogram fit (goodness of fit R 2 >0.95) by the Gamma distribution function [Eq. (1)]. Although depth-decay compensation of the signal likely caused changes in its statistical distribution due to amplification of the overall noise level, it affected both tumor and normal regions in a similar "offset" fashion thus causing a global relative change without influencing the delineation. Typical experimental histograms plotted from normal (skin) and tumor (pancreatic adenocarcinoma) VOIs within OCT tissue volumes are shown in Fig. 3(b), exhibiting differences in the distribution of voxel intensities and fitting parameter values k and θ. In total, 120,000 VOIs were analyzed from ten (n=10) BxPC-3 and ten (n=10) B16-F10 tumors. Using fluorescence images, photomicrographs and microvascular 3D images as guides, 300 × 300 × 115 voxels volumes were chosen in each OCT data set and categorized into 6,000 tumor VOIs and 6,000 normal tissue VOIs. Analysis showed that the values of the fitting parameter k overlapped for tumors and surrounding normal tissues, whereas θ values demonstrated excellent separation between normal and tumor regions; this is shown in the k versus θ plot for BxPC-3 tumor type in Fig. 3(c). Based on this analysis, only the θ fitting parameter was used for volumetric tumor delineation as shown in Fig. 3(d)-(f). θ threshold of 250 was selected to be high enough to ensure tumor boundary delineation (light blue color in Fig. 3(e)), yet low enough to segment out normal skin tissues (refer to θ intervals of normal and tumor tissue in Fig. 3(c)). Although the chosen VOI size resulted in reduced parametric image resolution, this potential limitation may be addressed in the future by dense OCT scanning to record more image voxels within the particular OCT system resolution volume. Histological evaluation (Fig. 3(l)) confirmed much better tumor-normal tissue separation with Gamma distribution textural analysis (Figs. 3(j),(k)), compared to conventional structural imaging (Fig. 3(g)) its direct intensity thresholding (Fig. 3(h)) or average-filtering (Fig. 3(i)) with a 24 × 3 sampling window.

Parametric imaging
Parametric images can substantially improve OCT contrast relative to regular structural images, enabling accurate 3D tumor delineation as demonstrated in Fig. 3; this is especially true in settings involving tumors that are highly heterogeneous, are not well circumscribed (presenting a diffuse rather than a sharp border), and are aggressively invading surrounding tissues. Figure 4 shows representative images on the two examined tumor models, culminating in 3D parametric visualization of melanoma (Fig. 4(c)) and pancreatic adenocarcinoma (Fig. 4(h)). The shape and location of the pigmented melanoma grown in the nude mouse thigh (Fig. 4(a)) is confirmed by the fused photoacoustic and ultrasound image that shows excellent spatial concordance with the OCT θ map that follows, (red color in Fig. 4(b), taken from the location labeled with dotted black line in (a)). In the corresponding OCT θ-parametric image (Fig. 4(c)), the skin appears in red, melanoma in blue and subcutaneous fat in yellow-green. Microphotograph of a pancreatic adenocarcinoma in a window chamber (Fig. 4(d)) presents a rapidly growing tumor several weeks after window chamber installation.
Location of the 3D OCT structural image acquisition is labeled with a white dotted rectangle. Fluorescently labeled cancer cells provide excellent tumor contrast in the epi-fluorescent microscope image in Fig. 4(e), exhibiting the visualization of 2D spatial extent of the live tumor cell mass. Major blood vessels seen in the image allow for its accurate (manual here for validation purposes) co-registration with the acquired OCT structural image through the simultaneouslyderived OCT depth-encoded microvasculature map (Fig. 4(f)). A two-dimensional average intensity projection (Fig. 4(g)) of a 3D Gamma fit parameter θ image (Fig. 4(h)) and its thresholded tumor volume (Fig. 4(i)), obtained from the same OCT dataset as Fig. 4(f), provides detailed information about tumor tissues growing within the skin. The sharp tumor-skin "interface" at θ = 250 allows for accurate tumor delineation through θ thresholding for further analysis.
Each tumor type possesses unique optical scattering characteristics due to structural and compositional differences of its cellular (nuclear and cellular size, cellularity) and surrounding tumor microenvironmental (stromal organization, interstitial fluids content, microvascular density) compartments, thus changing the statistical distribution of backscattered light in a unique and tumor-specific way. It is therefore unlikely that the product kθ of the Gamma fit parameters that we use in our approach and which represents the mean of the statistical distribution will be identical for various tumor types. Indeed, this is evident from Fig. 4(c) (melanoma) and 4(h) (pancreatic adenocarcinoma), where tumors in θ parametric space exhibit slightly different numerical parameter values.

Delineation method application to a pilot study of early tumor response to radiotherapy
Armed with a validated and objective method to delineate tumor extent in 3D, we applied this methodology to the challenging problem of radiotherapy response monitoring. We chose BxPC-3 tumor model because of the current emerging clinical interest in stereotactic body radiotherapy for treating pancreatic cancer [37,38]. Our previous work showed that detecting the subtle radiobiological changes in OCT images is prone to artefacts stemming from subjective tumor volume selection, and further exacerbated by the need for repeated imaging and corresponding tumor delineation throughout the post-treatment monitoring process. In addition to the challenges associated with substantial tissue heterogeneity that affects the intra-tumor voxel intensity distribution and thus resultant fitting parameters, the tumor may change its shape and size following irradiation. It may move out of the initial field of view, thus making longitudinal tissue volume selection / tracking very complicated. Its microenvironmental components (e.g., blood and lymph vasculature, collagenous stroma, interstitial fluids) may also substantially change depending on the treatment dose, time following treatment, tumor size, location and other parameters. In light of these various challenges, assessment of the whole tumor response with the developed delineation methodology aimed to reduce the human error and fully automate the objective evaluation process.
In this pilot study, we set up a novel scheme to irradiate only half of each tumor. The reasons for this choice were (1) the complexity of the longitudinal monitoring problem (e.g., are the observed changes truly due to radiation, or instead simply due to the dynamics of tumor evolution / development unrelated to RT?) and (2) the fact that the number of experimental animals may be reduced if RT-induced changes are assessed by having a non-treated "reference" from the same tumor. The alternative of having non-treated tumors in different mice act as controls is not very useful owing to large inter-animal variations in tumors [35]. However, our choice of intra-animal "reference" is not ideal either, as it may be argued that the by-stander effects manifest in the un-irradiated half may be present and thus introduce additional complexity (biological effects expressed after irradiation by cells that have not been received the primary given dose; these can include DNA damage, chromosomal instability, mutation, and apoptosis [39]). A related complication, the abscopal effect -a fascinating phenomenon of shrinkage of untreated tumors after irradiation of the primary tumor due to the activation of the immune system response [40] is unlikely to be a major artefact here due to the immunocompromised nature of experimental animals. With these caveats in mind, we performed the half-tumor irradiations in pancreatic adenocarcinoma (BxPC-3) model and monitored the two sides (irradiated and unirradiated) of each tumor for 5 weeks following dose deposition.
A total of eleven longitudinal observations of each BxPC-3 tumor for 5 weeks after single-dose 15Gy irradiation were recorded. Figure 5 shows six characteristic "snapshots" (one per week) for one of the animals. Photomicrographs (column (a) in Fig. 5) help visually observe tumor changes in response to the treatment. The treated part of the tumor was initially slightly larger than the untreated part in this particular mouse (seen clearly in the top-row white light image of Fig. 5(a) and fluorescence viability image of Fig. 5(b)). One can then see that throughout the monitoring period, the treated part gradually shrinks and the untreated part gradually increases in size over 5 weeks. This observation is confirmed in Fig. 5 column (b), where tumor shape changes similarly in fluorescence images. This irradiated tumor shrinkage was expected, in light of the ∼0.6% tumor cell survival anticipated in the first days after a 15Gy dose [41], and the avalanche of tumor cell deaths in the weeks following vascular network destruction [35,42]. Indeed, depth-encoded vascular maps in Fig. 5 column (c) demonstrate gradual vascular destruction of the treated part of the tumor starting from week 1 followed by revascularization closer to week 5; we have previously noted similar temporal damage trajectory [35]. In contrast, the untreated part of the tumor shows the opposite behavior: it continuously enlarges and its vascular network becomes more complex within this non-irradiated tumor region. In our previous BxPC-3 tumor monitoring study in the same window chamber / NRG mouse model, untreated tumors doubled their size in 5-6-weeks [35]. Here, similar trends are seen in the unirradiated part of the tumor as well. Obviously, VOI size and location selection to accurately segment tumor volumes in these 'changing target' conditions in tens of different image sets over time would be highly subjective and greatly prone to error, effectively masking the detection of any underlying radiobiological responses. Yet the 3D θ parametric images in Fig. 5 column (d) successfully overcome these complexities and potential artefacts, allowing for accurate tumor delineation and further quantitative analysis. Indeed, in this half-irradiation case, θ parameter thresholding and division of the whole tumor into two parts (treated vs untreated) was automated, objective and straightforward.
Moving forward, the same Gamma distribution fit parameters used for tumor delineation may also serve as early biomarkers of tumor response to the treatment, as radiation likely causes characteristic microstructural transformations at the cellular level that are known to change the optical properties of the cell [21]. Previous studies reported that OCT backscattered signal spectral [43] and speckle texture [26] analyses can be related to in vitro and in vivo apoptosis, necrosis and mitotic catastrophe of tumor cells. With these three types being major pathways of radiation-induced cell death [44], we hypothesize the changes in OCT image speckle texture, and hence the Gamma distribution fit parameters of speckle patterns may reflect tumor response.
Our pilot study supports this hypothesis. Figure 6(a),(b) shows average signal intensity distributions of non-irradiated and irradiated parts of the tumor shown in Fig. 5 for the whole 5-week period of observation. The non-irradiated tissue's slight deviations at different times after treatment around the central pre-treatment week 0 distribution (black dashed curve in Fig. 6(a)) are small, suggesting modest changes in the tumor microstructure, perhaps owing to natural disease progression. This is in stark contrast to the large changes seen over time in the irradiated tumor side (Fig. 6(b)). Plotted for simplicity in red and blue colors over time for the untreated and treated regions, respectively, two interesting features characterize the differences of these distributions. First, for approximately two weeks after RT, the non-irradiated region signals roughly follow the irradiated region trends (green upward arrows in (a) and (b)), although on a smaller scale. This may be an indicator of a by-stander effect in light of significant cell death expected on the treated side, radiation-induced short-term severe blood vessel occlusion [38] and subsequent feeding vasculature destruction, as seen in the Fig. 5(c) column. Second, after the week-2 time point, the similarity of the OCT signal distribution changes for the two regions essentially stops, and their dynamics begin to differ significantly (grey downward arrows in (a) and (b)). Fig. 6. Average signal intensity distributions of (a) non-irradiated and (b) irradiated parts of the BxPC-3 tumor plotted for the whole 5-week period of observations, using the mouse data whose images are presented in Fig. 5. The initial (Week 0) distribution is plotted in dashed black in both (a) and (b). Green arrows show the direction of the signal distribution changes during first two weeks after radiotherapy, grey arrows -after the two-week period; (c) and (d) Longitudinal changes in k and θ parameters for two experimental tumors, as represented by four symbols at each time point (mouse 1 irradiated and control halves, ditto for mouse 2). Non-irradiated parts of each tumor are represented with green circles, irradiated parts -with yellow triangles. In accord with the color scheme of panels (a) and (b), corresponding red and blue dashed lines connect average values for two tumors. Note the relative constancy of both k and θ in the non-irradiated regions of the tumor, and their respective temporal changes with irradiation (∼30% reduction in k, ∼500% increase in θ).
Distribution differences may be more effectively characterized by the parameters of the Gamma distribution fits. Figures 6(c) and 6(d) show longitudinal changes in k and θ parameters respectively, for two mice (44 volumetric data sets in total). Non-irradiated parts of each tumor are represented with green circles while radiation treated parts are shown with yellow triangles, with corresponding red and blue dashed lines (same color scheme as in (a) and (b)) displaying the average values for the two animals. Both parameters are normalized to the day of the treatment (week 0) individually for each animal. Together with the two abovementioned differencessimilar dynamics of changes on the two sides up to ∼ 2 weeks, and significant divergence in response afterwards -other important features may be singled out. What immediately attracts attention is the magnitude of significant changes in the Gamma fit parameter values of the treated regions, compared to their relatively steady behavior of the non-irradiated regions of each tumor. Further, parameter k is ∼20 times less sensitive than θ to the 15Gy treatment: its maximum change is at max ∼25% decrease from initial values, compared to a greater than 500% rise in θ after week 3. This may be an indicator of different sources of tissue changes influencing the k and θ parameters as discussed in the Introduction. Previous ultrasound studies suggest that k is related to the effective tissue scatterer number density [21]. In view of OCT as an optical analog to ultrasound [45] and in accord with the recent OCT study in a controlled phantom system of different-sized microspheres [25], our observed small k change after radiotherapy may also relate to the small changes in tissue components like interstitial fluids and porous structures (e.g., blood and lymphatic vessels). This may be an important conjecture to investigate further, given the delayed and dose-threshold-dependent response of, for example, blood vessels compared to tumor cells [35,46]. This also may partially explain the k parameter distribution overlap observed in Fig. 3(c): blood, lymph and other fluids are present in both tumors and non-malignant tissues, albeit likely in different amounts.
A significant θ increase for the treated part of the tumor three weeks after radiation correlates well with the tumor volume shrinkage (related to the ∼0.6% tumor cell survival rate mentioned above) seen in Fig. 5 (columns (a) and (b)). In the 1300 nm near-infrared spectral region of the OCT system used for imaging, elastic scattering is the dominant interaction of light with biological tissue [47]. Assuming that light absorption remains constant and low, and the θ parameter is proportional to attenuation of the backscattered light, these novel pilot results suggest that θ may indeed be related to cell death as this process does change the scattering characteristic of cells. This association has been previously suggested in the context of chemotherapy for AML [21,26] and for human fibroblast cells [48] in-vitro, and may also hold true for radiotherapeutic treatments in-vivo as in the current study.
Previous promising results obtained by several research groups, limited to well-controlled in-vitro conditions, demonstrated that backscattered signal characteristics and parametric fits nevertheless varied from one batch of cells to another [21]. Further, the follow-up in-vivo study was limited to only a few animals without longitudinal assessment [42], owing to even greater challenges of experimentation and analysis, and highlighting the fact that live animal experimental data, where the microenvironment is less well controlled, should be interpreted with caution. This is further compounded by the challenges of repeated imaging over weeks, and the need for accurate, reproducible and consistent tumor VOI delineation at each time point to enable elucidation of longitudinal trends. Thus, this pilot radiation response study is limited to a detailed temporal examination of only two animals and one dose level (yet spanning 6-7 weeks and containing >10 longitudinal high-quality 3D image sets for each). However, it is clear even at this early stage that Gamma distribution fit parameters can (1) robustly delineate tumor from normal tissue, enabling a variety of biomedical applications (e.g., early disease detection / diagnosis / assessment and hitherto impossible longitudinal response monitoring), and (2) reflect some of the tumor changes following irradiation.
The observed changes in the Gamma parameter fits post-RT may be indicative of tumor cell death, which must be confirmed by a thorough histological assessment at various post-irradiation stages. If confirmed, this OCT textural analysis methodology appears promising for a variety of interesting studies. For example, the proposed dorsal skin window chamber "half-irradiated" model may be used for direct investigation of abscopal effect by comparing the responses of immune-compromised and immune-competent animals developing tumors of the same type. An extensive OCT study involving a larger cohort of animals, various radiation schedules, and other independent readout metrics (tumor volume, vascular density, histological assessment) is currently ongoing in our laboratory, in an effort to further "shed light" on radiotherapy and possibly optimize / personalize this important cancer treatment modality.

Summary and conclusion
This study has demonstrated that it is possible to delineate the 3D tumor in vivo using the Gamma distribution fits of OCT image voxels: in the fitting parameter space, sharp tumor boundaries allowed for accurate tumor volume delineation through θ parameter thresholding. The developed methodology was then applied to the challenging problem of radiotherapeutic monitoring over a multi-week observation period, showing that the same fitting parameters are sensitive to early tumor response to radiation. Five-week tumor monitoring after irradiation demonstrated significant changes in Gamma fit parameter values compared to their relatively steady behavior on the non-irradiated regions of the tumor. Some putative biophysical interpretations of these trends have been offered, but need further thorough investigation and independent validation. Automatic and objective delineation of irradiated tumors and simultaneous detection of early signs of their radiobiological response in a larger cohort of animals at clinically relevant doses and schedules is currently ongoing and will be reported in a separate publication.
Funding. Canadian Institutes of Health Research (173971).