Drug homogeneity index in mass spectrometry imaging (MSI

6 days ago - ... homogeneity is usually manually assessed by an expert and this approach is biased by inter-observer variability and lacks reproducibi...
0 downloads 0 Views 398KB Size
Subscriber access provided by University of Winnipeg Library

Article

Drug homogeneity index in mass spectrometry imaging (MSI) Mridula Prasad, Geert Postma, Lavinia Morosi, Silvia Giordano, Raffaella Giavazzi, Maurizio D'Incalci, Francesca Falcetta, Enrico Davoli, Jeroen J. Jansen, and Pietro Franceschi Anal. Chem., Just Accepted Manuscript • DOI: 10.1021/acs.analchem.8b01870 • Publication Date (Web): 25 Oct 2018 Downloaded from http://pubs.acs.org on October 25, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 23 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Analytical Chemistry

Drug homogeneity index in mass spectrometry imaging (MSI)

Mridula Prasad,†,‡,¶ Geert Postma,¶ Lavinia Morosi,§ Silvia Giordano,k Raffaella Giavazzi,§ Maurizio D’Incalci,§ Francesca Falcetta,§ Enrico Davoli,k,# Jeroen Jansen,¶,# and Pietro Franceschi †,#,*

†Computational

Biology Unit, Research and Innovation Centre, Fondazione Edmund Mach, 38010 San

Michele all’Adige, Italy. ‡ DNanotechnology

in Medicinal Chemistry, Department of Molecular Biotechnology and Health Sciences,

Univeristà di Torino, 10124 Torino, Italy. ¶ IMM

§

/ Analytical Chemistry, Radboud University, Heyendaalseweg, 6525 AJ Nijmegen, The Netherlands.

Department of Oncology, IRCCS Istituto di Ricerche Farmacologiche Mario Negri, Via La Masa 19,

20156 Milan, Italy. k Department

of Environmental Health Science, Mass Spectrometry Laboratory, IRCCS Istituto di Ricerche

Farmacologiche Mario Negri, Via La Masa 19, 20156 Milan, Italy.

#Contributed

equally to this study

* Corresponding Author: Pietro Franceschi Computational Biology Unit, Research and Innovation Centre, Fondazione Edmund Mach, 38010 San Michele all’Adige, Italy. Tel: +39-0461-615556 E-mail: [email protected]

1 ACS Paragon Plus Environment

Analytical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Abstract Enhancing drug penetration in solid tumors is an interesting clinical issue, of considerable importance. In pre-clinical research mass spectrometry imaging is a promising technique to visualize drug distribution in tumors in different treatment conditions and its application in this field is rapidly increasing. However, in view of the huge variability among MSI datasets, drug homogeneity is usually manually assessed by an expert and this approach is biased by inter-observer variability and lacks reproducibility. We propose a new texture-based feature, the ‘Drug Homogeneity Index (DHI)’ which provides an objective, automated measure of drug homogeneity in MSI data. A simulation study on synthetic datasets showed that previously known texture features do not give an accurate picture of intra-tumor drug distribution patterns and are easily influenced by the tumor tissue morphology. The DHI has been used to study the distribution profile of the anti-cancer drug paclitaxel in various xenograft models, either not pretreated or pretreated with antiangiogenesis compounds. The conclusion is that drug homogeneity is better in the pretreated condition, which is in agreement with previous experimental findings published by our group. This study shows that DHI could be useful in preclinical studies as a new parameter for the evaluation of protocols for better drug penetration.

Introduction Mass spectrometry imaging (MSI) is a molecular imaging technique that can measure the distribution of metabolites and small molecules in tissue slices1–5 with high specificity and sensitivity, providing quantitative information on the diffusion of drugs.3 Since complete and efficient penetration of anticancer drugs in solid tumors is essential for optimal therapeutic effect, this technique would certainly help in understanding one of the causes of inefficient responses to chemotherapy or drug resistance, i.e. the heterogeneous and irregular tumor drug distribution because of abnormal tumor microenvironment. 6,7 Our group research is focused on improving the penetration of drugs by altering the tumor microenvironment.3 Tumor xenograft models are extremely valuable to monitor with imaging techniques the distribution of the drug in the tumor tissue, correlating it with the pharmacokinetic profiles and the antitumor activity8-10. In particular, we have recently shown the effect of pretreatment with the anti-angiogenic drug bevacizumab in the normalization of tumor vasculature with the consequent improvement of penetration in solid tumor and the antitumor activity of the chemotherapy agent paclitaxel8. Unfortunately, however, the assessment of the level of homogeneity of drug distribution is rather complex due to the lack of a suitable parameter to 2 ACS Paragon Plus Environment

Page 2 of 23

Page 3 of 23 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Analytical Chemistry

describe it quantitatively. We assessed it manually in our MSI dataset with a procedure that is labor-intensive but also poorly reproducible and is easily biased by human variability. In the manual approach applied in our previous publications8-9 drug and internal standard peaks are manually identified in the mass spectrum and extracted using a fixed mass range. The drug image was normalized dividing by the internal standard image obtaining the PTX/d5-PTX ratio. The distribution of the drug in these normalized images was analyzed by the software ImageJ (imagej.nih.gov/ij) 11. Regions of interest (ROIs) enclosing the tumor section were drawn. We have defined a threshold for enhancing pixels as a signal change greater than background noise in untreated samples. The number of pixels that exceeded this threshold was calculated and divided by the total number of pixels in the ROIs to obtain the percentage of pixels where the drug is present. In addition, we did a “particles analysis” on images after applying the threshold. This measures the characteristics (number, mean size, mean perimeter or circularity) of clusters formed by adjacent positive pixels (above the threshold) in the image in order to describe the drug distribution in detail. This image analysis approach, however, was really time consuming and error prone, relying entirely on manual procedures. The development of a Drug Homogeneity Index (DHI) that can objectively quantify drug homogeneity in MSI datasets is therefore highly desirable. This index, combined with a highly reproducible bioinformatics pipeline could pursuit automatic processing of series of MSI datasets. The medical literature frequently describes texture methods for the quantification of spatial intensity differences in clinical images, particularly for radiographic data.12–16 In this type of analysis, the characteristics of each image are translated into a set of texture features that contain information about the different spatial intensities,12 mostly using texture methods. The texture features are extracted from a series of matrices containing information on how the intensities are related across different pixels. The three commonly used matrices are: a) gray level co-occurrence matrix (GLCM),17 b) gray level run-length matrix (GLRLM),18 and c) gray level size-zone matrix (GLSZM).19 Details of the three matrices are given in the Supporting Information. These texture analysis methods represent the natural framework for the development of a DHI and have been also used to optimize the processing of untargeted MSI datasets inside the EXIMS pipeline20 and as a part of the pySM metabolite annotation framework21. In order to apply a homogeneity measure 3 ACS Paragon Plus Environment

