High-speed 2D light-sheet fluorescence microscopy enables quantification of spatially varying calcium dynamics in ventricular cardiomyocytes

Introduction: Reduced synchrony of calcium release and t-tubule structure organization in individual cardiomyocytes has been linked to loss of contractile strength and arrhythmia. Compared to confocal scanning techniques widely used for imaging calcium dynamics in cardiac muscle cells, light-sheet fluorescence microscopy enables fast acquisition of a 2D plane in the sample with low phototoxicity. Methods: A custom light-sheet fluorescence microscope was used to achieve dual-channel 2D timelapse imaging of calcium and the sarcolemma, enabling calcium sparks and transients in left and right ventricle cardiomyocytes to be correlated with the cell microstructure. Imaging electrically stimulated dual-labelled cardiomyocytes immobilized with para-nitroblebbistatin, a non-phototoxic, low fluorescence contraction uncoupler, with sub-micron resolution at 395 fps over a 38 μm × 170 µm FOV allowed characterization of calcium spark morphology and 2D mapping of the calcium transient time-to-half-maximum across the cell. Results: Blinded analysis of the data revealed sparks with greater amplitude in left ventricle myocytes. The time for the calcium transient to reach half-maximum amplitude in the central part of the cell was found to be, on average, 2 ms shorter than at the cell ends. Sparks co-localized with t-tubules were found to have significantly longer duration, larger area and spark mass than those further away from t-tubules. Conclusion: The high spatiotemporal resolution of the microscope and automated image-analysis enabled detailed 2D mapping and quantification of calcium dynamics of n = 60 myocytes, with the findings demonstrating multi-level spatial variation of calcium dynamics across the cell, supporting the dependence of synchrony and characteristics of calcium release on the underlying t-tubule structure.

For the quantification of calcium dynamics, the Fluo-4 channel background for each individual acquisition was determined by estimating the average signal across a rectangular area manually selected from outside of the cell area. The background was subtracted as a constant offset from each Fluo-4 channel image. Cell segmentation was achieved by intensity thresholding using the Otsu method (Otsu, 1979) of the mean image across 3,160 frames from the Fluo-4 channel, and a series of subsequent modifications of the resultant mask. To exclude any unwanted particles or parts of other cells appearing within the imaged FOV, connected components of the thresholded image containing fewer than 30,000 pixels (corresponding to a cell area of ~650 µm 2 ) were removed, and the thresholding step was repeated with that region masked out. Any holes within the masked region were filled and the remaining cell mask was smoothed by convolution with a 5×5-pixel x-y kernel. with the resultant cell mask shown as a white outline in the merged channel image in Supplementary Figure 2d.

Supplementary Figure 2
Pre-processing of the dual channel data illustrated for a single frame from a 12,000-frame acquisition of representative Cell X at 395 fps at the peak of the first transient. (a) Raw single frame with green and magenta dashed rectangles indicating the ROI used for the Fluo-4 and CMO channels respectively. (b) Separated and co-registered spectral channels with subtracted fixed pattern noise. (c) The Fluo-4 (green) and CMO (magenta) spectral channels overlayed with the segmentation mask outlined in white. Intensities are in digital numbers. Scalebar: 10 μm.

