Sinoatrial Beat to Beat Variability Assessed by Contraction Strength in Addition to the Interbeat Interval

Beat to beat variability of cardiac tissue or isolated cells is frequently investigated by determining time intervals from electrode measurements in order to compute scale dependent or scale independent parameters. In this study, we utilize high-speed video camera recordings to investigate the variability of intervals as well as mechanical contraction strengths and relative contraction strengths with nonlinear analyses. Additionally, the video setup allowed us simultaneous electrode registrations of extracellular potentials. Sinoatrial node tissue under control and acetylcholine treated conditions was used to perform variability analyses by computing sample entropies and Higuchi dimensions. Beat to beat interval variabilities measured by the two recording techniques correlated very well, and therefore, validated the video analyses for this purpose. Acetylcholine treatment induced a reduction of beating rate and contraction strength, but the impact on interval variability was negligible. Nevertheless, the variability analyses of contraction strengths revealed significant differences in sample entropies and Higuchi dimensions between control and acetylcholine treated tissue. Therefore, the proposed high-speed video camera technique might represent a non-invasive tool that allows long-lasting recordings for detecting variations in beating behavior over a large range of scales.

Introduction

Heart rate variability (HRV) refers to variations in the time intervals between two consecutive heart beats and serves as a diagnostic and prognostic tool for cardiac as well as non-cardiac diseases, e.g., heart failure, aging, Parkinson's disease, diabetes, and sepsis (Goldberger et al., 2002; Devos et al., 2003; Kudat et al., 2016; de Castilho et al., 2017; Elstad et al., 2018; Sessa et al., 2018). These variations are mainly attributed to dynamic changes of neuroendocrine inputs on ion channel activity in the sinoatrial node SAN, but a certain degree of beat to beat variability is inherently present at the level of the isolated heart, within the isolated SAN and also at the level of single sinoatrial pacemaker cells (Lombardi and Stein, 2011; Papaioannou et al., 2013; Yaniv et al., 2014a; Zaniboni et al., 2014).

There is a large number of quantitative algorithms to investigate these interval variations in autorhythmic cardiac tissue, cell clusters, or single cells, including spectral, linear, and nonlinear methods. Power-law behavior of beat to beat intervals BBIs analyzed by the power spectral method has been shown for neonatal rat cardiomyocytes in cultured tissue layers measured by microelectrode arrays (Ponard et al., 2007). Long-range correlations were also detected in extracellular electrograms of human embryonic stem cell-derived cardiomyocyte clusters by using again spectral methods (Mandel et al., 2012). Furthermore, fractal-like behavior has been reported for rabbit sinoatrial node tissue and for a small percentage of single sinoatrial node cells by using power law and detrended fluctuation analysis (Yaniv et al., 2014a) and in small clusters of chick embryonic cardiomyocytes (Ahammer et al., 2013). So far, investigations have focused on variabilities in the time domain of both, electrical and contraction signals. The underlying processes are tightly linked via the excitation-contraction-coupling (Eisner et al., 2017) and hence, the time structure of the electrical process substantially shapes not only the frequency of contraction but also its magnitude. Thus, it is reasonable to assume that also the variability of the contraction strength shows long-term correlations.

Therefore, in this study we propose the investigation of contraction strengths and their variabilities additionally to interval variabilities in SAN tissue. In detail, we evaluated beat to beat interval variabilities and beat to beat contraction strength CS variabilities of murine atrial preparations that contained the SAN region by means of high-speed camera video recordings. Each image of a video represented a time stamp and contractions of the tissue were recorded as changes in average gray values. Simultaneous measurements of extracellular potentials using a cardiac-near-field electrode validated beat to beat intervals of video recordings. Measurements of the spontaneous activity of tissue samples were performed before and after the administration of acetylcholine ACh, the predominant transmitter of the parasympathetic nervous system. Its effects on atrial tissue are already well investigated and include a decrease in beating rate and force of contraction (Kitazawa et al., 2009).

Our main objectives were to determine the suitability of video recordings to register BBIs and CSs and to analyze changes of nonlinear measures in the variabilities of these two parameters due to ACh treatment. Sample entropy and Higuchi dimension are popular estimators capturing intrinsic nonlinear patterns in time series of measured signals (Higuchi, 1988; Richman and Moorman, 2000). We hypothesized that ACh significantly affects sample entropies and Higuchi dimensions of BBI and CS variabilities. To distinguish actual values from white noise, surrogate data series were constructed and analyzed.

In summary, this high-speed camera video recording-technique provides a promising tool to thoroughly investigate beat to beat behavior regarding absolute values of beating rate and contraction strength as well as their variabilities in autorhythmic tissue.

Methods

Tissue Preparation

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).

Video Acquisition

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.

FIGURE 1

