2.1.

LSS System

A customized LSS setup was developed for application in cardiac tissue. A custom-built spectroscopy probe (Berkshire Photonics, Washington Depot, CT) was used to both illuminate the tissue with broad-spectrum light and gather light from the tissue. The setup shown in Fig. 1 includes a stabilized tungsten-halogen broad-spectrum light source (SLS201L/M, Thorlabs, Newton, New Jersey) providing light in the 350- to 2000-nm range, two Czerny-Turner type charge-coupled diode (CCD) spectrometers (CCS175/M, Thorlabs), a computer with spectrum acquisition software (OSA, Thorlabs), and the spectroscopy probe.

Fig. 1

Spectroscopy system. (a) System overview. Broad-spectrum light was applied to tissue using an illumination fiber. Backscattered light was collected by collection fibers and measured by Czerny-Turner type CCD spectrometers. (b) Schematic view of fiber arrangement showing central illumination fiber and orthogonal fiber pairs for collecting light. (c) Tip of prototype LSS probe.

The probe utilized a single $200\text{-}\mu \mathrm{m}$ core diameter central illumination fiber (FVP200220240, Molex, Phoenix, Arizona) accompanied by four $100\text{-}\mu \mathrm{m}$ core diameter light collection fibers (FVP100110125, Molex). These four light collection fibers were fixed into positions along two orthogonal radial arms adjacent to the central illumination fiber, with two of the four fibers positioned along each of the arms. Figures 1(b) and 1(c) show the arrangement of the fibers at the tip of the probe. The two fibers along each arm were positioned with center-to-center distances of 210 and $345\mu \mathrm{m}$ from the illumination fiber. The fibers positioned $210\mu \mathrm{m}$ from the central illumination fiber were connected to a single CCD spectrometer, and the fibers positioned at $345\mu \mathrm{m}$ were connected to a separate spectrometer. These two sets of fibers are referred to collectively as the R1 and R2 fibers, respectively. The arrangement of two sets of fibers equidistant yet orthogonal in relation to the central illumination fiber was chosen to mitigate the effects of optical anisotropy inherent in cardiac tissue.²⁵ The tip of the probe comprises a 1.7-mm outer diameter stainless steel tube and a 1-mm thick laser-cut quartz glass end cap into which the illumination and absorption fibers are fixed using epoxy. A prior study established the capabilities of this system regarding quantifying depth, arrangement, and composition of fibrotic cardiac tissues.²⁰

2.2.

Spectroscopy of Cardiac Tissue

The LSS system was used to gather spectra from ovine ventricular tissue according to the process outlined in Fig. 2. All animal usage was approved by the Institutional Animal Care and Use Committee (IACUC) at Boston Children’s Hospital. We obtained an $\sim 1{\mathrm{cm}}^3$ transmural tissue sample from the right ventricular free wall (RV) of formalin-fixed hearts from 18 animals with gestational ages ranging 4.3 to 56 months. Two additional samples from the left ventricular free wall and ventricular septum of the two youngest hearts brought the total sample size to 22. Samples were placed on a black foam pad and kept submerged in phosphate-buffered saline (PBS) before and during spectroscopy. Black open-cell foam prevented backscattering of photons after passing through the sample. Samples were examined using LSS to identify changes in ND. The LSS probe was positioned on the epicardial surface at multiple locations to gather 20 spectra for each sample. The spectra were gathered with an integration time of 200 ms. Spectra were recorded with full width at half maximum resolution of 0.6 nm in the wavelength range of 500 to 1100 nm. Following acquisition, the spectra were saved for analysis. Tissue samples were stored in PBS for histology. Each spectrum was normalized to its mean intensity. The spectra were then filtered using a 1D gaussian filter with a standard deviation of 20 measures ($\sim 3.3\mathrm{nm}$) and a kernel size of 80 measures.

Fig. 2

Experimental design. (a) Spectra were gathered from $\sim 1{\mathrm{cm}}^3$ transmural samples of fixed ovine hearts. (b) A $100\text{-}\mu \mathrm{m}$ thick section from the transverse center of the tissue sample at the same location as for spectroscopy. The section of tissue was stained with WGA and DAPI. The spectra from (a) and measured ND from (b) were applied for (c) machine learning. Also, spectra and measured ND were used to test the prediction of NDs.

2.3.

Histology