Analytical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

in an automated fashion across multiple MS imaging datasets collected on tissue sections of different sizes and orientation, the index must show have some specific characteristics. It should be invariant to the size and shape of a tumor section, and it should also be “gray level” invariant because equivalent sections may present different absolute intensities due to the intrinsic analytical variability of MSI. Here we propose a new texture-based DHI feature to provide a robust estimate of drug homogeneity in MSI data. We introduce the index and compare its efficacy and invariance properties with other known texture features on a set of synthetic images. Then we validated our approach, verifying and reproducing the manually obtained results on MSI datasets of two ovarian tumor models3 and an additional new MSI datasets of an ovarian and a colon cancer tumor model.

Experimental Section To describe our approach, we ran a simulation study with two synthetic imaging data sets and a series of MSI datasets from different tumor cell lines. All the synthetic images are created using OpenCV (https://github.com/opencv/opencv) library in Python free software. Synthetic imaging data set I: Tumor tissue morphology invariance There three main requirements for robust texture features, are: a) gray level, b) shape and c) size invariance. To investigate this, we designed a series of images constructing three images for each property. To investigate gray-level invariance we created similar images with different gray levels. To understand how shape affected texture features we constructed images with approximately the same frequency of each gray levels but slightly different spatial arrangements. To check size invariance, we designed the three images with similar texture patterns, increasing in size. The synthetic imaging data used for the analysis are shown in Figure 1. Synthetic imaging data set II: Drug homogeneity To check whether a specific texture feature is sensitive to drug homogeneity, we created five basic images under four different scenarios (Figure 2a). Each image has a fixed object size with different intra-tumoral drug distribution patterns: images B to E are derived from image A by gradually adding dark zones, mimicking an increasing amount of drug-inaccessible regions, hence reducing the drug-accessible regions. Thus in each category image A has the most homogeneous drug distribution, because of the absence of dark-zones and image E is the least homogeneous. In 4 ACS Paragon Plus Environment

Page 4 of 23

Page 5 of 23 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Analytical Chemistry

scenario I, all images have a single gray level, to give an idea of the behavior of a texture feature with the increase in the number of dark zones. In scenarios II and III, images are produced as in scenario I, after adding respectively 2 and 4 gray levels. This helps verify the change in feature values due to the simultaneous increase of drug inhomogeneity (dark zones) and intensity transitions. One more factor that could influence the feature performance is technical noise in some images so in scenario IV we added a fixed amount of salt-pepper noise in each image. Texture analysis Texture analysis refers to the various mathematical methods used to quantify structural complexity in an image, using texture features. Several texture methods are available in which second- and higher-order statistics based methods are the most commonly used, as they describe the distribution and inter-relationships of gray-level values in the image. These methods are sensitive not only to the intensity measured at each pixel, but also to the spatial arrangement of the intensities. For this reason, they are assumed to closely depict the human perception of image texture. All the texture features are constructed from a set of matrices constructed from the initial image. Among these matrices we will focus on: a) the gray-level co-occurrence matrix (GLCM) b) the gray-level run-length matrix (GLRLM), and c) the gray-level size zone matrix (GLSZM). An illustration of the characteristics of these matrices is included in the Supplementary Information.

All

the

texture

features22

were

calculated

using

the

radiomics

