All animal-handling procedures were concordant with the Canadian Council on Animal Care guidelines and were approved by the institutional Animal Research Ethics Committee. Isolation was performed according to the protocol of the neonatal cardiomyocyte isolation kit from Worthington. In summary, 1–3 days old rats (Sprague-Dawley, Charles River) were sacrificed by decapitation. Beating hearts were removed from the rats and immediately put in cold Ca²⁺ and Mg²⁺-free Hank’s Balanced Salt Solution. The ventricular muscle was selected by excision and the tissue was minced on ice into 1–3 mm³ pieces. The mixture was subjected to purified enzymatic digestion (trypsin and collagenase). Isolated cells (enriched cardiomyocytes) were counted and seeded at a density of 10⁶ cells/mL in the seeding area of the membrane pre-coated with 0.2% porcine-derived gelatin (G1890, Sigma) and 0.00125% fibronectin solution (F1141, Sigma). Cells were grown for 24 h in DMEM (319-050-CL, Wisent) with 5% fetal bovine serum (FBS, SH30396.03, Thermo Fisher Scientific Co., Ltd.) and 1% penicillin/streptomycin (450-201-EL, Wisent). Cardiomyocytes were then FBS starved with 1% penicillin/streptomycin in DMEM 24 h prior to the experiments. For the set of experiments used to study the effects of electrically pacing the cardiomyocytes at different pacing cycle length (PCL), the cells were washed and the medium changed to fresh DMEM without FBS just prior to record the activity (thus not being starved for 24 h prior to the experiments).

Acute effects of the β-adrenergic agonist, isoproterenol (ISO, I6504, Sigma-Aldrich) was studied by videomicroscopy with final concentration of 100 nmol/L at 1 min after injection.

Phase contrast videos of neonatal cardiomyocytes were acquired after 48 h post seeding in vitro with a Dalsa HM640 camera (N_y = 640 × N_x = 480 pixels) at rates of 30, 50, or 100 frame per second (fps) coupled to an inverted Nikon optical microscope (10× magnification). The field of view (FOV) covered by the camera was 0.44 mm by 0.33 mm.

The Morotti et al. (2014) model of mouse ventricular myocyte modified to integrate the Negroni et al. (2015) model of myofilament contraction as used by Surdo et al. (2017)¹ was simulated at PCL ranging from 200 to 1,000 ms to study the effects of period of activity. To account for the duration of transient solution when starting from the resting steady-state, simulation output at 30 and 60 s are considered in this study. The Negroni et al. (2015) of rabbit ventricular myocyte was also simulated at different PCL for a 30 s duration of simulation. A modified version of the rabbit model was also studied where two parameters of the myofilament section were multiplied by a factor of 5: Yb (binding rate of Ca²⁺) and f (transition rate between the 1st and 2nd state of the troponin system). The parameter Yv (part of the detachment function) was scaled down by a factor of 5.

The models were numerically integrated in Matlab (R2008, MathWorks Inc., Natick, MA) using the variable-step, variable-order solver based on the numerical differentiation formulas of orders 1–5 for stiff system (ode15s). Relative error tolerance was set to 10^–5.

All statistical analyses have been done in R (version 3.1.3). Linear regression was used to test for variable dependency on period of activity. The analysis of covariance (ANCOVA) was used to compare group and covariate effects. Means were compared with Student’s t-test.

The proposed initial approach is using the raw video signal and calculating for each pixel the difference between frames separated by a delay corresponding to a defined frame interval.

Let M(t) be the global camera frame recorded at time t. The composite signal ΔS is calculated using the following Equation:

where τ is a discrete time delay (or number of frames for delay multiplied by the time between frames Δt = 1/fps), N_y and N_x are, respectively, the number of horizontal and vertical pixels such that i and j are the pixels coordinates along the horizontal and vertical axis. An example of the composite signal ΔS for a spontaneously beating sample is shown in Figure 1 for different delay (τ = 1,4, and 8 multiplied by Δt). The shortest delay τ = Δt corresponding to a frame-by-frame difference shows the highest noise level (Figure 1A) while increasing τ to 4Δt (Figure 1B) and 8Δt (Figure 1C) decreased the amplitude difference between the high and low amplitude peaks.

In silico data and analysis reveal that the composite signal ΔS calculated from Eq. (1) can be interpreted as follow. The absolute derivative of the time-dependent cell length calculated from simulated sarcomere length (SL) obtained from the simulation is given by:

where SL is the sarcomere length. The change in SL calculated with the Surdo et al. (2017) model paced at a PCL of 500 ms is presented in Figure 2A. The time derivative of this contraction signal is depicted in Figure 2B which shows the initial contraction (negative derivative) followed by the positive relaxation signal. As shown in Figure 2C, there is a clear similarity between the rectified derivative given by Eq. (2) and the videomicroscopy signal obtained with Eq. (1) as shown in Figures 1A, 3A (highlighted in the inset of panel A).

