Pipeline_Validation Additional_Examples Pillar_Tracking 11affiliationtext:
Department of Mechanical Engineering, Boston University, Boston, MA 02215, USA
22affiliationtext: Center for Multiscale and Translational Mechanobiology,
Boston University, Boston, MA 02215, USA 33affiliationtext: Department of
Biomedical Engineering, University of Michigan, Ann Arbor, MI 48109, USA
44affiliationtext: Photonics Center, Boston University, Boston, MA 02215, USA
55affiliationtext: Department of Biomedical Engineering, Boston University,
Boston, MA 02215, USA 66affiliationtext: Harvard-MIT Program in Health
Sciences and Technology, Institute for Medical Engineering and Science,
Massachusetts Institute of Technology, Cambridge, MA 02139, USA
77affiliationtext: Wyss Institute for Biologically Inspired Engineering,
Harvard University, Boston, MA 02215, USA 88affiliationtext: Department of
Electrical and Computer Engineering, New York Institute of Technology, New
York, NY 10023, USA 99affiliationtext: Department of Cardiac Surgery,
University of Michigan, MI 48109, USA 1010affiliationtext: Department of
Biomedical Engineering, School of Imaging Sciences and Biomedical Engineering,
King’s College London, King’s Health Partners, London SE1 7EH, United Kingdom

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

Hiba Kobeissi  Javiera Jilberto  M. Çağatay Karakan  Xining Gao  Samuel J.
DePalma  Shoshana L. Das  Lani Quach  Jonathan Urquia  Brendon M. Baker
Christopher S. Chen  David Nordsletten  Emma Lejeune

###### Abstract

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,
requires little to no parameter tuning, and runs quickly on a personal
computer. 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.

_Keywords_ cardiac tissue engineering $\cdot$ mechanobiology $\cdot$ open
science $\cdot$ microbundle $\cdot$ hiPSC-CMs $\cdot$ bright field imaging
$\cdot$ phase contrast microscopy

00footnotetext: * Corresponding author<EMAIL_ADDRESS>
## Author Summary

