Hearts from 22 C57/BL6 wildtype mice (aged 12–20 weeks) of both sexes in equal number were used for this study. The preparation of atria including the intact SAN region was carried out as previously described (Torrente et al., 2015). Briefly, mice were heparinized and anesthetized with ketamine (100 mg/kg) and xylazine (10 mg/kg) and the hearts were quickly removed. The atria including the intact SAN region were dissected from the ventricles and fixed with needles on a silicone ground of an experimental chamber. For this study, an extracted SAN tissue of one mouse represented a single experiment. Therefore, the number of experiments corresponds to the number of mice.

The experimental procedure and number of used animals were approved by the ethics committee of the Federal Ministry of Science, Research and Economy of the Republic of Austria (BMWFW-66.010/0101-WF/V/3b/2016). The experiments were conducted according to the Directive of the European Parliament and of the Council of September 22, 2010 (2010/63/EU).

The experimental chamber containing the intact SAN tissue was mounted on the stage of an upright microscope (Olympus, BX51W1, 4x objective, light source TH4-200) and the tissue was superfused with oxygenated standard external solution (containing in mM: NaCl 137, KCl 5.4, CaCl₂ 1.8, MgCl₂ 1.1, NaHCO₃ 2.2, NaH₂PO₄ 0.4, HEPES/Na⁺ 10, D(+)-glucose 5.6, pH 7.4 adjusted with NaOH) which was kept at a constant temperature of 23°C. Recordings were started 20 min after the onset of superfusion in order to allow the tissue to establish and maintain a stable beating rhythm. Close to the primary pacemaking site of the tissue, a small image region of interest ROI showing distinct contractions was selected for recording. After recording of the first video (Con), acetylcholine (ACh, 3 μmol/L) was added to the perfusion solution and after 5 min superfusion time the second video was recorded. A number of nine tissue samples yielded 18 videos.

Tissue samples under investigation showed a beat to beat interval of about 500 ms (~2 Hz). In order to measure such intervals, it is necessary to sample the temporal course of the beating with enough data samples per second. The Nyquist-Shannon sampling theorem with a sample rate that is the twofold of the highest frequency in the signal is not applicable, because the content of harmonic frequencies are not of interest. More important is, that the sample rate determines a minimal jitter between consecutive time stamps. This jitter must be small, because it influences the nonlinear analysis including variability measures. For an accuracy of e.g., 1%, a number of 100 data samples is needed between two succeeding beats. Particularly, this would yield a sample rate of 200 Hz. Please note that video acquisitions using standard frame rates of 30 Hz would yield an accuracy of only ~6.67% for a beating rate of ~2 Hz. We decided to set the accuracy to 0.2% and consequently a sample rate of 1 kHz was used.

Video recordings of beating tissue samples were taken by using a high-speed camera system (MotionBlitz, GMCLTR1.3CL-SSL, Mikrotron, Germany) and a video camera (EoSens CL, MC1362, Mikrotron). This system implemented a hardware recording unit and therefore, avoided erroneous jitter effects of software trigger events such as e.g., USB camera solutions do usually show. The resolution of the camera was set to 1,280 × 1,024 pixels. A ROI with a pixel size of 160 × 160 was selected to get maximal gray level changes during beating of the tissue. Figure 1A shows a sample image.

The size of the ROI was empirically optimized by inspections of the temporal signals gained. Larger ROIs yielded too large and inconvenient video files and smaller ROIs yielded too low signal to noise ratios. Regions with thick tissue layers (trabeculae, crista terminalis) moving into the ROI turned out to yield the highest signal to noise ratios. With a sample rate of 1,000 fps and a recording duration of 5 min, a number of 300,000 single uncompressed images (each of them with 160 × 160 = 25,600 pixels) where taken and stored on hard disk in an avi container format. One avi file needed 21 GB of memory. The Mikrotron system saved the individual images in RGB format although the used camera was a gray value camera. Gray value cameras usually give higher signal to noise ratios than color cameras which is important for high-speed acquisitions with very small exposure times. Thus, images were converted in a first step to 8 bit, lowering the memory demand to 7, 6 GB per video. The whole measuring setup was tested against electrical and optical inferences coming from ambient light sources such as the laboratory light or the microscope light source itself. Fourier analyses of videos capturing a static scene revealed no residual frequency components of power supply frequencies or other additional noise components.

