Versatile open software to quantify cardiomyocyte and cardiac muscle contraction in vitro and in vivo

Contraction of muscle reflects its physiological state. Methods to quantify contraction are often complex, expensive and tailored to specific models or recording conditions, or require specialist knowledge for data extraction. Here we describe an automated, open-source software tool (MUSCLEMOTION) adaptable for use with standard laboratory and clinical imaging equipment that enables quantitative analysis of normal cardiac contraction, disease phenotypes and pharmacological responses. MUSCLEMOTION allowed rapid and easy measurement of contractility in (i) single cardiomyocytes from primary adult heart and human pluripotent stem cells, (ii) multicellular 2D-cardiomyocyte cultures, 3D engineered heart tissues and cardiac organoids/microtissues in vitro and (iii) intact hearts of zebrafish and humans in vivo. Good correlation was found with conventional measures of contraction in each system. Thus, using a single method for processing video recordings, we obtained reliable pharmacological data and measures of cardiac disease phenotype in experimental cell- and animal models and human echocardiograms.


40
The salient feature of cardiomyocytes (CMs) is their ability to undergo cyclic contraction and 41 relaxation, a feature critical for cardiac function. In many research laboratories and clinical settings it 42 is therefore essential that cardiac contraction can be quantified at multiple levels, from single cells to 43 multicellular or intact cardiac tissues. Measurement of contractility is relevant for analysis of disease 44 phenotypes, cardiac safety pharmacology, and longitudinal measures of cardiac function over time, 45 both in vitro and in vivo. In addition, human genotype-phenotype correlations, investigation of cardiac 46 disease mechanisms and the assessment of cardiotoxicity are increasingly performed on human 47

