2.1.

Study Population

Study participants were recruited between March 1, 2000, and September 30, 2004, from the Marvelle Koffler Breast Imaging Centre at Mount Sinai Hospital in Toronto. This study was approved by the Research Ethics Boards of the University of Toronto, Mount Sinai Hospital, and the University Health Network. Inclusion criteria were an analogue standard screening mammogram within approximately $12\text{months}$ prior to recruitment, but not exhibiting any radiological suspicious lesions $(N=232)$ . Exclusion criteria included prior fine-needle aspiration, core biopsies, or any other type of breast surgery including breast reduction or augmentation and any type of tattoos on the breast(s). Additionally, women showing left and right asymmetry in MD, based on classification of mammograms by an expert radiologist ( $>25\%$ difference), were excluded, thus limiting participants to women whose breast tissue retained symmetric tissue aging.

Information concerning participants’ age, menopausal status, height, and weight were collected by means of a self-administered questionnaire. Postmenopausal status was defined as having had no menstrual period for at least $12\text{months}$ . Height and weight were used to calculate BMI defined as weight in kilograms divided by the square of the height in meters.

2.2.

Quantification of MD from Mammograms

Mammographic breast tissue density (MD) was used as the gold standard to evaluate the potential for TiBS to estimate breast cancer risk. All film mammograms in cranial-caudal view (2 $\text{films}\times 232$ volunteers, $N=464$ ) were digitized using a Lumisys Digital Scanner (Kodak, Rochester, New York, USA) at $12\text{-}\text{bit}$ resolution and a pixel pitch of $260\mu \mathrm{m}$ . Digitized mammograms were examined using Cumulus,¹² an interactive density-threshold software. For each image, the trained reader interactively selects two pixel level values using a special outlining tool. The first level selected $\left(iEDGE\right)$ separates the outer edge of the breast from the image background. The second level $\left(iDY\right)$ defines the edges of x-ray dense tissue regions, where all pixels within the outlined region of interest are considered to represent mammographic densities. Additionally, the program enables the user to outline the pectoral muscle thus excluding contributions of the x-ray dense muscle from MD. Pixels between the delineated pectoral muscle and the edge threshold $\left(iEDGE\right)$ represent total breast area, while those above the density threshold $\left(iDY\right)$ represent dense tissue areas. The percent dense area (MD) is defined by the ratio of dense area(s) to the total breast area, multiplied by 100.

MD measurements on all mammograms were performed by two individual raters (KMB, LL) after being trained in the use of Cumulus by an expert rater (NB). Mammograms were presented in a randomized order over eight different sessions ( $N=58$ mammograms per session). To assess the reliability of MD measurements, a randomly selected subset of mammograms (15%) were interspersed throughout the eight sessions ( $N=72$ or nine films per session). For the trained raters, the reading of MD was repeated three times (reads 1 to 3) for the entire data set with a period of at least $1\text{month}$ separating each read. The expert rater only read the mammograms once, including repeats. Table 1 displays the breakdown of study volunteers, as well as MD, demographic, reproductive, and anthropometric information.

Table 1

MD, demographic, reproductive, and anthropometric information of all participants with analogue mammograms.

2.3.

Optical Setup and Procedure

The instrumentation used to gather transillumination spectra was previously described in detail.³⁰ A $50\text{-}\mathrm{W}$ halogen lamp served as broadband light source. UV, part of the visible spectrum, and mid-IR radiation were eliminated using a cut-on and a heat rejection filter, respectively. A total power of $0.25\mathrm{W}$ , covering the $550\text{to}1300\text{-}\mathrm{nm}$ bandwidth, was delivered to the skin via a $5\text{-}\mathrm{mm}$ light guide. Transmitted light was collected via a $7\text{-}\mathrm{mm}$ -diam optical fiber bundle (140 $\mathrm{Si}∕\mathrm{Si}$ fibers, $200\mu \mathrm{m}$ core diameter, numerical aperture: 0.36, P & P Optica, Kitchener, Canada), pointed coaxialles toward the light source. The interoptode distance was provided by a caliper holding the optodes. The source fiber was placed against the skin on the top surface of the breast with minimal compression. Wavelength-resolved detection was achieved using a spectrophotometer (Kaiser Optical Systems, Inc., Ann Arbor, MI, USA) with holographic transillumination grating ( $15.7\text{rules}∕\mathrm{mm}$ blazed at $850\mathrm{nm}$ ) and a CCD (Photometrics, Tucson, AZ, USA) with a spectral resolution of better than $3\mathrm{nm}$ from $625\text{to}1060\mathrm{nm}$ ( $N=436$ wavelengths). These wavelengths were selected as they include the absorption spectra of the primary tissue chromophores, namely, hemoglobin (deoxy- and oxy-), water, and lipid. Health Canada Investigational New Device Class II approval was obtained.