For comparison of video and electrical measurements we carried out 13 separate control experiments in order to simultaneously record video as well as extracellular electrical signals. After positioning of the microscope's objective and choosing a ROI, a cardiac-near-field CNF electrode (Hofer et al., 2006) was placed close to the ROI (Figure 1A). This ensured that the electrical and optical measurement sites corresponded in the spatial domain. For evaluation of beat to beat intervals, only one of the four CNF channels was used. Electrical signals were amplified (gain 100), anti aliasing lowpass filtered (4th order Bessel, cutoff frequency 20 kHz) and recorded with custom software (LabVIEW, National Instruments, Austin, Texas) at a sampling rate of 100 kHz (NI USB-6210, National Instruments, Austin, Texas). Signals were digitally filtered (Butterworth lowpass, 4th order, cutoff frequency 1.5 kHz and Butterworth highpass, 4th order, cutoff frequency 1.5 Hz). A sample electrical recording can be seen in Figure 1B. Subsequent time stamps of individual beats were computed by setting a threshold well above the noise level to the decreasing slopes of the signal (denoted by red plus signs in Figure 1B). For signals with lower signal to noise ratio the same threshold criterion was used in the first temporal derivative of the signal where the steep downslope during electrical activation was more pronounced.

Time signals of the beating tissues were reconstructed by computing the average gray value G¯i of each image i of a video which can be seen in Figure 1C.

with g_{i, p} the gray value in the range [0, 255] of pixel p in the image i, i = 1, 2, 3, …, N_I, N_I the number images in a video (N_I = 300, 000), and p = 1, 2, 3, …, N_P, N_P the number of pixels in an image (N_P = 25, 600).

This yielded a temporal time signal comprising 300,000 data points with a data compression of (25,600 to 1). The algorithm for finding the time stamps of contractions was designed around finding subsequent local minima (denoted by red plus signs in Figure 1C). The original signal was slightly smoothed by applying a moving average filter with 25 data points to improve the shape of the minima and to ensure that minima are right between the adjacent declining and rising slopes. A threshold was set between the doubled noise value of the baseline and the smallest minimum in the video. Further on, only values smaller than that threshold were investigated. Then, a minimum was computed by simply looking for the smallest value between two threshold points. Very rarely a minimum consisted of two neighboring points with exactly the same value in which case the second value was taken as the minimum. Finally, a beat to beat interval BBI was defined as the temporal interval between two succeeding minima (beats).