(https://github.com/cran/radiomics) R package. The co-occurrence matrix was only calculated for adjacent pixels. To check direction independence, the GLCM and GLRLM matrices were calculated in all possible directions, and the feature values derived from all these matrices were averaged. Tumor models and mass spectrometry imaging The distribution of paclitaxel was measured in several xenograft models of human ovarian and colon carcinoma cells. Two ovarian cancer lines (A2780, HOC84) and one colon cancer line (HCT116) were ectopically transplanted and one ovarian cell line (IGROV1) was orthotopically transplanted. Tumor-bearing mice (4-5 per time point) were given 60 mg/kg of paclitaxel intravenously (day 0) alone or 5 and 1 days after bevacizumab (two intraperitoneal injections 150 µg/mouse). After 6 hours, the tumor was collected, divided into two, and frozen in liquid nitrogen. 5 ACS Paragon Plus Environment

Analytical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The IRFMN adheres to the principles set out in the following laws, regulations, and policies governing the care and use of laboratory animals: Italian Governing Law (D.lgs 26/2014; Authorization no. 19/2008-A issued March 6, 2008 by Ministry of Health); Mario Negri Institutional Regulations and Policies providing internal authorization for persons conducting animal experiments (Quality Management System Certificate - UNI EN ISO 9001:2008 - Reg. no. 6121); the NIH Guide for the Care and Use of Laboratory Animals (2011 edition) and EU directives and guidelines (EEC Council Directive 2010/63/UE). The Statement of Compliance (Assurance) with the Public Health Service (PHS) Policy on Human Care and Use of Laboratory Animals was recently reviewed (9 September 2014) and expires on 30 September, 2019 (Animal Welfare Assurance A5023-01). The animals were regularly checked by a certified veterinarian responsible for health monitoring, animal welfare supervision, experimental protocols and procedure revision. All surgery was done under general anesthesia, and all efforts were made to minimize suffering. For imaging analysis as described in Cesca et al.8, the frozen tissues were cut into 10-μm-thick sections using a cryo-microtome (Leica Microsystems, Wetzler, Germany) at −20°C and mounted on a pre-cooled MALDI plate (Opti-TOF 384 Well insert) by standard thawmounting techniques. One section every 100 μm was cut starting from the central part of the tissue. Four tumors per group and three sections per each tumor were analyzed. The plate, dried in a vacuum drier at room temperature overnight, was sprayed with TiO2 matrix suspension containing d5-PTX 3 μg mL-1, as internal standard, using a BD 180 precision double-action trigger airbrush (FENGDA, Zhejiang, China) with 0.20 mm nozzle diameter, with nitrogen at 0.2 atm. A MALDI 4800 TOF-TOF (AB SCIEX, Old Connecticut Path, Framingham, MA 01701, USA) was used. It was equipped with a 355-nm Nd:YAG laser with a 200 Hz repetition rate, controlled by 4000 Series ExplorerTM software (AB SCIEX). MS spectra were acquired in scan mode with 20 laser shots at an intensity of 6000 arbitrary units, with a bin size of 1.0 or 0.5 ns, in reflectron negative-ion mode. Mass spectra were recorded in full scan profile mode over a limited mass range (m/z 250–300). Images of tissue sections were acquired using the 4800 Imaging Tool software [ www.maldi-msi.org, M. Stoeckli, Novartis Pharma, Basel, Switzerland], with an imaging raster of 100 × 100 μm (pixel dimension ca. 0.01 mm2). MSI data pre-processing

6 ACS Paragon Plus Environment

Page 6 of 23

Page 7 of 23 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Analytical Chemistry

The initial data files were stored in Analyze 7.5 format which was processed in R free software23 version 3.1 using the MALDIquant24 package. All subsequent analyses were done in multiple steps using

R

script

available

as

a

Github

repository

ProcessMSI

(https://github.com/

pietrofranceschi/ProcessMSI). The proposed pre-processing pipeline was optimized for the dataset under consideration and have to be adapted to the specific characteristics of each dataset and instrument set-up (in particular in the case of high resolution spectrometers). MSI data were unpacked in a two-dimensional matrix where each row contained a complete spectrum collected in sequential order. To correct for the expected misalignment between the different spectra, binning over the m/z scale was done. The bin size was set at 0.5 daltons, considering the characteristics of the instrument used for data acquisition and the low chemical noise associated to the use of TiO2 matrix. As this is a targeted study the mass-to-charge ions of interest are known. Thus in the desired bin, peaks were selected using a local maxima search method above a certain signal-to-noise threshold using msProcess25 R package. The optimal threshold was determined spectrum-wise, as complex phenomena occurring during ionization mean that a fixed threshold value results in either false peak selection or rejection of several true ion peaks. The intensity values in each spectrum were therefore sorted in ascending order, and the median value of the initial 200 was used as an adaptive measure of the noise. The rationale behind this choice is that the majority of the low-intensity values are in fact due to the noise; this was verified by visual inspection. In the MSI dataset under analysis, we focused on three specific m/z bins: the deuterated internal standard (m/z 289.1-289.3) dissolved in a TiO2 matrix suspension, a side chain with an amide-acyl group of paclitaxel (m/z 284.01-284.23), and a tissue ion, most likely corresponding to oleic acid (m/z 281.16-281.44). Bins for these three main ions were selected to create an ms ion image. Drug and tissue ion images were normalized with the internal standard ion image. A density-based thresholding approach was used to ensure an accurate delineation of the tumor tissue morphology on the MALDI plate. A density plot was created for a tissue landmark peak inside the m/z 281.16-281.44 bin. Usually this plot showed a single local minimum which was used to separate the background (MALDI plate) from the foreground (tissue). However, for some of the sections the profile of the density plot was more complex because of the presence of background noise from the plate. To account for that and retrieve a reliable mask of the tissue, when multiple minima identified more than two clusters in the image, we discarded low-intensity clusters accounting for less than 5% of the total pixels (again, this was decided on the basis of our data 7 ACS Paragon Plus Environment

Analytical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 23

characteristics). This thresholding step must be done properly as we use this threshold image to estimate the tumor area in DHI calculation. The inherent analytical variability of MSI means that the absolute intensity range for each individual drug image collected from different samples can be different. To make the texture matrices constructed on these images comparable, the intensity values are “quantized” into a limited set of gray levels (8, 16, or 32) within the intensity range of each image. The formula is given below: 𝐼𝑞𝑢𝑎𝑛𝑡𝑖𝑧𝑒𝑑 =

𝑄𝑢𝑎𝑛𝑡𝑖𝑧𝑎𝑡𝑖𝑜𝑛 𝐿𝑒𝑣𝑒𝑙 max (𝐼𝑑𝑟𝑢𝑔)

𝐼𝑑𝑟𝑢𝑔

Here the quantization levels are the different maximum allowed intensity bins (8, 16, 32) used to convert an original image (Idrug) into the quantized image (Iquantized). Finally, after the quantization step, to remove the spatial background noise, the individual drug ion image was denoised with a two-dimensional three-pixel median filter to account for strong pixel-to-pixel variability. The DHI was calculated on the drug ion image of the tissue. Statistical analyses To test the difference in DHI resulting from of the drug delivery protocol, taking into account the correlation between the tumor slices from a single tumor organ, we used a linear mixed effect (LME) model approach.26 P-values under 0.05 were considered significant after the application of Bonferroni correction for multiple comparisons. The LME model was fitted and corrected using the R packages nlme27 and multicomp28

Results and discussion DHI definition The proposed DHI is based on the gray level size-zone matrix (GLSZM),19 which is already rotation-invariant by construction. The DHI is defined as: 𝑁

𝐷𝐻𝐼 =

𝑁

∑𝑖 =𝑔 1∑𝑗 =𝑧 𝑁 𝑗𝑃(𝑖,𝑗) 𝑢

1

𝑁 𝑁 ∑𝑖 =𝑔 1∑𝑗 =𝑧 𝑁 𝑃(𝑖,𝑗) 𝑇𝑢𝑚𝑜𝑟 𝑢

𝐴𝑟𝑒𝑎

where, P(i,j) is the gray level size-zone matrix containing the gray levels (Ng) as rows and the size zones (Nz) as columns (see Methods). Each element of P(i,j) is the number of times a contiguous area of gray level i and size j is found in a given image. To favor big homogeneous areas, each 8 ACS Paragon Plus Environment

Page 9 of 23 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Analytical Chemistry

column of P(i,j) is multiplied by the index of column j which is also equal to the absolute area of each gray zone. To make different images comparison feasible normalization step is required. The 𝑁𝑧

∑𝑗 = 𝑁𝑢𝑃(𝑖,𝑗)term accounts for the overall “weight” of P(i,j), while 1/Tumor Area compensates for the size of the tumor. The final DHI value obtained from the overall formula lies between 0 and 1. DHI is expected to grow monotonically with the visual homogeneity of the image. However, the homogeneity measure may be biased by noise. The presence of many small zones with uniform intensity will produce a low value for the index because the first normalization factor becomes large on account of the high frequency of small regions. A noisy image can obviously be considered non-homogeneous, so to be useful in real-life applications the index has to be less sensitive to small variations which can be associated with pixel-to-pixel technical variability. In order to do that, small structures can be excluded from the DHI by summing only the values of the columns above a specific index Nu (Nu ∈ [1,Nz]). This is equivalent to setting the size of the smallest zone in the DHI calculation. An implementation of the proposed algorithm is available as R package on GitHub (https://github.com/pietrofranceschi/HomogenMSI). In order to allow the testing of the significance of the measured homogeneity the package also implements a permutation based approach to calculate the empirical null distribution of DHI. Calculation of DHI on synthetic and real MSI datasets We examined the performances of known texture features (see the Methods section) and our newly developed DHI in two steps using various synthetic images. Then we applied the DHI to MSI datasets from ovarian and colon cancer xenograft models under two treatment conditions. Testing invariance properties: synthetic data set I The prototype images used to investigate the robustness of the different texture features towards changes in tumor tissue morphology are set out in Figure 1, grouped in three main categories corresponding to the three required invariance properties. In the first row, the images are the same shape and size but have different absolute values for three gray levels in each image. In the second row, the number of pixels in different gray level areas is approximately the same, but the spatial arrangement is different. The last three images have the same shape and the same gray level intensities, but the absolute sizes are different. Even in this oversimplified scenario, a robust homogeneity index should return a constant value for each row. 9 ACS Paragon Plus Environment

Analytical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Including our DHI, a total of 44 features can be obtained from the different texture matrices. All texture feature values calculated for the nine images are listed in Supplementary Table 1, and Table 1 summarizes the ones showing good invariance. Out of 21 GLCM-based features, 11 are gray level invariant, two are shape invariant and none are size invariant; of the GLRLM matrix-based features five are gray level invariant, three are shape invariant, and none are size invariant. In the GLSZM-based features, intensity variability (IV), size zone variability (SZV), low-intensity emphasis (LIE), and high-intensity emphasis (LIE) are robust in relation to image object shape and size. As their name indicates, LIE and HIE features highlight high and low-intensity regions, so their dependency on image gray levels is to be expected. In summary, these results show that the values of GLCM and GLRLM based features very much depend on image object morphology so they cannot be used to compare the homogeneity of images of different tissue sections. Out of all the features, only GLSZM-based IV and SZV features show consistent invariant properties. The proposed DHI provides a constant value for each property Testing the sensitivity to drug homogeneity: synthetic dataset II Feature invariance to tumor tissue morphology is the main desirable property for a reliable homogeneity index, but a good index should also monotonically decrease with any decrease in drug homogeneity. We tested this aspect with another series of synthetic images. Two main factors indicate a decrease in drug homogeneity: 1) large areas with different intensity (i.e. a large number of gray levels, 2) discontinuities in drug distribution, i.e. the absence of drug in contiguous regions. In other words, homogeneous distributions have fewer intensity transitions, resulting in a smooth image with few gray levels. The set of synthetic images used to test this is described in the Methods section and shown in Figure 2a. The normalized feature values for these images are summarized in Figure 2b. An optimal homogeneity index should give a decreasing value from left to right and from top to bottom. In scenario I, where all the images have the same gray level, the decrease in drug homogeneity is due to an increasing number of dark spots, mimicking the inaccessibility of a drug in certain areas. Here both IV and SZV returned a constant value for all the images. This is because the GLSZM of all images contains a single size-zone of different dimensions; however, the IV and SZV features take into account only the count of size-zones, resulting in a constant value. The DHI worked better here because it is weighted by the absolute size-zone value, hence it gradually decreases from 10 ACS Paragon Plus Environment

Page 10 of 23

Page 11 of 23 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Analytical Chemistry

image A to E. Again, in scenarios II and III, the behavior of SZV and IV were not in accordance with the expected evolution of drug distribution from image A to E. In scenario II, for example, SZV was constant for all images, while IV showed incoherent results, giving similar values for image A and B, with a sudden increase in value from image D to E. Also here DHI showed a coherent behavior. Scenario IV, where the salt-pepper noise was added to all images, deserves special attention. The SZV value was indeed gradually decreasing from image A to E (that could be the sign of the gradual decrease in size-zone count values), but IV had the highest value for image B, which was again unexpected. With Nu set to one also the DHI value was almost constant, while it showed the expected decreasing trend only for larger values of. This behavior is not unexpected considering that salt and pepper noise impact the structure of the dataset at the single pixel scale (it actually adds random black or white pixels to the image) so it will dominate the measure of homogeneity when small size zones are included. This scenario clearly demonstrates the importance of a principled choice of in real life applications. Assessment of paclitaxel homogeneity in MSI data For a final validation of the proposed approach with a real-life dataset, we calculated the DHI was on a series of MSI images showing the drug (Paclitaxel) distribution in several tumor xenograft models. Manual findings on these datasets were partially already presented in our previous publications giving as a validation reference.3,4 Half the images came from animals that had received anti-angiogenesis treatment (with bevacizumab) before drug injection mimicking a clinical relevant setting; the other half were used as non-pretreated controls. An example of the paclitaxel drug ion images from MSI data of ovarian (A2780-1A9, HOC84, and IGROV1) and colon (HCT116) cancer models is shown in Figure 3 and the full set of images is included in the Supplementary Material. Pretreatment with bevacizumab in two of these models (A2780-1A9 and IGROV1) promoted more homogeneous drug distribution4, so here we assessed whether the DHIbased automatic pipeline could confirm this. The DHI for each ion image, with different numbers of gray levels (8, 16, and 32) at multiple Nu parameter settings (1, 5, and 10), are included in Supplementary Table 2. Figure 4 gives a subset of the DHI values for all the cell data (Iquantized = 8, Nu = 10). To check for the presence of significant differences in drug homogeneity between pretreated and not pretreated animals, the DHI from the various MSI datasets were analyzed using a linear mixed effects model, taking into account the study design. This analysis showed a 11 ACS Paragon Plus Environment

Analytical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

significant difference in homogeneity between pretreated and not pretreated animals in the A27801A9 and IGROV1 tumor models (means and corrected p-values are given in Supplementary Table 2), which is in accordance with previous results3. Bevacizumab pretreatments effect on the homogeneity of drug distribution also is evident in the additional colon cancer model HCT116, further reinforcing the conclusions of our previous paper and underlining the potentiality of DHI index. In the HOC84 model, instead, there was no real difference in homogeneity between pretreated and not pretreated animals, but this too agrees with our previous finding on HOC84’s histopathological characteristics, indicating a basal (even without pretreatment) homogeneous paclitaxel distribution.4 In conclusion, the proposed automatic approach confirmed the conclusions about the different distribution of paclitaxel without pretreatment in different tumor models. The DHI was then able to approximate the human visual perception of intra-tumoral drug distribution and can serve quantitatively to obtain statistically robust results.

Discussion MSI offers an unprecedented potential for studying drug access inside the different regions of tumors. To drew a quantitative, statistically robust conclusion from these experiments one has to deal with the huge morphological and experimental variability naturally occurring in MSI datasets from on animal models. For these reasons, visually comparing drug homogeneity is not only subjective but is also extremely challenging, even for a trained researcher. Here we propose a reproducible and standardized texture-based approach that can automatically provide a robust quantitative measure of drug homogeneity. Texture-based methods can highlight the local and global variation in spatial structures of an image and, for this reason, they have been also recently proposed for the discrimination of ions showing a “spatially relevant” distribution in untargeted MSI datasets either as part of the EXIMS pipeline20 or for the calculation of “Spatial Chaos”21 Even if these implementations, have been designed to process a single MSI dataset and not to analyze independent MSI datasets, an extensive comparison of these two alternative approaches with DHI on real and synthetic datasets is presented in the Supporting Information. The analysis of the curated imaging dataset described in21 (Figure S3 and figure S4) clearly shows that Spatial chaos and DHI are strongly associated and this is not unexpected since both methods are relying on size zone calculation and all images have exactly the same size. The correlation with EXIMS is instead low and negative, indicating that this method does not give coherent results when 12 ACS Paragon Plus Environment

Page 12 of 23

Page 13 of 23 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Analytical Chemistry

compared with the other two. Interestingly, the relation between Spatial Chaos and DHI is not linear, supporting the conclusion that DHI could provide a more contrasted quantification of the homogeneity of the ion images as shown in Figure S4. When the comparison of the three indices is performed on the data obtained in different imaging experiments (Figure S5b), the relation between DHI, EXIMS and Spatial Chaos change radically, but it is important to remember we are dealing with sections of different size and shape so the normalization by the total tumor area (included in DHI) is necessary make the different images comparable. The characteristics of the DHI approach, then, make it ideal to compare multiple MSI datasets, substantially enlarging the potential application of MSI in biomedical research. Our simulation study clearly shown that none of the previous texture features do this, as is also suggested in another study on radiomics data.29 The tumor morphology-based invariance properties are clearly achieved with our DHI formulation. Gray level invariance is obtained by not taking account of the absolute intensity of each gray level, but considering the size-zones occupied by the various gray levels. The shape invariance results from construction based on the characteristics of the GLSZM. Size invariance is achieved by normalization with the total tumor tissue area. Without this normalization, the DHI increases with the tumor tissue area (with constant tissue homogeneity). To obtain reliable results from the DHI, the tissue area must be determined exactly. In our case, on the basis of expert knowledge we identified a single ion peak at m/z = 281.2 which was present in all the tumors. Using this peak we separated the tissue from its background and calculated the tissue area exactly. In the absence of this information, one has first to identify ion peaks in the data and then create and compare the images of the various ions, to identify the optimal ion peak that outlines the tissue best (tissue morphology). The co-registered image made up of optical image of the tumor section and the MSI dataset can also be used to extract information for the tumor histological boundaries. Since some instrumental variability is inevitable in MSI datasets, pixel-to-pixel analytical variability can adversely influence the homogeneity calculation. As this type of variability is not related to drug distribution, we suggest neglecting the smallest homogeneous size zones by acting on the Nu parameter. This could be optimized taking into consideration the size of the laser spot, the rastering step and the histological characteristics of the tissue. Analyzing the real datasets, we noted a slight dependence of the homogeneity measure on the number of intensity quantization levels used. For example, the HCT116 cell line images created with 32 gray levels did not show any significant difference in therapy outcome, but with 8 and 16 gray levels, the difference was 13 ACS Paragon Plus Environment

Analytical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

significant. This is because if too many gray levels are included the homogeneity assessment becomes too dependent on very small gray-level differences. In any case, this was not seen with the other cell lines. The proposed DHI was calculated on real MSI datasets with the aim of automating manual and quantitative findings and test whether it can be applied reliably on reallife datasets, facilitating fast and efficient homogeneity assessment. A proper biological study should include biological replicates30, for this reason we studied an average of five mice per tumor model and treatment condition, and three tissue slices per mouse. In certain mouse models, the variability within the MSI data was vast (see Supplemental Figure 1) and was taken into consideration by applying the LME model for statistical analysis. However, DHI does seem to indicate an increase in drug homogeneity for the A2780-1A9, IGROV1, and HCT116 tumor models, confirming other reports and extending our knowledge on this topic. Thus, we encourage involving multiple tumor models and slices per tumor model regarding final conclusion over therapy performance based on our DHI.

Conclusions We propose a highly reproducible systematic approach for drug homogeneity assessment in imaging data, in this case hyperspectral MSI data. The DHI is a valuable tool to mimic human visual perception without the bias of inter-observer variability. More importantly, it allows efficient comparison of different imaging datasets and this could greatly boost the potential of MSI in tumor research. DHI can summarize intensity-area relationships in an optimal way, giving the oncologist an objective robust measure of drug homogeneity, facilitating the evaluation of drug distribution and the efficacy of new therapies or combined approaches. Apart from drug homogeneity assessment, our discussion of the requirements of texture features and our approach using synthetic images to validate possible indices could help researchers dealing with other texture-related problems.

Acknowledgement We thank J.D. Baggot for Editorial Assistance. This work was supported by grants from the Cariplo Foundation for the project “Nanostructuredinitiators for matrix-free, surface-based mass spectrometry imaging of antitumor drugs in tissues” (Project 2013-0692).

14 ACS Paragon Plus Environment

Page 14 of 23

Page 15 of 23 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Analytical Chemistry

Supporting Information Available The following files are available as Supporting Information: • SupplementaryInformation.pdf: 1) set of images showing PTX distribution in all tissue sections included in this study; 2) concise introduction to texture matrices and texture properties. 3) Comparison of DHI with Spatial Chaos and EXIMS. • SupplementaryTable1.xls: complete list of texture features for the images in Figure 1. • SupplementaryTable2.xls: DHI values for the full set of MSI tumor sections.

