2.1.

Tissue Excision, Sample Processing, and 3D Model Creation

Tissue samples were extracted, prepared, sectioned, and histologically stained following established protocols from Refs. 27 and 28, and reconstructed 3D models from Ref. 28 were used in this study. In summary of the methods described in Ref. 28, the tissue samples were excised, sectioned into $4\mu \mathrm{m}$ thick sections spaced approximately $25\mu \mathrm{m}$ apart and stained using an automated Masson’s Trichrome staining method. The sections were imaged on an automated slide scanner at a $0.44\mu \mathrm{m}/\text{pixel}$ resolution. The resulting images were registered and stacked to create serial images of the section tissue samples.

These serial images were then segmented using a random-forest segmentation network created in the software Ilastik.²⁹ Most of the technical aspects of random forest parameter optimization, such as the number of trees, nodes, and other hyperparameters, were managed by the Ilastik software. For the model used in this work, the default value of 100 trees was used. The ground truth segmentations used for training the segmentation classifier were created in the Ilastik software using the provided image annotation tools. Forty images sampled randomly from 7 different nodal tissue sections were included in the ground truth segmentations used to train the random forest network. The model was trained to segment four different classes including myocardium, connective tissue, and separate classes for nuclei within the myocardium and connective tissue regions. As described in Ref. 28, the segmentation results were inspected and verified by an external expert cardiac pathologist. The details of post-processing the segmentation results to refine the nuclear segmentation are described in a later section. These segmented images were then used to create 3D reconstructions of the excised tissue samples. The average distance between slides was calculated using methods outlined in Ref. 28, and this distance was used to interpolate values for the connective and myocardial tissue samples between histological sections in the direction of sectioning using methods outlined by Schenk et al.³⁰ The final voxel size of the 3D tissue composition model used for this work was $0.44\mu \mathrm{m}\times 0.44\mu \mathrm{m}\times 7\mu \mathrm{m}$, where the third dimension is the direction of sectioning and interpolation. Figure 1 provides an overview of this process.

Fig. 1

Pipeline used to create 3D reconstructions from excised cardiac tissue samples containing nodal tissue regions.

2.2.

Tissue Prep for Data Acquisition

Following extraction of the nodal tissue regions of interest and before tissue sectioning and histological staining, both the sinus and AV nodal region samples had three fiducial markers placed transmurally on the epicardial surface for the sinus nodal regions and through the membranous septum and ventricular septum for the AV nodal regions. These fiducial markers were created using a 1.0 mm diameter disposable biopsy punch (Rapid-Core 1.0, World Precision Instruments, Sarasota, Florida, United States).

In preparation for LSS spectra acquisition, each tissue sample was placed onto and pinned to a thin sheet of open-pore black foam, preventing the backscattering of photons after passing through the sample. The well dish was then filled with 1× phosphate buffered solution to the point of submersion of the tissue within the solution.

The custom-designed LSS probe was connected to a tungsten-halogen broad-spectrum light source (SLS201L/M, Thorlabs, Newton, New Jersey, United States), which emitted light in the 350 to 2000 nm range and two Czerny-Turner type charge-coupled diode (CCD) spectrometers (CCS175/M, Thorlabs). The spectrometers were also connected to a computer running custom spectrum analysis software (Cite LSS Nuclear paper). The tip of the LSS probe was then attached with a custom fixture, holding the probe normal to the workbench’s surface and pointing down. The fixture was then attached to a three-axis micromanipulator (MP-285, Sutter Instrument Company, Novato, California, United States) magnetically fixed to the workbench. The well dish containing the excised tissue sample was placed directly beneath the LSS probe tip and aligned with markings on the workbench for consistent placement and orientation. SN samples were oriented by placing the epicardial surface facing up, with the superior vena caval aspect toward the left of the acquisition setup and the atrial appendage aspect toward the right of the setup. AV nodal tissue samples were oriented with the right atrial aspect facing up, the membranous septum aspect toward the top of the acquisition setup, and the ventricular septal aspect toward the bottom of the setup. The acquisition setup’s top/bottom and left/right orientations correspond to the $x$ and $y$ axes of the micromanipulator fixture, respectively. An illustration of this setup can be found in Fig. 2(a).

