MicroBundleCompute: Automated segmentation, tracking, and analysis of subdomain deformation in cardiac microbundles

Advancing human induced pluripotent stem cell derived cardiomyocyte (hiPSC-CM) technology will lead to significant progress ranging from disease modeling, to drug discovery, to regenerative tissue engineering. Yet, alongside these potential opportunities comes a critical challenge: attaining mature hiPSC-CM tissues. At present, there are multiple techniques to promote maturity of hiPSC-CMs including physical platforms and cell culture protocols. However, when it comes to making quantitative comparisons of functional behavior, there are limited options for reliably and reproducibly computing functional metrics that are suitable for direct cross-system comparison. In addition, the current standard functional metrics obtained from time-lapse images of cardiac microbundle contraction reported in the field (i.e., post forces, average tissue stress) do not take full advantage of the available information present in these data (i.e., full-field tissue displacements and strains). Thus, we present “MicroBundleCompute,” a computational framework for automatic quantification of morphology-based mechanical metrics from movies of cardiac microbundles. Briefly, this computational framework offers tools for automatic tissue segmentation, tracking, and analysis of brightfield and phase contrast movies of beating cardiac microbundles. It is straightforward to implement, runs without user intervention, requires minimal input parameter setting selection, and is computationally inexpensive. In this paper, we describe the methods underlying this computational framework, show the results of our extensive validation studies, and demonstrate the utility of exploring heterogeneous tissue deformations and strains as functional metrics. With this manuscript, we disseminate “MicroBundleCompute” as an open-source computational tool with the aim of making automated quantitative analysis of beating cardiac microbundles more accessible to the community.


Introduction
Despite significant recent advances in cardiovascular disease prevention and diagnosis [1,2], heart disease remains the leading cause of death among adults worldwide [3]. This is due, in part, to the fact that the native heart has a poor regenerative ability [4,5], and thus damage to the heart muscle during an adverse medical event such as a myocardial infarction is irreversible [6]. Cardiac tissue engineering is a promising approach to address this unmet societal need [5]. In particular, cardiac tissue engineering with human induced pluripotent stem cell derived cardiomyocyte (hiPSC-CM) based technology [7] is a promising approach to disease modeling [8,9,10,11], drug discovery [5,9,11,12,13], and regenerative tissue engineering [14,15,16]. However, the development of viable hiPSC-CM technology is very much ongoing. In particular, one major challenge is that differentiated hiPSCs initially resemble fetal cardiomyocytes -they are morphologically and functionally different compared to adult cardiomyocytes [12]. Thus, developing technology to promote the maturation of hiPSC-CMs and, likewise, hiPSC-CM based tissue is an active area of research [16]. One impactful approach to promoting the maturation of hiPSC-CMs is the use of engineered tissue culture platform designs in both two [17,18] and three dimensions [19,20,21,22,23] (Fig. 1). Across these different platforms, there are multiple physical [17,19,22,24], electrical [19,20,23,24,25], and chemical [26,27,28] knobs to tune to promote maturation and explore different physiological and pathological conditions. Even if we restrict our focus to microbundles (i.e., aligned, electromechanically coupled, microscale cardiac tissue bundles formed with natural hydrogel materials suspended between pillars), there is massive variability across different experimental setups [17,19,20,22].
Driven by this diversity in experimental approaches and the rapid growth of the field, the mechanical behavior of cardiac microbundles is challenging to compare across studies. Fundamentally, this challenge is driven by multiple factors, ranging from the high volume of data collected with these testbeds [29], to challenges associated with reproducing results when software and data are not shared under open-source licenses, or when extracting quantities of interest from data requires significant manual processing. To date, there have been multiple non-destructive image-based methods for quantifying the contractile action of cardiac microbundles [21,30,31,32,33,34,35,36,37,38,39], often inspired by related approaches to assessing the contractile behavior of cardiomyocytes [30,40,41,42,43,44,45,46,47,48,49,50]. Broadly speaking, most of these tools can be grouped into four main categories: (1) edge detection systems [31,32], (2) pillar tracking-based methods [33,34,35,36,37,39], (3) inter-frame pixel disparity methods [21,38], and (4) optical flow-based tracking [30,39]. Each one of these methods has benefits and limitations that suit specific platforms, conditions, and research questions (some of these approaches are elaborated on in Section 2). Despite this wide range of computational tools, many of which are available under open-source licenses [21,33,34,38,39,51], few options can compare across multiple experimental testbeds and function on new datasets out of the box. In addition, the most popular approach to assessing microbundle contractile behavior -pillar tracking -does not necessarily capture the full richness of mechanical behavior in these systems. realistic synthetic data, and then show the results of implementing "MicroBundleCompute" on multiple testbeds, shown in Fig. 1. We not only portray the efficacy of our approach, but also show that examining full-field tissue deformation consistently reveals heterogeneous contractile behavior throughout the domain. Looking forward, our goal is to use this work as a starting point to move beyond previously developed ad hoc approaches to analyzing these data [53], and establish a computational foundation for performing cardiac microbundle assessment and quality control at scale. Type 1 Type 2 Type 3 Figure 1: Examples of 3 microbundle experimental testbeds of different types. We briefly note here that "Type 1" and "Type 3" represent 3D culture platforms while "Type 2" portrays a significantly thinner (almost 2D) platform. We elaborate more on each type in Section 2.1.2 and Supplementary Document B.
The remainder of this document is organized as follows. In Section 2, we describe our methods to generate synthetic movies of beating cardiac microbundles with known ground truth to validate our computational pipeline, the diverse set of real examples on which we test the software, and finally, our code pipeline along with the main output metrics. Then in Section 3, we present the main findings of our code validation on one synthetic example and show sample outputs on a few real movies. Finally, in Section 4, we share our final thoughts regarding "MicroBundleCompute." Along with this document, we provide three supplementary documents where we explain in more details our validation pipeline (Supplementary Document A), share our complete set of real examples (Supplementary Document B), and finally provide more information on basic pillar tracking (Supplementary Document C), an additional feature of our computational tool. Overall, the intention of this manuscript is to outline our method and approach to making a tool of broad utility for the cardiac tissue engineering research community.
With regards to applications to cardiac microbundles in specific, available image-based techniques for computing the contractile behavior of microbundles, as categorized in Section 1 into 4 broad approaches, focus mainly on quantifying temporal profiles for the entire construct and as such, provide averaged results that lump the spatially heterogeneous tissue behavior into a single value per frame (3 of the 4 approaches). However, tools developed based on optical flow (4 th approach) can extract full-field, as well as directional outputs, such as directional displacement fields and strains.
To elaborate more, edge-detection approaches rely on quantifying the microbundle shape change between a relaxed (reference) state and a contracted (deformed) state. For example, in Ronaldson-Bouchard et al. [31], tissue contractility was calculated by tracking the change in the tissue area while previously, in Hansen et al. [32], the difference between the ends of the tissue was used to measure contraction. Yet, these tools were custom-developed and are not readily available online for the broader research community. As for pillar tracking-based methods, there are currently a number of available tools [33,34,39] out of the identified implementations that were developed and kept in-house [35,36,37]. In these approaches, the pillar or cantilever head deflection is estimated to generate contraction waveforms and extract contraction kinetics including contraction frequency, force, as well as the time to achieve 10%, 50%, or 90% of the peak contraction or relaxation.
Available tissue tracking methods via pixel intensity disparity have been implemented based on different approaches. For example, in "MUSCLEMOTION" [51], which is offered as an ImageJ [78] plugin, the absolute difference in pixel intensity between a reference frame and a frame of interest is calculated, whereas the MATLAB-based [79] "CardiacContractileMotion" [21] identifies the tissue region in a relaxed baseline state and tracks, within this region of interest, changes in pixel motion with time. Another MATLAB-based [79] tool, "ContractQuant" [38], which was specifically developed for implementation with micron-scale 2D cardiac muscle bundles, uses cross-correlation to track pixel features and find the best match for a specified region of interest across consecutive frames. Outputs from these software are in general similar to those extracted with pillar tracking approaches and include contraction and relaxation profiles and velocities, as well as contraction and relaxation times.
Overall, these approaches are suitable when only basic information, such as the mean value and the time rate or velocity of tissue shortening, shrink or contraction, are enough. However, tools developed based on optical flow can provide richer outputs. For example, Huebsch et al. [30] implemented block-matching optical flow [80] methods in MATLAB [79] to estimate absolute as well as directional full-field contractility. And very recently, a particle image velocimetry toolbox in MATLAB [79] was utilized in [39] to calculate displacement vectors of sub-regions or patches identified within a manually selected region of interest in the "reference" frame using a cross-correlation approach. And from these displacement fields, strain maps were subsequently derived. To the best of our knowledge, these methods appear to be robust and relatively versatile, yet at present they lack automation.
Within this scope, we present here our high-throughput optical flow-based computational framework to extract full-field deformation metrics from lab-grown cardiac microbundles. In the Sections that follow, we describe both the data used for testing and validating the "MicroBundleCompute" software, and the details of our computational methods. First, in Section 2.1, we introduce the two general categories of cardiac microbundle data that we have used in developing the software: synthetic data with a known ground truth and experimental data. Next, in Section 2.2, we explain the details of the code pipeline and describe core functionalities and output metrics.