All measurements were taken in the dark, with the participant seated comfortably in an upright position and the breast resting on the support platform containing the attached caliper. A total of eight measurements in cranial-caudal projection were taken per individual, four per breast (center: midline close to the pectoral muscle; medial: $2\mathrm{cm}$ from the inner edge; distal: $2\mathrm{cm}$ behind the nipple; lateral: $2\mathrm{cm}$ from the outer edge). With the measurement of four positions on each breast, an ovoid shaped volume of approximately $23{\mathrm{cm}}^3$ for a $5\text{-}\mathrm{cm}$ -thick breast is sampled.^{30, 31} Numerical modeling demonstrated that 70% of the total breast tissue contributed to 98% of the optical signal. Temporal and spatial reproducibility of the optical measurements is good, as previously addressed.^{30, 31} The entire TiBS procedure takes approximately $15\min$ .

2.4.

Preparation of Spectra for Data Analysis

All spectra were corrected for variations in the wavelength-dependent signal transfer function of the optical system and the thickness of the interrogated tissue, such that all spectra used in further data processing are independent of the instrument and interoptode distance resulting in units of optical density per centimeter (OD ${\mathrm{cm}}^{-1}$ ) for the spectra. For correction of the signal transfer function, spectra were referenced to a $1\text{-}\mathrm{cm}$ -thick ultra-high-density polyurethane transmission standard (Gigahertz Optics, Munich, Germany). Its optical properties (OD $\sim 1.8$ to 2.3 over the wavelength range of interest) were measured separately using an integrating sphere diffuse reflectance setup.³⁰

2.5.

PLS Analysis

Prior to PLS analysis, the error in MD measurements for each rater and between raters (intra- and interrater error, respectively) was determined as this can directly affect the quality of PLS algorithm training and hence the strength of the attainable correlation between TiBS spectra and MD (see the following). Intraclass correlation coefficients (ICCs) were calculated to assess intrarater repeatability for each read (1 to 3) for the two trained raters (KMB and LL) and for the expert rater (NB, one read only). To assess interrater error, two methods were used. First, the mean absolute MD difference between each trained rater and that of the expert was calculated for each mammogram and each read using a mixed linear model (PROC MIXED in SAS 9.1). Second, an interclass correlation coefficient was calculated between each trained rater and the expert rater for each read.

Training of the PLS algorithm was executed on a subset of randomly selected spectra ( $N_T=348$ or 75%) and the predictive power of the algorithm tested by including the remaining 25% of spectra in a validation data set $(N_V=116)$ . The PLS function in MATLAB²⁹ extracts a common spectrum, called a $\mathbf{b}$ vector $(\%\mathrm{cm}{\mathrm{OD}}^{-1})$ , from the entire training spectral data set that when multiplied as a cross-product with an individual’s spectrum $\left(X_i\right)$ (OD ${\mathrm{cm}}^{-1}$ ) produces a scalar or the target, MD $\left(Y_i\right)$ where $Y_i=X_i\mathbf{b}$ . The $\mathbf{b}$ vector identifies those wavelengths, and indirectly also morphological and structural traits, that contribute positively, negatively, or not appreciably to MD $\left(Y_i\right)$ . The 75 versus 25% ratio for training and validation ensures that the training set covers a sufficient range of the variation within the population without over training the system, thereby permitting the PLS algorithm to retain validity on the validation set.²⁹ Because the PLS algorithm uses absolute intensities, spectra were not mean centered as in our previous analysis.²⁸