Fig. 2

Placement of the custom LSS probe in a fixture and used to gather spectral data in an orientation normal to the surface of the tissue. (a) Placement of the tissue onto a black open-pore foam within a well-dish. (b) Example grid pattern created and used to gather data across the surface of the tissue. Yellow spots indicate individual, unique spectral sampling locations.

Custom software written in MATLAB (R2022b, MathWorks, Natick, Massachusetts, United States) was used to control the micromanipulator and extract spectral measurements from the two spectrometers attached to the LSS probe. This software also recorded the location in 3D space of where each spectral measurement was taken.

2.3.

LSS Data Acquisition

The 3D locations of the three fiducial marks in the tissue were acquired by using the micromanipulator to position the LSS probe tip concentric with the fiducial mark above the tissue, and then lowering the tip of the probe to contact the surface of the tissue and recording the tip location. Each fiducial mark location was sampled three times, and the probe tip was moved away from the tissue between each measurement before performing a subsequent measurement.

After recording the locations of the fiducial marks, LSS data acquisition locations were mapped across the surface of the tissue sample in a diamond grid pattern, with $280\mu \mathrm{m}$ between every individual sample location across the tissue surface [Fig. 2(b)]. A trained operator confirmed sufficient contact of the probe tip with the tissue surface at each grid location. After confirming the contact of each grid position with the tissue surface, automated acquisition of spectral samples using the LSS probe was initiated within the custom MATLAB software.

Three hundred spectral samples were acquired at each grid location. Each spectrum was sampled with an integration time of 70 ms and a full width at half maximum resolution of 0.6 nm within the 500 to 1000 nm wavelength range. Upon measuring the 300 spectra, their values at each wavelength were averaged and stored as a single spectral measurement associated with the 3D location at the grid point.

After the topographical spectral measurements were completed, the 3D locations of each of the three fiducial markers were recorded following the same protocol as before grid sampling. The grid locations and spectral data were resampled if any drift or notable change was noted between the initial fiducial mark locations and the post-spectral acquisition locations.

Once the tissue samples were fully sectioned, stained, imaged, and digitally processed, a manual segmentation process was used to identify the location of the fiducial markers within the tissue samples using the image processing software FIJI.³¹ These segmentations were then processed following the same protocol as the myocardial and connective tissue segmentations to create 3D volumes for each fiducial marker in the tissue samples.

2.4.

Nuclear Segmentation Verification