15 ACS Paragon Plus Environment

Analytical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

References (1) Cornett, D. S.; Frappier, S. L.; Caprioli, R. M. Analytical Chemistry 2008, 80, 5648– 5653. (2) Prideaux, B.; Stoeckli, M. Journal of Proteomics 2012, 75, 4999–5013. (3) Morosi, L.; Zucchetti, M.; D’Incalci, M.; Davoli, E. Current Opinion in Pharmacology 2013, 13, 807–812. (4) Karlsson, O.; Hanrieder, J. Archives of Toxicology 2017, 91, 2283–2294. (5) Nilsson, A.; Goodwin, R. J. A.; Shariatgorji, M.; Vallianatou, T.; Webborn, P. J. H.; Andrén, P. E. Analytical Chemistry 2015, 87, 1437–1455. (6) Minchinton, A. I.; Tannock, I. F. Nature Reviews. Cancer 2006, 6, 583–592. (7) Trédan, O.; Galmarini, C. M.; Patel, K.; Tannock, I. F. Journal of the National Cancer Institute 2007, 99, 1441–1454. (8) Cesca, M.; Morosi, L.; Berndt, A.; Nerini, I. F.; Frapolli, R.; Richter, P.; Decio, A.; Dirsch, O.; Micotti, E.; Giordano, S.; D’Incalci, M.; Davoli, E.; Zucchetti, M.; Giavazzi, R. Molecular Cancer Therapeutics 2016, 15, 125–135. (9) Giordano, S.; Zucchetti, M.; Decio, A.; Cesca, M.; Nerini, I. F.; Maiezza, M.; Ferrari, M.; Licandro, S. A.; Frapolli, R.; Giavazzi, R.; Maurizio, D.; Davoli, E.; Morosi, L. Scientific Reports 2016, 6, 39284. (10) Morosi, L.; Spinelli, P.; Zucchetti, M.; Pretto, F.; Carrà, A.; D’Incalci, M.; Giavazzi, R.; Davoli, E. PLOS ONE 2013, 8, e72532. (11) I Schneider, C. A.; Rasband, W. S.; Eliceiri, K. W. Nat. Methods 2012, 9 (7), 671–675. (12) Davnall, F.; Yip, C. S. P.; Ljungqvist, G.; Selmi, M.; Ng, F.; Sanghera, B.; Ganeshan, B.; Miles, K. A.; Cook, G. J.; Goh, V. Insights into Imaging 2012, 3, 573–589. (13) Chicklore, S.; Goh, V.; Siddique, M.; Roy, A.; Marsden, P. K.; Cook, G. J. R. European Journal of Nuclear Medicine and Molecular Imaging 2013, 40, 133–140. (14) Tixier, F.; Rest, C. C. L.; Hatt, M.; Albarghach, N.; Pradier, O.; Metges, J.-P.; Corcos, L.; Visvikis, D. Journal of Nuclear Medicine 2011, 52, 369–378. (15) El Naqa, I.; Grigsby, P. W.; Apte, A.; Kidd, E.; Donnelly, E.; Khullar, D.; Chaudhari, S.; Yang, D.; Schmitt, M.; Laforest, R.; Thorstad, W. L.; Deasy, J. O. Pattern Recognition 2009, 42, 1162–1171. (16) Michoux, N.; Van den Broeck, S.; Lacoste, L.; Fellah, L.; Galant, C.; Berlière, M.; Leconte, I. BMC Cancer 2015, 15, 574. 16 ACS Paragon Plus Environment