Data
To validate our pipeline, we first invest significant effort into creating labeled data with a known ground truth. In Section 2.1.1, we outline our process for synthetic data generation. As a brief note, this is complemented with additional information in Supplementary Document A where we show comparisons to another form of synthetic data and manually labeled experimental data. In Section 2.1.2, we then show three distinct classes of microbundle experimental testbeds that we will use to showcase the function and versatility of our software.

Synthetic data generation
Here, we describe the steps to generate realistic synthetic brightfield movies of beating cardiac microbundles based on examples from "Type 1," as described in Section 2.1.2. In Fig. 2, we summarize the main steps of this pipeline. For each synthetic data example, we begin with a frame from an experimental movie (Fig. 2a). As the first step of this pipeline, we manually trace the tissue region in a relaxed valley frame, and use the traced region to obtain both tissue geometry and image texture. From the traced region, we extract the coordinates of the external contour to generate Finite Element (FE) simulation geometry, and isolate the tissue texture that will be warped following the FE simulation results (Fig. 2b).
To inform our FE simulations, we first generate a three-dimensional microbundle geometry based on the contour coordinates of a mask extracted from the single representative valley frame (Fig. 2c). Specifically, we extrude the 2D surface created by connecting these contour coordinates to a thickness of 400µm, a reasonable microbundle thickness given our target experimental setups [22]. To approximate the pillars, we implement the geometry and dimensions detailed in [22], which matches one of the main platforms used in our experimental dataset ("Type 1"). To create the FE mesh, we use Gmsh 4.10.5 [81], where the final mesh consists of 205, 524 tetrahedral elements which was deemed sufficient for our purpose, following a mesh refinement study. We provide a detailed schematic of the three-dimensional mesh geometry in Supplementary Document A.
In Fig. 2d, we briefly summarize the main components of the FE model as implemented in FEniCS 2019.1.0 [82,83]. Following popular recent work in the field of soft tissue biomechanics [84,85,86], we model the cardiac tissue as a nearly-incompressible transversely isotropic hyperelastic material where deformation is driven by periodic activation [84,86,87]. Of note, we create synthetic data with heterogeneous deformation fields (specifically, (1) a fully actively contracting tissue domain, (2) a passive circular inclusion at the center of the actively contracting tissue), and we vary the direction of contractile alignment (1) through the thickness or depth of the tissue or (2) along its length. We model  Figure 2: Illustrations of the main elements of the synthetic data generation pipeline in order of implementation: (a) an image of a microbundle movie frame, (b) a mask of the microbundle, extracted contour coordinates, and a segmented tissue texture; (c) a three-dimensional mesh for the FE model; (d) main variables within the FE simulations in order from left to right: profile of a linearly variable fiber direction with respect to depth, an illustration of the tissue depth direction, a uniform time series activation; (e) extracted surface displacement results and a warped image texture based on an estimated projective transformation. the poly(dimethylsiloxane) (PDMS) pillars as passive Neo-Hookean material, and treat the interface between the cardiac tissue and the pillars as perfectly bonded. For each FE simulation, we extract, with respect to time, the X, Y , and Z positions of the mesh cell centers at the top surface of the microbundle for each step and save the results as text files. Then, following the schematic illustration in Fig. 2e, we estimate a projective transformation based on the initial and deformed positions of the mesh cell centers and warp the image texture accordingly using the warp transform function in the scikit-image 0.19.3 Python library [88]. To enable a heterogeneous transformation, we subdivide the image domain and perform subdomain specific warping. Additional details are available in Supplementary Document A.
Overall, our main synthetic dataset consists of 60 generated movies of beating experimentally derived image textures. To obtain these 60 examples, we use 15 different base texture images extracted from 5 experimental movies of "Type 1" as described in 2.1.2. We then deform these extracted textures with FE results obtained from 4 different FE simulations run under the variable conditions specified above. To perform quantitative evaluation, we extract a 90 × 90 pixel region from each domain center and make all direct comparisons based on this domain.
For more information on implementing the Finite Element model, image warping, and the addition of Perlin noise [89], we refer the reader to Supplementary Document A. We also make the entire synthetic dataset prior to the addition of Perlin noise available along with all the Python code and files that are necessary to re-generate our dataset available on Github (https://github.com/HibaKob/SyntheticMicroBundle).

Experimental data
Our experimental dataset can be systematically categorized into 3 different types (Fig. 1). "Type 1" includes movies obtained from standard experimental microbundle strain gauge devices [90,91,92]. We refer to data collected from non-standard platforms termed FibroTUGs [93] as "Type 2" data. As for "Type 3," they represent data obtained from a highly versatile experimental platform [19,20] and as such, include the most diverse examples in this collection.
Specifically, "Type 1" examples were prepared as previously detailed in [53]. Briefly, PDMS (Dow Silicones Corporation, Midland, MI) microbundle devices were first cast from 3D printed molds (Protolabs). Each device contains 6 wells, each with two cylindrical pillars with rounded caps, where the cardiac microbundles are seeded. Up to 2 days before seeding, the devices were sequentially treated with 0.01% poly-l-lysine (ScienCell) and then with 0.1% glutaraldehyde (EMS) to promote cell attachment to the caps. On the day of seeding, devices were cleaned with 70% ethanol and ultraviolet (UV) sterilized. Next, the device wells were incubated with a small volume of 2% Pluronic F-127 (Sigma) to prevent cell attachment at the base of the well. hiPSC-CMs, differentiated and purified as described by Lian et al. [94], were seeded with human ventricular cardiac fibroblasts in a Matrigel (Corning) and fibrin (Sigma) extracellular matrix (ECM) solution. Microbundles were maintained in growth medium, with replacement every other day. Time-lapse videos of tissue contractions were acquired 5-7 days after seeding at 30 Hz using a 4× objective on a Nikon Eclipse Ti (Nikon Instruments Inc.) with an Evolve EMCCD camera (Photometrics), while maintaining a temperature of 37 • C and 5% CO 2 .
As for the second type, FibroTUG microbundles were fabricated as described previously [17,93]. First, arrays of PDMS cantilevers were fabricated by soft lithography as detailed in [93]. Then, fiber matrices, suspended between pairs of these cantilevers, were generated by selective photo-crosslinking of electrospun dextran vinyl sulfone (DVS) fiber matrices deposited onto the microfabricated PDMS cantilevers [17,93]. Matrix and cantilever stiffnesses were tuned by adjusting photoinitiator concentrations and cantilever height, respectively, while matrix alignment was controlled by altering the translation speed of collection substrates during fiber deposition [17,93]. Following functionalization of the electrospun fiber matrices with cell adhesive cRGD peptides, iPSC-CMs, differentiated and purified [17], were patterned onto matrices using microfabricated seeding masks cast from 3D-printed molds. Finally, time-lapse videos of the microtissue's spontaneous contractions were acquired at ∼ 65 Hz on Zeiss LSM800 equipped with an Axiocam 503 camera while maintaining a temperature of 37 • C and 5% CO 2 .
Examples from "Type 3" were generated using the protocol previously described in [19,20]. In brief, a combination of soft lithography and two-photon direct laser writing (DLW, Nanoscribe) was used to fabricate the seeding platforms. The process involves printing negative master molds using DLW, casting PDMS onto the molds, followed by sandwiching, curing and demolding. This results in 0.5 -0.6 mm-thick PDMS devices with embedded microfluidic channels and deformable seeding wells. As a final step, cage-like microstructures were printed using DLW on the sides of the wells of the demolded PDMS devices to facilitate cell attachment. After device fabrication, differentiated hiPSC-CMs as per the procedures described in [19], were seeded into the wells with human mesenchymal stem cells in a collagen ECM solution, with the growth medium changed every other day. Time-lapse videos of the tissue contractions were acquired 4-9 days after seeding at 30 Hz using 4× or 10× objectives on a Nikon Eclipse Ti (Nikon Instruments Inc.) with an Evolve EMCCD camera (Photometrics) equipped with a temperature and CO 2 equilibrated environmental chamber.
In total, we include in this framework 24 real experimental data, 11 examples from "Type 1," 7 from "Type 2," and 6 from "Type 3." This diverse pool of examples allows us to not only demonstrate the adaptability of our computational pipeline to different input examples, but also gain valuable insight about the heterogeneous contractile action of cardiac tissue by extracting and observing relevant mechanical metrics, such as full-field displacements, subdomain-averaged strains, and displacement and strain-derived outputs, as shown in Section 3.2 and Supplementary Document B. We note that details about each specific experimental example are provided in Supplementary Document B as well as on Dryad [95] (https://doi.org/10.5061/dryad.5x69p8d8g) where the whole dataset, "Microbundle Time-lapse Dataset", is made available.

Code
In this Section, we describe the main working components of our "MicroBundleCompute" software for the automatic analysis of deformation in brightfield and phase contrast movies of cardiac microbundles. Because our goal is to implement an approach with simplicity, versatility, and adaptability in mind, our pipeline is structured with four modular components: (1) image pre-processing and mask creation, (2) deformation tracking, (3) post-processing (e.g., rotation, interpolation, strain analysis), and (4) visualization. In Fig. 3, we provide a graphical summary of the major functionalities included in this pipeline and the computational workflow. As a brief note, the software GitHub repository (https://github.com/HibaKob/MicroBundleCompute) contains instructions on how to install and run the code, detailed explanations of each main function, and a more thorough description of the formatting of output files.
The implementation of these methods is divided across four Python files: create_tissue_mask, image_analysis, strain_analysis, and optional_preprocessing. The bulk of core functionality is contained in image_ analysis, and all functions of the code are designed to be modular when possible such that they can be replaced in the future if a need arises. In addition, many of the specific post-processing steps are optional, and adaptation to compute additional quantities of interest should be straightforward. As we describe our pipeline, we will specify which Python file a given function is located in. Essential functions included in our pipeline are as follows, following the order illustrated in Fig. 3.

Tracking region mask generation
Within create_tissue_mask, we provide multiple options for creating a binary mask of the tissue region (Fig. 3 I). At present, we provide 3 basic segmentation functions for microtissue mask creation: 1) a straightforward threshold-based mask, 2) a threshold-based mask that is applied after a Sobel filter, and 3) a threshold-based mask that is applied to either the minimum or the maximum (specified by the user as an input) of all movie frames. We implement our tissue segmentation pipeline mainly based on the threshold_otsu and sobel functions provided within the filters module in scikit-image 0.19.3 Python library [88].
In addition to programming alternative automated mask generation functions, software users can provide an externally generated tissue mask by including a file named "tissue_mask.txt" where the mask is a two-dimensional array in which the tissue domain is denoted by "1" and the background domain is denoted by "0" to allow for the analysis of domains that may fall outside the original scope of this endeavor.