It has been well established that human cardiac tissue is not capable of
regenerating after an injury. For example, after a heart attack, the damaged
cardiac muscle is replaced by scar tissue that is not able to contribute to
the heart’s contractility and blood pumping action, subjecting the intact
heart muscle to strenuous loads. As such, this heart damage is not only
irreversible, but can also become worse with time. One promising approach to
help address this societal challenge is cardiac tissue engineering.
Specifically, researchers are working on human induced pluripotent stem cell
derived cardiomyocytes (hiPSC-CMs) technology as an in vitro platform to both
understand causes of cardiac disease, and develop viable pharmaceutical and
tissue engineering-based treatments. In this article, we introduce a
versatile, validated, and open-source software for analyzing the contractile
functional behavior of in vitro cardiac tissue. In addition to describing our
methodology, we showcase the results of applying our method across several
different cardiac tissue engineering testbeds. Looking forward, we anticipate
that this tool will directly enable advances in cardiac tissue engineering by
making it easier to quantitatively assess the outcomes of laboratory
experiments. 00footnotetext: Email addresses<EMAIL_ADDRESS>(Hiba Kobeissi),
<EMAIL_ADDRESS>(Javiera Jilberto<EMAIL_ADDRESS>(M. Çağatay Karakan),
<EMAIL_ADDRESS>(Xining Gao<EMAIL_ADDRESS>(Samuel J. DePalma<EMAIL_ADDRESS>(Shoshana L. Das<EMAIL_ADDRESS>(Lani Quach<EMAIL_ADDRESS>(Jonathan Urquia<EMAIL_ADDRESS>(Brendon M. Baker<EMAIL_ADDRESS>(Christopher S. Chen<EMAIL_ADDRESS>(David Nordsletten<EMAIL_ADDRESS>(Emma Lejeune)

##  1 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.

Given this research landscape, there is a clear need for an open-source
computational tool to extract functional metrics from time-lapse images of
microbundle contraction. To address this need, we have developed the
“MicroBundleCompute” computational framework that we will present in this
paper. In brief, “MicroBundleCompute” is disseminated as a Python package and
is based on the Lucas-Kanade optical flow algorithm [52] for computing full-
field displacements, subdomain-averaged strains, and displacement and strain
derived metrics. In this manuscript, we introduce the methods underlying the
“MicroBundleCompute” framework, discuss our approach to validating the
pipeline primarily via 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.

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.

##  2 Materials and methods

In a broader sense, capturing full-field soft tissue deformation is critical
for a number of research applications ranging from inverse material
characterization [54, 55, 56, 57], to high-fidelity biomechanical modeling of
in vivo mechanisms [58, 59, 60], and patient-specific modeling and procedure
planning [61, 62]. And recently, there has been keen interest in implementing
image registration-based techniques widely used in the field of computer
vision [54, 63] for full-field measurements, including digital image
correlation (DIC) [64, 65, 66, 67, 68], 3D-DIC [69, 70], digital volume
correlation (DVC) [71, 72, 73, 74], and optical flow algorithms [75, 76, 77].

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.

###  2.1 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.

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

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.

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\mu$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 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\times 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).

####  2.1.2 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\times$ objective on a Nikon Eclipse Ti (Nikon Instruments Inc.)
with an Evolve EMCCD camera (Photometrics), while maintaining a temperature of
$37^{\circ}$C and $5\%$ CO2.

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 $\sim$ $65$ Hz on Zeiss LSM800 equipped with an
Axiocam 503 camera while maintaining a temperature of $37^{\circ}$C and $5\%$
CO2.

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\times$ or $10\times$
objectives on a Nikon Eclipse Ti (Nikon Instruments Inc.) with an Evolve EMCCD
camera (Photometrics) equipped with a temperature and CO2 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.

###  2.2 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.

Figure 3: Schematic illustration of the “MicroBundleCompute” computational
framework. Two main inputs are required: I) a binary mask generated either
automatically using the software or manually and consecutive movie frames of
the beating microbundle. II) For tracking, (a) marker points are identified on
the first frame and (b) tracked across all frames to identify individual beats
and perform time segmentation. (c) This allows us to perform the analysis per
beat and correct for the observed drift as discussed in Section 2.2.3. (d)
Finally, we save the row and column positions of the tracked markers per a
single beat and use these saved outputs to compute full-field displacements
and derive strain results. III) Post-processing functionalities include (a)
rotation of the images and tracking results, (b) interpolation of the results
at query points, and (c) segmentation into subdomains for which (d) average
deformation gradients and subsequently, subdomain-averaged strains are
calculated. IV) Finally, to visualize the results, the software outputs (a)
time series plots per beat and movies of (b) full-field and (c) subdomain-
averaged results.

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.

####  2.2.1 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.

####  2.2.2 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\times 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.

####  2.2.3 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\times$ 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.

####  2.2.4 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.

####  2.2.5 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.

####  2.2.6 Strain computation

With these defined subdomain regions, we then compute the average deformation
gradient $\mathbf{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 $\mathbf{F}$ as follows:

$\displaystyle\mathbf{F}d\mathbf{X}=d\mathbf{x}$ (1)

where $\mathbf{F}$ maps a vector $d\mathbf{X}$ in the initial or reference
configuration to its deformed configuration $d\mathbf{x}$ [101]. To apply this
within the context of our tracking pipeline, we define a set of $n$ vectors,
$\mathbf{\Lambda_{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, $\mathbf{\Lambda_{0}}$
is defined in frames that represent the most relaxed tissue state. We then
compute $\mathbf{\Lambda}$ following the same structure as
$\mathbf{\Lambda_{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:

$\displaystyle\mathbf{F_{avg}}\mathbf{\Lambda_{0}}=\mathbf{\Lambda}\quad\text{where}\quad\mathbf{\Lambda_{0}}=[\mathbf{v_{01}},\mathbf{v_{02}},...,\mathbf{v_{0n}}]\quad\text{and}\quad\mathbf{\Lambda}=[\mathbf{v_{1}},\mathbf{v_{2}},...,\mathbf{v_{n}}]$
(2)

where $\mathbf{F_{avg}}$ is a $2\times 2$ matrix, and $\mathbf{\Lambda_{0}}$
and $\mathbf{\Lambda}$ are $2\times n$ matrices of vectors in the initial
(reference) and current (deformed) configurations, respectively. Of note, when
the initial and current frames are identical, $\mathbf{\Lambda_{0}}$ =
$\mathbf{\Lambda}$, $\mathbf{F_{avg}}$ = $\mathbf{I}$, a $2\times 2$ identity
matrix. To solve this over-determined system, we can use the normal equation
to find the best fit average deformation gradient as:

$\displaystyle\mathbf{F_{avg}}=\mathbf{\Lambda}\mathbf{\Lambda^{T}_{0}}\left[\mathbf{\Lambda_{0}}\mathbf{\Lambda^{T}_{0}}\right]^{-1}\,.$
(3)

We schematically illustrate our method to compute the mean deformation
gradient in Fig. 3 III.d. With the computed $\mathbf{F_{avg}}$, finding the
average Green-Lagrange strain tensor is straightforward:

$\displaystyle\mathbf{E_{avg}}=\frac{1}{2}\left(\mathbf{C_{avg}}-\mathbf{I}\right)\quad\text{where}\quad\mathbf{C_{avg}}=\mathbf{F^{T}_{avg}}\mathbf{F_{avg}}$
(4)

and we can then compute strain on a per subdomain per beat basis to obtain
subdomain time series results for $\mathbf{F_{avg}}$ and subsequently
$\mathbf{E_{avg}}$.

####  2.2.7 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\times 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\times 2$ mean deformation gradient per subdomain
per beat.

####  2.2.8 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 $\mu$m/pixel.

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

####  2.2.10 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.

##  3 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.

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

Figure 4: Validation results against synthetic example ST_1 in its original
(no noise) state as well as after the addition of Perlin noise (magnitude
ratio of $12\%$ and octaves of $40$): (a) a schematic representation of the
synthetic data domain and subdomain divisions for strain calculations, (b)
tracked and FEA-extracted (ground truth) mean absolute displacement of the
full-field data for a single beat (beat $4$), where maximum error for
synthetic data with noise is observed, (c-i)–(c-iv) tracked and FEA-extracted
(ground truth) subdomain-averaged $E_{cc}$ strain for beat $4$ for each
subdomain. Overall, the $R^{2}$ values indicate relatively good agreement
between the tracking output and ground truth data. For more information on
validation with synthetic data, refer to Supplementary Document A.

In Fig. 4a, we specify the synthetic data domain for which the validation
studies were performed, a $90\times 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. And overall, considering all validation examples in Supplementary
Document A prior to the addition of Perlin noise, peak mean displacement
errors were found to be less than $7\%$ for super-pixel displacements
(Supplementary Document A Table A_1 and Fig. A_5) whereas sub-pixel
displacements had higher errors that reach up to $15\%$ (Supplementary
Document A Table A_1 and Fig. A_6) with a minimum $R^{2}$ value of $0.939$.
Furthermore, it is evident that the addition of Perlin noise produces higher
shifts within the tracked outputs, with higher errors reported for higher
magnitude ratios and lower octaves (Supplementary Document A Figs. A_23 and
A_24). Specifically, as shown in Supplementary Document A, Table A_3 and Fig.
A_23, the percentage error at the peak mean absolute displacement increases to
around $14.4\%$ for super-pixel displacement examples. For the sub-pixel
displacement examples, mean absolute displacement errors become unreasonably
high in the majority of the synthetic examples and are on the order of
$10^{3}-10^{4}$ (Supplementary Document A Fig. A_24). And, in some of these
sub-pixel displacement examples where the added Perlin noise is
unrealistically extreme (higher 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.

As described in this Section and in Supplementary Document A, we have
performed extensive validation studies for the “MicroBundleCompute”
computational framework against a total of $400$ synthetic examples ($16$
examples generated under baseline conditions and the remaining $384$ created
by adversely altering the original $16$ examples via the addition of Perlin
noise) with known ground truth. These validation studies: 1) corroborate the
output of the software against FEA-labelled data, 2) test its performance
against realistic synthetic examples, and 3) evaluate its robustness against
challenging examples with excessive noise artifacts. 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.

###  3.2 Experimental data examples

We provide here a summary of implementing “MicroBundleCompute” on the
experimental dataset described in Section 2.1.2 and Supplementary Document B,
and shared on Dryad (https://doi.org/10.5061/dryad.5x69p8d8g) under a Creative
Commons CC0 $1.0$ Universal Public Domain Dedication. Specifically, the
“Microbundle Time-lapse Dataset” [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) time series plot of $E_{cc}$ strains with respect to beat $0$ frames.
These are representative examples from the comprehensive list of outputs
described in the software GitHub repository
(https://github.com/HibaKob/MicroBundleCompute). While displacement and strain
visualizations reveal spatial information on tissue heterogeneous behavior and
spatial contraction patterns, time series strain plots highlight the synchrony
(or lack of synchrony) of subdomain contraction across each beat.

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 $\sim$ \- 6 kPa peak amplitude (equivalent to $\sim$ 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, we show, in Supplementary Document B,
representative results obtained for all $24$ real examples following the same
format. Critically, all of the results presented in the Supplementary
Document, for all $3$ of the data types, were obtained without the need for
any parameter tuning. As is evident in Supplementary Document B, Fig. B_1, we
can group data from “Type 1” into $2$ different categories: 1) examples
$1$–$6$ which constitute relatively challenging examples to the code since the
imaging brightness/contrast is not optimized to fully accentuate the natural
tissue texture, and 2) examples $7$–$11$ where a texture is clearly visible in
the microbundle region, making them a set of optimal examples for running the
code. Despite the fact that “Type 1” data comprises the most consistent and
well studied experimental platform, the full-field absolute displacement plots
reveal wide variation in spatially distributed behavior. However, magnitude
ranges are relatively consistent across samples with mean absolute
displacement values at the peaks varying between $1.2$ and $2.3$ pixels ($4.8$
and $9.2$ $\mu$m). We anticipate that this observation might be attributed to
variations in pillar-tissue attachment behavior across different examples (see
Supplementary Document C Fig. 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$ $\mu$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.

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.

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
“MicroBundleCompute” 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.

Figure 7: Outputs of “MicroBundleCompute” run on Example $6$ from “Type 3”
reveal regions of tissue extension at the center. We briefly note that for
displacement outputs, rotating the results is optional, while for strain
outputs, the frames are rotated such that the major axis of the microbundle
automatically aligns with the column (i.e., horizontal) axis. To view the
original orientation of this example, refer to Supplementary Document B, Fig.
B_3.

##  4 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.

##  Appendix A Pipeline validation

aa_template addtofilelistaa_template.tex addtofilelistaa_template.aux

Supplementary Document 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 Methods

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

Figure A_1:  Schematic representation of the “Type 1” microbundle mesh
implemented in our Finite Element simulations.

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 $\Psi$ written
as:

$\displaystyle\Psi=\Psi(\mathbf{F})-p(J-1)$ (5)

where $\mathbf{F}$ is the deformation gradient tensor, $p$ is a Lagrange
multiplier to impose the incompressibility constraint, and $J=det(\mathbf{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 $\mathbf{F_{iso}}$
and purely volumetric $\mathbf{F_{vol}}$ components where
$\mathbf{F}=\mathbf{F_{iso}}\mathbf{F_{vol}}$. As such we can write
$\mathbf{F_{iso}}=J^{-1/3}\mathbf{I}$ and
$\mathbf{F_{vol}}=J^{1/3}\mathbf{I}$.

The two components of the deformation gradient tensor contribute in an
additive manner to the strain energy density function, where we can write:

$\displaystyle\Psi(\mathbf{F})=\Psi_{iso}(\mathbf{F_{iso}})+\Psi_{vol}(\mathbf{F_{vol}})$
(6)

The isochoric component of the strain energy density function can be defined
as any hyperelastic model. Here, we consider a Neo-Hookean material model:

$\displaystyle\Psi_{iso}=\frac{1}{2}\mu\big{[}tr\left(\mathbf{C_{iso}}\right)-3\big{]}$
(7)

where $\mu$ is the shear modulus and
$\mathbf{C_{iso}}=\mathbf{F^{T}_{iso}}\mathbf{F_{iso}}=J^{-2/3}\mathbf{F^{T}}\mathbf{F}$
with $\mathbf{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 $\Psi_{vol}=\frac{\kappa}{2}ln(J)^{2}$ where $\kappa$ 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 $\Psi_{vol}$ where $p=-d\Psi_{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 ($P2$)
while the pressure field is represented by a continuous piecewise linear
Lagrange shape function ($P1$).

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 $\mathbf{F}$
into an active component $\mathbf{F_{a}}$ and an elastic component
$\mathbf{F_{e}}$ [84, 86]. For a transversely isotropic activation,
$\mathbf{F_{a}}$ is defined as:

$\displaystyle\
\mathbf{F_{a}}=\left(1-\gamma\right)\mathbf{f_{0}}\otimes\mathbf{f_{0}}+\frac{1}{\sqrt{1-\gamma}}\left(\mathbf{I}-\mathbf{f_{0}}\otimes\mathbf{f_{0}}\right)$
(8)

where $\gamma$ is the activation in the fiber direction and $\mathbf{f_{0}}$
is the fiber axis. The elastic part $\mathbf{F_{e}}$ is defined as
$\mathbf{F_{e}}=\mathbf{F}\mathbf{F^{-1}_{a}}$. Applying this to the isochoric
deformation gradient, we get
$\mathbf{(F_{iso})_{e}}=\mathbf{F_{iso}}\mathbf{F^{-1}_{a}}$. Finally, the
strain energy density will be a function of $\mathbf{F_{e}}$ only, and
following our previous decomposition, the Neo-Hookean material model can be
written as:

$\displaystyle\Psi_{iso}=\frac{1}{2}\mu\left[tr\left(\mathbf{C_{iso}}\right)_{\mathbf{e}}-3\right]$
(9)

where
$\left(\mathbf{C_{iso}}\right)_{\mathbf{e}}=\mathbf{(F^{T}_{iso})_{e}}\mathbf{(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:

$\gamma=\begin{cases}-A\cos(\omega t)\qquad&\gamma\geq 0.008\\\
0&\gamma<0.008\end{cases}\ \qquad\text{with}\quad 0\leq\omega\leq 2\pi$ (10)

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 $\alpha$ defined with the global
$X$-axis (horizontal) from a lower value of $9.33^{\circ}$ to a maximum of
$15.33^{\circ}$ in two different tissue dimensions: 1) with respect to tissue
depth (global $Z$-axis, into the page), where the fiber angle is
$9.33^{\circ}$ on the top surface of the tissue ($Z=0$) and becomes
$15.33^{\circ}$ at the tissue’s bottom ($Z=-400\mu$m) and 2) with respect to
tissue length (global $X$-axis, horizontal), where the fiber angle is
$9.33^{\circ}$ on the left side of the tissue ($X=0$) and becomes
$15.33^{\circ}$ at the tissue’s right end.

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:

$r=\begin{cases}1\qquad&d<r_{min}\\\ \frac{r^{2}_{min}}{\beta
r^{2}_{min}+(1-\beta)d}\qquad&d\geq r_{min}\end{cases}\ \qquad\text{with}\quad
d=S_{\alpha}(\|P-p_{1}\|,...,\|P-p_{n}\|)$ (11)

where $r_{min}$ is the radius of the circular inclusion, $\alpha=-80$,
$\beta=0.9$, $P$ is the center of the inclusion, $p_{1},...,p_{n}$ are the
tissue spatial mesh coordinates, $\|\quad\|$ denotes the Euclidean distance,
and $S_{\alpha}$ is the smooth maximum approximation function given by:

$S_{\alpha}(x_{1},...x_{n})=\frac{\sum^{n}_{i=1}x_{i}e^{\alpha
x_{i}}}{\sum^{n}_{i=1}e^{\alpha x_{i}}}\qquad\text{where}\quad
S_{\alpha}\rightarrow min\quad\text{as}\quad\alpha\rightarrow-\infty$ (12)

and calculate the heterogeneous activation as $\gamma_{h}=\gamma(1-r)$.

As for the pillars, we define a passive compressible Neo-Hookean constitutive
model as:

$\Psi=\frac{1}{2}\,\mu\,\big{[}\textbf{F}:\textbf{F}-3-2\ln(\det\textbf{F})\big{]}+\frac{1}{2}\,\lambda\,\big{[}\frac{1}{2}[(\det\textbf{F})^{2}-1]-\ln(\det\textbf{F})\big{]}$
(13)

where F is the deformation gradient, and $\mu$ and $\lambda$ are the Lamé
parameters equivalent to Young’s modulus $E$ and Poisson’s ratio $\nu$ as
$E=\mu\,(3\lambda+2\mu)/(\lambda+\mu)$ and $\nu=\lambda/(2(\lambda+\mu))$. 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 $\mathbf{F_{FEA}}$ and its equivalent that is obtained
from approximating a projective transformation $\mathbf{F_{projective}}$ (Fig.
A_2). In our final implementation, we divide each synthetic texture domain
into $16\times 16$ subdomains.

Figure A_2: Convergence of the percentage error of the $x$ and $y$ positions
of the mesh cell centers in comparison to a ground truth with respect to grid
size: (a) The peak (maximum) percentage error in all the $s\times s$ grids for
each respective grid size; (b) The peak (maximum) of the mean percentage error
in each grid of the total $s\times s$ grids.

Overall, our synthetic dataset consists of $60$ generated movies of beating
microbundle textures. We extract $90\times 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\times 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.

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.

Figure A_4: To make our synthetic dataset more challenging and more
representative of real experimental data, we add Perlin noise to the synthetic
microbundle textures: (a) the upper triangle shows an illustration of Perlin
noise (a function of magnitude and octave), and the lower triangle shows the
Perlin noise superimposed on the synthetic microbundle texture; (b) time
series plots of (i) an example of tracked real data showing the drift in
valleys along with the value of the slope of the linear fit of the valleys,
(ii) an example of synthetic data without the addition of Perlin noise where
it 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 Results

####  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 (MAD) (Figs. A_5 and A_6) and subdomain-averaged
$E_{cc}$ strain (strain in the column-column or horizontal-horizontal
direction) (Figs. A_7 – A_22) obtained by comparing the tracking software
output to the known ground truth. For the $16$ original validation examples,
we base our error analysis on $4$ basic metrics: 1) percentage error at the
peak MAD, 2) coefficient of determination ($R^{2}$) between tracked and ground
truth MAD for each beat, 3) coefficient of determination ($R^{2}$) between
tracked and ground truth subdomain-averaged $E_{cc}$ for each beat per
subdomain, and 4) the mean absolute error (MAE) in $E_{cc}$ per subdomain per
beat.

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

Table A_1: A summary of the mean absolute displacement validation errors for
synthetic data without any added Perlin noise.  Table A_2: A summary of the
mean absolute $E_{cc}$ strain validation errors for synthetic data without any
added Perlin noise.

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.

Table A_3: A summary of the mean absolute displacement validation errors for
synthetic data with added Perlin noise.  Table A_4: A summary of the mean
absolute $E_{cc}$ strain validation errors for synthetic data with added
Perlin noise. The “Worst finding” example corresponds to the largest obtained
error whereas the “Less extreme” finding corresponds to the second largest
error. For the “General” case, we consider the mean of the maximum MAE values
for the remaining synthetic examples based on homogeneous activation excluding
examples ST_1, ST_18, and ST_42, as well as the minimum absolute ground truth
$E_{cc}$.

####  A.2.2 Validation against manual tracking

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 Green-Lagrange strain results,
with the peak errors for both outputs falling below $10\%$ for noiseless data
exhibiting microbundle contractions that exceed a single pixel. For these
examples with relatively large deformation, our software proves to be robust
even against significant added noise. However, for sub-pixel deformation
examples, the results are less reliable, with added Perlin noise drastically
affecting their accuracy. This is important context for highlighting the
applicability and limitations of our pipeline.

Figure A_5: Error in mean absolute displacement for synthetic data of “Type 1”
based on FE simulations with homogeneous activation. Figure A_6: Error in mean
absolute displacement for synthetic data of “Type 1” based on FE simulations
with heterogeneous activation. Figure A_7: Error in $E_{cc}$ strain for
‘‘SyntheticTextures_1’’. The numbers $0,1,2,3$ shown on beat 0 MAE plot
indicate the subdomain number. Figure A_8: Error in $E_{cc}$ strain for
‘‘SyntheticTextures_3’’. The numbers $0,1,2,3$ shown on beat 0 MAE plot
indicate the subdomain number. Figure A_9: Error in $E_{cc}$ strain for
‘‘SyntheticTextures_5’’. The numbers $0,1,2,3$ shown on beat 0 MAE plot
indicate the subdomain number. Figure A_10: Error in $E_{cc}$ strain for
‘‘SyntheticTextures_7’’. The numbers $0,1,2,3$ shown on beat 0 MAE plot
indicate the subdomain number. Figure A_11: Error in $E_{cc}$ strain for
‘‘SyntheticTextures_14’’. The numbers $0,1,2,3$ shown on beat 0 MAE plot
indicate the subdomain number. Figure A_12: Error in $E_{cc}$ strain for
‘‘SyntheticTextures_16’’. The numbers $0,1,2,3$ shown on beat 0 MAE plot
indicate the subdomain number. Figure A_13: Error in $E_{cc}$ strain for
‘‘SyntheticTextures_18’’. The numbers $0,1,2,3$ shown on beat 0 MAE plot
indicate the subdomain number. Figure A_14: Error in $E_{cc}$ strain for
‘‘SyntheticTextures_20’’. The numbers $0,1,2,3$ shown on beat 0 MAE plot
indicate the subdomain number. Figure A_15: Error in $E_{cc}$ strain for
‘‘SyntheticTextures_29’’. The numbers $0,1,2,3$ shown on beat 0 MAE plot
indicate the subdomain number. Figure A_16: Error in $E_{cc}$ strain for
‘‘SyntheticTextures_33’’. The numbers $0,1,2,3$ shown on beat 0 MAE plot
indicate the subdomain number. Figure A_17: Error in $E_{cc}$ strain for
‘‘SyntheticTextures_42’’. The numbers $0,1,2,3$ shown on beat 0 MAE plot
indicate the subdomain number. Figure A_18: Error in $E_{cc}$ strain for
‘‘SyntheticTextures_46’’. The numbers $0,1,2,3$ shown on beat 0 MAE plot
indicate the subdomain number. Figure A_19: Error in $E_{cc}$ strain for
‘‘SyntheticTextures_51’’. The numbers $0,1,2,3$ shown on beat 0 MAE plot
indicate the subdomain number. Figure A_20: Error in $E_{cc}$ strain for
‘‘SyntheticTextures_52’’. The numbers $0,1,2,3$ shown on beat 0 MAE plot
indicate the subdomain number. Figure A_21: Error in $E_{cc}$ strain for
‘‘SyntheticTextures_55’’. The numbers $0,1,2,3$ shown on beat 0 MAE plot
indicate the subdomain number. Figure A_22: Error in $E_{cc}$ strain for
‘‘SyntheticTextures_56’’. The numbers $0,1,2,3$ shown on beat 0 MAE plot
indicate the subdomain number. Figure A_23:  Results of our pipeline
validation against noisy synthetic data of “Type 1” based on FE simulations
with homogeneous activation for different Perlin noise octaves and magnitude
ratios (MR). 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).
Figure A_25:  Results of our pipeline validation against synthetic data of
“Type 2”. Figure A_26:  Results of our pipeline validation against manually
tracked points for $2$ different cases of “Type 2” data.

##  Appendix B Additional examples

aa_template addtofilelistaa_template.tex addtofilelistaa_template.aux

Supplementary Document B: Additional Experimental 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_1: A summary of the experimental conditions associated with each
example movie. We note that the symbol † for “Type 3” “Imaging modality” is
used to indicate that Example 1 is obtained via phase contrast microscopy ,
whereas the remainder of the examples are obtained via brightfield microscopy.

Table B_2: A summary of the example-specific information, code implementation
details and subdomain segmentation parameters for each example of “Type 1”
data.

Table B_3: A summary of the example-specific information, code implementation
details, and subdomain segmentation parameters for each example of “Type 2”
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 $\sim$ \- 6 kPa peak amplitude (equivalent to $\sim$ 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.

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.

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.

##  Appendix C Pillar tracking

aa_template addtofilelistaa_template.tex addtofilelistaa_template.aux

Supplementary Document 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:

$\displaystyle F=k\delta$ (14)

where the deflection $\delta$ 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:

$\displaystyle k=\frac{6EI}{a^{2}(3L-a)}$ (15)

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=\frac{wt^{3}}{12}$ where $w$ is the pillar width and
$t$ is the pillar thickness while for circular cross section beams,
$I=\frac{\pi 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\mu$N$/\mu$m and
length scale ($ls$) $4\mu$m/pixel. For the “Type 2” examples, we specify
$k=0.41\mu$N$/\mu$m and $ls=0.908\mu$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 ($\mu$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”. Figure C_2:  Pillar
absolute force ($\mu$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”.

###  C.3 Future Work

In this document, we have provided a description of the pillar tracking
functionality included within the “MicroBundleCompute” 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.

## References

  * Brush et al. [2022] John E Brush, Jonathan Sherbino, and Geoffrey R Norman.  Diagnostic reasoning in cardiovascular medicine.  _bmj_ , 376, 2022.  doi:https://doi.org/10.1136/bmj-2021-064389.
  * Vardas et al. [2022] Panos E Vardas, Folkert W Asselbergs, Maarten van Smeden, and Paul Friedman.  The year in cardiovascular medicine 2021: digital health and innovation.  _European heart journal_ , 43(4):271–279, 2022\.  doi:https://doi.org/10.1093/eurheartj/ehab874.
  * Kaptoge et al. [2019] Stephen Kaptoge, Lisa Pennells, Dirk De Bacquer, Marie Therese Cooney, Maryam Kavousi, Gretchen Stevens, Leanne Margaret Riley, Stefan Savin, Taskeen Khan, Servet Altay, et al.  World health organization cardiovascular disease risk charts: revised models to estimate risk in 21 global regions.  _The Lancet Global Health_ , 7(10):e1332–e1345, 2019.  doi:https://doi.org/10.1016/S2214-109X(19)30318-3.
  * Nguyen et al. [2021] Phong D Nguyen, Dennis EM de Bakker, and Jeroen Bakkers.  Cardiac regenerative capacity: an evolutionary afterthought?  _Cellular and molecular life sciences_ , 78(12):5107–5122, 2021.  doi:10.1007/s00018-021-03831-9.
  * Bursill et al. [2022] Christina A Bursill, Nicola J Smith, Nathan Palpant, Isabella Tan, Margaret Sunde, Richard P Harvey, Benjamin Lewis, Gemma A Figtree, Jamie I Vandenberg, and Australian Cardiovascular Alliance.  Don’t turn off the tap! the importance of discovery science to the australian cardiovascular sector and improving clinical outcomes into the future.  _Heart, Lung and Circulation_ , 2022.  doi:https://doi.org/10.1016/j.hlc.2022.06.669.
  * Thygesen et al. [2007] Kristian Thygesen, Joseph S. Alpert, and Harvey D. White.  Universal definition of myocardial infarction.  _Circulation_ , 116(22):2634–2653, 2007.  doi:10.1161/CIRCULATIONAHA.107.187397.  URL https://www.ahajournals.org/doi/abs/10.1161/CIRCULATIONAHA.107.187397.
  * Takahashi et al. [2007] Kazutoshi Takahashi, Koji Tanabe, Mari Ohnuki, Megumi Narita, Tomoko Ichisaka, Kiichiro Tomoda, and Shinya Yamanaka.  Induction of pluripotent stem cells from adult human fibroblasts by defined factors.  _cell_ , 131(5):861–872, 2007.  doi:https://doi.org/10.1016/j.cell.2007.11.019.
  * Brandão et al. [2017] Karina O Brandão, Viola A Tabel, Douwe E Atsma, Christine L Mummery, and Richard P Davis.  Human pluripotent stem cell models of cardiac disease: from mechanisms to therapies.  _Disease models & mechanisms_, 10(9):1039–1059, 2017.  doi:https://doi.org/10.1242/dmm.030320.
  * Nakao et al. [2020] Shu Nakao, Dai Ihara, Koji Hasegawa, and Teruhisa Kawamura.  Applications for induced pluripotent stem cells in disease modelling and drug development for heart diseases.  _European Cardiology Review_ , 15, 2020.  doi:doi:10.15420/ecr.2019.03.
  * Daly et al. [2021] Andrew C Daly, Matthew D Davidson, and Jason A Burdick.  3d bioprinting of high cell-density heterogeneous tissue models through spheroid fusion within self-healing hydrogels.  _Nature communications_ , 12(1):753, 2021.  doi:https://doi.org/10.1038/s41467-021-21029-2.
  * Dou et al. [2022] Wenkun Dou, Manpreet Malhi, Qili Zhao, Li Wang, Zongjie Huang, Junhui Law, Na Liu, Craig A Simmons, Jason T Maynes, and Yu Sun.  Microengineered platforms for characterizing the contractile function of in vitro cardiac models.  _Microsystems & Nanoengineering_, 8(1):26, 2022\.  doi:https://doi.org/10.1038/s41378-021-00344-0.
  * Hnatiuk et al. [2021] Anna P Hnatiuk, Francesca Briganti, David W Staudt, and Mark Mercola.  Human iPSC modeling of heart disease for drug development.  _Cell chemical biology_ , 28(3):271–282, 2021\.  doi:https://doi.org/10.1016/j.chembiol.2021.02.016.
  * Ronaldson-Bouchard and Vunjak-Novakovic [2018] Kacey Ronaldson-Bouchard and Gordana Vunjak-Novakovic.  Organs-on-a-chip: a fast track for engineered human tissues in drug development.  _Cell stem cell_ , 22(3):310–324, 2018.  doi:https://doi.org/10.1016/j.stem.2018.02.011.
  * Hirt et al. [2014] Marc N Hirt, Arne Hansen, and Thomas Eschenhagen.  Cardiac tissue engineering: state of the art.  _Circulation research_ , 114(2):354–367, 2014\.  doi:https://doi.org/10.1161/CIRCRESAHA.114.300522.
  * Masumoto et al. [2014] Hidetoshi Masumoto, Takeshi Ikuno, Masafumi Takeda, Hiroyuki Fukushima, Akira Marui, Shiori Katayama, Tatsuya Shimizu, Tadashi Ikeda, Teruo Okano, Ryuzo Sakata, et al.  Human ips cell-engineered cardiac tissue sheets with cardiomyocytes and vascular cells for cardiac regeneration.  _Scientific reports_ , 4(1):1–7, 2014.  doi:https://doi.org/10.1038/srep06716.
  * Karbassi et al. [2020] Elaheh Karbassi, Aidan Fenix, Silvia Marchiano, Naoto Muraoka, Kenta Nakamura, Xiulan Yang, and Charles E Murry.  Cardiomyocyte maturation: advances in knowledge and implications for regenerative medicine.  _Nature Reviews Cardiology_ , 17(6):341–359, 2020\.  doi:https://doi.org/10.1038/s41569-019-0331-x.
  * DePalma et al. [2021] Samuel J. DePalma, Christopher D. Davidson, Austin E. Stis, Adam S. Helms, and Brendon M. Baker.  Microenvironmental determinants of organized iPSC-cardiomyocyte tissues on synthetic fibrous matrices.  _Biomater. Sci._ , 9:93–107, 2021.  doi:10.1039/D0BM01247E.  URL http://dx.doi.org/10.1039/D0BM01247E.
  * Batalov et al. [2021] Ivan Batalov, Quentin Jallerat, Sean Kim, Jacqueline Bliley, and Adam W Feinberg.  Engineering aligned human cardiac muscle using developmentally inspired fibronectin micropatterns.  _Scientific Reports_ , 11(1):1–14, 2021.  doi:https://doi.org/10.1038/s41598-021-87550-y.
  * Jayne et al. [2021] Rachael K Jayne, M Çağatay Karakan, Kehan Zhang, Noelle Pierce, Christos Michas, David J Bishop, Christopher S Chen, Kamil L Ekinci, and Alice E White.  Direct laser writing for cardiac tissue engineering: a microfluidic heart on a chip with integrated transducers.  _Lab on a Chip_ , 21(9):1724–1737, 2021.  doi:10.1039/D0LC01078B.
  * Karakan [2023] M Çağatay Karakan.  _A Direct-Laser-Written Heart-on-a-Chip Platform for Generation and Stimulation of Engineered Heart Tissues_.  PhD thesis, Boston University, 2023.
  * Ronaldson-Bouchard et al. [2019] Kacey Ronaldson-Bouchard, Keith Yeager, Diogo Teles, Timothy Chen, Stephen Ma, LouJin Song, Kumi Morikawa, Holly M Wobma, Alessandro Vasciaveo, Edward C Ruiz, et al.  Engineering of human cardiac muscle electromechanically matured to an adult-like phenotype.  _Nature protocols_ , 14(10):2781–2817, 2019.  doi:https://doi.org/10.1038/s41596-019-0189-8.
  * Javor et al. [2021] Josh Javor, Subramanian Sundaram, Christopher S. Chen, and David J. Bishop.  A microtissue platform to simultaneously actuate and detect mechanical forces via non-contact magnetic approach.  _Journal of Microelectromechanical Systems_ , 30(1):96–104, 2021.  doi:10.1109/JMEMS.2020.3036978.
  * Zhao et al. [2019] Yimu Zhao, Naimeh Rafatian, Nicole T Feric, Brian J Cox, Roozbeh Aschar-Sobbi, Erika Yan Wang, Praful Aggarwal, Boyang Zhang, Genevieve Conant, Kacey Ronaldson-Bouchard, et al.  A platform for generation of chamber-specific cardiac tissues and disease modeling.  _Cell_ , 176(4):913–927, 2019.  doi:https://doi.org/10.1016/j.cell.2018.11.042.
  * Ruan et al. [2016] Jia-Ling Ruan, Nathaniel L Tulloch, Maria V Razumova, Mark Saiget, Veronica Muskheli, Lil Pabon, Hans Reinecke, Michael Regnier, and Charles E Murry.  Mechanical stress conditioning and electrical stimulation promote contractility and force maturation of induced pluripotent stem cell-derived human cardiac tissue.  _Circulation_ , 134(20):1557–1567, 2016.  doi:https://doi.org/10.1161/CIRCULATIONAHA.114.014998.
  * Feric et al. [2019] Nicole T Feric, Isabella Pallotta, Rishabh Singh, Danielle R Bogdanowicz, Marietta M Gustilo, Khuram W Chaudhary, Robert N Willette, Tim P Chendrimada, Xiaoping Xu, Michael P Graziano, et al.  Engineered cardiac tissues generated in the biowire ii: a platform for human-based drug discovery.  _Toxicological sciences_ , 172(1):89–97, 2019\.  doi:https://doi.org/10.1093/toxsci/kfz168.
  * Huang et al. [2020] Chen Yu Huang, Rebeca Peres Moreno Maia-Joca, Chin Siang Ong, Ijala Wilson, Deborah DiSilvestre, Gordon F Tomaselli, and Daniel H Reich.  Enhancement of human iPSC-derived cardiomyocyte maturation by chemical conditioning in a 3D environment.  _Journal of molecular and cellular cardiology_ , 138:1–11, 2020.  doi:https://doi.org/10.1016/j.yjmcc.2019.10.001.
  * Lee et al. [2019] Junmin Lee, Vijayan Manoharan, Louis Cheung, Seungkyu Lee, Byung-Hyun Cha, Peter Newman, Razieh Farzad, Shreya Mehrotra, Kaizhen Zhang, Fazal Khan, et al.  Nanoparticle-based hybrid scaffolds for deciphering the role of multimodal cues in cardiac tissue engineering.  _ACS nano_ , 13(11):12525–12539, 2019.  doi:https://doi.org/10.1021/acsnano.9b03050.
  * Huebsch et al. [2022] Nathaniel Huebsch, Berenice Charrez, Gabriel Neiman, Brian Siemons, Steven C Boggess, Samuel Wall, Verena Charwat, Karoline H Jæger, David Cleres, Åshild Telle, et al.  Metabolically driven maturation of human-induced-pluripotent-stem-cell-derived cardiac microtissues on microfluidic chips.  _Nature Biomedical Engineering_ , 6(4):372–388, 2022.  doi:https://doi.org/10.1038/s41551-022-00884-4.
  * Pointon et al. [2017] Amy Pointon, James Pilling, Thierry Dorval, Yinhai Wang, Caroline Archer, and Christopher Pollard.  From the cover: high-throughput imaging of cardiac microtissues for the assessment of cardiac contraction during drug discovery.  _Toxicological Sciences_ , 155(2):444–457, 2017\.  doi:https://doi.org/10.1093/toxsci/kfw227.
  * Huebsch et al. [2015] Nathaniel Huebsch, Peter Loskill, Mohammad A Mandegar, Natalie C Marks, Alice S Sheehan, Zhen Ma, Anurag Mathur, Trieu N Nguyen, Jennie C Yoo, Luke M Judge, et al.  Automated video-based analysis of contractility and calcium flux in human-induced pluripotent stem cell-derived cardiomyocytes cultured over different spatial scales.  _Tissue Engineering Part C: Methods_ , 21(5):467–479, 2015.  doi:https://doi.org/10.1089/ten.tec.2014.0283.
  * Ronaldson-Bouchard et al. [2018] Kacey Ronaldson-Bouchard, Stephen P Ma, Keith Yeager, Timothy Chen, LouJin Song, Dario Sirabella, Kumi Morikawa, Diogo Teles, Masayuki Yazawa, and Gordana Vunjak-Novakovic.  Advanced maturation of human cardiac tissue grown from pluripotent stem cells.  _Nature_ , 556(7700):239–243, 2018.  doi:https://doi.org/10.1038/s41586-018-0016-3.
  * Hansen et al. [2010] Arne Hansen, Alexandra Eder, Marlene Bönstrup, Marianne Flato, Marco Mewe, Sebastian Schaaf, Bülent Aksehirlioglu, Alexander Schwörer, June Uebeler, and Thomas Eschenhagen.  Development of a drug screening platform based on engineered heart tissue.  _Circulation research_ , 107(1):35–44, 2010.  doi:https://doi.org/10.1161/CIRCRESAHA.109.211458.
  * Tamargo et al. [2021] Manuel Alejandro Tamargo, Trevor Ray Nash, Sharon Fleischer, Youngbin Kim, Olaia Fernandez Vila, Keith Yeager, Max Summers, Yimu Zhao, Roberta Lock, Miguel Chavez, et al.  millipillar: a platform for the generation and real-time assessment of human engineered cardiac tissues.  _ACS biomaterials science & engineering_, 7(11):5215–5229, 2021.  doi:https://doi.org/10.1021/acsbiomaterials.1c01006.
  * Rivera-Arbeláez et al. [2022] José M Rivera-Arbeláez, Carla Cofiño-Fabres, Verena Schwach, Tom Boonen, Simone A Ten Den, Kim Vermeul, Albert van den Berg, Loes I Segerink, Marcelo C Ribeiro, and Robert Passier.  Contractility analysis of human engineered 3d heart tissues by an automatic tracking technique using a standalone application.  _Plos one_ , 17(4):e0266834, 2022.  doi:https://doi.org/10.1371/journal.pone.0266834.
  * Oyunbaatar et al. [2016] Nomin-Erdene Oyunbaatar, Deok-Hyu Lee, Swati J Patil, Eung-Sam Kim, and Dong-Weon Lee.  Biomechanical characterization of cardiomyocyte using pdms pillar with microgrooves.  _Sensors_ , 16(8):1258, 2016.  doi:10.3390/s16081258.
  * Dostanić et al. [2020] Milica Dostanić, Laura M. Windt, Jeroen M. Stein, Berend J. van Meer, Milena Bellin, Valeria Orlova, Massimo Mastrangeli, Christine L. Mummery, and Pasqualina M. Sarro.  A miniaturized eht platform for accurate measurements of tissue contractile properties.  _Journal of Microelectromechanical Systems_ , 29(5):881–887, 2020.  doi:10.1109/JMEMS.2020.3011196.
  * Thavandiran et al. [2020] Nimalan Thavandiran, Christopher Hale, Patrick Blit, Mark L Sandberg, Michele E McElvain, Mark Gagliardi, Bo Sun, Alec Witty, George Graham, Van TH Do, et al.  Functional arrays of human pluripotent stem cell-derived cardiac microtissues.  _Scientific reports_ , 10(1):6919, 2020.  doi:https://doi.org/10.1038/s41598-020-62955-3.
  * Tsan et al. [2021] Yao-Chang Tsan, Samuel J DePalma, Yan-Ting Zhao, Adela Capilnasiu, Yu-Wei Wu, Brynn Elder, Isabella Panse, Kathryn Ufford, Daniel L Matera, Sabrina Friedline, et al.  Physiologic biomechanics enhance reproducible contractile development in a stem cell derived cardiac muscle platform.  _Nature Communications_ , 12(1):6167, 2021.  doi:https://doi.org/10.1038/s41467-021-26496-1.
  * Méry et al. [2023] Adrien Méry, Artur Ruppel, Jean Revilloud, Martial Balland, Giovanni Cappello, and Thomas Boudou.  Light-driven biological actuators to probe the rheology of 3d microtissues.  _Nature Communications_ , 14(1):717, 2023.  doi:https://doi.org/10.1038/s41467-023-36371-w.
  * Steadman et al. [1988] BW Steadman, KB Moore, KW Spitzer, and JHB Bridge.  A video system for measuring motion in contracting heart cells.  _IEEE Transactions on Biomedical Engineering_ , 35(4):264–272, 1988.  doi:doi: 10.1109/10.1375.
  * Lim et al. [2008] Chee Chew Lim, Haijun Yang, Mingfeng Yang, Chien-Kao Wang, Jianru Shi, Eric A Berg, David R Pimentel, Judith K Gwathmey, Roger J Hajjar, Michiel Helmes, et al.  A novel mutant cardiac troponin c disrupts molecular motions critical for calcium binding affinity and cardiomyocyte contractility.  _Biophysical journal_ , 94(9):3577–3589, 2008\.  doi:https://doi.org/10.1529/biophysj.107.112896.
  * Hossain et al. [2010] Mohammad Mosharraf Hossain, Eiichi Shimizu, Masato Saito, Sathuluri Ramachandra Rao, Yoshinori Yamaguchi, and Eiichi Tamiya.  Non-invasive characterization of mouse embryonic stem cell derived cardiomyocytes based on the intensity variation in digital beating video.  _Analyst_ , 135(7):1624–1630, 2010.  doi:10.1039/C0AN00208A.
  * Hayakawa et al. [2012] Tomohiro Hayakawa, Takeshi Kunihiro, Suguru Dowaki, Hatsume Uno, Eriko Matsui, Masashi Uchida, Seiji Kobayashi, Akio Yasuda, Tatsuya Shimizu, and Teruo Okano.  Noninvasive evaluation of contractile behavior of cardiomyocyte monolayers based on motion vector analysis.  _Tissue Engineering Part C: Methods_ , 18(1):21–32, 2012.  doi:https://doi-org.ezproxy.bu.edu/10.1089/ten.tec.2011.0273.
  * Ahola et al. [2014] Antti Ahola, Anna L Kiviaho, Kim Larsson, Markus Honkanen, Katriina Aalto-Setälä, and Jari Hyttinen.  Video image-based analysis of single human induced pluripotent stem cell derived cardiomyocyte beating dynamics using digital image correlation.  _Biomedical engineering online_ , 13(1):1–18, 2014.  doi:https://doi.org/10.1186/1475-925X-13-39.
  * Hayakawa et al. [2014] Tomohiro Hayakawa, Takeshi Kunihiro, Tomoko Ando, Seiji Kobayashi, Eriko Matsui, Hiroaki Yada, Yasunari Kanda, Junko Kurokawa, and Tetsushi Furukawa.  Image-based evaluation of contraction–relaxation kinetics of human-induced pluripotent stem cell-derived cardiomyocytes: Correlation and complementarity with extracellular electrophysiology.  _Journal of Molecular and Cellular Cardiology_ , 77:178–191, 2014.  ISSN 0022-2828.  doi:https://doi.org/10.1016/j.yjmcc.2014.09.010.  URL https://www.sciencedirect.com/science/article/pii/S0022282814002910.
  * Rajasingh et al. [2015] Sheeja Rajasingh, Jayakumar Thangavel, Andras Czirok, Saheli Samanta, Katherine F Roby, Buddhadeb Dawn, and Johnson Rajasingh.  Generation of functional cardiomyocytes from efficiently generated human iPSCs and a novel method of measuring contractility.  _PLoS One_ , 10(8):e0134093, 2015.  doi:https://doi.org/10.1371/journal.pone.0134093.
  * Czirok et al. [2017] Andras Czirok, Dona Greta Isai, Edina Kosa, Sheeja Rajasingh, William Kinsey, Zoltan Neufeld, and Johnson Rajasingh.  Optical-flow based non-invasive analysis of cardiomyocyte contractility.  _Scientific reports_ , 7(1):10404, 2017.  doi:10.1038/s41598-017-10094-7.
  * Shradhanjali et al. [2019] Akankshya Shradhanjali, Brandon D Riehl, Bin Duan, Ruiguo Yang, and Jung Yul Lim.  Spatiotemporal characterizations of spontaneously beating cardiomyocytes with adaptive reference digital image correlation.  _Scientific reports_ , 9(1):1–10, 2019.  doi:https://doi.org/10.1038/s41598-019-54768-w.
  * Scalzo et al. [2021] Sérgio Scalzo, Marcelo QL Afonso, Néli J da Fonseca, Itamar CG Jesus, Ana Paula Alves, Carolina ATF Mendonça, Vanessa P Teixeira, Diogo Biagi, Estela Cruvinel, Anderson K Santos, et al.  Dense optical flow software to quantify cellular contractility.  _Cell Reports Methods_ , 1(4), 2021.  doi:https://doi.org/10.1016/j.crmeth.2021.100044.
  * Cheng et al. [2023] Zhiyang Cheng, Yuxin Yang, Kai Jiang, Hongyi Nie, Xingbo Yang, Zizhuo Tu, Jiayi Liang, and Yaozu Xiang.  Quantification of cardiomyocyte contraction in vitro and drug screening by myocytobeats.  _Journal of Cardiovascular Translational Research_ , pages 1–10, 2023\.  doi:10.1007/s12265-023-10357-x.
  * Sala et al. [2018] Luca Sala, Berend J Van Meer, Leon GJ Tertoolen, Jeroen Bakkers, Milena Bellin, Richard P Davis, Chris Denning, Michel AE Dieben, Thomas Eschenhagen, Elisa Giacomelli, et al.  Musclemotion: a versatile open software tool to quantify cardiomyocyte and cardiac muscle contraction in vitro and in vivo.  _Circulation research_ , 122(3):e5–e16, 2018\.  doi:https://doi.org/10.1161/CIRCRESAHA.117.312067.
  * Bouguet et al. [2001] Jean-Yves Bouguet et al.  Pyramidal implementation of the affine lucas kanade feature tracker description of the algorithm.  _Intel corporation_ , 5(1-10):4, 2001.
  * Das et al. [2022] Shoshana L Das, Bryan P Sutherland, Emma Lejeune, Jeroen Eyckmans, and Christopher S Chen.  Mechanical response of cardiac microtissues to acute localized injury.  _American Journal of Physiology-Heart and Circulatory Physiology_ , 323(4):H738–H748, 2022.  doi:https://doi.org/10.1152/ajpheart.00305.2022.
  * Raghupathy [2011] Ramesh Raghupathy.  _Form from function: generalized anisotropic inverse mechanics for soft tissues_.  University of Minnesota, 2011.
  * Genovese et al. [2014] Katia Genovese, Luciana Casaletto, Jay D Humphrey, and Jia Lu.  Digital image correlation-based point-wise inverse characterization of heterogeneous material properties of gallbladder in vitro.  _Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences_ , 470(2167):20140152, 2014.  doi:https://doi.org/10.1098/rspa.2014.0152.
  * Narayanan et al. [2021] Bharath Narayanan, Max L Olender, David Marlevi, Elazer R Edelman, and Farhad R Nezami.  An inverse method for mechanical characterization of heterogeneous diseased arteries using intravascular imaging.  _Scientific Reports_ , 11(1):22540, 2021.  doi:https://doi.org/10.1038/s41598-021-01874-3.
  * Wang et al. [2021] Zhenlin Wang, Jonathan B Estrada, Ellen M Arruda, and Krishna Garikipati.  Inference of deformation mechanisms and constitutive response of soft material surrogates of biological tissue by full-field characterization and data-driven variational system identification.  _Journal of the Mechanics and Physics of Solids_ , 153:104474, 2021.  doi:10.1016/j.jmps.2021.104474.
  * Weiss et al. [2020] Dar Weiss, Cristina Cavinato, Authia Gray, Abhay B Ramachandra, Stephane Avril, Jay D Humphrey, and Marcos Latorre.  Mechanics-driven mechanobiological mechanisms of arterial tortuosity.  _Science advances_ , 6(49):eabd3574, 2020.  doi:10.1126/sciadv.abd3574.
  * Sree et al. [2019] Vivek D Sree, Manuel K Rausch, and Adrian B Tepole.  Linking microvascular collapse to tissue hypoxia in a multiscale model of pressure ulcer initiation.  _Biomechanics and modeling in mechanobiology_ , 18(6):1947–1964, 2019.  doi:https://doi.org/10.1007/s10237-019-01187-5.
  * Kong et al. [2018] Fanwei Kong, Thuy Pham, Caitlin Martin, Raymond McKay, Charles Primiano, Sabet Hashim, Susheel Kodali, and Wei Sun.  Finite element analysis of tricuspid valve deformation from multi-slice computed tomography images.  _Annals of biomedical engineering_ , 46(8):1112–1127, 2018.  doi:https://doi.org/10.1007/s10439-018-2024-8.
  * Joldes et al. [2019] Grand Joldes, George Bourantas, Benjamin Zwick, Habib Chowdhury, Adam Wittek, Sudip Agrawal, Konstantinos Mountris, Damon Hyde, Simon K. Warfield, and Karol Miller.  Suite of meshless algorithms for accurate computation of soft tissue deformation for surgical simulation.  _Medical Image Analysis_ , 56:152–171, 2019.  ISSN 1361-8415.  doi:https://doi.org/10.1016/j.media.2019.06.004.
  * de Zélicourt et al. [2010] Diane A de Zélicourt, Alison Marsden, Mark A Fogel, and Ajit P Yoganathan.  Imaging and patient-specific simulations for the fontan surgery: current methodologies and clinical applications.  _Progress in Pediatric Cardiology_ , 30(1-2):31–44, 2010.  doi:https://doi.org/10.1016/j.ppedcard.2010.09.005.
  * Wang et al. [2015] Zhaoyang Wang, Hien Kieu, Hieu Nguyen, and Minh Le.  Digital image correlation in experimental mechanics and image registration in computer vision: Similarities, differences and complements.  _Optics and Lasers in Engineering_ , 65:18–27, 2015.  doi:https://doi.org/10.1016/j.optlaseng.2014.04.002.
  * Chu et al. [1985] TC Chu, WF Ranson, and Michael A Sutton.  Applications of digital-image-correlation techniques to experimental mechanics.  _Experimental mechanics_ , 25:232–244, 1985.  doi:https://doi.org/10.1007/BF02325092.
  * Zhang and Arola [2004] Dongsheng Zhang and Dwayne D Arola.  Applications of digital image correlation to biological tissues.  _Journal of Biomedical Optics_ , 9(4):691–699, 2004.  doi:https://doi.org/10.1117/1.1753270.
  * Annaidh et al. [2012] Aisling Ní Annaidh, Karine Bruyère, Michel Destrade, Michael D Gilchrist, and Mélanie Otténio.  Characterization of the anisotropic mechanical properties of excised human skin.  _Journal of the mechanical behavior of biomedical materials_ , 5(1):139–148, 2012.  doi:https://doi.org/10.1016/j.jmbbm.2011.08.016.
  * Blaber et al. [2015] J Blaber, B Adair, and Antonia Antoniou.  Ncorr: open-source 2d digital image correlation matlab software.  _Experimental Mechanics_ , 55(6):1105–1122, 2015\.  doi:https://doi.org/10.1007/s11340-015-0009-1.
  * Yang and Bhattacharya [2021] Jin Yang and Kaushik Bhattacharya.  Fast adaptive mesh augmented lagrangian digital image correlation.  _Experimental Mechanics_ , 61(4):719–735, 2021\.  doi:https://doi.org/10.1007/s11340-021-00695-9.
  * Solav et al. [2018] Dana Solav, Kevin M. Moerman, Aaron M. Jaeger, Katia Genovese, and Hugh M. Herr.  Multidic: An open-source toolbox for multi-view 3d digital image correlation.  _IEEE Access_ , 6:30520–30535, 2018.  doi:10.1109/ACCESS.2018.2843725.
  * Murienne and Nguyen [2016] Barbara J Murienne and Thao D Nguyen.  A comparison of 2d and 3d digital image correlation for a membrane under inflation.  _Optics and lasers in engineering_ , 77:92–99, 2016.  doi:https://doi.org/10.1016/j.optlaseng.2015.07.013.