Two specific measures are proposed that are presented in Figure 2D which corresponds to the amplitude ratio between the contraction and relaxation peaks (Ratio_C/R) and the time between the contraction and relaxation peaks (Δ_R–C). The effects of the PCL on contraction signal and specific measures are shown in Figure 2D which highlights the sensitivity of the measures on rhythm and calcium dynamics. Decreasing the PCL from 500 to 220 ms lead to a decrease of 6% of Δ_R–C and of 33% of Ratio_C/R.

Automatic discrimination between contraction and relaxation peaks is important in the perspective of user-friendliness and for the approach to have a clear potential in a high throughput screening/testing system. The approach is based on the observations of experimental data showing differences in amplitude and rate of decrease around the peaks between the contraction and relaxation. Starting from a ΔS signal as shown in Figure 3A, three conditions are being used to detect and validate classification of the peaks: 1- there is alternance between contraction and relaxation peaks, 2- clusters are usually separated in the variable space given by the time difference between peaks (Δt_markers) and amplitude of the peaks (ΔS_peaks) as depicted in Figure 3B, 3- experimental data has sharpest peaks for contraction and widest peaks for relaxation resulting in more rapid amplitude loss around the maximum ΔS peak amplitude. The last condition can be easily evaluated by taking the amplitude for 1 sample (first sample before and after the position of the peak) and 2 samples (second sample before and after the peak) around the peak divided by the peak amplitude (Ratio_drop: the ratios given by the value of the first and second samples around the peak divided by the value of the signal for the peak define). As such, Ratio_drop serves to evaluate how fast the amplitude decreases around the peak distinguishing fast changes (narrow contraction peaks)_and slower changes (wide relaxation peaks). The obtained data are presented in Figure 3C where a lower average Ratio_drop is found for the contraction peaks (left circles; 0.77 ± 0.03 and 0.40 ± 0.06 n.u.) compared to the relaxation peaks (right dots; 0.92 ± 0.02 and 0.78 ± 0.04 n.u.). Using the second neighboring points (what we labeled drop by two samples) helps separating the contraction and relaxation groups of points (the red points behind more clearly separated as a function of Ratio_drop).

An example of an analyzed acquisition is presented in Figure 4. The signal ΔS is shown in Figure 4A with detected peaks (red circles) and resulting classification highlighted for contraction (red dotted lines) and relaxation peaks (red dashed lines). The resulting series of spontaneous activity shows a constant period T (Figure 4B: 0.64 ± 0.01 s) almost no variation in Δ_R–C (Figure 4C: 0.22 ± 0.01 s). However, Ratio_C/R shows greater variability between samples (Figure 4C: 1.89 ± 0.06 n.u.) mainly due to variations in the maximum amplitude of the contraction peaks as seen on the signal in Figure 4A.

A set of 8 samples data obtained at PCL ranging from min. 500 to max. 2,000 ms are presented in Figure 5 (each color corresponding to a given cell culture sample except for the black line that corresponds to the mean value of the samples). <Δ_R–C> (Figure 5A) is the lowest for shortest cycle length and is augmenting with PCL increasing except for one sample. The average curve (black line) shows a slight biphasic shape with intermediate PCL values being slightly higher than for the highest PCL (6/8 samples clearly showing a biphasic shape).

Data for <Ratio_C/R> are presented in Figure 5B. The lowest value is usually found for lower PCL and is increasing with PCL increasing (7/8 samples). There is again a tendency for a biphasic shape of the curves which is more visible for the average curve (black line).

Long-term simulations of the Surdo et al. (2017) mouse ventricular model ionic model have been done for a comparison purpose with the experimental data. Simulation results for <Δ_R–C> are shown in Figure 6A. The general shape of the curves differs between the experiments and simulations. Similar to experimental data, simulated <Δ_R–C> is the lowest for shortest cycle length and is augmenting with PCL increasing for shortest PCL. There is, however, a particular biphasic shape (more like a narrow peak) as the maximum value is found around 240–245 ms before rapidly decreasing followed by a slow increase with PCL increasing. This narrow peak comes from the morphology of the |dSL/dt| relaxation phase which has a double hump shape. The peak in the <Δ_R–C> as a function of PCL comes from a change from the first hump to be higher to the second bump thus resulting from an added delay in the marker for relaxation (∼0.02 s). The maximum (experiments: 0.54 s; simulation: 0.21 s) and range (experiments: 0.14 s; simulation: 0.05 s) of <Δ_R–C> is also less than in experiments. We compared the results for two duration of simulation (30 and 60 s of simulation for different constant PCL) to highlight the possible transient behavior of the model. The rate main difference between shorter (30 s) and longer (60 s) simulations are a higher <Δ_R–C> right after the peak resulting in a slower increase with increasing PCL for PCL > 300 ms.