Page 16 of 23

Page 17 of 23 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Analytical Chemistry

(17) Haralick, R. M.; Shanmugam, K.; Dinstein, I. IEEE Transactions on Systems, Man, and Cybernetics 1973, SMC-3, 610–621. (18) Galloway, M. M. Computer Graphics and Image Processing 1975, 4, 172–179. (19) Thibault, G.; Fertil, B.; Navarro, C. L.; Pereira, S.; Cau, P.; Lévy, N.; Sequiera, J.; Mari, J.L. Texture indexes and gray level size zone matrix. Application to cell nuclei classification. 10th International Conference on Pattern Recognition and Information Processing, PRIP 2009. Minsk, Belarus, 2009; pp 140–145. (20) Wijetunge, C. D.; Saeed, I.; Boughton, B. A.; Spraggins, J. M.; Caprioli, R. M.; Bacic, A.; Roessner, U.; Halgamuge, S. K. Bioinformatics 2018, 31 (June 2015), 3198–3206. (21) Palmer, A.; Phapale, P.; Chernyavsky, I.; Lavigne, R.; Fay, D.; Tarasov, A.; Kovalev, V.; Fuchser, J.; Nikolenko, S.; Pineau, C.; Nat. Methods 2016, 14 (November), 1–8. (22) Parmar, C.; Velazquez, E. R.; Leijenaar, R.; Jermoumi, M.; Carvalho, S.; Mak, R. H.; Mitra, S.; Shankar, B. U.; Kikinis, R.; Haibe-Kains, B.; Lambin, P.; Aerts, H. J. W. L. PLOS ONE 2014, 9, e102107. (23) R Core Team, R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing: Vienna, Austria, 2017. (24) Gibb, S.; Strimmer, K. Bioinformatics 2012, 28, 2270–2271. (25) Lixin, G.; William, C.; Yu Alex, C. msProcess: protein mass spectra processing. 2011. (26) Laajala, T. D.; Corander, J.; Saarinen, N. M.; Mäkelä, K.; Savolainen, S.; Suominen, M. I.; Alhoniemi, E.; Mäkelä, S.; Poutanen, M.; Aittokallio, T. Clinical Cancer Research: An Official Journal of the American Association for Cancer Research 2012, 18, 4385–4396. (27) Pinheiro, J.; Bates, D.; DebRoy, S.; Sarkar, D.; R Core Team, nlme: Linear and Nonlinear Mixed Effects Models. 2017; R package version 3.1-131. (28) Hothorn, T.; Bretz, F.; Westfall, P. Biometrical Journal 2008, 50, 346–363. (29) Tixier, F.; Hatt, M.; Le Rest, C. C.; Le Pogam, A.; Corcos, L.; Visvikis, D. Journal of Nuclear Medicine: Official Publication, Society of Nuclear Medicine 2012, 53, 693–700. (30) Williams, J. A.; Lalonde, R.; Koup, J. R.; Christ, D. D. Predictive Approaches in Drug Discovery and Development: Biomarkers and In Vitro / In Vivo Correlations; John Wiley & Sons, 2012; Google-Books-ID: ipPrgYmUukYC.