Figure 1. Microscopic sample image and sample recordings. (A) Sample image of a tissue preparation including the intact atria and the region of the sinoatrial node SAN. The left atrium cannot be seen in the image. RA, right atrium; CT, crista terminalis; IAS, interatrial septum; IVC, inferior vena cava; SVC, superior vena cava; CNF, cardiac-near-field electrode. The region of interest ROI (red rectangle) with a size of 160 × 160 pixels was set to a region yielding high changes of average gray values. This was usually the case when regions of a thick tissue layer (e.g., trabeculae) moved into the ROI during contractions. (B) Sample electrical recording lasting about 4 s (down-sampled to 5 kHz for the graphical representation). Four subsequent beats are depicted and the red plus symbols indicate data points for computing beat to beat intervals BBIs. A single beat is zoomed out. (C) Sample average gray values $\overline{G}$ according to Equation (1) from about 2,000 (out of 300,000) images. Red plus symbols mark data points that defined the BBIs and in addition with the baseline points (red dots) the contraction strengths CSs. Interval as well as contraction strength variations are clearly visible for these four contractions. (D) Finding correct baseline points is crucial for the determination of the CSs. This signal sample shows some beats (red plus symbols) and baseline points (red dots). The first two baseline points were found very well, in contrast to the last two baseline points. Accordingly, CSs for these two points cannot be accurately computed.

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.

Electrical Recordings

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 Signal Generation From Videos

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

$\begin{array}{ll}{\overline{G}}_i=\frac{1}{N_P}{{\displaystyle \sum }}_{p=1}^{N_P}{g}_{i,p},& (1)\end{array}$

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).

$\begin{array}{ll}BB{I}_b=t_b-t_{b-1},& (2)\end{array}$

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.

FIGURE 2

Figure 2. Beat to beat intervals BBIs and contraction strengths CSs, ΔCSs, rCSs computed from optical (video) recordings. (A) BBIs computed according to Equation (2). This tissue sample yielded actually 660 data points (number of beats N_B = 661). (B) Contraction strengths CSs computed according to Equation (3) for the same tissue sample with 660 data points. (C) Differences of contraction strengths ΔCS computed as the differences of the average gray values $\overline{G}s$ of the beating signal according to Equation (4). (D) Relative contraction strengths rCSs computed by integration of the ΔCS values according to Equation (5).

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

$\begin{array}{ll}C{S}_b=|{\overline{G}}_b-baselin{e}_b|,& (3)\end{array}$

with CS_b the contraction strength of the beat b, ${\overline{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.

$\begin{array}{ll}\Delta {CS}_b={\overline{G}}_b-{\overline{G}}_{b-1},& (4)\end{array}$

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.

$\begin{array}{ll}rC{S}_b=\Delta {CS}_b+rC{S}_{b-1},& (5)\end{array}$

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 and Higuchi Dimension

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:

$\begin{array}{ll}{X}_m(i)=\left[x(i),x(i+1),\dots ,x(i+m-1)\right].& (6)\end{array}$

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

$\begin{array}{ll}d\left[X_m(i),X_m(j)\right]={\displaystyle \underset{k=1,2,..,m}{max}}\left\{|X_m(i+k-1)-X_m(j+k-1)|\right\}& (7)\end{array}$

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:

$\begin{array}{ll}{C}_i^m(r)=\frac{numberof{X}_m(j)whered\left[X_m(i),X_m(j)\right]\le r}{N_B-m+1}.& (8)\end{array}$

The normalized number of sums can be computed using

$\begin{array}{ll}{B}^m(r)=\frac{1}{N_B-m}{{\displaystyle \sum }}_{i=1}^{N_B-m}{C}_i^m(r).& (9)\end{array}$

Finally, SampEn was computed with

$\begin{array}{ll}SampEn(m,r,N_B)=-ln\left(\frac{B^{m+1}(r)}{B^m(r)}\right).& (10)\end{array}$

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:

$\begin{array}{ll}{S}_m(k):x(m),x(m+k),x(m+2k),\dots ,x(m+\lfloor \frac{N_B-m}{k}\rfloor k)& (11)\end{array}$

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

$\begin{array}{ll}{L}_m(k)=\frac{1}{k}\left\{\left({\displaystyle {\sum }_{i=1}^{\lfloor \frac{N_B-m}{k}\rfloor }\left|x\left(m+ik\right)-x(m+\left(i-1\right)k)\right|}\right)\frac{N_B-1}{\lfloor \frac{N_B-m}{k}\rfloor k}\right\}& (12)\end{array}$

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

$\begin{array}{ll}L(k)=\frac{1}{k}{{\displaystyle \sum }}_{m=1}^k{L}_m(k).& (13)\end{array}$

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

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.

Results

Linear variance measures of the beat to beat interval are dependent on absolute values and are only well suited for linear stochastic processes. Nonlinear signals with random correlations or Random walk like signals can be well investigated with scale independent measures such as the sample entropy SampEn or the Higuchi dimension D_H.

Double Log Plots for Higuchi Dimensions

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).

FIGURE 3

Figure 3. Higuchi dimension D_H double log plots. Data from optical recordings. (A) Typical double logarithmic plots of D_H showing very linear regressions for the control and the acetylcholine treated cases. (B) Another sample of double logarithmic plots showing occasional deviations from the linear regression, in this particular case for acetylcholine (blue dots).

Correlation of Optically and Electrically Recorded Beat to Beat Intervals

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].