Sparse optical flow algorithm for tracking
Within image_analysis, we provide all of the essential functions to automatically run the tracking algorithm illustrated in Fig. 3 II. Our tracking pipeline is built on OpenCV's [96] pyramidal implementation of the Lucas-Kanade sparse optical flow algorithm [52] where we identify markers as Shi-Tomasi "good features to track" corner points [97] that fall within the specified tissue mask. In brief, these corner points are points where a slight shift in location leads to "large" changes in intensity along both the horizontal and vertical axes. There are two key parameters required to tune OpenCV's goodFeaturesToTrack function: qualityLevel, a minimum score to measure if a feature can be tracked well, and minDistance, the minimum permitted Euclidean distance between two identified corners. We initialize minDistance = 3 and qualityLevel= 0.1. Then, we iteratively decrease qualityLevel until coverage > 40, where coverage is defined as the average number of pixels associated with each tracking point, for up to 15 iterations. As for minDistance, we define a local_coverage measure, which is the coverage computed on 20 × 20 pixel subdivisions. The minDistance is automatically incremented by 1 as long as the largest 3 local_coverage values are less than or equal to 50 and the number of iterations does not exceed 2. We note that we adjust qualityLevel and minDistance simultaneously.
For running OpenCV's [96] Lucas-Kanade optical flow [52] function (calcOpticalFlowPyrLK), we automatically tune the parameter winSize, which dictates the size of the integration window. Crucially, the window size in both horizontal and vertical directions, w x and w y , should be larger than the maximum tracked pixel motion between frames.
To specify winSize, we adopt a pragmatic approach where we initialize winSize = 5, perform a preliminary tracking step, calculate the maximum absolute displacement, and compare its magnitude to the initial window size. If the calculated displacement is larger, we increase winSize by 5, and continue to iterate until the condition is met or the number of iterations exceeds 15. Critically, keeping winSize from being larger than necessary reduces error during tracking. In the remainder of this Section, we will describe our methods for leveraging this basic sparse optical flow algorithm to effectively and automatically analyze microbundle domains.

Temporal segmentation
After automatic identification of features and tracking algorithm parameters, we run a preliminary tracking step. Representative results of preliminary tracking are shown in Fig. 3 II.b as a plot of mean absolute displacement vs. frame number. We use these preliminary tracking results to perform temporal segmentation where we delineate individual tissue beats. To accomplish this, we use the SciPy signal processing library function find_peaks [98]. The find_peaks input parameters, distance, the minimum horizontal distance between neighbouring peaks, and prominence, minimum prominence for a perturbation to be recognized as a peak, are identified automatically. Specifically, we initialize distance and prominence to values of 20 and 0.1 respectively. Then distance is updated to take a value equal to 1.5× the horizontal distance separating two consecutive intersection points between the time series and its mean. For the prominence parameter, we keep it constant at a value of 0.1, which we found to be suitable for all the example videos that we have analyzed to date. After peaks are identified, we define valleys as the midpoints of two consecutive peaks. By design, our computational framework requires a minimum of 3 beats to be present in any movie for automatic analysis as this is required for an accurate approximation of beat period. Temporal segmentation into individual beats is then performed based on the temporal location of each valley.
Beyond the determination of microbundle time-related properties (e.g., beat period), temporal segmentation is an essential part of our pipeline as we use it to work around the tracked feature drift observed over the duration of the movie (see Fig. 3 II.b for an illustration of drift).
After identifying individual beats, we split, based on the first and last beat frame numbers, the main two arrays storing column (horizontal) and row (vertical) locations of marker points obtained during preliminary tracking into multiple arrays corresponding to each segmented beat. Likewise, instead of frame 0 being the reference configuration for the whole tracking duration, the fiducial marker positions in the first frame of each beat become the baseline for all future output calculations within the beat.

Optional rotation and interpolation
After the tracking step is complete, we include two optional features for post-processing: sample rotation and fiducial marker displacement interpolation. First, we include an option to rotate both the images and the tracking results based on a specified center of rotation and desired horizontal axis vector. The center of rotation and desired horizontal axis vector can be either specified manually or identified automatically based on the geometry of the tissue mask. As a brief note, rotation is performed after tracking as the process involves interpolation which can lead to loss of image resolution. Also, we automatically rotate the tissue domain before performing strain subdomain calculations to match the global row (vertical) and column (horizontal) coordinate system. The second optional feature, interpolation, allows the user to interpolate the tracking results returned at the automatically identified fiducial marker points to user-specified locations, on a structured grid for example. This step can be performed after tracking and optional rotation, and will output the interpolated displacement fields to specified sampling points for either visualization or downstream analysis.