Algorithm development 135
The principle underlying the algorithm of MUSCLEMOTION is the assessment of contraction using 136 an intuitive approach quantifying absolute changes in pixel intensity between a reference frame and 137 the frame of interest, which can be described as values, while pixels that are highly changed result in high (white) values (Fig. 1a). Next, the mean 145 pixel intensity of the resulting image is measured. This is a quantitative measure of how much the 146 pixels have moved compared to the reference frame: more white pixels indicate more changing pixels 147 and, thus, more displacement. When a series of images is analysed relative to the same reference 148 image, the output describes the accumulated displacement over time (measure of displacement, Fig.  149   1b). 150 However, if a series of images is analysed with a reference frame that depends on the frame of interest 151 (e.g. '() = $/0 ), this results in a measure of the relative displacement per interframe interval. 152 We defined this parameter as contraction velocity (measure of velocity, Fig. 1b). 153 Since velocity is the first derivative of displacement in time, the first derivative of the measure 154 of displacement should resemble the measure of velocity derived from image calculations. To test the 155 linearity of the method, three movies of moving blocks were analysed. The block moved back and 156 forth at two different speeds in each direction (where 2 = 2 • 0 ): i) along the x-axis, ii) along the y-157 axis and iii) along both axes (Movie S1). As expected, the measure of displacement and velocity 158 showed a linear correlation (Fig. S1). This does not hold when the position of the block in $ does 159 not overlap the position of the block in '() , with a consequent saturation in the measure of displacement (i.e. max pixel white value, Fig. S2). Therefore, comparison of the differentially derived 161 velocities should approximately overlap in the absence of pixel saturation. This was used as a 162 qualitative parameter to determine whether the algorithm outputs were reliable. 163 164

Algorithm implementation 165
MUSCLEMOTION was then modified to handle typical experimental recordings by (i) improving the 166 signal-to-noise ratio (SNR), (ii) automating reference frame selection and (iii) programming built-in 167 checks to validate the generated output data (Fig. 1c). The SNR was increased by isolating the pixels 168 of interest in a three-step process: i) maximum projection of pixel intensity in the complete 169 displacement stack, ii) creation of a binary image of this maximum projection with a threshold level 170 equal to the mean grey value plus standard deviation and iii) multiplication of the pixel values in this 171 image by the original displacement and speed of the displacement image stack (Fig. S3). This process 172 allowed the algorithm to work on a region of interest with movement above the noise level only. 173 Next, a method was developed to identify the correct '() from the speed of displacement image 174 stack by comparing values obtained from the frame-to-frame calculation with their direct neighbouring 175 values, while also checking for the lowest absolute value (Fig. S4). 176 The reliability of MUSCLEMOTION for structures with complex movements was validated using a 177 custom-made contracting 3D "synthetic CM" model ( Fig. 1d,f,g) that was adapted to produce 178 contractions with known amplitude and duration. Linearity was preserved during the analysis of the 179 contraction and velocity; other output parameters of the analysis matched the input parameters ( Fig.  180   1e). A second 3D model (Fig. 1g), with a repetitive pattern aimed to create out-of-bounds problems 181 was also generated. As expected, contraction amplitude information here was not linear (Fig. 1e), 182 although contraction velocity and temporal parameters did remain linear (Fig. 1e,g). To mitigate this 183 problem, we implemented an option for a 10-sigma Gaussian blur filter that can be applied on demand 184 to biological samples that presented highly repetitive patterns (e.g. sarcomeres in adult CMs). shortening for adult CMs. Remarkably, standard methods currently used measure only contraction or 193 contraction velocity. Linearity was preserved in all cases during the analyses, demonstrating the 194 reliability of the results (Fig. S5). 195 First, single hPSC-CMs (Fig. 2a, Movie S2) exhibited concentric contraction (Fig. 2a ii) and 196 contraction velocity amplitudes correlated well with the amplitudes obtained by optical flow analysis 197 (R 2 = 0.916) (Fig. 2a v). In contrast to single cells, the area of displacement for hPSC-CM monolayers 198 was distributed heterogeneously throughout the whole field ( Fig. 2b ii, Movie S3). Optical flow 199 analysis was compared with our measure of velocity (Fig. 2b iv); this showed a good linear correlation 200 (R 2 = 0.803) (Fig. 2b v). Complex (mixed, multicellular) 3D configurations were also investigated by 201 analyzing hPSC-derived cardiac organoids 17 (Movie S4) and EHTs 14 (Movie S5). Cardiac organoids 202 showed moderate levels of displacement throughout the tissue (Fig. 2c ii), while the EHTs showed 203 high deflection throughout the bundle (Fig. 2d ii). The contraction velocity of the organoids correlated 204 well with the output of optical flow analysis (R 2 = 0.747, Fig. 2c v). Similarly, contraction amplitudes 205 in EHTs showed high linear correlation (R 2 = 0.819) with the absolute force values derived from 206 measurement of pole deflection (Fig. 2d v). Finally, single adult rabbit ventricular CMs were analyzed 207 ( Fig. 2e, Movie S6). Large displacements were evident around the long edges of the CM (Fig. 2e ii). 208 These cells were analyzed with a 10-sigma Gaussian blur filter, which also minimized (unwanted) 209 effects of transverse movements on contraction patterns. Linearity was preserved (Fig. S5) despite the 210 repetitive pattern of the sarcomeres and this resulted in accurate measures of both contraction ( Fig. 2e  211 iii) and speed of contraction (Fig. 2e iv). The contraction amplitude of the adult CMs stimulated at 1 212 Hz correlated well with the output of sarcomeric shortening using fast Fourier transform analysis 20 (R 2 213 = 0.871, Fig. 2e v). Thus, the MUSCLEMOTION algorithm yielded data in these initial studies 214 comparable with methods of analysis tailored for the individual platforms.