Simulated data for <Ratio_C/R> are presented in Figure 6B. Similar to experimental data, the lowest value is found for lower PCL and is increasing with PCL increasing. There is a clear biphasic shape of the curve that seems also to be found in some experimental samples in Figure 5B. The highest values in simulated <Ratio_C/R> is found near PCL = ∼375 ms. The biphasic shape comes from the difference in rate of change between the contraction and relaxation peaks with increasing PCL. Starting at PCL = 360 ms, the rate of increase in the amplitude of the relaxation peak as a function of increasing PCL is greater than the rate of increase of the contraction peak. As such, the Ratio_C/R for PCL > 360 ms is decreasing with PCL increasing. There is, however, a decrease in highest <Ratio_C/R> for longer simulations (60 s compared to 30 s) indicating that it is, at least in part, linked to transient behavior of the model. More importantly is the range of values which differs from the experimental data. <Ratio_C/R> is always less than one for simulated data with the mouse model while experimental data are mostly greater than one (except for some cases at lowest PCL). The differences are important because value lower than one means that the rate of relaxation is faster than the rate of contraction.

The differences between the experimental and simulated data lead us to test if it was because of the ionic model. Thus, the rabbit ventricular model by Negroni et al. (2015) was simulated with normal parameter values (ctl) and with two contraction parameters (modified, see “Materials and Methods” section for details). The change in parameters (between ctl and modified models) did not significantly impact the action potential duration when paced at 500 ms as shown in Supplementary Figure 1A. While the ctl model has a relaxation rate higher than the contraction rate (black curve in Supplementary Figure 1C) resulting in a Ratio_C/R < 1, the modified model (green curve) shows a faster contraction than relaxation (Ratio_C/R> 1). A visual comparison of the curves clearly shows that the modified model has longer Δ_R–C compared to the ctl model. <Δ_R–C> and <Ratio_C/R> as a function of PCL for the Negroni et al. (2015) model can be found in Supplementary Figures 2A,B, respectively. The ctl rabbit ventricular model shows an increase in <Δ_R–C> compared to Morotti et al. mouse model (maximum value around 0.25 s vs. ∼0.21 s). Interestingly, simulating the modified Negroni et al. (2015) rabbit model resulted in even longer <Δ_R–C> (∼0.35 s) getting closer to the experimental values. While <Ratio_C/R> for the ctl Negroni et al. (2015) rabbit ventricular model is always less than 1 the modified version resulted in <Ratio_C/R>> over almost all the PCL range.

A set of samples (n = 29) has been analyzed. For each video of 30 s duration, the average values <T>, <Δ_R–C>, and <Ratio_C/R> were calculated as the average of the temporal values obtained from the signal analysis of individual acquisition. A clear monotone increase in <Δ_R–C> with increasing <T> is found as depicted in Figure 7A with ∼40% change between the minimum and maximum values. The linear regression (p < 0.001) has a slope of 0.072 s/s (represented as a red line) and intercept of 0.166 s. As expected, the variability is greater for <Ratio_C/R> but a trend to an increasing ratio as <T> augments is found. However, the variation with the period is less with ∼23% change between the minimum and maximum values in the dataset (see Figure 7B). The fitted regression line (p < 0.005, red line on the panel) has a slope of 0.276 s^–1 and an intercept of 1.93.

The variation in the contraction measures that can be evaluated by our videomicroscopy approach has been tested with isoproterenol (ISO), a β–adrenergic agonist. Results are presented in Figure 8. As expected, the period of activity is significantly decreased with ISO compared to CTL (0.9 ± 0.6 s vs. 1.9 ± 2.4 s in CTL, p < 0.05). Similar variation of <Δ_R–C> is found between CTL and ISO groups as a function of <T> although <T> is in average smaller with ISO as expected (Figure 8A). Interestingly the <Ratio_C/R> showed statistically significant <T> and group effects with p < 0.001. The slope of the linear regression being larger for the ISO group compared to CTL (0.146 for CTL vs. 0.490 for ISO).

In order to see if Δ_R–C and Ratio_C/R may be interesting to evaluate cardiotoxicity, we used the change in sarcomere length data presented in Figure 5B in the study by Timolati et al. (2009) comparing the rat ventricular cardiomyocyte contraction paced with a period of 500 ms in control and treated with 1 μmol/L of doxorubicin for 48 h. The data was digitized and sampled with 10 ms between samples (shown in Supplementary Figure 3). Using the same approach as depicted in Figure 2 for the Surdo et al. (2017) ionic model, we calculated the mean of the last two contractions (skipping the first one) yielding Δ_R–C = 0.11 s in control compared to 0.07 s for the doxorubicin-treated cardiomyocyte. Ratio_C/R was increased in cardiomyocyte treated with doxorubicin (1.3 in control vs. 1.5 for the doxorubicin-treated cardiomyocyte). Although this result is very limited, it points to the interest of the measures in detecting contraction changes and evaluating cardiotoxicity.

All the previous analyses presented in this study are based on a global composite signal calculated using the entire FOV. The same approach can be used for sub-regions of the FOV by calculating a composite signal for each section. We present here two additional approaches aiming to study spatial-temporal differences in videomicroscopy signals.