FIGURE 4

Figure 4. Opto-electrical correlation. (A) Some representative beat to beat intervals BBIs computed from simultaneously measured electrical and optical (video) recordings. (B) Regression of electrically and optically gained BBIs from one experiment (i.e., one mouse). The linear slope is very close to one and the coefficient of determination R² is very high.

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.

FIGURE 5

Figure 5. Sample entropy SampEn (left column) and Higuchi dimension D_H (right column) of beat to beat intervals BBIs determined optically and electrically. (A) Net plots showing 13 SampEn and D_H values for BBIs. (B) Box plots of these 13 experiments, p-values with df = 12, median test. (C) Scatterplots of these 13 experiments. Dashed lines represent theoretical correlations with the slope of one.

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.

Correlation of Contraction Strengths and Beat to Beat Intervals

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.

FIGURE 6

Figure 6. Scatterplots of contraction strengths CSs vs. beat to beat intervals BBIs in milliseconds. Data from optical recordings. Red dots correspond to control tissue, blue dots to ACh treated tissue. For each individual plot the experiment (mouse) number #, in brackets the actual Spearman's rank correlation coefficient r_s and the degrees of freedom df are depicted.

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.

Acetylcholine Induced Changes of Beat to Beat Intervals and Contraction Strengths

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).

FIGURE 7

Figure 7. Medians of beat to beat intervals BBIs (left column) and contraction strengths CSs (right column). Data from optical recordings. The control group is depicted in red and the values for the ACh treated group in blue. (A) Net plots of all experiments, showing individual and pairwise values. (B) Box plots of all experiments show statistically significant differences between Con and ACh groups, p < 0,001, df = 8, median test.

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).

Acetylcholine Induced Changes to Sample Entropies and Higuchi Dimensions

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.

FIGURE 8

Figure 8. Sample entropy SampEn (left column) and Higuchi dimension D_H (right column) of beat to beat intervals BBIs. Data from optical recordings. The control group is depicted in red and the ACh treated group in blue. (A) Net plots of all experiments, showing individual and pairwise values. (B) Box plots of all experiments. Differences are not statistically significant (SampEn p = 0.11, df = 8, D_H p = 0.97, df = 8, median test). (C) Mean SampEn and D_H values of shuffled (50x) data series.

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.

TABLE 1

Table 1. Median test of sample entropy SampEn and Higuchi dimension D_H values from optical recordings against means of 50x shuffled data series according to SurrEval-2.

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).

FIGURE 9

Figure 9. Sample entropy SampEn (left column) and Higuchi dimension D_H (right column) of contraction strengths CSs. Data from optical recordings. The control group is depicted in red and the ACh treated group in blue. (A) Net plots of all experiments, showing individual and pairwise values. (B) Box plots of all experiments. The difference in SampEn is statistically significant, p < 0.001, df = 8, median test. The difference in D_H is statistically borderline (p = 0.05, df = 8). (C) Mean SampEn and D_H values of shuffled (50x) data series.

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).

FIGURE 10

Figure 10. Sample entropy SampEn (left column) and Higuchi dimension D_H (right column) of relative contraction strengths rCSs. Data from optical recordings. The control group is depicted in red and the ACh treated group in blue. (A) Net plots of all experiments, showing individual and pairwise values. (B) Box plots of all experiments. The differences in SampEn and D_H are statistically significant, p < 0.001, df = 8 and p = 0.03, df = 8 respectively, median test. (C) Mean SampEn and D_H values of shuffled (50x) data series.

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).

Discussion

Beat to beat intervals are commonly investigated in order to detect nonlinear correlations in time signals. This study proposes an optical method, particularly, a high-speed video technique to detect mechanical contractions of the heart tissue. Usual video recordings with frame rates of about 30 fps or software triggered acquisitions are a convenient way for spectral analyses or computing beating frequencies (Kojima et al., 2006; Chan et al., 2009; Fassina et al., 2011; Hsiao et al., 2013; Ahola et al., 2014). But obviously, video frame rates must be higher for high beating rates (De Luca et al., 2014) or accurate detections of beating events (Stummann et al., 2008). Our high-speed video recordings allowed us to extract beat to beat intervals BBIs as well as the contraction strengths CSs and the relative contraction strengths rCSs, because the average gray value of an image was directly proportional to the mechanical contraction. Variation analyses with two distinguished nonlinear measures, the sample entropy SampEn and the Higuchi dimension D_H, revealed that this video technique is able to produce consistent results for BBIs as well as for CSs and rCSs.