Application of MUSCLEMOTION to multiple imaging and recording platforms 216
To examine whether MUSCLEMOTION could potentially be used in applications that measure other 217 aspects of CMs functionality in parallel, we first determined the electrophysiological properties of 218 hPSC-CMs using patch clamp whilst recording their contractile properties through video imaging. 219 This allowed simultaneous quantitative measurement of action potentials (APs) and contraction ( Fig.  220   3a), for in-depth investigation of their interdependence. We observed a typical 21 profile of AP 221 followed by its delayed contraction. 222 To measure contractile force in combination with contractile velocity in single CMs, we integrated 223 fluorescent beads into polyacrylamide substrates patterned with gelatin ( Fig. 3b) Having shown that MUSCLEMOTION was fit-for-purpose in analyzing contraction over a variety of 232 platforms, we next sought to demonstrate its ability to detect the effects of positive and negative 233 inotropes. This is essential for ensuring the scalability of the tool over multiple platforms, particularly 234 in the context of hiPSC-CMs where regulatory authorities and pharmaceutical companies are 235 interested in using these cells as human heart models for drug discovery, target validation or safety 236 what has been reported 27 , contraction amplitude decreased at doses higher than 1 nM. NIFE treatment 242 decreased both contraction amplitude and duration starting from 3 nM, respectively (Fig. 4a). In paced 243 (1.5 Hz) hPSC-CMs monolayers, no significant effects were measured after addition of ISO on either 244 relaxation time or contraction amplitude. NIFE caused a progressive decrease in contraction duration 245 and amplitude in a concentration-dependent manner starting at 100 nM (Fig. 4b). Similarly, cardiac 246 organoids paced at 1.5 Hz showed no significant effects on both relaxation time and contraction 247 amplitude with ISO, while both parameters decreased after NIFE, starting from 100 nM and 300 nM, 248 respectively (Fig. 4c). EHTs paced at 1.5 times baseline frequency and analyzed with 249 MUSCLEMOTION showed a positive inotropic effect starting from 1 nM ISO and a negative 250 inotropic effect starting at 30 nM NIFE as previously reported 14 (Fig. 4d). 251 Paced (1 Hz) adult rabbit CMs exhibited no significant increase in relaxation time and contraction 252 amplitude at any ISO concentration. At concentrations higher than 3 nM, adult CMs exhibited after-253 contractions and triggered activity during diastole, which hampered their ability to be paced at a fixed 254 frequency. No significant effects were observed on contraction duration with NIFE, while contraction 255 amplitude significantly decreased in a dose-dependent manner starting from 100 nM (Fig. 4e). Data 256 generated by post deflection and sarcomere fractional shortening are available for comparison 257 purposes in Fig. S6. 258 259

Analysis of disease phenotypes in vivo 260
To extend analysis to hearts in vivo, we took advantage of the transparency of zebrafish, which allows 261 recording of contracting cardiac tissue in vivo (Fig. 5a, Movie S9). It was previously shown that 262 mutations in G protein β subunit 5 (GNB5) are associated with a multisystem syndrome in human, 263 with severe bradycardia at rest. Zebrafish with loss of function mutations in gnb5a and gnb5b were 264 generated. Consistent with the syndrome manifestation in patients, zebrafish gnb5a/gnb5b double 265 mutant embryos showed severe bradycardia in response to parasympathetic activation 23 . Irregularities 266 in heart rate were visually evident and were clearly distinguishable from the wild type counterpart 267 after analysis with MUSCLEMOTION (Fig. 5b). Quantification of the heart rate of these zebrafishes 268 with MUSCLEMOTION highly correlated (R 2 = 0.98) with the results of the published manual 269 analyses 23 (Fig. 5c). There was however, a striking time-saving for operators in carrying out the 270 analysis using the algorithm (5-10 times faster than manual analysis; 150 recordings were analysed in 271 5 hours versus 4 days) without compromising accuracy of the outcome. Qualitative analysis of 272 contraction patterns allowed rapid discrimination between arrhythmic vs non-arrhythmic responses to 273 carbachol treatment (Fig. 5c). 274 Finally, we examined human echocardiograms from five healthy and cardiomyopathic individuals 275 (Fig. 5d). To assess ventricular function, videos were cropped to exclude movement contributions of 276 the atria and valves. MUSCLEMOTION enabled rapid quantification of temporal parameters from 277 standard ultrasound echography (Fig. 5e) such as time-to-peak, relaxation time, RR interval and the 278 contraction duration (Fig. 5f). 279 280