Subdomain spatial segmentation
In contrast to standard approaches to analyzing cardiac microbundles [51], our approach is unique in that we compute full-field quantities of interest over the tissue domain. In order to better analyze these full-field results and reliably post-process displacement fields to compute strain, we perform spatial segmentation to define tissue subdomains over which we can report average strain quantities [99,100]. This subdomain spatial segmentation is implemented in the strain_analysis file within "MicroBundleCompute." We provide two options to specify the subdomain extents defined as a rectangle: 1) automatic subdomain generation via clipping the input tissue mask or 2) manually providing subdomain extent coordinates. Given rectangular subdomain extents, we then delineate individual subdomain tiles by specifying either the target number of tiles in each column and row or by specifying the target tile dimensions in pixels.
Representative subdomain segmentation results are illustrated in Fig. 3 III.c.
gth. As for the second method, specified via tile_style = 2, it fixes the input values for the number of columns and rows and adjusts the side length accordingly.

Strain computation
With these defined subdomain regions, we then compute the average deformation gradient F avg of each subdomain and use it to compute relevant strain metrics. As stated previously, we compute average subdomain strain rather than full-field strain due to: (1) the desire to reduce the influence of imaging artifacts and noise, and (2) increased ease of comparison between samples. Within each subdomain, we define the standard continuum mechanics deformation gradient F as follows: where F maps a vector dX in the initial or reference configuration to its deformed configuration dx [101]. To apply this within the context of our tracking pipeline, we define a set of n vectors, Λ 0 , that connect each potential pair of fiducial markers that lie within the extents of each subdomain in the reference configuration. We define the reference configuration with respect to each cardiac tissue beat as the first frame in the segmented beat frame. Thus, Λ 0 is defined in frames that represent the most relaxed tissue state. We then compute Λ following the same structure as Λ 0 for each subsequent movie frame, where the updated fiducial marker positions capture the subdomain deformed configuration.
With this definition, we can set up the over-determined system of equations: where F avg is a 2 × 2 matrix, and Λ 0 and Λ are 2 × n matrices of vectors in the initial (reference) and current (deformed) configurations, respectively. Of note, when the initial and current frames are identical, Λ 0 = Λ, F avg = I, a 2 × 2 identity matrix. To solve this over-determined system, we can use the normal equation to find the best fit average deformation gradient as: We schematically illustrate our method to compute the mean deformation gradient in Fig. 3 III.d. With the computed F avg , finding the average Green-Lagrange strain tensor is straightforward: and we can then compute strain on a per subdomain per beat basis to obtain subdomain time series results for F avg and subsequently E avg .

Data structure preparation
Tracking all identified fiducial markers for an extended number of frames produces a large quantity of output data for each movie. Thus, we selectively save output results such that they are both comprehensible and easily accessible for downstream data analysis. First, we save information regarding the column (horizontal) and row (vertical) positions of the tracked marker points per beat. Specifically, we store one row-position text file and one col-position text file for each beat formatted as a M × N array where M is the number of markers that were tracked and N is the number of frames in the beat. Similarly, we output mean deformation gradient results as text files saving the column and row positions of the center for each subdomain and the 4 components of the 2 × 2 mean deformation gradient per subdomain per beat.

Data analysis, key metrics, and visualization
In Fig. 3 IV, we show key output metrics and their visualizations. In brief, we provide tools to visualize full-field displacement and average subdomain strain, and provide key quantities of interest such as maximum strain, beat period, and synchrony. In all cases, we build on the popular matplotlib package [102] for producing all visualizations. In addition to the metrics directly enabled by our novel pipeline, we provide other commonly pursued relevant outputs including beat frequency, beat mean amplitude, and tissue width at the domain center. To convert the numerical outputs from dimensions of pixels and frame numbers to physical units, the user can specify: 1) frames per second and 2) the length scale in units of µm/pixel.

Quality checking and rejection of unsuitable examples
In our extensive experience testing our code on different synthetic and real examples, we have identified three main instances that negatively influence the fidelity of our outputs and decrease our confidence in analysis results: 1. blurred input movie frames that prevent effective identification of corner points for the tracking described in Section 2.2.2.
2. movies that start from a contracted tissue position that confounds the temporal segmentation step described in Section 2.2.3.
3. movies where all displacement is on sub-pixel length scales lead to large relative errors given our choice of tracking algorithm.
The specific influence of these conditions are explored in Supplementary Document A. To address these conditions, we provide both warnings in the code when we detect that these conditions arise and functions to correct these scenarios (e.g., removing initial movie frames to address case 2).

Note on pillar tracking
In the broader cardiac microbundle and microtissue literature [33,39], tissue force and tissue stress are two commonly computed metrics. Broadly speaking, pillar force is computed by tracking pillar displacement and converting displacement to force via cantilever beam equations where the pillars are treated as beams with known elastic modulus and geometry [103]. Stress is then computed from pillar force by dividing force by tissue cross sectional area where tissue width is measured from the in-plane images [34] and tissue depth is assumed based on typical values observed via three dimensional imaging modalities [104]. We are able to adapt our framework to readily track pillar displacement by simply specifying pillar regions as the area to track rather than the tissue domain. Due to lack of novelty, we do not emphasize this functionality in this paper. However, we do provide additional details on this topic in Supplementary Document C.

Results and discussion
In this Section, we show a summary of our main results. In Section 3.1, we present representative validation examples from validating our pipeline on an example from the synthetic dataset with known ground truth behavior, and in Section 3.2, we share examples of implementing our computational framework on 3 different experimental platforms (i.e., 3 different data types as explained in Section 2.1.2). As a brief note, we provide a more comprehensive set of results in Supplementary Document A on validation, Supplementary Document B on additional real data examples, and Supplementary Document C on pillar tracking.

Synthetic data examples
Here, we briefly summarize the results of our validation studies. In Fig. 4, we show the performance of our pipeline on a single representative synthetic example ST_1. In Supplementary Document A, we provide a comprehensive summary of all 400 synthetic validation examples. In our validation studies, we compare "MicroBundleCompute" displacement and strain outputs against the known ground truth for synthetic data as originally generated (without noise) and for the same examples with added Perlin noise of different magnitude ratios and octaves. The code required to reproduce all synthetic data examples can be found on Github (https://github.com/HibaKob/SyntheticMicroBundle).
In Fig. 4a, we specify the synthetic data domain for which the validation studies were performed, a 90 × 90 domain warped with FEA-informed displacements as briefly described in Section 2.1.1 and presented in more details in Supplementary Document A. In Fig. 4b, we plot the mean absolute displacement with respect to frame number for beat 4, the beat with the maximum error at the peak displacement (full contraction) for the synthetic example with Perlin noise of magnitude ratio of 12% and octaves of 40. We specifically highlight these Perlin noise parameters because they mimic, to a great extent, noise artifacts found in the real experimental data. A direct comparison between both tracked and ground truth mean absolute displacements per beat revealed that the percentage error at peak displacement is less than 2.5% for the example shown in Fig. 4b, with the addition of Perlin noise slightly raising this error to 3.2%.
In addition, we assessed the ability of our computational framework to "predict" the mean absolute displacement by calculating the coefficient of determination (R 2 ) and found good agreement.  magnitude ratios and lower octaves), the code fails to produce meaningful outputs as indicated by missing data points in Supplementary Document A, Fig. A_24.
In Fig. 4 c-i-iv, we show the error on the column-column direction (i.e., horizontal) Green-Lagrange strain (E cc ) per subdomain for beat 4. We refer to this direction as the "column-column" direction to be consistent with the row (vertical) and column (horizontal) directions defined by our input images. Here, the reported R 2 values reveal that the minimum value of 0.943 occurs in subdomain B1. And, in general, these plots indicate that strain magnitudes tend to be underestimated. Furthermore, errors on E cc follow the observed trend for displacement errors, where synthetic examples of sub-pixel displacements exhibit higher errors with some cases having negative R 2 values (Supplementary Document A Table A_2 and Fig. A_8).
In Supplementary Document A, we share the validation results for all synthetic examples. We note that strain outputs are sensitive to subdomain divisions, specifically subdomain size. Determining a suitable subdomain size is a delicate process that is governed by two main opposing factors. The subdomain size should be large enough such that each subdivision contains an appropriate number of automatically identified fiducial markers to ensure that the computations are less sensitive to noise. On the contrary, the subdomain size should be small enough to avoid the loss or reduction of information due to averaging the heterogeneous deformation, or put more explicitly, obtaining attenuated or zero strain values due to lumping regions that are experiencing opposing deformations, for example extension versus compression.
In general, based on our comprehensive experience with implementing the code on a number of synthetic and real examples, we recommend a subdomain side length that is between 30 and 40 pixels. From a methodological perspective, we propose that the user observes the displacement field and avoids having subdomains that span regions where the change in displacement flips sign. In the future, we plan to investigate the approach adopted in [76] and [77] where strains are directly informed from affine warping functions optimized via the Lucas-Kanade algorithm [75] without computing displacement fields. The obtained results reveal that for super-pixel displacements, our software is quite robust to all tested cases with relatively low magnitudes and high octaves, or in other words, when the Perlin noise patterns resemble speckle noise rather than pronounced textures (see Supplementary Document A Fig. A_4). However, for examples with entirely sub-pixel displacements, the performance of the software degrades. In summary, "MicroBundleCompute" breaks down when the synthetic data has noise artifacts that appear similar to the original texture, and for sub-pixel displacements. Yet, given that real experimental examples of beating microbundles generally produce displacements that exceed a single pixel and that Perlin noise examples with higher octaves more faithfully represent naturally occurring noise in real microbundle data than lower values, we anticipate that "MicroBundleCompute" will output reliable mechanical metrics in real experimental settings on condition that the natural contrast of the microbundle textures is visibly present and in focus in the time-lapse videos.