The detection of the baseline (relaxed tissue) may be prone for errors such as measurement noise or optical drifts during the recording and thus, we proposed the second contraction parameter rCS. This is basically the varying contraction content of the signal, without the absolute value, drift or offset. SampEn, D_H and other scale independent nonlinear measures or fractal dimensions are not dependent on absolute values and consequently, rCS is an appropriate and very promising parameter for variance analyses.

To our knowledge, studies of BBV using isolated SAN tissue are very limited. Since no consensus exists to classify possible physiological artifacts (e.g., ectopic beats) in this in vitro preparation, we analyzed the original signals without any editing that could lead to a loss of valuable information.

We observed that ACh changed the beating behavior of the sinus node tissue by significantly reducing beating frequency as well as contraction strength, which is in accordance to previously published data (Kitazawa et al., 2009). Application of SampEn and D_H, two frequently used nonlinear measures for time signal variations, revealed a significant change of variabilities in the contraction strength but not in the beat to beat interval. The observed reduction of nonlinear measures indicates that the contraction process estimated by CS and rCS becomes more regular in the SAN tissue after ACh application. The spontaneous activity of pacemaker cells in the SAN tissue is based on two tightly linked clocks referred to as calcium and membrane clock (Lakatta and DiFrancesco, 2009). Both clocks exhibit inherent random components which arise from stochastic opening and closing of transmembrane ion channels (Krogh-Madsen et al., 2017) in the case of the membrane clock and from spontaneous stochastic calcium release via sarcoplasmic ryanodine receptors (Yaniv et al., 2014b) in the case of the calcium clock. The spontaneous calcium release in turn activates the sodium-calcium exchanger thereby triggering the action potential upstroke and subsequently a massive calcium release from sarcoplasmic reticulum, thus coupling excitation with contraction. ACh is an important modulator of SAN beating frequency as well as of contraction strength, particularly in the adjacent atrial tissue (Okada et al., 2013). Activation of muscarinic receptors by ACh causes multiple effects on the membrane and the calcium clock via G-protein coupled signaling ultimately reducing beating frequency and contractility (Harvey and Belevych, 2003). This is in line with our results. The physiological mechanisms underlying the observed increase in contraction strength regularity by ACh in our study are currently unknown. Theoretically, a reduced randomness in membrane and/or calcium clock as well as in the contraction process itself could account for our observation. It is noteworthy that in our study ACh increases CS regularity but not BBI regularity. This may be due to the fact that the beating behavior in the time domain is determined solely by sinus node pacemaking, whereas CS regularity may also depend on the effect of ACh on atrial tissue present in our preparations.

Studies on ACh effects on SAN cells/tissue using nonlinear measures are very scarce. Yaniv et al (Yaniv et al., 2014a) investigated the beating rate variability at different levels of integration from the heart in vivo to single pacemaker cells by linear (coefficient of variation) and nonlinear (approximate entropy, power law and detrended fluctuation analysis) measures. Their results show that beating interval regularity increased in the order in vivo, denervated heart, isolated SAN tissue, but decreased again in single pacemaker cells. However, single SAN cells showed fractal-like behavior only to a small percentage. Carbachol, a parasympathomimetic drug, decreased regularity of beating intervals of single SAN cells. Since this group analyzed the effect of parasympathetic stimulation on beating behavior only in the time domain and at the single cell level, a direct comparison to our results does not seem to be reasonable. Clearly, further studies are needed to elucidate the underlying physiological mechanisms of muscarinic stimulation on nonlinear measures in SAN cells/tissue.

Compared to commonly used mechanical force transducer measurements (Kihara and Morgan, 1991; Torres and Janssen, 2011; Koyani et al., 2017), the high-speed video technique used in this study seems to be an appropriate, contact-free tool to quantify changes in contraction strength variability. The SAN preparation represents a very sensitive and fragile tissue which could be easily damaged by hooks or threads of mechanical transducers. Moreover, SAN tissue may not be very suitable for mechanical force measurements because forces developed by the low tissue mass are very small leading to low signal amplitudes and hence low signal to noise ratios. The suggested video method does not provide absolute values of contraction forces. Absolute values are certainly a prerequisite for linear analyses but not for nonlinear investigations of variabilities that are per se independent of the absolute value. The electrical and the contractile processes underlying the measured BBIs and CSs are not reducible to each other, although tightly linked. Thus, the measured video signal provides information about these two distinct processes and consequently allows a more comprehensive characterization compared to frequently used electrode techniques.