with t_b, the absolute time of the beat b (minimum) in a video, b = 2, 3, …, N_B, and N_B the number of beats in a video. This algorithm (named PointFinder) was implemented in the software IQM (Kainz et al., 2015) and is available from the authors or from the IQM project page (https://sourceforge.net/projects/iqm/). A sample result of BBIs is shown in Figure 2A.

The number of data points in such a graph was equal to the number of beats N_B − 1 during the recording interval of 5 min and therefore, was not constant from video to video. The sample graph in Figure 2A consists of actually 660 data points because the specific tissue sample contracted 661 times in 5 min. Roughly, this reflects another data compression of 450 to 1. Overall the method has a compression rate of about 11.6 million to 1.

On further inspection of Figure 1C, it is obvious that not only interbeat intervals can be computed from these local minima. Additionally, the height of a minimum measured from the baseline, which is the horizontal line containing the points of relaxed tissue only, can be computed. Such a height reflects the mechanical contraction strength CS of the tissue. Stronger contraction of cells in the tissue yielded a higher light absorption and hence darker images in the video and consequently lower minima in the average graphs such as in Figure 1C. We created correlation plots of subsequent height changes vs. subsequent interbeat intervals (actual plots are described in section Correlation of Contraction Strengths and Beat to Beat Intervals). A correlation between these two variables was given in most cases, which is in accordance with previously published data (Torres and Janssen, 2011). Nevertheless, correlations were not perfect and some experiments showed only weak or negative correlations. Consequently, we decided to additionally evaluate variations of contraction strengths.

First, in order to obtain the CS, it was necessary to determine the baseline, although it drifted during the time course of 5 min. To avoid drifting and offset errors we computed a separate baseline value for each minimum. The actual baseline value was computed by the median of all points in between the actual minimum and the preceding minimum and thus we obtained moving baseline values (for every minimum a separate baseline value, some sample baseline values are depicted graphically in Figure 1C with red points). This ensured that baseline drifts, unavoidable during a recording time of 5 min, did not contribute. Consequently, CS was computed by

with CS_b the contraction strength of the beat b, G¯b  the average gray value of the image detected as the “beat image,” baseline_b the corresponding baseline value, b = 1, 2, 3, …, N_B, and N_B the number of beats in a video. A sample signal can be seen in Figure 2B.

Although the signals seemed to be reliable, it turned out that this moving median algorithm produced some variation errors due to residual noise components of the baseline. Additionally, but only occasionally, for one mouse treated with acetylcholine (#5 in Figure 6), subsequent contractions showed some consecutively fast repeating bursts instead of full contractions. For such short burst intervals, the median algorithm yielded erroneous baseline values leading to too small CS values. A graphical representation of such errors is shown in Figure 1D. An alternative algorithm for finding the baseline may be feasible, but we present a convenient way that does not need baseline detection at all.

In this approach we eliminated the baseline (offset), since we actually were interested in variabilities of these temporal signals and not in absolute values. Discrete differentiation of the time signal (average gray values according to Equation (1) and exemplarily depicted in Figure 1C) revealed differences of contraction strengths ΔCSs and eliminated the hassle of finding baseline points.

with b = 2,3,…, N_B.

Subsequent discrete integration generated back the changing content of the contraction strength signal but without the baseline and was termed relative contraction strength rCS.

with b = 2, 3, …, N_B and rCS₁ = 0.

Integration usually gives the anti-derivative plus an unknown constant, which was in our case the baseline (offset). Discrete differentiation followed by integration was actually carried out with software IQM (Kainz et al., 2015) using the mathematics feature for one dimensional signals.

Sample graphs of ΔCS and rCS can be seen in Figures 2C,D. BBI, CS and rCS data series for each experiment (Con and ACh treated) are provided as “Data Sheet 1” csv file in the supplement.

Sample entropy SampEn and Higuchi dimension D_H are two well-known and successfully applied nonlinear descriptors for time signal variations (Higuchi, 1988; Richman and Moorman, 2000). Approximate entropy ApEn is also widely used, but is not suitable for this particular study because the number of beats changed from video to video. SampEn is proportional to the conditional probability that a sequence which is similar for m points remains similar for m+1 points. A tolerance distance r is defined so that repetitions must not be exact. Usually, r is defined as a multiple of the standard deviation SD of the signal and therefore, SampEn is a scale invariant measure (Richman and Moorman, 2000). Self matches are not included.

The discrete time signals {x(1), x(2), …, x(N_B)} (x stands for values from CS or rCS signals) with length N_B were taken and (N_B – m + 1) sequences were created:

The parameter m was set to two (m = 2). Distances in between these data series were computed using the maximum metric:

The normalized sums of distances smaller than the tolerance distance r = 0.15SD were computed for each i, j with 1 ≤ i, j ≤ N_B − m + 1 and i ≠ j:

The normalized number of sums can be computed using

Finally, SampEn was computed with

Higuchi proposed a method to compute the fractal dimension of a signal by using sums of differences with varying inter data point intervals (delays) (Higuchi, 1988). The Higuchi dimension is frequently applied in contemporary neurophysiology and neuropathology (Kesić and Spasić, 2016) and is well known for its accuracy, speed and robustness including high linearities of the double log plots. Phase space reconstructions are not involved and therefore, the number of data points available can be restricted. Initial data points are set to m = 1,2,…, k with a delay interval k = 1,2,…, 30. Following data point series were constructed:

The lengths L_m(k) of these series, depending on the initial data points m and k were computed according to:

The symbol ⌊ ⌋ stands for the floor function. For each k, the mean length was determined by

Finally, a double logarithmic plot of L(k) vs. k was constructed and the slope of a linear regression was used to compute D_H. Values of k above 30 were not used because they introduced noticeable deviations of data points from the linear regression. Signals were processed as they were recorded without editing. Algorithms were implemented in the Software IQM (Kainz et al., 2015) and are available from the authors or from the IQM project page.

Statistics was computed with public domain software R, version 3.3.3 and RStudio software version 1.0.136 (RStudio, 2016; R Core Team, 2017). Due to small sample sizes, differences of paired samples were statistically analyzed with a two-sided median test using the R function sintv2 (Wilcox and Rousselet, 2018). This method performs very well in terms of controlling the probability of a Type I error (Wilcox, 2016). Data acquisition via video and electrode setup did not start synchronously (lag of 1–3 beats). Cross correlation (unbiased estimate, MATLAB® R2017b) was used to remove this start dependent asynchronism between the optical and electrical signal. Correlation of CSs with BBIs was computed using Spearman's rank correlation coefficient r_s. Coefficients of determination R² were computed for double log plots to estimate Higuchi dimensions.

Surrogate analysis was performed to provide further evidence of long-range nonlinear correlations in the optically measured signals and to demask possible white noise components indicated by some relatively high SampEn and D_H values. Each optically recorded signal was shuffled 50 times using IQM (Kainz et al., 2015). A total number of 5,400 surrogate data series were constructed considering nine SAN tissues, two nonlinear measures (SampEn, D_H), two treatments (Con, ACh), three signal types (BBI, CS, rCS), and 50x shuffling. Following evaluation types were carried out:

SurrEval-1: Each individual experimentally gained value was tested against the normally distributed surrogate values applying a two-sided one sample Student's t-test.

SurrEval-2: The experimentally gained values were tested against the respective means of the shuffled signals by a two-sided median test using the R function sintv2.

SurrEval-3: The respective means of shuffled control against means of shuffled ACh signals were tested by a two-sided median test using the R function sintv2.

Double log plot linear regressions for estimating the Higuchi dimensions revealed very high coefficients of determination R² within a range of [0.959–0.999]. Sample linear regressions can be seen in Figures 3A,B for control and ACh treated cases. With this high linearity, the application of the fractal concept seems to be very appropriate and robust. Nevertheless, we found a marginally lower R² for some signals (4 out of 18) which was visible as a slight wobbling of data points around the straight line (see a sample graph in Figure 3B, blue dots). This occurred for control as well as ACh treated cases. The reason was that the specific signals showed subsequent alternating values which could be interpreted as binary oscillations or negative correlations. Fractal dimensions of period doubling signals cannot be directly calculated, but it is known that periodic components in time series yield distinct and periodic differences in the double log plot (Galvez Coyt et al., 2013).

Optically (video) recorded values of BBIs correlated very well with electrically recorded values. Representative BBI value pairs showing just negligible differences can be seen in Figure 4A. A sample regression plot of a whole 5 min recording can be seen in Figure 4B. The slopes for all 13 experiments were in the range of [0.977–1.067] and the coefficients of determination R² were in the range of [0.984–0.999].

Additionally, we computed sample entropies and Higuchi dimensions for these 13 control samples. Results can be seen in net plots for each individual experiment in Figure 5A. Traditional box plots including the data values can be seen in Figure 5B.

SampEn values are very close and statistically not different (p = 0.46, df = 12). For D_H values a p-value of 0.05 (df = 12) indicates a possible effect. Since the corresponding median difference was very small (only in the third decimal place), we show scatterplots of data point pairs in Figure 5C. The minimal deviations from the straight line (no effect) suggest no practical relevance.

Scatterplots of CSs vs. BBIs in milliseconds are shown in Figure 6. Control tissue (red) and ACh treated tissue (blue) showed mostly positive correlations. A few correlations are weak and/or negative.

Additionally, scatterplots and correlations of relative contraction strengths rCSs vs. BBIs were computed. Actual plots are not shown, because the correlations were quite similar compared to Figure 6.

The median beat to beat interval of the control group was 512 ms and increased to 614 ms after applying ACh. Net and box plots can be seen in Figure 7 (left column). The significant increase in the beating rate of ~20% (3 μM ACh) is in line with previously published data (Glukhov et al., 2010).

The median beat to beat contraction strength of the control group was 20.85 and decreased to 13.24 after treatment with ACh. Net- and box plots can be seen in Figure 7 (right column). The significant decrease of CS of ~37% is in accordance with previously published data (Kitazawa et al., 2009).

Figure 8 depicts net- and box plots of nonlinear measures of BBIs. Median values of SampEn decrease from 1.58 (Con) to 0.92 (ACh) and are not significantly different (p = 0.11, df = 8, Figure 8B, left column). This is also the case for D_H (p = 0.97, df = 8, Figure 8B, right column) which decreased from 1.98 to 1.97.

SampEn surrogate evaluation according to SurrEval-1 revealed that all experimentally gained values (BBIs for Con and ACh) were significantly lower than the shuffled ones with p < 0.001, df = 49. This agrees with SurrEval-2 showing that all experimental values were also significantly lower than the means of the corresponding shuffled ones with p < 0.001, df = 8, see Table 1. Furthermore, according to SurrEval-3 the means of SampEn values for shuffled signals showed no indication for statistical significance between Con and ACh, p = 0.82, df = 8, see Figure 8C.

SurrEval-1 for D_H yielded similar results with p < 0.001, df = 49, except for some experimental values >2 in the third decimal (two cases for Con and three cases for ACh) where obviously the shuffled values could not further increase. This is well reflected by SurrEval-2 with a borderline p = 0.05, df = 8 for Con and a non-significant p-value for ACh (Table 1). As expected, SurrEval-3 reveals no significant difference between the groups Con and ACh, p = 0.77, df = 8, see Figure 8C.

Figure 9 shows net and box plots of nonlinear measures of the proposed CSs. Median values for SampEn show a significant difference between Con and ACh (p < 0.001, df = 8), namely a decrease from 1.60 (Con) to 0.72 (ACh) which can been seen in Figures 9A,B, left column. Median values for D_H tend to decrease slightly (Figure 9B, right column) as indicated by a borderline p-value of p = 0.05 (df = 8). In detail median values decrease from 1.97 (Con) to 1.91 (ACh).

For SampEn, all three CS surrogate evaluations yielded consistent results regarding nonlinear patterns in the measured signals, since shuffled values were always statistically higher (SurrEval-1: p < 0.001, df = 49, and SurrEval-2: Table 1) and the difference between Con and ACh vanished compared to the experimental case (SurrEval-3: Figure 9C). Now, SurrEval-1 for D_H yielded no exception with p < 0.001, df = 49 and all three surrogate evaluations are again consistent (SurrEval-2: Table 1 and SurrEval-3: Figure 9C).

Finally, Figure 10 depicts net and box plots of nonlinear measures of the proposed rCSs. For SampEn, the decrease and the significance of rCS is similar to CS with median values from 1.76 (Con) to 0.49 (ACh) and with p < 0.001, df = 8 (Figures 10A,B, left column). The decrease of D_H values is slightly more pronounced for rCS compared to CS with median values from 1.98 (Con) to 1.91 (ACh) and now exceeds the 95% statistical significance level with p = 0.03, df = 8 (Figure 10B, right column).

Surrogate analyses for rCs are again consistent for both, SampEn and D_H, indicating nonlinear long-range correlations. Shuffled values are always statistically higher (SurrEval-1: p < 0.001, df = 49, and SurrEval-2: Table 1) and the difference between Con and ACh vanishes compared to the experimental case (SurrEval-3: Figure 10C).

MM-R was employed by KML vision OG, Graz, Austria. KML vision OG did not sponsor this study. The other 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 handling Editor declared a past co-authorship with one of the authors HA.