Experimental data examples
We provide here a summary of implementing "MicroBundleCompute" on the experimental dataset described in Section 2.  [95] contains all raw videos in ".tif" format for the 24 experimental time-lapse images of beating cardiac microbundles, 23 of which are brightfield videos, while the remaining single example is a phase contrast video. Besides the raw videos and the experimental metadata describing the conditions under which they were obtained, we include the tissue mask used for each example, whether generated automatically via our computational pipeline or manually via tracing in ImageJ [78]. These time-lapse videos and masks were used to generate the results shown here and in Supplementary Document B. We note briefly that it is only possible to develop automatic mask segmentation functions for examples where there is imaging consistency and when we have an ample number of examples to identify a pattern. Future extensions of this framework will include automatic mask functions tailored to specific experimental needs.
In Fig. 5, we show visualizations of output results generated by running "MicroBundleCompute" on 3 experimental examples, each from a different data type. Of note, we only visualize 3 of many potential software outputs: (1) full-field absolute displacement, (2) spatially distributed subdomain-averaged Green-Lagrange strain E cc (automatically rotated to align with the column-column horizontal direction and plotted at beat 0 strain peak), and (3) Figure 5: We show sample outputs of the code when run on 3 examples from the 3 different data types, specifically Example 9 from "Type 1," Example 5 from "Type 2," and Example 1 from "Type 3" as shared on Dryad [95]. Note that the strain output for "Type 2" data is automatically rotated such that E cc aligns with the column-column direction. Additionally, the "Type 3" example shows an actuated microbundle at 0.5 Hz by applying sawtooth pressure waves with ∼ -6 kPa peak amplitude (equivalent to ∼ 2.5% strain) on the right side using a microfluidic pump (Elveflow OB1), and hence, the considerable discrepancy of the E cc time series plot from 0 at the end of beat 0.
To complement the results in Fig. 5 C_1). However, further investigations are required to develop a systematic way to assess the effect of this pillar-tissue interaction on microbundle contractility and analyze the extracted metrics in this context. Similar spatial heterogeneity is observed for the "Type 2" samples (Supplementary Document B Fig. B_2), where examples 1-5, which were prepared under the same experimental conditions (soft, aligned matrix), exhibit discrepant displacement distributions but fairly consistent magnitude ranges (mean absolute displacement values at the maximum contraction vary between 5.0 and 7.9 pixels or equivalently, 4.54 and 7.17 µm). For examples 6 and 7, which are prepared with stiff matrices, the results reveal lower contraction magnitudes where an aligned stiff matrix (example 6) shows higher contractions than the randomly distributed one (example 7) [93]. Of note, by examining all results of cardiac microbundle"Type 1" and "Type 2" data, a noticeable diagonal contraction pattern is apparent, especially in examples 2, 7, 9, and 10 from "Type 1" (Supplementary Document B Fig. B_1) and examples 2, 3, and 5 from "Type 2" (Supplementary Document B Fig. B_2). This observation, enabled by full-field tracking, indicates that future study to investigate the association between microbundle contraction, fiber alignment, and emergent load paths between the pillars, would be meaningful future work.
The diversity of "Type 3" experimental data prevents direct comparison between samples. However, the time series strain visualizations (Supplementary Document B Fig. B_3) reveal that cardiac microbundles grown on these experimental constructs typically contract more synchronously than microbundles of "Type 1" and "Type 2," where the strain time series plots in Supplementary Document B, Figs. B_1 and B_2 respectively, show aspects of temporal heterogeneity for which peak contractions do not occur at the same frame within all subdomains per beat. Finally, the visualized average subdomain E cc strains, as well as the remaining two strain components (row-row direction Green-Lagrange strain E rr and column-row direction Green-Lagrange strain E cr ) that are computed and saved but not included in the set of representative outputs that we visualize here, give insight about the regions within the microbundle that are contracting or bulging in a given direction across each beat.
As we mention in Section 2.2.10 and describe in more detail in Supplementary Document C, it is straightforward to implement our computational framework to track pillar displacements, adding to the versatility of the "MicroBundle-  Figure 6: Direct comparison of pillar tracking and tissue tracking: implementing the pillar tracking functionality of "MicroBundleCompute" on (a) Example 10 from "Type 1" data, (a-i) we obtain the pillar displacement based tissue force to describe the microbundle beating behavior, whereas full-field tissue tracking reveals the heterogeneous (a-ii) full displacement field as well as (a-iii) subdomain-averaged strains. In (b) we show the same outputs for Example 3 from "Type 2." We note that, for this example, we show the rotated displacement output to be consistent with the subdomain segmentation orientation. To view the original non-rotated displacement results, refer to Supplementary Document B, Fig. B_2.
Compute" software framework. In Fig. 6, we show, side by side, the outputs obtained via pillar tracking (Figs. 6a-i and b-i) and via tracking the entire tissue domain (Figs. 6a-ii-iii and b-ii-iii) for an example from "Type 1" and "Type 2" data respectively. Pillar tracking enables the calculation of an average absolute or directional pillar displacement based tissue force, which can be used to infer an average tissue stress given approximate tissue width and thickness as described in more detail in Supplementary Document C. On the other hand, tissue tracking reveals abundant information regarding the inherent heterogeneous nature of microbundle beating. For example, full-field displacement fields (Figs. 6a-ii and b-ii) show regions where maximum or minimum tissue contractions are taking place. Furthermore, while subdomain-averaged strains also underline the spatial heterogeneity of the tissue contractions, visualizing them with respect to time (or frame number) reveals the nature of temporal synchrony across the different subdivisions of the beating microbundle (Figs. 6a-iii and b-iii). We note that a similar comparison can be carried out for all "Type 1" and "Type 2" data for which pillar tracking results are included in Supplementary Document C, while tissue tracking results are shared in Supplementary Document B.
Finally, we include in Fig. 7 an example of "Type 3" data which clearly indicates that the microbundle is experiencing positive E cc strains at the center, specifically in subdomains A2 and A3 as indicated in the inset legend on the lower left corner of the strain time series plot. According to Wang et al. [105], this observation suggests that a necking instability is forming on this tissue over time, which is also supported by the thinning tissue width at the center. Within this context, being able to extract and examine strain distributions and magnitudes enables further investigations into understanding and determining the factors, such as pillar stiffness and ECM density [105], that produce microbundles that are more stable against necking.
Based on the results shown in this Section and in Supplementary Document B, conventional metrics comprising tissue force and tissue stress obtained via basic pillar tracking offer insightful yet lumped information about the microbundle behavior. Complemented with reliable and reproducible full-field data, such as displacement distributions and subdomain-averaged strains, as in [53] for example, these metrics become more useful to assessing the highly heterogeneous cardiac tissue behavior and understanding the complex underlying mechanisms driving this behavior. Furthermore, this spatial information would allow us to study injury models of tissue [53], as well as further investigate mechanical observations from pillar tracking, for example, the unbalanced pillar forces noted in Supplementary Document C, Figs. C_1 and C_2. Finally, with sufficient collected data, this approach would enable the development of appropriate statistical models of the whole tissue. These endeavors are planned as part of our future work.