17 ACS Paragon Plus Environment

Analytical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure Captions Figure 1: Synthetic imaging data set I: a tumor tissue morphology invariance assessment. Each image measures 200 x 200 pixels. In gray level invariance, all circular objects in a given image have a diameter of 140 pixels. In shape invariance a circular object has a diameter of 140 pixels, the middle square measures 120x130 pixels, and ellipse largest and smallest axes are respectively 80 and 60 pixels. In the bottom under size invariance, the circular objects have diameters of 60, 140, and 180 pixels. Figure 2: Synthetic imaging data set II for drug homogeneity assessment. (a) There are five basic images and four main scenarios from top to bottom, namely: I) single gray-level image showing a gradual decrease in the drug accessible region, II) and III) after adding different gray levels, IV) after adding salt and pepper noise. The dark zones (small circles inside the main image) simulate the part of the tumor where a drug is absent. In a real MS image they have a value similar to the background but here, for illustration, we have given it a different gray level. Both background and dark regions have not been included in any texture matrix calculation. (b) Plot of DHI and GLSZM-based texture features (SZV, IV) for the above images. DHI was calculated at Nu = 1, except for scenario IV, where Nu = 5 and Nu = 10 was also used. Figure 3: Paclitaxel drug ion images in tumor sections from mice pretreated with bevacizumab (Bev+PTX) and not pretreated (PTX). One representative section is shown for each tumor model and each treatment schedule. The complete results are available in Supplementary Material. Figure 4: DHI plot for Paclitaxel MS images (Iquantized = 8, Nu = 10) from ovarian (A2780- 1A9, IRGOV1, HOC84) and colon (HCT116) cancer cell lines. Blue and red indicate respectively tumor-bearing mice that had received bevacizumab pretreatment or not before paclitaxel.