The Ilastik-based nuclear segmentation was validated on serial histological sections of ovine cardiac tissue. The samples were acquired, excised, sectioned, and stained following the protocols outlined by Johnson et al.²⁷ Four serial section pairs of tissue, all acquired within $400\mu \mathrm{m}$ of each other from within the tissue block, were identified for nuclear segmentation validation. One section from each pair was stained, imaged, and segmented using the Masson’s trichrome-based segmentation described previously. The other tissue section from each pair was stained using a DAPI staining medium (#D1306, ThermoFisher Scientific, Waltham, Massachusetts, United States) and imaged using the same automated slide scanner (Carl Zeiss AG, Oberkochen, Germany). DAPI nuclear segmentation and post-processing image refinement involved thresholding the DAPI image at two standard deviations above the median intensity value, applying a water-shedding algorithm to separate the nuclei, and filtering out segmented objects that were smaller than $20{\mu \mathrm{m}}^2$ and larger than $100{\mu \mathrm{m}}^2$. Nuclear densities within and without the segmented nodal regions were calculated for both sinus and AV nodal regions. Nodal tissue segmentations were performed according to methods outlined in Ref. 28. A trained cardiac pathologist validated segmentations of other tissues, including connective, myocardial, and nodal tissues.

2.5.

3D Model Exclusion Criteria

After the stained tissue slides were prepared, tissue samples were excluded from the analysis if they contained at least five slides that were stained insufficiently. Insufficient staining is defined in this context as including inconsistent and abnormally light or heavy staining across the tissue section. Additionally, tissue sections were excluded from subsequent analysis if more than five slides contained extensive damage to the tissue or excessive inflammation or calcification.

2.6.

LSS Point Registration Methods

Registration of the LSS probe grid locations to the surface of the 3D tissue model was performed by registering a minimum of 10 paired points corresponding to the fiducial marks and other landmarks within the 3D model and the probe grid. The locations of the fiducial markers were identified as the centroids of the segmented 3D volumes associated with each fiducial mark placed in the tissue samples. A transformation matrix was found by optimizing the Procrustes distances between the paired points.^32,33 The average distance between serial tissue sections was determined by altering the distance between sections and optimizing the registration between the probe grid fiducial marker locations and their corresponding paired points along the surface of the 3D tissue model. These methods resulted in applying a similarity transform to the LSS data points, associating the spectra to a registered location on the tissue surface of a 3D reconstructed tissue model.

2.7.

Volume Composition Information Extraction

After registering the probe grid locations onto the surface of the 3D tissue models, tissue composition information was extracted from the 3D models at each point within the LSS probe grid. A conical frustum volume was used to determine the 3D volume of interest associated with each point in the LSS probe location grid. The top diameter of the frustum was set to $540\mu \mathrm{m}$, which is the maximum distance between sensing fibers at the tip of the LSS probe. The angle of the sides of the frustum was set to the full acceptance angle, also called the divergence angle, of the sensing fibers, 12.71 deg, which was calculated using Eq. (1) using a numerical aperture (NA) of 0.22 for the silica optical fiber (FVP100110125, Molex)

Figure 3(a) shows how the frustum relates to the tip of the custom LSS probe. The frustum was then oriented in 3D space such that the top surface was positioned at a given 3D location in the probe grid concentric with the orientation of the LSS probe and the Z axis. The frustum volume extended transmurally through the tissue, containing tissue spanning the entire depth of the tissue sample directly below the given location. An example of how this extends through segmented tissue can be seen in Fig. 3(b). All 3D model voxels contained within the frustum were included in the tissue composition measurements, and statistics, such as nuclei per μm and muscle volume fraction, were associated with the spectral measurements acquired at the given point. All image processing and association of LSS spectra with tissue composition information were performed using custom software written in MATLAB.

Fig. 3

Relationship between the tip of the LSS probe, calculated frustum, and segmented tissue models. (a) Probe tip describing how the sensing fibers relate to the initial diameter of the frustum, as well as how the angle calculated from the numerical aperture of the fibers relates to the increasing diameter of the frustum as a function of distance from the probe tip. (b) Cross-section of segmented tissue. Dark blue indicates segmented connective tissue, green indicates myocardium, light blue and yellow indicate segmented nuclei; dark red indicates the center of the location where the probe was placed; and the light red overlay indicates the region, based on the frustum that would be included in the analysis for the given point on the tissue surface.

2.8.

Preprocessing the LSS Signals

All LSS spectra were calibrated using a spectrum from a 99% reflectance factor white diffuse standard reference (Spectralon, Labsphere Inc., North Sutton, New Hampshire, United States). Calibration spectra were re-acquired before measuring LSS data from each tissue sample. These individual calibration spectra were used only with the LSS spectra from the associated tissue sample. The spectra were then individually $z$-scored and grouped according to the tissue sample. Mean and standard deviations used in $z$-scoring were calculated within each tissue sample group of spectra. After preprocessing, the spectra from the 210 and $315\mu \mathrm{m}$ sensing fibers for each measurement location were concatenated to facilitate one-dimensional input into the linear regression methods used in this study. All preprocessing of spectra and training of regression models was performed in custom code written in the Python programming language.

A principal component analysis was performed on the LSS spectra by extracting the principal components from the spectra and visualizing the first two components and their association with the nuclear density and muscle volume fractions. Additionally, dimensionality reduction was also performed using the non-linear uniform manifold approximation and projection (UMAP) method.³⁴

2.9.

Train-Test Split and Leave-One-Out Cross-Validation

The calibrated, $z$-scored, and concatenated spectra were then used for training both linear regression models and multi-layer perceptron regressors to predict nuclear density in nuclei per $\mu {\mathrm{m}}^3$ and muscle or connective tissue volume fraction as a percent of the volume of the extracted frustum. After preprocessing, the input for each model was a collection of calibrated, $z$-scored spectra acquired from the custom LSS probe. The output or ground truth associated with each individual spectra was either a nuclear density value or a fraction. To facilitate training and model convergence, both ground truth measurements were scaled to a range between 0 and 1000 to improve convergence while training linear and neural network-based regression methods.

The data were split into training and testing datasets following an LOOCV method.³⁵ This method entailed creating a training dataset that included all the spectra minus the spectra and data associated with a single tissue sample. This “hold out” set of data, which consisted of all the spectra associated with the given tissue sample, was then used to test both the linear and non-linear regressor models trained on the given training dataset. A total of 14 train-test splits were created using this method, where each test set contained all the spectra associated with a given, unique tissue sample that was excluded from the training set in the given test-training split.

2.10.

Linear Approximation Methods

The processed spectra were then used, in conjunction with the scaled muscle volume fraction and nuclear density ground truth values, to train a linear regressor model. The model was trained in Python using the scikit-learn package. A linear regressor using $L1$ and $L2$ data priors as regularization parameters, also called an Elastic Net,³⁶ was trained and evaluated using the LOOCV approach described earlier. A small hyperparameter search was carried out and determined that an alpha regularization coefficient of 1.0, with a maximum number of iterations of 100 produced the most stable, consistent results throughout the cross-validation. A new linear regressor model was trained for each of the 14 train-test splits in the cross-validation. Default values were kept for the $L1$ ratio, and alpha parameters were used for training. The results from the individual test predictions were aggregated and evaluated using r-squared and concordance correlation coefficients.

2.11.

Weight Analysis of Linear Models

The input weights of the linear model were collected from each training session within the cross-validation. The weights were then averaged across the training sessions and visualized for inspection to determine their relative importance in the outcome of the regression analysis.

2.12.

Neural Network Analysis

A small, fully connected neural network was created and trained as a regressor model to predict muscle volume fraction and nuclear density from the processed spectra. This network consisted of an input layer and two fully connected layers, each of which was subsequently followed by a batch normalization and dropout layer. A hyperparameter search that included the dropout percentage and number of nodes in each fully connected layer was carried out, resulting in a final network topology with values of 35 and 0.2 for the number of nodes and dropout percentage, respectively.

The neural networks were trained and evaluated following the same LOOCV method as the linear model. A new model was trained for each of the 13 train-test splits in the cross-validation. The networks were trained for a total of 40 epochs with a learning rate of 0.1.

3.1.

Nuclear Segmentation Validation

The average area of nuclei segmented from DAPI-stained reference tissue images was $43.55\pm 18.96{\mu \mathrm{m}}^2$. The average area of nuclei segmented using the random-forest segmentation algorithm on Masson’s trichrome-stained slides was $36.03\pm 14.1{\mu \mathrm{m}}^2$. The mean number of nuclei counted on the DAPI reference slides was $\mathrm{91,807}\pm \mathrm{12,499}$ nuclei, with the mean number of nuclei counted on the Masson’s slides being $\mathrm{91,116.63}\pm 9877$ nuclei. The difference between the two as a percentage of the count from the DAPI reference slides was less than 1%.

The nuclear density within the sinus nodal regions of the tissue samples was found to be $5.3•{10}^{-5}\text{nuclei}/{\mu \mathrm{m}}^3$. This was found to be greater than the surrounding tissues ($P\ll 0.05$), which had an average nuclear density of $3.7•{10}^{-5}\text{nuclei}/{\mu \mathrm{m}}^3$. The nuclear density within the AV nodal regions was $6.24•{10}^{-5}\text{nuclei}/{\mu \mathrm{m}}^3$, which was found to be less than the density of the surrounding tissues at $7.87•{10}^{-5}\text{nuclei}/{\mu \mathrm{m}}^3$ ($P<0.05$).

3.2.

Dimensionality Reduction Analysis

Visualizations of the first two principal components of the spectral data can be seen in Fig. 6. This principal component analysis demonstrates the correlation between the nuclear density and the spectral data [Fig. 6(a)] and the muscle volume fraction and the spectral data [Fig. 6(b)] as a gradient increasing toward the bottom right of the plot. The results from the UMAP analysis demonstrated greater sensitivity to the batch effect, clustered the spectra more according to tissue sample, and demonstrated less correlation with the nuclear density and the muscle volume fraction.

Fig. 6

Principal components and associated tissue composition features. (a) Correlation of the nuclear density with the first two principal components. (b) Correlation muscle volume fraction with the first two principal components.

3.3.

Linear Regression Analysis

A total of 14 elastic net regressor models were trained following the LOOCV approach. The nuclear density regressor achieved a prediction r-squared value of 0.64 and a concordance correlation coefficient of 0.78. The muscle volume fraction regressor achieved a prediction r-squared of 0.41 and a concordance correlation coefficient of 0.59. Figure 7 shows the aggregated test-set predictions of the individual regressor networks across the cross-validation.

Fig. 7

Relation between predicted values and ground truth resulting from the elasticnet regressor via the aggregated LOOCV prediction results. (a) Results from the prediction of nuclear density. (b) Results from the prediction of muscle volume fraction.

Figure 8 shows the weights of the nuclear density elastic net regressors for both the $215\mu \mathrm{m}$ sensing fibers and the $310\mu \mathrm{m}$ sensing fibers. This analysis indicates that the nuclear density regressors weighted wavelengths below 600 nm and wavelengths between 650 and 700 nm highly for spectra acquired from the 215 nm sensing fiber. The nuclear density regressors also weighted wavelengths below 600 nm highly for measurements from the 310 nm sensing fiber but weighted other frequency bands differently, including between 800 and 900 nm as well as 950 to 1000 nm.

Fig. 8

Elasticnet linear model weights aggregated and averaged across all LOOCV training sessions associated with the nuclear density. The light blue trace shows raw weight values, and the red trace indicates the weights after a smoothing function was applied. (a) Weight values per wavelength for the sensing fiber $210\mu \mathrm{m}$ from the illumination fiber. (b) Weight values per wavelength for the fiber that is $315\mu \mathrm{m}$ from the illumination fiber.

Figure 9 shows the weights associated with the muscle volume fraction elasticnet regressors. This analysis indicated that, unlike the nuclear density regressors, the muscle volume fraction regressors applied a more similar weighting to both the 215 and $310\mu \mathrm{m}$ sensing fibers, placing higher weights on wavelengths between 600 and about 700 nm. The volume fraction regressors also placed higher weights on wavelengths greater than 850 nm for the $215\mu \mathrm{m}$ sensing fiber and greater than 800 nm for the $310\mu \mathrm{m}$ sensing fiber.

Fig. 9

Elasticnet linear model weights aggregated and averaged across all LOOCV training sessions associated with the nuclear density. (a) Weight values per wavelength for the sensing fiber $210\mu \mathrm{m}$ from the illumination fiber. (b) Weight values per wavelength for the fiber that is $315\mu \mathrm{m}$ from the illumination fiber.

3.4.

Neural Network Regression Analysis

Similar to the linear regression analysis, 14 fully connected neural networks were trained and tested following the LOOCV approach. The nuclear density regressor achieved an r-squared value of 0.62 and a concordance correlation coefficient of 0.79. The muscle volume fraction regressor achieves an r-squared of 0.41 and a concordance correlation coefficient of 0.63, producing results nearly identical to those from the linear regressor.