Image analysis -identifying t-tubule structure and nuclei
As the cell was immobilized and there were minimal motion artefacts during acquisition, an average across 12,000 frames in the CMO-channel was taken to maximize the SNR (Supplementary Figure 3a). To account for local background spatially varying across the cell, the average image was blurred by a 15 × 15-pixel median filter, with the resultant image (Supplementary Figure 3b) subtracted from the original (Supplementary Figure 3c). The tubule microstructure was extracted from the fluorescence of the CMO membrane stain using an intensity threshold. To compensate for the CMO intensity variation between different cells, the threshold used for CMO tubule segmentation was set to 1% of the membrane intensity determined for that cell. The cell membrane intensity was calculated by taking a one-dimensional maximum intensity projection (MIP) along the vertical axis of the CMO image and calculating the average in the horizontal direction across the central third of the FOV (Supplementary Figure 3d), which allowed exclusion of the peaks in the vertical MIP at the ends of the cells. The intensity threshold was used to produce a binary mask (Supplementary Figure  3e), and connected components with areas less than 20 pixels were excluded, as they were predominantly found to be due to noise artefacts. The tubule mask was multiplied by the cell segmentation mask to exclude any regions outside of the cell and smoothed by a 3 x 3 pixel median filter (Supplementary Figure 3f).
Nuclei were segmented by visual inspection of the stained cell membrane in the CMO-channel averaged the across the acquisition. The nuclear envelope was outlined manually using the drawassisted MATLAB tool, creating a ROI. The selected ROI were combined to create a binary mask, within which connected components were identified and characterized using the bwconncomp and regionprops MATLAB tools.   Table 3 Comparison of average time-to-half-maximum (T50) for different cellular regions: the cell (excluding the nucleus), the nucleus, the central IQR and the outer quarters. The mean, standard deviation (SD) and standard errors (SE = SD/√n) calculated across the medians for each T50 distribution for each of the n = 40 cells with identified nuclei for the first category, and all n = 60 cells for the second category. The Δ column indicates the difference between the T 0 for each region, averaged over all cells. Based on statistical considerations outlined in Section 2.4, the p-values are calculated using the unpaired Mann-Whitney (MW), the paired Wilcoxon signed rank test (WSRT), and one-sample nested t-test with cell-by-cell differences averaged for each heart. For all three statistical analyses (unpaired, paired by cell, and averaged for each heart) the nuclei and central regions of the cell had a significantly (p-values < 0.05 are highlighted in green) longer time-to half than the rest of the cell and the central region respectively.  Table 5 omparison of correlation slopes (ΔT 0/ΔDNT, including 95% confidence intervals (CI) of the linear fit) in tubulated (DNT < 5) and detubulated (DNT > 5) cell areas, calculated as an average across the slopes for each cell. The p-values are calculated using the unpaired Mann-Whitney (MW), the paired Wilcoxon signed rank test (WSRT), and one-sample nested t-test with cellby-cell differences nested by heart. For all three statistical analyses (unpaired, paired by cell, and nested by heart) the T50/DNT slope was not significantly different between the two DNT categories.  Table 6 Comparison of the spark and transient parameters calculated for LV and RV CM the full n = 60 dataset quoted as the mean, standard deviation (SD) and standard errors (SE = SD/√n) calculated from the medians of the cell parameters in each of the two cell populations.

Comparison of left and right ventricle cardiomyocytes LV (n = 32 Cells) RV (n = 28 Cells) Δ (%) p-value t-test
The differences (Δ) are calculated as the difference between corresponding parameters in left and right ventricles, divided by the former value, and converted to a percentage. The p-values are calculated using the unpaired t-test (UT) and the unpaired Mann-Whitney (MW). For the full n = 60 dataset there were no significant differences across any of the spark parameters.
To evaluate the potential systematic bias in imaging time due to un-randomized blinding of the two ventricles, the correlation of each parameter with the relative imaging time point (to the first imaged cell in that batch) was considered. Correlation was assessed by linear fits to the data for cells from each ventricle and considering the coefficient of determination R 2 . The slopes (Δ), vertical axis intercepts (Y0) and R 2 coefficients for the linear fits are summarized in Supplementary Table 7. A weak correlation with the time of imaging was found for some spark and transient parameters (R 2 = 0.12-0.49), with a positive time correlation for spark rate and FDHM in both ventricles and LV respectively, and a negative time correlation for spark amplitude and transient dyssynchrony for LV CM and RV CM respectively.
The scatterplots for parameters with a weak correlation (R 2 > 0.1) present are shown in Supplementary  Figure 7. Two of the selected left ventricle acquisitions had a significantly earlier relative imaging time point compared to the rest (0 and 7 min respectively). To evaluate the influence of these outliers on the time correlation of spark amplitude, the linear correlation was calculated for the dataset excluding those two points (Supplementary Figure 8). As a result, the slope and intercept of the new linear correlation fit remained nearly unchanged (Δ = -0.003, Y0 = 0.93), while the coefficient of determination decreased from R 2 = 0.49 to R 2 = 0.31, indicating a lower negative correlation level, but with a similar time dependence. Supplementary Figure 10 Correlation of spark amplitude with the relative time of imaging for LV (orange) and RV (blue) CM, with two earliest LV CM datapoints excluded. The equation and R-squared (R 2 ) coefficients of determination for the linear fits are provided next to the datapoints. , quoted as the mean, standard deviation (SD) and standard errors (SE = SD/√n) calculated from the medians of the cells in each of the two populations. The differences (Δ) are calculated as the average for each parameter for LV myocytes minus the average for RV myocytes, divided by the former value, and converted to a percentage. Positive Δ indicates larger value for LV myocytes. ll datasets were normally distributed, and each dataset pair passed the F-test for equal variance, and hence p-values were calculated using two sided, two sample, homoscedastic unpaired t-test (UT). For the calculated parameters, only the spark amplitude was significantly different between the two ventricles, with LV CM having 25% larger spark amplitude (p = 0.0016).