18 ACS Paragon Plus Environment

Page 18 of 23

Page 19 of 23 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Analytical Chemistry

Tables Table 1: Texture features showing robust performance in different invariance categories for the images given in Figure 1. The features with invariant behavior in all categories are in bold style. Texture Matrix

Size Invariance

Shape Invariance

Gray-level Invariance

GLCM

difference, entropy,

autocorrelation,

maxProbability

contrast, dissimilarity, energy, entropy, homogeneity(1,2), IDMN, IDN, inverse variance, maxProbability

GLRLM

GLSZM

LGLRE, HGLRE,

GLN, LRE, RLN,

RP

RP, SRE

IV, SZV, LIE, HIE,

IV, SZV, ZP, LIE,

SAE, LAE, IV,

DHI

HIE, LISAE,

SZV, ZP, DHI

HISAE, DHI

19 ACS Paragon Plus Environment

Grlv

img3 Page 20 of 23

Shap

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26

img2 Chemistry Analytical

Size

img1

ACS Paragon Plus Environment

B

C Chemistry D Analytical

E

II III IV I

II

value

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 (b) 20 1.00 21 22 0.95 23 24 0.90 25 26 0.85 27 28 29 30 31 32

A

I

(a) 21 of 23 Page

III

1.50

1.50

1.25

1.25

1.00

1.00

0.75

0.75

IV 1.2 1.0 0.8 0.6

0.50 A B C D E

A B C D E

A B C D E

A B C D E

id Texture Feature

ACS Paragon Plus Environment DHI_1

DHI_10

DHI_5

IV

SZV

BevPtx

Page 22 of 23

A2780 HCT116 HOC84

ACS Paragon Plus Environment

IGROV1

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38

Analytical Chemistry Ptx

Page 23 of 23

0.04

0.03

Treatment

DHI

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32

Analytical Chemistry

Bev+PTX PTX 0.02

0.01

A2780-1A9

HCT116 HOC84 IGROV1 ACS Paragon Plus Environment

Cell.line