Conclusion
In this work, we describe our approach to creating a computational tool for analyzing brightfield and phase contrast movies of beating microbundles. In brief, we describe our process for converting movies of beating microbundles to meaningful quantitative metrics, validating our approach against synthetically generated data, and testing it on a diverse pool of real experimental examples. In addition to ease of use, limited requirements for user intervention, relatively short run time, and no parameter tuning, "MicroBundleCompute" is easy to implement out of the box on new experimental datasets. Because it relies solely on the natural contrast of brightfield and/or phase contrast microbundle images and the resulting intensity gradients, it is also straightforward to integrate with existing experimental workflows.
To enable broad adoption, we share the software under open-source license and look forward to receiving feedback from different users on how to adapt the code to tailor to their specific needs and enhance the overall user experience.
Looking forward, we aim to constantly improve the code and benchmark it against available cardiac microbundle tissue analysis tools. In addition, we plan to continue developing our pipeline to address alternative quantitative metrics and different imaging modalities such as calcium imaging [106] and integrate these results with structural information such as sarcomere geometry and alignment [100,107]. From the results shown in Section 3.2 and Supplementary Documents B and C, it is also clear that there is significant variation across individual microbundle behavior both within and across testbeds. One key motivation for applying the "MicroBundleCompute" framework to these data moving forward is that it will make it possible to better understand and analyze this heterogeneity. Overall, our intention is for other researchers to directly benefit from disseminating this work. As a final note, here we demonstrate the utility and functionality of "MicroBundleCompute" for a particular highly used engineered cardiac tissue format: cardiac microbundles. In concurrent work [108], we leverage the fundamental core of this computational framework and make some minor modifications and extensions to extract mechanical metrics from actuated 2D muscle sheets. Looking forward, we will continue to generalize our framework to automatically analyze other engineered contractile tissue platforms to maximize the utility to the tissue engineering research community.

A Pipeline validation
In this Supplementary Document, we elaborate on Section 2.1.1: "Synthetic data generation" and provide further details on our approach to validate the output of our cardiac microbundle tracking software. In Section A.1, we describe the steps to generate Finite Element Analysis (FEA)-informed synthetic movies of beating microbundles with a known ground truth. We also provide details on our approach to collect ground truth data via manual tracking performed on real examples. Then in Section A.2, we present the results obtained by comparing the tracked output generated by the software and the ground truth data, and comment on the accuracy and robustness of our computational framework.

A.1.1 Synthetic dataset
In this Section, we describe additional details of our approach to creating realistic brightfield movies of beating cardiac microbundles, as depicted in Fig. 2 of the main document. We note that our synthetic dataset is primarily derived from "Type 1" examples, since "Type 1" cardiac microbundle data is a more common approach in the literature.
For the "Type 1" synthetic data, we first start by identifying frames from which we can extract tissue textures. We systematically choose frames that represent valley and peak positions of the beating microbundle (i.e., the most relaxed or most contracted states) in order to sample multiple distinct textures (Fig. 2a). We then manually trace the microbundle element in a valley frame. From the obtained mask, we extract the coordinates of the mask contour and crop out the background of the movie frame to obtain an image that contains just the tissue texture (Fig. 2b). For the Finite Element (FE) simulations, we generate a microbundle model based on the contour coordinates of a mask extracted from a single representative valley frame using Gmsh 4.10.5 [81] (Fig. 2c & Fig. A_1). By connecting the points representing the mask contour, we obtain a 2D tissue surface, which we then extrude to a representative thickness for our real "Type 1" data. To simulate pillars, we follow the design and dimensions detailed in [22], which corresponds to a commonly implemented platform design. The generated mesh consists of 205, 524 tetrahedral elements.
In Fig. 2d, we briefly summarize the main components of the FE model as implemented in FEniCS 2019.1.0 [82,83]. Specifically, we model the cardiac microbundle as a nearly-incompressible hyperelastic material, following recent work in the literature [84,85,86]. We describe the incompressible hyperelastic material with the strain energy density function Ψ written as: where F is the deformation gradient tensor, p is a Lagrange multiplier to impose the incompressibility constraint, and J = det(F) is the volume ratio and should be equal to unity for full incompressibility [101]. The scalar p acts as a reaction hydrostatic pressure resisting any volume change and can only be determined from the equilibrium equations and the boundary conditions.
To relax the incompressibility constraint, we adopt the standard volumetric-isochoric decomposition of the deformation gradient tensor to obtain a more robust mathematical formulation of the constitutive model equations and avoid numerical difficulties such as mesh locking that arise when finite element methods are implemented in the analysis of incompressible materials [84,109]. This representation can be derived following the multiplicative decomposition of the deformation gradient tensor into purely isochoric F iso and purely volumetric F vol components where F = F iso F vol . As such we can write F iso = J −1/3 I and F vol = J 1/3 I.
The two components of the deformation gradient tensor contribute in an additive manner to the strain energy density function, where we can write: The isochoric component of the strain energy density function can be defined as any hyperelastic model. Here, we consider a Neo-Hookean material model: where µ is the shear modulus and C iso = F T iso F iso = J −2/3 F T F with C being the right Cauchy-Green tensor. As for the volumetric strain energy, several forms are cited in the literature. In our simulation, we define Ψ vol = κ 2 ln(J) 2 where κ is the bulk modulus of the material, which is orders of magnitude higher than the nearly-incompressible material's shear modulus. Finally, for nearly-incompressible formulation, incompressibility is implemented by the constraint relating the Lagrange multiplier p and Ψ vol where p = −dΨ vol (J)/dJ.
In the numerical setting, we employ a mixed formulation finite element method with Taylor-Hood tetrahedral elements [110], where the displacement field is represented by a continuous piecewise quadratic Lagrange shape function (P 2) while the pressure field is represented by a continuous piecewise linear Lagrange shape function (P 1). Finally, to model the active properties of the cardiac tissue and its ability to actively contract and generate force without external loads, we implement the active strain approach [111]. In brief, this approach is based on a multiplicative decomposition of the deformation gradient tensor F into an active component F a and an elastic component F e [84,86]. For a transversely isotropic activation, F a is defined as: where γ is the activation in the fiber direction and f 0 is the fiber axis. The elastic part F e is defined as F e = FF −1 a . Applying this to the isochoric deformation gradient, we get (F iso ) e = F iso F −1 a . Finally, the strain energy density will be a function of F e only, and following our previous decomposition, the Neo-Hookean material model can be written as: where (C iso ) e = (F T iso ) e (F iso ) e . We define a uniform periodic time series activation ranging between 0 and a maximum activation value (see Fig. 2d lower right inset) along the fiber axis as: To introduce slight variability to the periodic function, we add correlated Brownian (red) noise [112]. The lower left and middle insets in Fig. 2d provide a general depiction of the defined fiber direction. In our simulations, we linearly vary the fiber angle α defined with the global X-axis (horizontal) from a lower value of 9.33 • to a maximum of 15. We also simulate spatially heterogeneous activation by including a circular passive region in the middle of the tissue where the activation is zero and increases gradually with distance from the periphery of the circle to reach the maximum defined activation value at the current step. We implement this by defining an inclusion ratio as: where r min is the radius of the circular inclusion, α = −80, β = 0.9, P is the center of the inclusion, p 1 , ..., p n are the tissue spatial mesh coordinates, ∥ ∥ denotes the Euclidean distance, and S α is the smooth maximum approximation function given by: where S α → min as α → −∞ (12) and calculate the heterogeneous activation as γ h = γ(1 − r).
As for the pillars, we define a passive compressible Neo-Hookean constitutive model as: where F is the deformation gradient, and µ and λ are the Lamé parameters equivalent to Young's modulus E and Poisson's ratio ν as E = µ (3λ + 2µ)/(λ + µ) and ν = λ/(2(λ + µ)). We specify a value of 0.47 for the Poisson's ratio, and a value of 1.6 MPa for the Young's modulus.
From these Finite Element simulations, we extract the X, Y , and Z, positions of the mesh cell centers at the top surface of the microbundle for each step. We then estimate a projective transformation based on the initial and deformed positions of the mesh cell centers and warp the image texture accordingly using the "warp" transform function in the scikit-image 0.19.3 Python library [88]. To deform the synthetic texture with a heterogeneous transformation, we divide the domain into smaller slightly overlapping subdomains and deform each one according to a projective transformation calculated specifically for the subdomain. We perform a convergence study to decide on an appropriate number of subdomains and choose a grid size that results in the lowest error between the deformation gradient mapping obtained from FE simulations F FEA and its equivalent that is obtained from approximating a projective transformation F projective (Fig. A_2). In our final implementation, we divide each synthetic texture domain into 16 × 16 subdomains.
Overall, our synthetic dataset consists of 60 generated movies of beating microbundle textures. We extract 90 × 90 pixel regions for the tissue textures from 3 different frames per movie (2 valleys and 1 peak) and employ for this purpose, 5 experimental movies from "Type 1." Thus, we obtain 3 × 5 = 15 different base texture images. We then deform these extracted textures with FE results from 4 different FE simulations: 1) homogeneous activation across the whole microbundle domain with the fiber direction varying linearly in the X direction, 2) homogeneous activation across the whole microbundle domain with the fiber direction varying linearly in the Z direction, 3) heterogeneous activation where the active microbundle domain has a passive inclusion in the middle with the fiber direction varying linearly in the X direction, and finally, 4) heterogeneous activation where the active microbundle domain has a passive inclusion in the middle with the fiber direction varying linearly in the Z direction. We briefly note that all the files required to reproduce our complete synthetic dataset of "Type 1" along with the dataset files are made available on Github (https://github.com/HibaKob/SyntheticMicroBundle).
As a final step, we add spatially correlated Perlin noise [89], to both make the dataset more realistic and to create more challenging examples to test the robustness of our computational framework. We add Perlin noise at varying levels of magnitude defined as a ratio of the maximum image intensity, and with a varying number of octaves as shown in Fig.  A_4. From these examples, we see that the addition of noise to the synthetic images leads to accrued drift that is similar to what we observe when tracking displacements in real data.
For code testing and validation, we select a subset of the entire synthetic dataset consisting of 16 original synthetic movies, with half of this testing data corresponding to homogeneous activation and the other half to heterogeneous activation, as depicted in Fig. A_3. For each of these 16 examples, we add Perlin noise of 5 different magnitude ratios and 5 different octaves (Fig. A_3). In total, we run our code on 400 test examples with known ground truth. We provide the results in Section A.2.

Heterogeneous activation
Name ST_3 ST_7 ST_16 ST_20 ST_51 ST_52 ST_55 ST_56  Figure A_3: A tabulated summary of the synthetic data along with the added Perlin noise parameters implemented for code validation. Here, "ST" refers to "synthetic texture." Finally, we include a single synthetic example with known ground truth based on "Type 2" data. The main steps to generate this example are similar to the ones followed when generating "Type 1" synthetic data; that is: 1) extracting the microbundle texture, 2) running a FE simulation mimicking the microbundle behavior, and finally, 3) warping the extracted texture with displacement results obtained from the FE simulations. These FE simulations are based on a tissue-specific model that has a detailed representation of fibers and cardiomyocytes as described in Jilberto et al. [113], which makes the kinematic behavior similar to that of real tissues. We provide this synthetic video with the supplementary material as "Synthetic_Type2.tif" and the performance of our framework on this example in Section A.2. is clear that the drift in valleys is negligible, (iii) noisy synthetic data generated with 32 octaves and 10% magnitude ratio for which the slope of the drift is a good match to real data; (iv) an example of excessive distortion to the time series caused by the addition of Perlin noise generated with 18 octaves and 8% magnitude ratio; (c) a quantitative comparison of the drift in valleys (measured by the base 10 logarithm of the slope of the drift) with respect to the octave and magnitude ratio. From this analysis, we observe that there are a range of Perlin noise parameters (highlighted with a yellow background) that lead to drift behavior that is consistent with the average observed in 5 real examples of "Type 1" data. From this investigation, we: 1) show that adding Perlin noise to our synthetic data leads to synthetic data that better recapitulates the challenges associated with the real experimental data, and 2) identify the characteristics of the Perlin noise that we will use in our validation dataset.