Since only MD of the entire breast (global assessment) is an established marker for breast cancer risk, ^{9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19} the four spatially collected spectra per breast were considered as a single vector $\left({\mathbf{X}}_i\right)$ for PLS training with the corresponding Cumulus MD as the target $\left(Y_i\right)$ . Therefore, spectra from all quadrants on each breast (center, medial, distal, lateral) were appended to create a single spectrum to approximate global information (Fig. 1 ). Hence, for each individual two spectra with corresponding MD (left and right breasts were treated as separate events) were used as the input $\left({\mathbf{X}}_i\right)$ and target $\left(Y_i\right)$ variables, respectively.

Fig. 1

Example of an appended spectrum from a participant where Cumulus MD $\left(Y_i\right)=1.6\%$ . Order of spectra is center (C), medial (M), distal (D), lateral (L). Each spectrum is measured in OD ${\mathrm{cm}}^{-1}$ over 436 wavelengths $\left(625\text{to}1060\mathrm{nm}\right)$ .

We considered multiple PLS models (Table 2 ) to determine which spectral processing technique on the input data $\left({\mathbf{X}}_i\right)$ and/or the target data $\left(Y_i\right)$ best predicts individual global MD $\left({\widehat{Y}}_i\right)$ . For the target variable $\left(Y_i\right)$ , we examined a model that included only Cumulus MD of the expert rater (NB), as well as models that used the average Cumulus MD of all three raters, to approximate the “gold standard” or best available read (see the following results). We also considered models that excluded mammograms (and associate appended spectra) with greater than 10% absolute disagreement in Cumulus MD $(n=32)$ between any two raters and models where the Cumulus MD of each rater was not averaged and hence each appended spectrum was used up to three times in the training of our PLS algorithm. In some models, detection shot or CCD readout noise seen in the TiBS spectra $\left({\mathbf{X}}_i\right)$ was reduced using either one of two MATLAB defined functions: Boxcar (window $\delta \lambda =7\mathrm{nm}$ ) or Savitsky-Golay (Savgol) smoothing and differentiation with a filter width $\delta \lambda =9\mathrm{nm}$ or a third-order polynomial.²⁹

Table 2

Results of linear regression analysis [Cumulus MD (Yi)=independent variable (x) and PLS predicted MD (Ŷi)=dependent variable (y) ] and Spearman’s rho correlation coefficients for different PLS models for the training and validation data sets.

Spectral processing and PLS algorithm training was executed in MATLAB version 12.0 and statistical analyses were conducted using SAS version 9.1. For each PLS model, the association between the PLS predicted MD $\left({\widehat{Y}}_i\right)$ and the target, Cumulus MD $\left(Y_i\right)$ , was established using Spearman’s rank correlation and linear regression analysis [where Cumulus MD $\left(Y_i\right)=\text{independent}$ variable $\left(x\right)$ and PLS MD $\left({\widehat{Y}}_i\right)=\text{dependent}$ variable $\left(y\right)$ ] for both the training and validation data sets. For each model, we also examined the residuals [i.e., Cumulus MD $\left(Y_i\right)$ minus PLS predicted MD $\left({\widehat{Y}}_i\right)$ ] plotted against Cumulus MD $\left(Y_i\right)$ . For statistical testing, $p$ values $<0.05$ were considered significant.

3.1.

Intra- and Interrater Error in Cumulus MD $\left(Y_i\right)$

The intrarater ICC for the expert (NB) was 0.97. For each trained rater (KMB and LL), the intrarater ICC improved with each read and was highest for the third read (0.96 and 0.93, respectively). The mean absolute difference in Cumulus MD for all mammograms between each rater and the expert also decreased with each read (i.e., from read 1 to read 3). The mean absolute difference between rater 1 and the expert and rater 2 and the expert for the third read was 6.0% (95% CI: 5.4%–6.4%) and 5.4% (95% CI: 4.9%–6.4%), respectively, and both were significantly different from zero $(p<0.0001)$ . The corresponding interrater ICC for rater 1 was 0.92, while for rater 2 it was 0.93.

