Edited by: Jichao Zhao, The University of Auckland, New Zealand
Reviewed by: Zhong Jian, University of California, Davis, United States; Arun V. Holden, University of Leeds, United Kingdom; Michael Alan Colman, University of Leeds, United Kingdom
This article was submitted to Computational Physiology and Medicine, a section of the journal Frontiers in Physiology
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
Calcium signaling plays a pivotal role in cardiomyocytes, coupling electrical excitation to mechanical contraction of the heart. Determining locations of active calcium release sites, and how their recruitment changes in response to stimuli and in disease states is therefore of central interest in cardiac physiology. Current algorithms for detecting release sites from live cell imaging data are however not easily validated against a known “ground truth,” which makes interpretation of the output of such algorithms, in particular the degree of confidence in site detection, a challenging task. Computational models are capable of integrating findings from multiple sources into a consistent, predictive framework. In cellular physiology, such models have the potential to reveal structure and function beyond the temporal and spatial resolution limitations of individual experimental measurements. Here, we create a spatially detailed computational model of calcium release in an eight sarcomere section of a ventricular cardiomyocyte, using electron tomography reconstruction of cardiac ultrastructure and confocal imaging of protein localization. This provides a high-resolution model of calcium diffusion from intracellular stores, which can be used as a platform to simulate confocal fluorescence imaging in the context of known ground truth structures from the higher resolution model. We use this capability to evaluate the performance of a recently proposed method for detecting the functional response of calcium release sites in live cells. Model permutations reveal how calcium release site density and mitochondria acting as diffusion barriers impact the detection performance of the algorithm. We demonstrate that site density has the greatest impact on detection precision and recall, in particular affecting the effective detectable depth of sites in confocal data. Our findings provide guidance on how such detection algorithms may best be applied to experimental data and give insights into limitations when using two-dimensional microscopy images to analyse three-dimensional cellular structures.
Each heartbeat is induced by an efflux of calcium (Ca2+) from the sarcoplasmic reticulum (SR) into the cytosol. This release of Ca2+ through clusters of ryanodine receptors (RyRs) raises the bulk cytosolic Ca2+ concentration from 0.1 μM to ≈1 μM within 30 ms and results in exposure of cross-bridge binding sites on the actin filaments, facilitating cellular contraction (Bers,
Functional imaging of Ca2+ release in cardiomyocytes is possible using fluorescence-labeled confocal microscopy. Unfortunately, these images suffer from low signal-to-noise ratio and their two-dimensional nature collapses the dynamic three-dimensional system. Recently, Tian et al. (
However, questions remain regarding the suitability and performance of this proposed approach. Foremost is the inability to define ground truth locations of active RyR clusters to compare against those detected when applying the method to experimental data. Further complicating detection is that, compared to the near-vacuum between astronomical objects, the diffusive volume between clusters of RyRs in cardiomyocytes is heterogeneous, consisting primarily of myofibrils and mitochondria. This issue is relevant to the local distribution of diffusing Ca2+: while homologs of cell membrane Ca2+ transport channels exist on mitochondria and exhibit a modest buffering effect, Ca2+ flux between the cytosol and intra-mitochondrial space is negligible compared to other cytosolic Ca2+ pathways under normal physiological conditions (Williams et al.,
To address these issues, we develop a spatially detailed computational model of an eight sarcomere section of a cardiomyocyte by extruding an electron tomography image. Within this three dimensional domain, we simulate reaction-diffusion of Ca2+ emanating from RyR clusters. Our model captures these mechanics during the rising phase of the Ca2+ transient- the first 30 ms following the membrane depolarization before contraction begins- upon which CaCLEAN detection operates (Tian et al.,
Finite element model of Ca2+ reaction-diffusion in an eight sarcomere section of a cardiomyocyte.
The results of this model are subsequently used to simulate confocal fluorescence images. These images are then analyzed with CaCLEAN to assess RyR cluster detection performance in the context of known ground truth locations to quantify true positives (hits), false negatives (misses), and false positives (algorithmic artifacts). Our findings indicate that the presence of mitochondria only has a marginal negative impact on detection at a typical experimental imaging resolution. We quantify the impact of inter-cluster spacing on detection performance, as well as how far from the imaging plane clusters are most accurately detected. We estimate the recall and precision of the algorithm as between 69–82%, depending on the density of cluster locations and their distance from the imaging plane. Our analysis therefore serves as a reference for future applications or extensions of CaCLEAN and similar release-site detection algorithms, providing quantitative analysis of performance using a physics-based modeling framework with known ground truths.
A spatially detailed finite element (FE) computational model of an eight sarcomere section of a cardiomyocyte was constructed to evaluate RyR cluster detection. The algorithms used to generate the model (Rajagopal et al.,
CASE 1 (high cluster density, no mitochondria): A case with high cluster density (
CASE 2 (high cluster density, with mitochondria): Mitochondrial regions (
CASE 3 (low cluster density, no mitochondria): A case with a relatively low cluster density (
CASE 4 (low cluster density, with mitochondria): Mitochondrial obstacles were included as in case 2; RyR cluster distributions were defined as in case 3.
In the three-dimensional FE model, the fluorescence-bound Ca2+ (FCa) field was calculated based on the time-dependent reaction-diffusion emanating from the surrounding cluster sources (see
Simulation of confocal fluorescence microscopy images from FE model results.
The resulting time-dependent, simulated confocal fluorescence microscopy images at each slice (
Example of CaCLEAN detection and classification against modeled locations.
Illustration of the admissible window parameter. The admissible window defines the distance tolerance from the simulated imaging plane for RyR cluster centers to be considered ground true positives.
For each simulated imaging plane, modeled RyR clusters were considered TP (ground truth) if their centers were within a distance threshold from the imaging plane we referred to as the “admissible window” (see
Clusters located very near the z-depth
We quantified algorithmic performance using well-established measures to capture the impact of false positives and false negatives on performance: recall, precision, and f1-score (Nisbet et al.,
Recall and precision of CaCLEAN applied to simulated microscopy data. Recall, precision and f1-score evaluated based on classification results. Shaded regions indicate one standard deviation of recall and precision values. The f1-score shown is evaluated based on the mean values for precision and recall (solid red and blue lines). In both cases shown, regions representing mitochondria act as barriers to diffusion.
Ca2+ release site detection performance results at maximum f1-score in each model.
Low RyR density, |
0.83 | 0.83 ± 0.06 | 0.83 ± 0.08 | 630 |
Low RyR density, |
0.82 | 0.82 ± 0.06 | 0.82 ± 0.07 | 620 |
High RyR density, |
0.70 | 0.67 ± 0.08 | 0.74 ± 0.10 | 290 |
High RyR density, |
0.69 | 0.68 ± 0.07 | 0.70 ± 0.09 | 290 |
Whereas there may be applications where optimizing for recall or precision may be more appropriate, we used the f1-score as a single general measure, equally weighting the influence of misses and false positives in interpretation of performance. Peak values of the f1-score in the cases including mitochondria barriers ranged between 0.69 at an admissible window of 620 nm in the high cluster density case and 0.82 at an admissible window of 290 nm in the low cluster density case (see
Ca2+ release site detection performance results at the precision-recall break-even point, where Precision ≈ Recall and the number of false positives ≈ the number of false negatives.
Low RyR density, |
0.83 | 0.83 ± 0.06 | 0.83 ± 0.08 | 630 |
Low RyR density, |
0.82 | 0.82 ± 0.06 | 0.82 ± 0.07 | 620 |
High RyR density, |
0.69 | 0.69 ± 0.08 | 0.69 ± 0.12 | 270 |
High RyR density, |
0.69 | 0.69 ± 0.08 | 0.68 ± 0.09 | 280 |
The effect of mitochondria acting as barriers to diffusion had a slight consistently negative impact on the evaluated performance when compared against cases where these regions were modeled as diffusing homogeneously with the cytosol. This effect was most evident in the precision-recall curves in
Precision-recall curve for CaCLEAN in four models. Recall and precision mean values from results shown in
Finally, we evaluated the fraction of clusters detected by CaCLEAN as a function of z-distance from the focal imaging plane to identify how detection of individual clusters decayed with increasing distance from the imaging plane (see
“Differential recall”: axial dependence of detection of true positive sites. Each curve shows the fraction of detected sites as a function of the distance z from the nominal focal plane. Differential recall at distance z is the fraction of all detected to modeled sites within 5 nm of that z depth i.e., within a 10 nm band centered about the distance z above and below the focal plane. A value of 1 is equivalent to the detection of all sites at a given z-depth. The various curves are calculated for different models, having a high or low density of sites and either using homogeneous diffusion (“no mito”) or obstacles to diffusion wherever mitochondria are (“with mito”). Data points shown are the result of applying a smoothing filter (see
We developed a computational model of the complex environment of Ca2+ diffusing into the intracellular space of a cardiomyocyte. Processing the reaction-diffusion model results to simulate confocal fluorescence microscopy data allowed for quantitative assessment of the performance of detection of Ca2+ release sites against known ground truth values in the context of realistic cellular physics. Statistical classification identified true positives, false positives, and false negatives
A key variable in the performance analysis was the definition of which modeled clusters were considered ground truths in the statistical classification at each imaging plane. We introduced the “admissible window” variable for admitting clusters as ground truth values based on the through-imaging-plane distance of cluster centers, as illustrated in
The admissible window is also of practical interest to those seeking to use CaCLEAN (or another detection algorithm) on their experimental data. This parameter can be used as an indicator of the maximum relevant depth of clusters detected by the algorithm in the image. As evident in
For example, if a user is confident that clusters in a sample are likely at least 1 μm apart and is satisfied with 80% precision and recall, the maximal relevant depth for clusters diffusing into the image would be ≈620 nm based on
Analysis of performance in terms of the admissible window provides a basis for assessing “cumulative recall,” i.e., the detection of all clusters within a given tolerance of the focal plane. We also evaluate “differential recall” in
The shape of the point spread function (PSF, see
In addition to distance from the imaging plane, further factors complicate whether the signal from a given cluster will reach the imaging plane and whether it is detectable in both our model and actual experimental data. Individual cluster firing time and strength variability biases detection of early and stronger events. Some false positives are produced as algorithmic artifacts. Proximity to other clusters can cause signals to merge or cover each-other- especially with increased cluster density. For instance, two separate clusters the same z-distance above and below the imaging plane but at the same x-y location in the imaging plane can only register as a single site in the imaged space. This identifies an inherent drawback resulting from using two-dimensional images as a basis for describing a three-dimensional system.
RyR cluster distribution spacing varies across species, with nearest-neighbor spacings reported as 0.66 ± 0.06 μm in rat and 0.78 ± 0.07 μm in human (Soeller et al.,
We investigated the impact of cluster spacing on detection performance by analysing two cluster distribution types. In the “high cluster density” cases, cluster locations at each z-disk were defined based on statistical analysis of cluster distributions from immuno-labeled confocal images and applied to admissible locations bordering mitochondria and myofibrils segmented from electron tomography data (Rajagopal et al.,
Our results suggest that CaCLEAN correctly detects the majority of Ca2+ release sites, with approximately one miss and one false positive out of every four or five sites (depending on site density). In our analysis, site spacing had a significant impact on both detection performance and the admissible window size associated with optimal performance. At the peak f1-score for cases including mitochondria acting as diffusion barriers, detection recall and precision were both 0.82 at 610 nm in the low cluster density case vs. 0.68 and 0.70 at 290 nm in the high cluster density case. In contrast, for the high cluster density results evaluated at the low cluster density peak f1-score admissible window of 610 nm, recall ≈0.43 thus indicating more false negative misses than true positive detection hits.
For useful detection performance, it is therefore important to consider the likely density of the events being detected in order to determine how far from the imaging plane such events are likely located. Those interested in using CaCLEAN to reconstruct three-dimensional maps of clusters should also be aware of how far from the imaging plane detected clusters are likely to be when choosing a slicing depth for reconstruction. In this case, the choice of slicing depth should be approximately twice the optimum performance admissible window. For example, in a case with high cluster density where precision and recall are equally important, a spacing of 580 nm would be recommended based on the values reported in
The presence of heterogeneous diffusion (in the form of mitochondrial obstacles) in the investigated models was found to have a marginal but consistent negative impact on both recall and precision in CaCLEAN. While mitochondrial diffusion barriers introduced some additional false positive events, localized increases from Ca2+ reflecting against these obstacles did not significantly impair overall detection performance. Resolution limits may actually help to mitigate this factor, as such local increases in [Ca2+] at mitochondria-cytosol boundaries may be offset by lower values within the mitochondrial regions during interpolation into images. We hypothesize that structural heterogeneity due to the presence of subcellular structures (e.g., mitochondria, nucleus, transverse tubules, z-disks, etc.) is therefore unlikely to have a major impact on detection of point sources using CaCLEAN or related algorithms under imaging resolution conditions similar to those simulated here. Subsequent models similarly focused on Ca2+ release site detection performance might therefore benefit from significantly simplifying the computational domain by avoiding distinguishing mitochondrial regions.
We evaluated the performance of a published algorithm for extracting calcium release site distributions from live confocal images using ground truth calcium release data that was generated using a spatially detailed model of a rat left ventricular cardiomyocyte. Several simplifying assumptions were made to create and analyse the model with this application in mind. Here we discuss these assumptions and also highlight ways in which the model could be extended and applied beyond the current study.
The modeled eight sarcomere domain was created by extruding a two-dimensional image from an electron tomography image stack. It therefore does not capture subtle structural variations along the longitudinal direction and assumes such changes are minimal over the relatively short (8 sarcomere, 16 μm) section modeled. Furthermore, sarcomere length in our model was fixed to 2 μm and we did not examine the effect of sarcomere length on detection results, as may occur in disease states (e.g., sarcomere lengths have been reported to shorten from 1.84 to 1.79 μm in a mouse model of lipotoxic diabetic cardiomyopathy (Flagg et al.,
Our model simulated the rising phase (first 30 ms) of the calcium transient that the detection algorithm (CaCLEAN) is designed to operate over (Tian et al.,
We chose to use a fixed Ca2+ release approach for the present study, since our goal was the evaluation of release site detection in simulated images rather than recapitulating inter-cluster feedback dynamics. This approach allowed greater control over the number of simulated releasing clusters in the model permutations investigated. Cases with mitochondrial diffusion barriers reduced the cytosolic volume by removing these regions, resulting in higher [Ca2+] throughout the cytosol. High cluster density permutations also increased global cytosolic [Ca2+]. These cases would therefore result in additional cluster activation from a CICR model compared to the low cluster density and homogeneously diffusing cases. We therefore found a fixed Ca2+ release model more appropriate for the control of this study given our focus. However, we acknowledge these mechanisms are important considerations for more general models seeking to explore inter-cluster Ca2+ signaling.
Ca2+ release in our model focused on sparks produced by active clusters of RyRs, the primary propagator of Ca2+ signals in healthy cardiomyocytes (Cheng et al.,
Sarco/endoplasmic reticulum Ca2+-ATPase (SERCA) re-uptakes Ca2+ from the cytosol into intracellular stores and has been shown to colocalize with RyRs at z-disks (Drago et al.,
Our model was based on experimental data from healthy male adult Wistar rat ventricular myocytes. We hypothesize our findings may extend to other cell types where RyRs are equally well organized (e.g., Purkinje, Hirose et al.,
The results presented in this work (particularly
The computational model results include the Ca2+ and fluorescence-bound Ca2+ data discussed here, along with solved fields for unbound fluorophore, calmodulin, Ca2+-bound calmodulin, ATP, Ca2+-bound ATP, and troponin C. These results may be useful to others interested in a high-resolution representation of the mechanics of these fields during the rising phase of the Ca2+ signal transient in a rat ventricular cardiomyocyte. The capability to simulate imaging data from the model results also means it could serve as a testing and validation platform for other analysis tools operating on confocal fluorescence imaging of Ca2+ release in cardiomyocytes that would benefit from using the high resolution continuum model as reference. For instance, this could be used to further improve CaCLEAN or other approaches seeking to identify RyR clusters by serving as a training set for improved detection or segmentation. The code and supporting datasets for the model, confocal data simulator, and detection performance analysis are therefore freely available for subsequent research (see Data Availability below). Finally, this study highlights how computational methods may be used to combine diverse experimental data into a consistent physical framework to establish ground truth values that may not otherwise be experimentally available.
The reaction-diffusion finite element (FE) model here builds on a previous FE model of a half-sarcomere (Rajagopal et al.,
From this half-sarcomere image stack, the central slice representing the level of the z-disk was extracted and regions representing myofibrils and mitochondria were manually segmented. This two-dimensional slice was then extruded 16 μm to create a three dimensional volume.
Two configurations were considered to assess the impact of mitochondria: one in which the interior of the modeled domain is a homogeneous material continuum and one in which the regions representing mitochondria were subtracted from the myofibrillar and cytoplasmic domain. The latter case was based on the assumption that the calcium buffering activity of mitochondria is negligible.
Ryanodine receptor (RyR) locations were defined algorithmically, using a spatial statistics method based on nearest neighbor distances of experimentally-derived RyR locations (Rajagopal et al.,
To evaluate the influence of spacing on the RyR cluster detection performance of CaCLEAN, two sets of constraints on the RyR distributions were considered: (1) RyR clusters were assigned with locations based directly on statistical analysis of experimental data and with the total number of RyR clusters
Spherical regions 100 nm in radius were defined around nodes nearest to the determined location for each RyR cluster. Nodes lying within this sphere were prescribed density amplitudes exponentially decreasing as a function of the square of their radial position from the central node. RyR cluster release times were sampled from an exponential distribution with a characteristic decay constant of 6.7 ms, following the findings of Wang et al. (
In our previous model (Rajagopal et al.,
We used an effective diffusion coefficient for diffusion of Ca2+ released in the cytosol. This diffusion is coupled with reaction terms that mediate Ca2+ binding and unbinding to buffers within the cytosol. Further details on the reaction-diffusion equations for the buffers modeled and the ordinary differential equation model describing release of Ca2+ from RyR clusters have been previously reported (Rajagopal et al.,
To simulate confocal fluorescence results, the irregularly distributed node-based fluorescence-bound Ca2+ field from the FE model was first interpolated onto a regular grid at a resolution of 53.75 nm in each direction. The data were temporally sampled at 5 ms intervals over the simulated time period of 30 ms, producing simulated imaging data for 7 timesteps. Nodal positions within mitochondria in the “with mitochondria” permutations were ascribed the initial and background value FCa = 2.08. Discrete natural neighbor (Sibson) interpolation (Park et al.,
A point spread function (PSF) was generated as a normalized function of a multivariate Gaussian distribution applied in three dimensions (see
The interpolated model imaging data was then convolved with the PSF in three dimensions using the SciPy convolve algorithm from the signal processing module. The resulting grid was then downsampled to a pixel resolution of 215 nm, following the resolution of the original CaCLEAN paper (Tian et al.,
The CaCLEAN algorithm was obtained from the author's GitHub repository:
A statistical classification approach was used to assess the performance of RyR cluster detection. Modeled cluster centers within the admissible window were considered the actual/ground truth class: TP (ground truth). Detection results were considered the predicted class. Detected RyR cluster sites were defined as determined by the CaCLEAN CRUProps function, which segments cluster regions using Matlab's built-in watershed algorithm and identifies centroids of segmented regions.
For each modeled cluster location, a TP (detected) classification was assigned if a TP (ground truth) cluster center lied within an available segmented CaCLEAN cluster region (or within a 1 pixel tolerance). When more than one TP (ground truth) fell within a CaCLEAN-detected cluster region, the detected cluster with the nearest centroid to the TP (ground truth) location was marked TP (detected). After classification as TP (detected), the associated CaCLEAN-detected site would be removed from the list of available matches. After iterating through the TP (ground truth) clusters, remaining TP (ground truth) unmatched with detected clusters were classified as false negatives (FN). Remaining CaCLEAN-detected sites unmatched with TP (ground truth) were classified as false positives (FP). Note that in this case
Three statistical binary classification performance measures were considered: recall, precision and f1-score (Nisbet et al.,
Recall therefore gives the fraction of actual modeled clusters within an admissible window that were correctly detected by CaCLEAN. Precision (also known as positive predictive value) was defined such that
Precision therefore identifies the fraction of the detected clusters within an admissible window that were correct (not false positives). Another useful parameter indicating the combined effect of both precision and recall is the f1-score (also known as f-measure or f-score), which measures the harmonic mean of these two variables i.e.,
This provides a single combined performance metric that equally weights the impact of false positive and false negative events. In all three performance metrics, higher values indicate better performance.
The above definition of recall may be considered “cumulative recall” in our application, identifying the detection fraction of all clusters within a given admissible window. To determine the detection fraction of clusters at a given z distance from the simulated imaging plane, we also defined an alternative “differential recall” such that
This measured the fraction of modeled clusters detected by CaCLEAN in 10 nm spaced bands above and below the imaging plane. Only bands with at least one TP (ground truth) were considered. Mean values for this detection fraction were acquired over the 22 simulated imaging planes. A Savitzky-Golay filter (polynomial order 3, frame length 21) was then applied to smooth the results as shown in
The code used to solve the FE model is available in the repository:
DL developed, analyzed, and curated the computational model. DL, CS, and VR contributed to the initial conceptualization of this study. CS, EC, and VR supervised the project. EC and VR directed administration of the project and funding acquisition. DL, AT, HR, CS, EC, and VR contributed to the development of the methodology and preparation of original and final drafts.
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest. The reviewer MAC declared a past collaboration with one of the authors CS.
This manuscript has been released as a pre-print at BioRxiv (Ladd et al.,
The Supplementary Material for this article can be found online at:
1Where the two-dimensional simulated image plane represented x-y axes.
2True negatives were not considered in this application, as they would represent the set of all remaining locations that were neither modeled nor detected.
3Immuno-labeled confocal data populations may still slightly underestimate cluster density owing to diffraction-limited resolution.
4See Figure 3Cb in Tian et al. (