A.1.2 Manual tracking
Another approach to validate our computational pipeline is to track the displacements of manually selected points and compare the results to those obtained by our tracking software. Specifically, we track the positions of 30 points across a single beat. We perform the manual tracking on two different examples of "Type 2" data, where each example is tracked by 2 different users ( Fig. A_26 "Case 1" & "Case 2"). We include the results of this manual tracking validation in Section A.2.

A.2.1 Validation against synthetic data
Here, we provide the results of validating "MicroBundleCompute" against synthetic data of beating microbundles. For the generated data based on "Type 1" data without any added Perlin noise, we show the computed errors for both mean absolute displacement ( For synthetic examples that are based on homogeneous activation simulations, where the resulting microbundle contractions exceed a single pixel, the percentage errors at peak MAD fall between 2% and 7% for all examples (Table A_1). Overall, as shown in Table A_1, the R 2 values indicate good agreement between the tracked and ground truth mean absolute displacement profiles, with 0.992 being the lowest R 2 value at the maximum error. For sub-pixel displacements, however, the peak MAD error range rises to 5% -15% and the lowest R2 value drops to 0.939 (Table A_1). As such, we consider tracked sub-pixel displacements to be of lower fidelity and advise the user of "MicroBundleCompute" software to interpret results derived from sub-pixel displacements with caution.
For the E cc strain outputs, the computed errors are subdomain dependent. Overall, consistent with the behavior observed for MAD errors, synthetic examples based on homogeneous activation simulations exhibit significantly lower strain errors, with a maximum MAE of 9% of the peak ground truth E cc in the corresponding subdomain versus an extreme value of 98% for sub-pixel displacement simulations (Table A_2). As anticipated, the introduction of Perlin noise to "Type 1" synthetic dataset results in higher tracking errors as shown in Figs. A_23 and A_24 and summarized in Tables A_3 and A_4 for the homogeneous activation cases. Again, movies with sub-pixel displacements are adversely affected by noise addition to a point where the code fails to identify beats to track. Such cases are indicated by missing data points in Fig. A_24. Overall, the errors that arise from the addition of Perlin noise decrease as the number of octaves increases, and the magnitude ratio decreases. For homogeneously activated data, the maximum MAD remains less than 14.5% (Table A_3) for all tested cases, even those marked as extremely noisy compared to realistic data, as shown in Fig. A_4. However, this is not the case with heterogeneously activated examples where the MAD errors become unreasonably high in most cases (on the order of 10 3 − 10 4 ).
In general, for E cc , the MAE is below 36% but higher than 15% of the peak ground truth E cc for all synthetic cases corresponding to homogeneous activation with added Perlin noise, except for the extreme case where the errors exceed 100% as shown in Table A_4. For "Type 2" synthetic data, we plot the analytical displacements in both X (horizontal) and Y (vertical) against their tracked equivalents obtained by running our optical flow pipeline (Fig. A_25). The results reveal good agreement between the two with R 2 values of 0.998 and 0.984 for displacements in X and Y , respectively.  In Fig. A_26, we show the comparison between the displacements in the X (horizontal) and Y (vertical) directions obtained via manual tracking to those obtained by "MicroBundleCompute" at the highlighted marker points. In both cases, the R 2 values indicate a good correlation between manual tracking results and the results from our computational tracking pipeline. Notably, the difference in tracking results between users 1 and 2 not only provides a range of values for comparison to our pipeline, but also indicates the fallibility of the manual tracking and the need for automated tools that lead to reproducible results.

A.3 Final remarks
Here, we have presented our two approaches to software validation: 1) comparison against synthetic data with known ground truth, and 2) comparison against manually labeled data (Section A.1). From the results shown in Section A.2, we see that our computational pipeline produces fairly accurate displacement results and slightly less accurate  Figure A_6: Error in mean absolute displacement for synthetic data of "Type 1" based on FE simulations with heterogeneous activation.           Figure A_24: Results of our pipeline validation against noisy synthetic data of "Type 1" based on FE simulations with heterogeneous activation for different Perlin noise octaves and magnitude ratios (MR).
Warped frame X Y Figure A_25: Results of our pipeline validation against synthetic data of "Type 2".
Case 2: Figure A_26: Results of our pipeline validation against manually tracked points for 2 different cases of "Type 2" data.