3.2.

PLS Predicted MD $\left({\widehat{Y}}_i\right)$

Because the third read showed the highest intra- and interrater ICC and the smallest absolute difference in Cumulus MD between each trained rater and the expert, the target $\left(Y_i\right)$ used in all PLS modeling was MD from Cumulus read 3 for the trained raters.

Table 2 displays the results of linear regression analysis and Spearman’s rank correlation coefficients for the association between Cumulus MD $\left(Y_i\right)$ (the independent variable, $x$ ) and PLS predicted MD $\left({\widehat{Y}}_i\right)$ (the dependent variable, $y$ ) for both the training and validation data sets.

Training the PLS algorithm using Cumulus MD of the expert rater only (model 1) or Cumulus MD of all three raters (not averaged) (models 7 to 9) as the target $\left(Y_i\right)$ , yielded the lowest $R^2$ values (0.61 to 0.74) and Spearman’s rank correlation coefficients (rho: 0.80 to 0.84) for the validation data set. Furthermore, for models 7 to 9 the slope of the regression line $\left(\beta \right)$ in the validation set was significantly less than 1.0, suggesting a bias in model development. Averaging Cumulus MD of all three raters (model 2) improved both the $R^2$ value and Spearman’s rho correlation coefficient of the validation set (i.e., versus models 1, 7 to 9) and the 95% CI of the regression slope included 1.0. Conversely, excluding mammograms (and associated spectra) that displayed greater than 10% disagreement between any two raters (model 3) or reducing detection or CCD readout shot noise using either the Boxcar or Savgol MATLAB functions (models 4 to 6) only improved the results for the training data set, not the prediction of MD $\left({\widehat{Y}}_i\right)$ (i.e., versus model 2). Because model 2 had the least assumptions regarding spectral and target processing and because it demonstrated both a large $R^2$ and Spearman’s correlation coefficient in the validation data set, it was considered the most parsimonious model. The PLS $\mathbf{b}$ vector associated with this model is shown in Fig. 2 . Figure 3a is a scatter plot of PLS predicted MD $\left({\widehat{Y}}_i\right)$ versus Cumulus MD $\left(Y_i\right)$ (model 2, validation set only) and Fig. 3b is a plot of the residuals [i.e., Cumulus MD $\left(Y_i\right)$ minus PLS predicted MD $\left({\widehat{Y}}_i\right)$ ] versus Cumulus MD $\left(Y_i\right)$ . Overall, for the majority of women (80%) the estimation of individual MD $\left({\widehat{Y}}_i\right)$ was within 10% of Cumulus MD $\left(Y_i\right)$ , although for the entire data set it ranged between $-26$ and $+22\%$ [Fig. 3b]. The slope of the residuals $(Y_i-{\widehat{Y}}_i)$ for model 2 was $\beta =0.16\%$ (95% CI: 0.08% to 0.24%), significantly larger than zero $(p<0.0001)$ , suggesting that as Cumulus MD $\left(Y_i\right)$ increases so does the error associated with the prediction of MD $(Y_i-{\widehat{Y}}_i)$ [Fig. 3b].

Fig. 2

Resulting PLS $\mathbf{b}$ vector for model 2 showing relative weights in% OD ${\mathrm{cm}}^{-1}$ versus wavelength number for each measurement position.

Fig. 3

(a) Plot of PLS predicted MD $\left({\widehat{Y}}_i\right)$ versus Cumulus MD $\left(Y_i\right)$ (model 2, validation set) showing regression line and 95% CI for individual PLS predicted MD $\left({\widehat{Y}}_i\right)$ ; (b) plot of residuals ( $Y_i$ minus ${\widehat{Y}}_i$ ) versus Cumulus MD $\left(Y_i\right)$ (model 2, validation set).