For optically computed SampEn values, all surrogate evaluations strengthen the experimental findings that long-range correlations are present in BBI, CS, and rCS signals. Shuffled values are always statistically higher compared to experimentally obtained values and the difference between control and ACh treated tissue vanishes compared to the experimental case. Hence, the investigated physiological signals contain inherent nonlinear patterns in the interval and the contraction strength domains, justifying the application of the chosen nonlinear measures. Furthermore, the increase of regularity due to ACh (lower SampEn and D_H values) seems to be caused by deterministic and not by random processes.

The surrogate analyses for D_H values agree very well to SampEn evaluations, except for BBI signals (see Table 1). Particularly, some experimental values were already close to two, implicating a high degree of underlying random processes. Obviously, data time series shuffling did not reveal any significant changes. Distinct D_H and SampEn surrogate results concerning BBIs imply different sensitivities to underlying random and deterministic physiological mechanisms. This may indicate that the Higuchi dimension and not the sample entropy is able to discriminate differences between interval and contraction strength signals, but further investigations are needed to corroborate this assumption. In order to rule out that high D_H values close to two were method specific, we additionally performed a detrended fluctuation analysis DFA (Peng et al., 1994; Goldberger et al., 2002). DFA characterizes white noise with α = 0.5 and Brownian noise with α = 1.5.The maximal window length was set to 30, comparable to k = 30 for D_H. The medians of all control cases including BBI, CS, and rCS (n = 27) were 1.98 for D_H and 0.61 for α (DFA) confirming the high degrees of randomness in the signals. Thus, the decreased D_H values under ACh treatment indicate a change from very low correlations to increased long-range correlations and self-affine processes.

Complexity can be defined as the presence of long-range correlations, arising from nonlinear interaction dominated dynamical processes being neither totally regular nor totally irregular (Van Orden et al., 2011). This concept has been successfully applied to discriminate healthy and pathological conditions, where a breakdown of long-term correlations and an according change in fractal dimension has been observed (Goldberger et al., 2002). In our case, SAN tissue preparations show a high degree of irregularity near white noise indicating a low complexity without long-range correlations or self-organizing mechanisms. This may be due to a loss of multiple and interwoven communication pathways and nonlinear dependencies present in the intact heart but not in the tissue preparation. This is supported by ACh as a relevant external stimulus that changes the interaction dominated dynamical system by reducing the degrees of freedom and by introducing long-range correlations, multiplicative interactions and feedback. ACh may be interpreted as a control parameter for the system. Phase transitions or bifurcations, dependent on control parameters exceeding critical values, may also play an important role for our tissue preparation, but at this stage further investigations are necessary to justify such interpretations.

In conclusion, the described technique represents a reliable, easy handling and long-lasting recording method, from which beating rate variabilities and contraction strength variabilities can be assessed.

Author Contributions

HA and SS are joint first authors. HA, SS, KZ-P, and BP designed the study. SS and PL dissected SAN tissues from mouse hearts. SS and ÁD carried out video experiments. SS and RA carried out electrical recordings. HA, MM-R, and ÁD developed mathematical algorithms. HA, SS, ÁD, and PL performed video analyses. RA analyzed electrical recordings. KZ-P and HA computed nonlinear parameters and statistical tests. HA, SS, and KZ-P wrote the manuscript. SS, RA, BP, KZ-P, HA, and MM-R revised the manuscript. HA and SS prepared the figures. All authors read and approved the manuscript.

Funding