B Additional examples
In this Supplementary Document, we (1) provide the details necessary to reproduce our analysis on all examples in the "Microbundle Time-lapse Dataset," and (2) show the results of running "MicroBundleCompute" on all of these data. In brief, the "Microbundle Time-lapse Dataset" contains 24 experimental time-lapse images of cardiac microbundles. The dataset is hosted under a CC0 open-source license on the Dryad Digital Repository [95]. Consistent with our description of the experimental dataset in Section 2.1.2: "Experimental data" in the main paper document, we categorize these data as "Type 1" (11 examples), "Type 2" (7 examples), and "Type 3" (6 examples). A brief metadata summary for each type is given in Table (Table B_1). These details are also elaborated on in Section 2.1.2: "Experimental data", and documented as metadata for the "Microbundle Time-lapse Dataset." In Tables B_2, B_3, and B_4, we provide the implementation details (mask type, first frame adjustment, and subdomain segmentation parameters) necessary to reproduce all results. For the subdomain segmentation parameters, we only include those that were changed from the default values provided within the "run_code.py" file located in the "MicroBundleCompute" GitHub repository https://github.com/HibaKob/MicroBundleCompute. Finally, in Figs. B_1, B_2, and B_3 we show representative outputs from running our code: full-field mean absolute displacement at the first tracked peak, subdomain averaged Green-Lagrange E cc (horizontal) strain at the first tracked peak, and time series plots of E cc strain for the first tracked beat. These results not only show the versatility of our framework, but also showcase the information that can be obtained from these data.    Table B_4: A summary of the example-specific information and subdomain segmentation parameters for each example of "Type 3" data. We note that example 1 (marked by an asterisk) is actuated by applying sawtooth pressure waves with ∼ -6 kPa peak amplitude (equivalent to ∼ 2.5% strain) using a microfluidic pump (Elveflow OB1) to stretch and release the tissue from one side at 0.5 Hz. All other examples are not subjected to any form of actuation. We also note that all masks for this data type were manually generated.  Figure B_1: Example outputs of "MicroBundleCompute" run on "Type 1" experimental data. Here we show, in order from left to right, full-field mean absolute displacement and subdomain-averaged Green-Lagrange E cc strain at the first tracked peak, as well as a time series plot of E cc strain for the first tracked beat.

Example_3
Example_4 Example_5 Example_6 Example_7 Figure B_2: Example outputs of "MicroBundleCompute" run on "Type 2" experimental data. Here we show, in order from left to right, full-field mean absolute displacement and subdomain-averaged Green-Lagrange E cc strain at the first tracked peak, as well as a time series plot of E cc strain for the first tracked beat.

Example_4
A PREPRINT

Example_5
Example_6 Figure B_3: Example outputs of "MicroBundleCompute" run on "Type 3" experimental data. Here we show, in order from left to right, full-field mean absolute displacement and subdomain-averaged Green-Lagrange E cc strain at the first tracked peak, as well as a time series plot of E cc strain for the first tracked beat.

C Pillar tracking
In this Supplementary Document, we elaborate on the basic pillar tracking feature that is included within our computational pipeline. In brief, we demonstrate how the functionality to track tissue deformation can be readily adapted and applied to compute the standard metrics of pillar displacement and subsequently cardiac microbundle twitch force. First, we describe our pipeline to extract pillar displacements, and then we describe our approach to convert these pillar displacements into twitch forces following standard equations from the literature (Section C.1) [103]. In Section C.2, we show the results obtained from running pillar tracking on our example data. Finally, in Section C.3, we briefly describe our future plans to advance our pillar tracking functionality beyond the scope of this project.

C.1 Methods
Here we describe our straightforward approach to tracking microbundle pillars and outputting twitch forces. Briefly, the pillar tracking functionality requires either manually or externally generating separate masks of each pillar, and follows a similar pipeline to the microbundle tissue tracking described in Section 2.2: "Code" in the main paper document with one main difference: temporal segmentation is by default skipped. Instead, we directly compute the mean position of all tracked points at every time step, and subsequently derive the mean absolute displacement relative to the first frame, considered to be a valley frame. However, we retain the temporal segmentation as an optional step to remove any drift present in the tracked results, as explained in Section 2.2.3: "Temporal segmentation". We briefly note that this approach is possible because we do not output full-field pillar displacement results.
Following standard approaches in the literature [103], the pillar directional and absolute forces can be computed from the obtained mean directional and absolute displacement or deflection results. Specifically, an approximation of the pillar force F can be found by applying Hooke's law: where the deflection δ refers to the mean tracked pillar displacement and the combined geometric and material pillar stiffness k is provided as an experimentally derived quantity.
Alternatively, if experimental testing data is not available, the poly(dimethylsiloxane) (PDMS) molded pillars can be approximated as cantilever beams [53,103], and the stiffness k can be determined by the pillar geometry and material properties as: where E is the material's elastic modulus, a is the location of force application, L is the cantilever length, and finally I is the moment of inertia defined as a function of the cantilever's geometry. For rectangular cross section beams, I = wt 3

12
where w is the pillar width and t is the pillar thickness while for circular cross section beams, I = πD 4 64 where D is the cylindrical pillar diameter. To compute mean tissue stress, pillar force is divided by tissue cross sectional area, where tissue cross sectional area is a function of tissue width computed from the tissue mask geometry, and experimentally measured tissue depth.
Experimental platforms that are based on pillar constructs (i.e., the "Type 1" and "Type 2" data introduced prior) are amenable to pillar tracking. We share our example implementation as well as general guidelines to use this feature on the "MicroBundleCompute" GitHub page (https://github.com/HibaKob/MicroBundleCompute). For our examples of "Type 1" data, we run the code with a pillar stiffness of k = 2.677µN/µm and length scale (ls) 4µm/pixel. For the "Type 2" examples, we specify k = 0.41µN/µm and ls = 0.908µm/pixel.

C.2 Results
Here, we present the results of implementing the pillar tracking pipeline with all 14 experimental examples of "Type 1" and "Type 2" data. Results are shown in Fig. C_1 and Fig. C_2 and reveal that the absolute force is not always consistent among the two pillars. While imaging artifacts, present in both types of data, might marginally contribute to the discrepancy in measured force outputs, we believe that there are a number of mechanistic reasons that we have yet to investigate, as we explain in Section C.3. Of note, "Type 1" data was run without temporal segmentation while "Type 2" data displayed non-trivial drift and required the temporal segmentation step. We attribute this difference to the out-of-plane motion that is more significantly visible in the substantially thinner "Type 2" tissues.  Figure C_1: Pillar absolute force (µN) obtained by running the pillar tracking pipeline with "Type 1" data. We note that, in these cases,"pillar 1" and "pillar 2" in the legends consistently refer to the left and right pillars, respectively. However, in general, "pillar 1" would refer to the user-defined pillar mask saved as "pillar_mask_1.txt" and "pillar 2" to the pillar whose mask is defined by "pillar_mask_2.txt".

C.3 Future Work
In this document, we have provided a description of the pillar tracking functionality included within the "MicroBundle-Compute" computational framework and have shared representative examples of implementing this functionality. We believe that, in its current state, the pillar tracking option would be useful to a number of researchers, given its ease of use. However, we have multiple plans to further improve our pipeline including: • automatic pillar mask extraction • validated outputs against synthetic data and against other tools currently available in the literature [33] • generalized pillar tracking functionality that can be implemented on experimental setups beyond the scope of the "Type 1" and "Type 2" data examples shown here In addition, by examining the time series plots of the pillar absolute force for both pillars in Fig. C_1 and Fig. C_2, we can see that the deflections of two pillars contained within one experimental sample do not always match. There are multiple potential mechanistic explanations of this observation. For example, the tissue may be at different vertical positions on each pillar, pillars may have variable material or geometric properties, or there may be an asymmetry in tissue-pillar interaction that we do not currently understand. Further investigations are needed to both formulate solid explanations of this finding, and develop the appropriate statistical tools to address the associated uncertainty as a companion to our software.  Figure C_2: Pillar absolute force (µN) obtained by running the pillar tracking pipeline with "Type 2" data. We note that, in these cases, "pillar 1" refers to the user-defined pillar mask saved as "pillar_mask_1.txt" and "pillar 2" to the pillar whose mask is defined by "pillar_mask_2.txt".