After spectral measurements, tissue samples were sectioned within the probed region at the center of each sample. A vibratome (Leica Model VT1200S, Wetzlar, Germany) was used to section at $100\text{-}\mu \mathrm{m}$ thickness perpendicular to the epicardial surface and parallel to the transmural plane of the ventricular wall. This orientation was chosen to include light scattering structures throughout the thickness of the ventricular wall. After sectioning, glycoconjugates of the extracellular matrix and glycoproteins of cell membranes were labeled using wheat germ agglutinin (WGA) conjugated to a fluorophore (WGA CF488A Conjugate, Biotium, Fremont, California). Cell nuclei were stained with $1\mu \mathrm{M}$ 6-diamidino-2-phenylindole (DAPI, D3571, Life Technologies, Carlsbad, California) according to a previously established protocol.²⁶ All sections were washed in PBS and mounted on a glass slide with Fluoromount-G (#17984-25, Electron Microscopy Science, Hatfield, Pennsylvania).

A Leica SP8 confocal microscope with a $40\times$ oil immersion lens (numerical aperture 1.3) was used to scan each section. Two-dimensional tile scans with a pixel size of $0.2\mu \mathrm{m}$ spanning $695\mu \mathrm{m}\times 2467\mu \mathrm{m}$ were performed for the capture of cardiac tissue, including the epicardium. DAPI and the WGA CF488A Conjugate were excited with a 405 nm and 488 laser, respectively.

Two or more different tile scans from each sample were used to calculate the ND. Manual segmentation was performed using the Fiji/ImageJ image processing software package.²⁷ Histogram-based thresholding (mode + 2 standard deviations) isolated nuclei in the DAPI images. Gaussian blurring with a sigma of 2 followed by erosion and dilation removed small segments. Connected nuclei were separated using a watershed function. Visual inspection of segmented images validated segmentation methodology. Automated segmentation and counting were performed to validate the ND values gathered from the manual segmentation. This was performed using Python 3.8.2 while following the same operands as the manual segmentation. The image processing was primarily done using the package scikit-image (0.17.2).²⁸ For both the manual and automated methods, ND was calculated by dividing the number of detected nuclei by the tissue area. The average ND produced by the manual segmentation for each tissue sample provided the ground truth for machine learning.

2.4.

Assessment of the Spectra–ND Relationship by Cluster Analysis

Unsupervised and supervised machine learning approaches were used to explore whether spectra gathered from the LSS system allowed for differentiation between the ND of the preparations. We employed MATLAB (2019b, MathWorks, Natick, Massachusetts) and its Deep Learning Toolbox for unsupervised learning evaluation of the correlation between LSS spectra and the measured ND.

We explored normalized, filtered, and concatenated spectra from the paired collection fibers. Each spectrum was annotated with the ND of the sample from which it was measured. We applied the MATLAB function spectralcluster to identify clusters of spectra from preparations of different NDs using the Euclidean distance metric. From this spectral clustering analysis, we generated a similarity graph between nodes of the first and second principal components of the spectra.²⁹ Analysis of the eigenvalues from the spectral cluster function indicated the spectra fell into five distinct groups. We utilized the MATLAB function gscatter to visualize the spread of the spectral set, with the resulting groups identified by color. Groups were associated with the NDs of the spectra contained in each group. The principal components analysis (PCA) was performed, and the first two principal components were used to visualize any potential groups within the dataset. A one-way ANOVA and Tukey-Kramer post hoc test using a significance level of 0.05 was used to identify differences between the clustered spectra. A box and whisker plot was created from these data, which compared the NDs associated with the spectrum within each group produced from the spectral clustering. Significant differences between clusters were indicated by brackets.

2.5.

Assessment of the Spectra–ND Relationship by CNN

We applied CNNs to differentiate NDs based on spectra from LSS. This work was built upon our previous work on identifying fibrosis from in vivo confocal imaging.¹⁹ Our software framework was based on Python scripts and the Keras neural network API (2.3.1)³⁰ with Tensorflow (2.2.0).³¹ We designed a CNN to classify the spectra into three groups of ND 0 to 2000, 2000 to 3800, and above $3800\text{nuclei}/{\mathrm{mm}}^2$.

The CNNs were systematically trained and validated using a form of k-fold cross-validation called leave-one-out cross-validation (LOOCV).³² This approach entailed creating 22 different unique datasets, each of which removed the spectra gathered from one of the 22 different samples. This removed sample of spectra was separated for use as a testing dataset for each fold. Another sample was also randomly selected from the then remaining 21 samples in each fold to be removed and used as a validation dataset to guide hyperparameter optimization. In all, for each LOOCV performed, 22 networks were trained, each on a training dataset that consisted of the spectra from 20 different samples, less the validation and training samples. To address and reduce the likelihood of the networks overfitting the datasets, initial network hyperparameters were selected that minimized the number of trainable parameters.

Two tests were performed to determine the potential sensitivity of the networks to specific wavelength ranges and to also determine the sensitivity of the networks to the resolution of the spectra. An initial network configuration was selected (Table 1). Hyperparameter variations were later performed using this network structure as a base. To determine the sensitivity of the networks to wavelength ranges, three training sessions of LOOCV were performed on spectra consisting of wavelengths 500 to 700 nm, 700 to 900 nm, and 900 to 1100 nm, respectively. One-way ANOVA and posthoc test using Fisher’s least significant difference method with a level of 0.05 were applied to determine any significant differences between the accuracies reported for each wavelength range. To determine the sensitivity of the networks to the resolution of the data, 15 LOOCV training sessions were performed using spectra with resolutions reduced by 15 different scales, K, distributed logarithmically between 10 and 1000. The spectra were downsampled using a 1D linear interpolation between every Kth point in the original spectra. When using the network structure outlined in Table 1 in conjunction with the downsampled spectra, the number of trainable parameters within the networks varied between 1755 and 135, for scales of 10 and 1000, respectively. Beyond the hyperparameters described in Table 1, both the wavelength-dependency and the resolution-reduction LOOCV training sessions were performed using a learning rate of 0.005 for all networks. Training for these sessions was terminated early if the loss did not decrease by more than 0.001, a value referred to as the “minimum delta” during parameter variations, after 100 epochs; otherwise, the training was permitted to run up to 1000 epochs. The categorical cross-entropy was calculated from the logits produced from the final layer of the network for use during the training process and post-analysis. This loss function was selected because each spectrum belongs to only one of the three ND ranges.

Table 1

CNN Configuration.

Although hyperparameter variation was not performed during these training sessions, a validation sample was still removed from the training dataset in addition to the testing sample. Network weights were saved during the training process when the loss value decreased. The network weights resulting in the lowest loss value were restored at the end of the training session. The predicted classes from the testing datasets for each fold were then compared to the assigned classes for each spectrum based on ND values. From this comparison of aggregated predictions, we calculated the accuracy for the CNN topology on the entire LOOCV as follows:

To improve the accuracy of the predictions resulting from the LOOCV, the above training process was altered to introduce hyperparameter variations within each training fold. Multiple networks were trained using the base topology outlined in Table 1 to maintain a low number of trainable parameters while varying the learning rate. The resolution of the spectra was also reduced by a factor of 50 using the methods described previously. This reduction decreased the number of trainable parameters while still maintaining most of the variations within the spectra that remained after filtering and produced a network topology with 441 trainable parameters. Ten different learning rates, logarithmically distributed between 0.0001 and 0.01, were used as well as four different minimum delta values: 0.001, 0.003, 0.005, and 0.007. The randomly selected validation sample associated with each network was then used to determine the validation accuracy and loss values of the network. The optimal learning rate and minimum delta for each network were selected based on the highest accuracy value reported. In cases where more than one hyperparameter set resulted in the same highest accuracy value, the hyperparameter set with the lowest validation loss score between them was selected. Once optimal hyperparameter values were selected for each of the 22 folds of the LOOCV, the optimal networks for each fold were used to predict the classes of their respective testing datasets. The average accuracy value of these predictions was then calculated using Eq. (1).

3.1.

Spectroscopy of Tissue Samples

We gathered 440 spectra from two collection fibers from 22 samples. Raw spectra displayed a moderate degree of variability in intensity within samples and in-between samples. Spectra from R1 and R2 were different with respect to intensity range and shape [Figs. 3(a) and 3(d)]. After normalization to their respective means, spectra from an individual sample exhibited a high correlation with the other spectra from the same sample ($R^2\ge 0.999$). Differences between spectra from different samples were subtle ($R^2=0.996\pm 0.002$). Figures 3(b) and 3(e) show mean spectra with standard deviations from four representative age groups. Figures 3(c) and 3(f) highlight regions within the spectra that present differences across the dataset associated with age.

Fig. 3

Comparison of spectra for selected gestational ages within the ranges of 0-5, 5-8, 8-13, and 13+ months. (a) Raw spectra from R1. (b) Averaged spectra for each age range from R1 after normalization and filtering. The black box indicates the region highlighted in 3C. (c) Enlarged view of the region indicated with a black box in 3B. This region was selected as it shows where the four age ranges differ significantly from each other. (d) Raw spectra from R2. (e) Averaged spectra for each age range from R2 after normalization and filtering. The black box indicates the region highlighted in 3F. (f) Enlarged view of the region indicated with a black box in 3E. This region was selected as it shows that spectra from tissues of different age ranges differ from each other.

3.2.

Histology and Determination of ND

Representative transmural confocal microscopy tile scans of RV free wall myocardium are shown in Fig. 4. Images of DAPI revealed higher ND and smaller myocytes in samples for younger animals vs. samples from older animals (Fig. 5). Decreases of ND in samples of 4.3, 6.5, and 56 months gestational age are visible in $250\mu \mathrm{m}\times 250\mu \mathrm{m}$ confocal images shown in Figs. 5(a)–5(c). Increased size of cardiac myocytes is visible from WGA images of the same samples is evident in Figs. 5(d)–5(f). Composite images are presented in Figs. 5(g)–5(i).

Fig. 4

Transmural section of heart wall from a 4.3-month old animal from confocal microscopy. (a) Nuclei stained with DAPI. (b) Extracellular matrix staining with WGA.

Fig. 5

Confocal images of cardiac tissue at gestational ages of (a) 4.3, (b) 6.5, and (c) 56 months. (a)–(c) DAPI labeled cell nuclei. (d)–(f) WGA labeled extracellular matrix and capillaries. (g)–(i) Composite DAPI and WGA images.

Manual segmentation using ImageJ indicated a 72.5% difference in ND between samples of 4.3 months and 10 months gestational age and a 78.2% difference in ND from 4.3 to 56 months. Automated segmentation yielded a 65.4% and 77.8% decrease in ND for the same age differences. The reduction of ND with age in our samples is summarized in Fig. 6(a) and consistent with prior studies.²³ ND differences between segmentation methods varied most for samples from the youngest hearts. Figure 6(b) shows the measured values for each segmentation method. The correlation between segmentation methods was high ($R^2=0.94$). Although the results from both the manual and automated segmentation techniques were similar, the values resulting from the manual segmentation were selected as ground truth for the duration of the study as they more closely reflected the reduction in ND over time in the literature.²²

Fig. 6

Counting of cell nuclei. (a) Scatterplot and regression fit of measured NDs relative to age. (b) Comparison of NDs from different segmentation methods.

3.3.

Cluster Analysis of Spectra

Dimensionality reduction via PCA revealed that the first two principal components accounted for $\sim 95\%$ of the variance between spectra. The cluster analysis revealed five clusters [Fig. 7(a)]. Two of the clusters contained only 20 spectra, and each contained spectra from only a single sample, whereas the other three clusters contained over 100 spectra with at least five samples in each cluster. Only a small percentage of spectra were not assigned to the same cluster as the majority of the spectra from the same tissue sample. In sum, a total of 23 of the 440 spectra ($\sim 5\%$) were miscategorized. Clusters were labeled 1 to 5 with increasing ND. Differences in ND were significant for many of the clusters [Fig. 7(B)]. For instance, NDs of clusters 3, 4, and 5 were different from clusters 1 and 2. Differences of ND for clusters 1 and 2, corresponding to the samples of lowest ND, were not significant.

Fig. 7

Spectral clustering of spectra measured from cardiac tissues. (a) The clustering identified 5 clusters shown in a principal component plot. (b) NDs of several identified clusters were different. Bracket identified clusters that are statistically different from one another at a 5% confidence level.

3.4.

CNN-Based Classification of ND Grouping

We explored the sensitivity of the CNN to different wavelength ranges. Accuracies of $59\pm 46\%$, $40\pm 48\%$, and $76\pm 33\%$ resulted from the wavelength ranges 500 to 700 nm, 700 to 900 nm, and 900-1100 nm, respectively. Differences between the 500 and 700 nm and 700 and 900 nm ranges were not significant. Also, differences between the 500 and 700 nm and 900 and 1100 nm ranges also were not significant. Accuracies for the 900 and 1100 nm were higher than for 700 to 900 nm ($p$-value: 0.0075). The exploration of the sensitivity of the network to the resolution of the spectra resulted in accuracies ranging from 67% to 95%. There were no statistically significant differences between the LOOCV sets based on the reduction in data resolution. Two reduced resolution datasets resulted in accuracies above 90%, i.e., the datasets reduced by factors of 37 and 268. The dataset reduced by a factor of 1000 resulted in an accuracy of $83\pm 31\%$. The accuracy from the optimal hyperparameters for each fold of the training process was $95\pm 12\%$. Classifications of spectra for this final network are summarized in the confusion matrix in Fig. 8. Among the 22 optimized networks, all the minimum delta values of 0.007, 0.005, 0.003, and 0.001 were selected 2, 6, 7, and 7 times each. Of the possible learning rates only five, 0.01, 0.00615, 0.00378, 0.00233, and 0.00143, were selected. These learning rates occurred 5, 9, 5, 2, and 1 time, respectively.

Fig. 8

Confusion matrix of aggregated CNN predictions of NDs from spectra. The resulting classifier accuracy was higher than 90%.