281
A reliable and easy-to-use method to quantify cardiac muscle contraction would be of significant 282 benefit to many basic and clinical science laboratories to characterize cardiac disease phenotypes, 283 understand underlying disease mechanisms and predicting cardiotoxic effects of drugs 14,24 . 284 Quantification of frame-to-frame differences in pixel intensity has been used in recent reports with 285 success 10 ; however, the full spectrum of applications for which these algorithms are relevant, how their 286 output data correlates with gold standards in each system and software performance, specifications, 287 license and software availability, have remained unclear. 288 Here we developed and tested a user-friendly, inexpensive, open source software platform that serves 289 this purpose in a variety of biological systems of heart tissue. Its integration into current research 290 practices would benefit data sharing, reproducibility, comparison and translation in many clinically 291 relevant contexts 25 . 292 The linearity and reliability of MUSCLEMOTION were validated using a 3D reconstructed artificial 293 CM which gave the expected linear correlations between known inputs and the outputs (Fig. 1d-f). 294 When random repetitive patterns were applied, amplitude outputs differed from inputs, suggesting a 295 potential limitation to measuring contraction amplitudes in highly repetitive biological samples (such 296 as when sarcomere patterns are well-organized), while temporal parameters remained valid (Fig.  297   1d,e,g). However, conditions such as these would be unlikely in standard biological samples, where 298 camera noise significantly reduces the possibility of saturating pixel movement. We partially 299 attenuated this problem by applying, on user demand, a 10-sigma Gaussian blur filter which 300 significantly increased the accuracy of MUSCLEMOTION with highly repetitive structures. Also, to 301 increase reliability, we built in additional controls to detect any mismatches and errors. 302 MUSCLEMOTION can automatically identify and select the reference frame and increase the signal-303 to-noise-ratio, features which were particularly relevant in reducing user bias and interaction while 304 improving user experience. MUSCLEMOTION is valid in a wide range of illumination conditions 305 without changing temporal parameters; however, exposure time was linearly correlated with 306 contraction amplitude (Fig. S7). Batch mode analyses and data storage in custom folders were also 307 incorporated to support overnight automated analyses. For accurate quantification of amplitude, time-308 to-peak and relaxation time, an appropriate sampling rate should be chosen. For applications similar to 309 those described here, we recommend recording rates higher than 70 frames per second to sample 310 correctly the fast upstroke of the time-to-peak typical of cardiac tissue. This recording rate is easily 311 achievable even using smartphone slow motion video options (~120/240 frames per second), obviating 312 the need for dedicated cameras and recording equipment if necessary. 313 We demonstrated excellent linear correlations between our software tool and multiple other standard 314 methods independent of substrate, cell configuration and technology platform and showed that 315 MUSCLEMOTION is able to capture contraction in a wide range of in vivo and in vitro applications 316 ( Fig. 2 and Fig. 3). Specifically, we identified several advantages compared to optical flow algorithms 317 in terms of speed and the absence of arbitrary binning factors or thresholds which, when modified, 318 profoundly affect the results. One limitation compared to optical flow or EHT standard algorithm is 319 that the tool lacks qualitative vector orientation, making it more difficult to assess contraction 320 direction. Particularly important was the correlation with force data calculated from the displacement We proposed and validated practical application in pharmacological challenges using multiple 329 biological preparations recorded in different laboratories; this means that immediate use in multiple 330 independent high-throughput drug-screening pipelines is possible without further software 331 development being required, as recently applied for a drug screening protocol on cardiac organoids 332 from hPSCs 17 . Intuitively, the possibility of having inter-assay comparisons will also be of particular 333 relevance where comparisons of contraction data across multiple platforms are required by regulatory 334 agencies or consortia (e.g. CiPA, CSAHi) 5,6,22,27 . Moreover, this might offer a quantitative approach to 335 investigating how genetic or acquired diseases of the heart (e.g. cardiomyopathies 7 , Long QT 336 Syndrome 28 ), heart failure resulting from anticancer treatments 29,30 or maturation strategies 18,31,32 affect 337 cardiac contraction. The possibility of linking in vitro with in vivo assays, with low cost technologies 338 applicable with existing hardware certainly represents an advantage as demonstrated by automatic 339 quantification of zebrafish heartbeats and human echocardiograms (Fig. 5) state is similar to that of the relaxed (e.g. approaching sinusoidal). We have implemented a "fast 357 mode" option that captures reliable baseline values even at high contraction rates. Furthermore, 358 recordings must be free of moving objects (e.g. debris moved by flow, air bubbles) other than those of 359 interest. 360