This work was partly supported by research grants from the Austrian National Bank (16435, Anniversary Fund) and the Franz-Lanyar-Foundation (#405) to BP.

Data Availability Statement

BBI, CS, and rCS data series for each experiment (Con and ACh treated) are provided as csv file in the supplement.

Conflict of Interest Statement

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.

Acknowledgments

We gratefully thank Dr. Matteo Mangoni and Dr. Pietro Mesirca from the Institut de Génomique Fonctionelle, Montpellier, France, for teaching the preparation technique of sinoatrial node tissue to our group. We thank the reviewers for their valuable suggestions.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphys.2018.00546/full#supplementary-material

References

Ahammer, H., Scherubel, S., Arnold, R., Zorn-Pauly, K., and Pelzmann, B. (2013). “Beat to beat variability of embryonic chick heart cells under septic conditions: application and evaluation of entropy as well as fractal measures,” in 35th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (Osaka: EMBC), 5566–5569.

Google Scholar

Ahola, A., Kiviaho, A. L., Larsson, K., Honkanen, M., Aalto-Setäl,ä, K., and Hyttinen, J. (2014). Video image-based analysis of single human induced pluripotent stem cell derived cardiomyocyte beating dynamics using digital image correlation. Biomed. Eng. Online 13:39. doi: 10.1186/1475-925X-13-39

PubMed Abstract | CrossRef Full Text | Google Scholar

Chan, P. K., Lin, C. C., and Cheng, S. H. (2009). Noninvasive technique for measurement of heartbeat regularity in zebrafish (Danio rerio) embryos. BMC Biotechnol. 9:11. doi: 10.1186/1472-6750-9-11

PubMed Abstract | CrossRef Full Text | Google Scholar

de Castilho, F. M., Ribeiro, A. L. P., da Silva, J. L. P., Nobre, V., and de Sousa, M. R. (2017). Heart rate variability as predictor of mortality in sepsis: a prospective cohort study. PLoS ONE 12:e0180060. doi: 10.1371/journal.pone.0180060

PubMed Abstract | CrossRef Full Text | Google Scholar

De Luca, E., Zaccaria, G. M., Hadhoud, M., Rizzo, G., Ponzini, R., Morbiducci, U., et al. (2014). ZebraBeat: a flexible platform for the analysis of the cardiac rate in zebrafish embryos. Sci. Rep. 4:4898. doi: 10.1038/srep04898

CrossRef Full Text | Google Scholar

Devos, D., Kroumova, M., Bordet, R., Vodougnon, H., Guieu, J. D., Libersa, C., et al. (2003). Heart rate variability and Parkinson's disease severity. J. Neural Transm. 110, 997–1011. doi: 10.1007/s00702-003-0016-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Eisner, D. A., Caldwell, J. L., Kistamás, K., and Trafford, A. W. (2017). Calcium and excitation-contraction coupling in the heart. Circ. Res. 121, 181–195. doi: 10.1161/CIRCRESAHA.117.310230

PubMed Abstract | CrossRef Full Text | Google Scholar

Elstad, M., O'Callaghan, E. L., Smith, A. J., Ben-Tal, A., and Ramchandra, R. (2018). Cardiorespiratory interactions in humans and animals: rhythms for life. Am. J. Physiol. Heart Circ. Physiol. doi: 10.1152/ajpheart.00701.2017. [Epub ahead of print].

PubMed Abstract | CrossRef Full Text | Google Scholar

Fassina, L., Di Grazia, A., Naro, F., Monaco, L., De Angelis, M. G., and Magenes, G. (2011). Video evaluation of the kinematics and dynamics of the beating cardiac syncytium: an alternative to the Langendorff method. Int. J. Art. Organs 34, 546–558. doi: 10.5301/IJAO.2011.8510

PubMed Abstract | CrossRef Full Text | Google Scholar

Galvez Coyt, G., Munoz Diosdado, A., Balderas Lopez, J. A., del Rio Correa, J. L., and Angulo Brown, F. (2013). Higuchi's method applied to the detection of periodic components in time series and its application to seismograms. Revista Mexic. Fisica S 59, 1–6.

Google Scholar

Glukhov, A. V., Fedorov, V. V., Anderson, M. E., Mohler, P. J., and Efimov, I. R. (2010). Functional anatomy of the murine sinus node: high-resolution optical mapping of ankyrin-B heterozygous mice. Am. J. Physiol. Heart Circ. Physiol. 299, H482–H491. doi: 10.1152/ajpheart.00756.2009

PubMed Abstract | CrossRef Full Text | Google Scholar

Goldberger, A. L., Amaral, L. A., Hausdorff, J. M., Ivanov, P. C. H., Peng, C. K., and Stanley, H. E. (2002). Fractal dynamics in physiology: alterations with disease and aging. Proc. Natl. Acad. Sci. U.S.A. 99, 2466–2472. doi: 10.1073/pnas.012579499

PubMed Abstract | CrossRef Full Text | Google Scholar

Harvey, R. D., and Belevych, A. E. (2003). Muscarinic regulation of cardiac ion channels. Br. J. Pharmacol. 139, 1074–1084. doi: 10.1038/sj.bjp.0705338

PubMed Abstract | CrossRef Full Text | Google Scholar

Higuchi, T. (1988). Approach to an irregular time-series on the basis of the fractal theory. Physica D 31, 277–283. doi: 10.1016/0167-2789(88)90081-4

CrossRef Full Text | Google Scholar

Hofer, E., Keplinger, F., Thurner, T., Wiener, T., Sanchez-Quintana, D., Climent, V., et al. (2006). A new floating sensor array to detect electric near fields of beating heart preparations. Biosens. Bioelectron. 21, 2232–2239. doi: 10.1016/j.bios.2005.11.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Hsiao, C.-W., Bai, M.-Y., Chang, Y., Chung, M.-F., Lee, T.-Y., Wu, C.-T., et al. (2013). Electrical coupling of isolated cardiomyocyte clusters grown on aligned conductive nanofibrous meshes for their synchronized beating. Biomaterials 34, 1063–1072. doi: 10.1016/j.biomaterials.2012.10.065

PubMed Abstract | CrossRef Full Text | Google Scholar

Kainz, P., Mayrhofer-Reinhartshuber, M., and Ahammer, H. (2015). IQM: an extensible and portable open source application for image and signal analysis in java. PLoS ONE 10:e0116329. doi: 10.1371/journal.pone.0116329

PubMed Abstract | CrossRef Full Text | Google Scholar

Kesić, S., and Spasić, S. Z. (2016). Application of Higuchi's fractal dimension from basic to clinical neurophysiology: a review. Comp. Methods Progr. Biomed. 133, 55–70. doi: 10.1016/j.cmpb.2016.05.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Kihara, Y., and Morgan, J. P. (1991). Abnormal Cai2+ handling is the primary cause of mechanical alternans: study in ferret ventricular muscles. Am. J. Physiol. Heart Circ. Physiol. 261, H1746–H1755. doi: 10.1152/ajpheart.1991.261.6.H1746

PubMed Abstract | CrossRef Full Text | Google Scholar

Kitazawa, T., Asakawa, K., Nakamura, T., Teraoka, H., Unno, T., Komori, S., et al. (2009). M3 muscarinic receptors mediate positive inotropic responses in mouse atria: a study with muscarinic receptor knockout mice. J. Pharmacol. Exp. Ther. 330, 487–493. doi: 10.1124/jpet.109.153304

PubMed Abstract | CrossRef Full Text | Google Scholar

Kojima, K., Kaneko, T., and Yasuda, K. (2006). Role of the community effect of cardiomyocyte in the entrainment and reestablishment of stable beating rhythms. Biochem. Biophys. Res. Commun. 351, 209–215. doi: 10.1016/j.bbrc.2006.10.037

PubMed Abstract | CrossRef Full Text | Google Scholar

Koyani, C. N., Kolesnik, E., Wölkart, G., Shrestha, N., Scheruebel, S., Trummer, C., et al. (2017). Dipeptidyl peptidase-4 independent cardiac dysfunction links saxagliptin to heart failure. Biochem. Pharmacol. 145, 64–80. doi: 10.1016/j.bcp.2017.08.021

PubMed Abstract | CrossRef Full Text | Google Scholar

Krogh-Madsen, T., Kold Taylor, L., Skriver, A. D., Schaffer, P., and Guevara, M. R. (2017). Regularity of beating of small clusters of embryonic chick ventricular heart-cells: experiment vs. stochastic single-channel population model. Chaos 27:093929. doi: 10.1063/1.5001200

PubMed Abstract | CrossRef Full Text | Google Scholar

Kudat, H., Akkaya, V., Sozen, A. B., Salman, S., Demirel, S., Ozcan, M., et al. (2016). Heart rate variability in diabetes patients. J. Int. Med. Res. 34, 291–296. doi: 10.1177/147323000603400308

PubMed Abstract | CrossRef Full Text | Google Scholar

Lakatta, E. G., and DiFrancesco, D. (2009). What keeps us ticking: a funny current, a calcium clock, or both? J. Mol. Cell. Cardiol. 47, 157–170. doi: 10.1016/j.yjmcc.2009.03.022

PubMed Abstract | CrossRef Full Text | Google Scholar

Lombardi, F., and Stein, P. K. (2011). Origin of heart rate variability and turbulence: an appraisal of autonomic modulation of cardiovascular function. Front. Physiol. 2:95. doi: 10.3389/fphys.2011.00095

PubMed Abstract | CrossRef Full Text | Google Scholar

Mandel, Y., Weissman, A., Schick, R., Barad, L., Novak, A., Meiry, G., et al. (2012). Human embryonic and induced pluripotent stem cells-derived cardiomyocytes exhibit beat rate variability and power-law behavior. Circulation 125, 883–893. doi: 10.1161/CIRCULATIONAHA.111.045146

CrossRef Full Text | Google Scholar

Okada, M., Noma, C., Yamawaki, H., and Hara, Y. (2013). Negative Inotropic effect of carbachol and interaction between acetylcholine receptor-operated potassium channel (K.ACh Channel) and GTP binding protein in mouse isolated atrium - a novel methodological trial. J. Vet. Med. Sci. 75, 377–380. doi: 10.1292/jvms.12-0369

PubMed Abstract | CrossRef Full Text | Google Scholar

Papaioannou, V. E., Verkerk, A. O., Amin, A. S., and de Bakker, J. M. (2013). Intracardiac origin of heart rate variability, pacemaker funny current and their possible association with critical illness. Curr. Cardiol. Rev. 9, 82–96. doi: 10.2174/157340313805076359

PubMed Abstract | CrossRef Full Text | Google Scholar

Peng, C.-K., Buldyrev, S. V., Havlin, S., Simons, M., Stanley, H. E., and Goldberger, A. L. (1994). Mosaic organization of DNA nucleotides. Phys. Rev. E Stat. Phys. Plasmas Fluids Relat. Interdisci.p Top. 49, 1685–1689. doi: 10.1103/PhysRevE.49.1685

PubMed Abstract | CrossRef Full Text | Google Scholar

Ponard, J. G., Kondratyev, A. A., and Kucera, J. P. (2007). Mechanisms of intrinsic beating variability in cardiac cell cultures and model pacemaker networks. Biophys. J. 92, 3734–3752. doi: 10.1529/biophysj.106.091892

PubMed Abstract | CrossRef Full Text | Google Scholar

R Core Team (2017). R: A Language and Environment for Statistical Computing. Vienna: R Core Team.

Richman, J. S., and Moorman, J. R. (2000). Physiological time-series analysis using approximate entropy and sample entropy. Am. J. Physiol. Heart Circ. Physiol. 278, H2039–H2049. doi: 10.1152/ajpheart.2000.278.6.H2039

PubMed Abstract | CrossRef Full Text | Google Scholar

RStudio (2016). RStudio: Integrated Development for R (Boston, MA: RStudio Inc.).

Sessa, F., Anna, V., Messina, G., Cibelli, G., Monda, V., Marsala, G., et al. (2018). Heart rate variability as predictive factor for sudden cardiac death. Aging 10, 166–177. doi: 10.18632/aging.101386

PubMed Abstract | CrossRef Full Text | Google Scholar

Stummann, T. C., Wronski, M., Sobanski, T., Kumpfmueller, B., Hareng, L., Bremer, S., et al. (2008). Digital movie analysis for quantification of beating frequencies, chronotropic effects, and beating areas in cardiomyocyte cultures. ASSAY Drug Dev. Technol. 6, 375–385. doi: 10.1089/adt.2008.129

PubMed Abstract | CrossRef Full Text | Google Scholar

Torrente, A. G., Zhang, R., Zaini, A., Giani, J. F., Kang, J., Lamp, S. T., et al. (2015). Burst pacemaker activity of the sinoatrial node in sodium–calcium exchanger knockout mice. Proc. Natl. Acad. Sci. U.S.A. 112, 9769–9774. doi: 10.1073/pnas.1505670112

PubMed Abstract | CrossRef Full Text | Google Scholar

Torres, C. A., and Janssen, P. M. (2011). Contractile strength during variable heart duration is species and preload dependent. Biomed. Res. Int. 2011:294204. doi: 10.1155/2011/294204

PubMed Abstract | CrossRef Full Text | Google Scholar

Van Orden, G., Kloos, H., and Wallot, S. (2011). “Living in the pink: Intentionality, wellbeing, and complexity,” in Philosophy of Complex Systems, Handbook of the Philosophy of Science, ed Hooker (Amsterdam: Elsevier), 629–672.

Google Scholar

Wilcox, R. R., and Rousselet, G. A. (2018). A guide to robust statistical methods in neuroscience. Curr. Protoc. Neurosci. 82, 8.42.1–8.42.30. doi: 10.1002/cpns.41

PubMed Abstract | CrossRef Full Text | Google Scholar

Wilcox, R. R. (2016). Understanding and Applying Basic Statistical Methods Using R. Hoboken, NJ: Joh Wiley & Sons.

Google Scholar

Yaniv, Y., Ahmet, I., Liu, J., Lyashkov, A. E., Guiriba, T.-R., Okamoto, Y., et al. (2014a). Synchronization of sinoatrial node pacemaker cell clocks and its autonomic modulation impart complexity to heart beating intervals. Heart Rhythm 11, 1210–1219. doi: 10.1016/j.hrthm.2014.03.049

PubMed Abstract | CrossRef Full Text | Google Scholar

Yaniv, Y., Lyashkov, A. E., Sirenko, S., Okamoto, Y., Guiriba, T.-R., Ziman, B. D., et al. (2014b). Stochasticity intrinsic to coupled-clock mechanisms underlies beat-to-beat variability of spontaneous action potential firing in sinoatrial node pacemaker cells. J. Mol. Cell. Cardiol. 77, 1–10. doi: 10.1016/j.yjmcc.2014.09.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Zaniboni, M., Cacciani, F., and Lux, R. L. (2014). Beat-to-beat cycle length variability of spontaneously beating guinea pig sinoatrial cells: relative contributions of the membrane and calcium clocks. PLoS ONE 9:e100242. doi: 10.1371/journal.pone.0100242

PubMed Abstract | CrossRef Full Text | Google Scholar