# Good, Better, Best: Choosing Word Embedding Context

James Cross1 Bing Xiang2
1School of EECS
Oregon State University
Corvallis, OR
<EMAIL_ADDRESS>& Bowen Zhou2
2IBM Watson
T. J. Watson Research Center
Yorktown Heights, NY
<EMAIL_ADDRESS>
###### Abstract

We propose two methods of learning vector representations of words and phrases
that each combine sentence context with structural features extracted from
dependency trees. Using several variations of neural network classifier, we
show that these combined methods lead to improved performance when used as
input features for supervised term-matching.

##  1 Introduction

Continuous-space embeddings of words and phrases are very appealing for many
applications in NLP. When used as a feature to represent words, they are much
more compact than “one-hot” vectors. In addition, their geometry can encode a
wealth of information about relationships between words, and thus their
meanings, which can be exploited by downstream applications.

The critical question is how to best achieve that geometry for a given task.
There are two types of training data which are often used to learn word
vectors. Approaches based on the “distributional hypothesis” assume that the
semantics of a word can be determined through the contexts in which it
appears. They model vectors through some form of dimension compression on word
co-occurrence statistics taken from an example text. Other methods learn
vectors from pre-existing information about specific relationships between
words, such as those found in a structured knowledge base.

We propose two vector-learning methods which combine sentence-context input
with structural features extracted automatically from dependency trees for the
same underlying corpus. We compare them with vectors trained from only one of
these types of input, and show that they lead to improved performance on
simple classification tasks.

##  2 Related Work

Many different approaches to learning vector representations of words have
been taken. Methods trained from large text examples by relying directly or
indirectly on word co-occurrence counts date back to latent semantic analysis
[Deerwester, 1990]. Recent examples include the window-based approach of
Mikolov et al. (2013c) and the GloVe method [Pennington et al., 2014]. There
has also been work incorporating prior knowledge with bag-of-words context [Yu
and Dredze, 2014].

It is also possible to perform unsupervised learning on structured information
about the underlying entities or concepts being represented, such as the
relationships encoded in an online knowledge base. Approaches for utilizing
this type of training data include the use of neural networks and/or kernel
density estimation [Bordes et al., 2011] and encoding relations as
translations on hyperplanes [Wang et al., 2014].

The most direct precursor to our work is the well-known Word2Vec algorithm
[Mikolov et al., 2013a, Mikolov et al., 2013b, Mikolov et al., 2013c]. It
consists of a simple single-layer neural model where the goal is to maximize
the ability of the word representations to predict neighbors of the word in a
large raw-text corpus. It was observed by Levy and Goldberg (2014) that this
algorithm could be generalized to use any type of discrete features with which
words could be paired as the “contexts” to be predicted.

##  3 Combined Objective

In this work, we combine sentence-window context with structural features
described in Section 4. The joint objective consists of using both types of
training data in parallel so that the training objective is to maximize the
combined log probability:

$\begin{split}\frac{1}{T}\sum_{t=1}^{T}\sum_{-c\leq j\leq c,j\neq 0}\log
p(w_{t+j}|w_{t})\\\ +\hskip 7.11317pt\frac{\alpha}{|P|}\sum_{(w,f)\in P}\log
p(f|w)\end{split}$ (1)

where $T$ is the number of words in the training corpus, $c$ is the context
window size, $P$ is the set of all pairs of words and structural features
$(w,f)$, and $\alpha$ is a weighting parameter as between the two types of
training data.

The conditional probabilities are modeled with vectors as in Word2Vec, and two
sets of output vectors, internal to the algorithm, are learned to represent
words and structural features as contexts. This objective is implemented with
a joint algorithm which performs updates using the sentence context and the
word/feature pairs in parallel. Both types of updates are based on the
negative sampling training method [Mikolov et al., 2013c].

We also implemented a sequential algorithm where vectors trained with just
sentence-window context are used as initial representations which are
subsequently refined by updates using the structural features. The relative
success of vectors trained in this way suggests that proceeding from a more
general representation toward one which is more directly tailored for a given
task is a promising avenue.

In the sections that follow, we describe a number of simple classification
experiments designed to compare sets of word embeddings trained using these
various methods. For the term-matching tasks we considered, the structural
features alone led to better-performing vectors than those trained with just
linear word context, while the combined-objective methods provided further
improvement.

Feature |  Dependency Pattern |  Example
---|---|---
such_as |  $mwe(\texttt{as}_{i},\texttt{such})\land\newline prep(X,\texttt{as}_{i})\land\newline pobj(\texttt{as}_{i},Y)$ |  … anarchists such as Frederico Urales …
known_as |  $prep(\texttt{known}_{i},\texttt{as}_{j})\land\newline nsubjpass(\texttt{known}_{i},X)\land pobj(\texttt{as}_{j},Y)$ |  The incident became known as the Haymarket affair.
name_for |  $nsubj(X,\texttt{name}_{i})\land\newline prep(\texttt{name}_{i},\texttt{for}_{j})\land\newline pobj(\texttt{for}_{j},Y)$ |  The trivial (non-systematic) name for alkanes is paraffins.

Table 1: Examples of noun-relationship pattern features extracted from
dependency structures. In each case, one feature is created for word $X$ which
includes the feature type and $Y$, and another is created for $Y$ which
includes the feature type and $X$. Subscripts indicate sentence position.

##  4 Structural Features

As a representative means of training word vectors on functional relationships
between words in sentences, rather than just word proximity, we began with
first-order dependency arcs of the type described in Levy and Goldberg (2013).
For every dependency arc, each word is annotated with a feature which consists
of the arc label, whether the word is the head or tail of the arc, and the
adjoined word.111Also included from that work and in the “Dependency” test set
are features coming from the “flattening” of prepositional phrases to directly
link the modified head and the object of the preposition (with the preposition
itself included in the feature name). We also extracted two types of higher-
order features from the dependency trees: noun-relationship patterns, and
subject-object relationships.

Noun-relationship features capture certain taxonomic relationships between
nouns, as reflected in the structure of a sentence, such as when one word
names an instance of a category defined by another word, or when they are
different names for the same thing. This methodology is based on the syntactic
patterns used to recognize is-a relationships in [Fan et al., 2011]. Examples
of these dependency patterns and the relationships they identify can be seen
in Table 1.

In addition, second-order features directly linking subjects and objects of
transitive verbs (and including the verb) are also extracted. Thus for the
sentence “The woman ate the paella” the word “woman” will be annotated with a
feature indicating that it is the subject of the verb “ate” with the direct
object “paella”, and “paella” will receive a feature marking it as the object
of “ate” with the subject “woman”.

The rationale for using both of these types of features is that nouns contain
much of the semantic nexus of a sentence, and they have direct dependency
relationships with essentially all other types of words. By refining word
vectors based on specific relationships between nouns, the quality of the
vector representations as a whole can be improved.

The features used to train vectors for our experiments were extracted from
dependency trees produced by the Stanford Neural Network Dependency Parser
[Chen and Manning, 2014] using the provided pre-trained model with Stanford
dependency labels. The underlying corpus for all methods was the full text of
English Wikipedia as of June 15, 2014.222As maintained by the KOPI project
[Pataki et al., 2012].

##  5 Qualitative Comparison

Target Word |  Raw Text
(Word2Vec) |  Structural
Features |  Joint
---|---|---|---
part-of-speech |  xml-like meta-language sxml script-based canonicalization |  file-type non-template single-word xml-style single-letter |  verb-noun html-like multiword script-based now-deprecated
dynamic |  decouples user-centric gesture-based multi-path process-based |  spacial [sic] timbral user-centric task-focused content-driven |  nonhierarchical composition-based non-virtual ever-changing human/machine
normalize |  normalise normalizing normalising rebalance normalization |  normalise regularize exponentiate recalculate quantify |  normalise regularize quantify define homogenize
Table 2: Most similar 5 words (by cosine similarity) to several target words
under different vector-training regimes.

To evaluate the relative qualities captured by the various embedding methods,
we looked at the closest words to a number of different target words, just as
in Levy and Goldberg (2014). Unsurprisingly, the nearest words were very
similar across the different vector sets, as all were trained with the same
base corpus.

In particular, jointly-trained vectors lead to similar words which are
subjectively very much like those produced by the vectors trained with
structural features only. There is reason to believe, however, that the joint-
objective vectors capture some additional semantic nuance stemming from the
intersection of (and interplay between) topical and functional similarity.

For the word “part-of-speech” the most similar words under all vector sets
skew toward hyphenated or compound modifiers (though the set induced by bag-
of-words training includes nominals). While all of the vector sets seem to
reflect the fact that this modifier is frequently used with the word “tag(s),”
the jointly-trained vectors also capture its affinity with linguistics in such
similar words as “verb-noun” and “multiword.” Likewise with other words such
as “dynamic” and “normalize” the jointly-trained vectors show properties which
are similar to those using structural features alone, but also seem to encode
meaning in a somewhat more targeted way (“ever-changing”; “homogenize”).

##  6 Experiments

In order to compare their predictive properties quantitatively, we tested word
embeddings produced with each method as inputs to binary term-matching
classifiers using two different data sets, described below. We used shallow
neural network classifiers with three alternative architectures, with a
uniform hidden-layer size of 200. The results are presented in Table 3.

For the standard multilayer perceptron (“MLP”) algorithm, the vector
representations of each word in a pair were concatenated to make the input to
a single hidden layer with a leaky rectifier activation function [Glorot et
al., 2011, Maas et al., 2013]. We also used a version where the two vectors
were each used as input to the same hidden layer, i.e., the weights were
shared, and the hidden-layer output for each, together with their element-wise
absolute-value difference, were concatenated as input to the logistic
regression classifier layer (“Shared MLP”).

The third classifier similarly used shared weights to process each vector
input (with a hyperbolic tangent activation function), but classification was
made based on the cosine similarity of the hidden-layer output, with a
threshold of 0.5. In this way, the network learns how to project the vectors
into a space where their angular separation is a good measure of their
similarity for the task at hand. We observed that this network converged very
fast during training, but had performance nearly on par with more complicated
models.

|  MLP |  Shared MLP |  Shared Cosine
---|---|---|---
WordNet
Raw Text (Word2Vec) |  75.88 |  73.76 |  77.10
Dependency |  82.71 |  79.56 |  80.74
\+ Stuctural Features |  82.86 |  80.70 |  81.61
Joint |  82.86 |  80.81 |  81.65
Sequential |  83.12 |  81.61 |  82.75
PPDB
Raw Text (Word2Vec) |  93.43 |  92.84 |  93.20
Dependency |  94.88 |  94.05 |  94.06
\+ Structural Features |  94.82 |  94.00 |  94.09
Joint |  95.52 |  94.63 |  94.51
Sequential |  95.48 |  95.03 |  94.53
Table 3: Accuracy results (%) for classification of WordNet synonym pairs and
PPDB lexical matches for different vector sets

###  6.1 WordNet Synonyms

The first source of term-matching pairs we considered were synonyms from the
WordNet lexical database [Miller, 1995]. WordNet groups words into synsets
which are groups of words (and multi-word phrases) with the same meaning. We
created match data by taking all distinct pairs of words within synsets of any
part of speech, and generated negative examples by random shuffling.333The
WordNet data, as filtered for coverage by these vector sets, consists of
21,089 training pairs, 2,636 development pairs, and 2,637 test pairs, split
evenly between matches and non-matches.

This turns out to be a relatively difficult task with unsupervised vectors
alone. A possible reason is that the vectors do not distinguish between word
senses, essentially folding together all usage cases into a single
representation. The WordNet data is characterized by relatively common words
where the synonymy may rely on a secondary or figurative definition. Examples
of this in our data are the pairs (fastball, bullet) and (arrest, collar). In
particular, since the raw-text data tends to group words by domain of
discourse rather than function, it makes sense that a figurative usage as for
“bullet” and “collar” in those pairs might cause problems.

In this case, the structural features were much more useful than bag-of-words
context. The pattern-based features also led to a clear improvement over using
just the subset of first-order dependency relations, which we hypothesize is
largely due to the subject-object relations. The joint vectors showed little
gain over using structural features only. This seems linked to the large gap
in utility between the two types of training data, and points to the
importance of application domain in determining embedding methods. The
sequentially trained vectors performed better across the board, however, which
suggests the usefulness of “priming” in this type of vector learning
algorithm.

###  6.2 Paraphrase Database Lexical Pairs

As another source of term-matching data, we used pairs from the small
(highest-confidence) lexical paraphrase dataset of the Paraphrase Database
[Ganitkevitch et al., 2013]. They consist of instances where two distinct
words were used to mean the same thing in bi-texts. In contrast to the WordNet
pairs, this dataset skews toward less common words. It includes many instances
of things like alternate spellings and alternate names (or transliterations)
for proper nouns.444The PPDB data consists of 79,696 training pairs, 9,962
validation pairs, and 9,962 test pairs, split 20%/80% between matches and non-
matches.

Here, the performance difference between the raw-text and structural-feature
vectors was much smaller. The combined-objective methods were better than
either. This suggests that when the two types of training data have comparable
degrees of utility, but work through different avenues, a combined training
method can create embeddings which capture the advantages of both approaches.

##  7 Conclusion and Future Work

It is clear that the best method of training continuous word embeddings
depends on the domain of application. When comparably useful information is
encoded in both sentence context and structured representations, combining
both as training data can lead to superior vector representations.

For future work, we plan to try using richer feature sets, such as knowledge-
base relations, as training inputs, and explore the mathematics behind the
sequential training mechanism.

## References

  * [Bergstra et al., 2010]  James Bergstra, O. Breuleux, F. Bastien, P. Lamblin, R. Pascanu, G. Desjardins, J. Turian, D. Warde-Farley and Y. Bengio.  2010\.  Theano: A CPU and GPU Math Expression Compiler.  Proceedings of the Python for Scientific Computing Conference (SciPy).
  * [Bordes et al., 2011]  Antoine Bordes, Jason Weston, Ronan Collobert, and Yoshua Bengio.  2011\.  Structured Embeddings of Knowledge Bases.  AAAI.
  * [Chen and Manning, 2014]  Danqi Chen and Christopher D. Manning.  2014\.  A fast and accurate dependency parser using neural networks.  In Empirical Methods in Natural Language Processing (EMNLP).
  * [Deerwester, 1990]  Scott C. Deerwester, Susan T. Dumais, Thomas K. Landauer, George W. Furnas, and Richard A. Harshman.  1990\.  Indexing by latent semantic analysis.  Journal of the American Society for Information Science, 41.
  * [Fan et al., 2011]  James Fan, Aditya Kalyanpur, J. William Murdock, and Branimir K. Boguraev.  2011\.  Mining Knowledge from Large Corpora for Type Coercion in Question Answering.  Web Scale Knowledge Extraction (WEKEX) Workshop at International Semantic Web Conference.
  * [Fan et al., 2012]  James Fan, Aditya Kalyanpur, D. C. Gondek, and David A. Ferrucci.  2012\.  Automatic knowledge extraction from documents.  IBM Journal of Research and Development 56, no. 3.4: 5-1.
  * [Ganitkevitch et al., 2013]  Juri Ganitkevitch, Benjamin Van Durme, and Chris Callison-Burch.  2013\.  PPDB: The Paraphrase Database.  In North American Chapter of the Association for Computational Linguistics (HLT-NAACL)
  * [Glorot et al., 2011]  Xavier Glorot, Antoine Bordes, and Yoshua Bengio.  2011\.  Deep Sparse Rectifier Neural Networks.  In Proceedings of the 14th International Conference on Artificial Intelligence and Statistics. JMLR W&CP Volume Vol. 15.
  * [Levy and Goldberg, 2014]  Omer Levy and Yoav Goldberg.  2014\.  Dependency-Based Word Embeddings.  In Association for Computational Linguistics (ACL).
  * [Maas et al., 2013]  Andrew L. Maas, Awni Y. Hannun, and Andrew Y. Ng.  2013\.  Rectifier nonlinearities improve neural network acoustic models.  In Workshop on Deep Learning for Audio, Speech, and Language Processing..
  * [McCord et al., 2012]  Michael C. McCord, J. William Murdock, and Branimir K. Boguraev.  2012\.  Deep parsing in Watson.  IBM Journal of Research and Development 56, no. 3.4: 3-1.
  * [Mikolov et al., 2013a]  Tomas Mikolov, Kai Chen, Greg Corrado, and Jeffrey Dean.  2013a.  Efficient estimation of word representations in vector space.  arXiv preprint arXiv:1301.3781.
  * [Mikolov et al., 2013b]  Tomas Mikolov, lya Sutskever, Kai Chen, Greg S. Corrado, and Jeff Dean.  2013b.  Continuous Space Word Representations.  In Advances in Neural Information Processing Systems 26: 27th Annual Conference on Neural Information Processing Systems 2013, pages 3111-3119.
  * [Mikolov et al., 2013c]  Tomas Mikolov, Wen-tau Yih, and Geoffrey Zweig.  2013c.  Distributed representations of words and phrases and their compositionality.  In North American Chapter of the Association for Computational Linguistics (HLT-NAACL)
  * [Miller, 1995]  George A. Miller.  1995\.  WordNet: a lexical database for English.  Communications of the ACM 38.11: 39-41.
  * [Pataki et al., 2012]  Máté Pataki, Miklós Vajna, and Attila Marosi.  2012\.  Wikipedia as Text.  ERCIM NEWS 89: 48-49.
  * [Pennington et al., 2014]  Jeffrey Pennington, Richard Socher, and Christopher D. Manning.  2014\.  GloVe: Global vectors for word representation.  In Empirical Methods in Natural Language Processing (EMNLP).
  * [Wang et al., 2014]  Zhen Wang, Jianwen Zhang, Jianlin Feng, and Zheng Chen.  2014\.  Knowledge graph embedding by translating on hyperplanes.  In Proceedings of the Twenty-Eighth AAAI Conference on Artificial Intelligence.
  * [Yu and Dredze, 2014]  Mo Yu and Mark Dredze.  2014\.  Improving Lexical Embeddings with Semantic Knowledge.  In Association for Computational Linguistics (ACL).


# Sarc-Graph: Automated segmentation, tracking, and analysis of sarcomeres in
hiPSC-derived cardiomyocytes

Bill Zhao
Department of Mechanical Engineering, Boston University
Boston, MA 02215, USA
<EMAIL_ADDRESS>
&Kehan Zhang
Department of Biomedical Engineering, Boston University
The Wyss Institute for Biologically Inspired Engineering, Harvard University
Boston, MA, 02115, USA
<EMAIL_ADDRESS>
&Christopher S. Chen
Department of Biomedical Engineering, Boston University
The Wyss Institute for Biologically Inspired Engineering, Harvard University
Boston, MA, 02115, USA
<EMAIL_ADDRESS>
&Emma Lejeune
Department of Mechanical Engineering, Boston University
Boston, MA 02215, USA
<EMAIL_ADDRESS>

###### Abstract

A better fundamental understanding of human induced pluripotent stem cell-
derived cardiomyocytes (hiPSC-CMs) has the potential to advance applications
ranging from drug discovery to cardiac repair. Automated quantitative analysis
of beating hiPSC-CMs is an important and fast developing component of the
hiPSC-CM research pipeline. Here we introduce “Sarc-Graph,” a computational
framework to segment, track, and analyze sarcomeres in fluorescently tagged
hiPSC-CMs. Our framework includes functions to segment z-discs and sarcomeres,
track z-discs and sarcomeres in beating cells, and perform automated
spatiotemporal analysis and data visualization. In addition to reporting good
performance for sarcomere segmentation and tracking with little to no
parameter tuning and a short runtime, we introduce two novel analysis
approaches. First, we construct spatial graphs where z-discs correspond to
nodes and sarcomeres correspond to edges. This makes measuring the network
distance between each sarcomere (i.e., the number of connecting sarcomeres
separating each sarcomere pair) straightforward. Second, we treat tracked and
segmented components as fiducial markers and use them to compute the
approximate deformation gradient of the entire tracked population. This
represents a new quantitative descriptor of hiPSC-CM function. We showcase and
validate our approach with both synthetic and experimental movies of beating
hiPSC-CMs. By publishing Sarc-Graph, we aim to make automated quantitative
analysis of hiPSC-CM behavior more accessible to the broader research
community.

Author summary (general audience)

Heart disease is the leading cause of death worldwide. Because of this, many
researchers are studying heart cells in the lab and trying to create
artificial heart tissue. Recently, there has been a growing focus on human
induced pluripotent stem cell-derived cardiomyocytes (hiPSC-CMs). These are
cells that are safely sampled from living humans, for example from the blood
or skin, that are then transformed into human heart muscle cells. One active
research goal is to use these cells to repair the damaged heart. Another
active research goal is to test new drugs on these cells before testing them
in animals and humans. However, one major challenge is that hiPSC-CMs often
have an irregular internal structure that is difficult to analyze. At present,
their behavior is far from fully understood. To address this, we have created
software to automatically analyze movies of beating hiPSC-CMs. With our
software, it is possible to quantify properties such as the amount and
direction of beating cell contraction, and the variation in behavior across
different parts of each cell. These tools will enable further quantitative
analysis of hiPSC-CMs. With these tools, it will be easier to understand,
control, and optimize artificial heart tissue created with hiPSC-CMs, and
quantify the effects of drugs on hiPSC-CM behavior.

##  1 Introduction

Quantitative analysis of movies of beating cardiomyocytes is a compelling
approach for connecting cell morphology to dynamic cell function [1, 2]. In
particular, connecting structure and function is a crucial step towards a
better fundamental understanding of human induced pluripotent stem cell-
derived cardiomyocytes (hiPSC-CMs) [3, 4]. In stark contrast to sarcomeres in
mature cardiomyocytes which have a highly ordered regular structure [5],
sarcomeres in hiPSC-CMs are typically immature and disordered [6]. This
irregular structure, combined with large variation between cells, makes
quantitative analysis both more difficult and more pressing [7]. Robust
quantitative analysis frameworks for analyzing beating hiPSC-CMs will help
enable technological advances in drug discovery, genetic cardiac disease, and
cardiac repair [8, 9, 10]. In this work, we aim to make automated quantitative
analysis of hiPSC-CM contractile behavior more accessible to the broader
research community.

Figure 1:  Analyzing movies of contracting human induced pluripotent stem
cell-derived cardiomyoctyes (hiPSC-CMs): (a) Still image of a beating hiPSC-CM
movie; (b) Movie still with a schematic illustration of deformation gradient
$\mathbf{F}_{\mathrm{avg}}$ and tracked sarcomeres overlayed, sarcomere color
corresponds to contraction level with red indicating the highest level of
contraction; (c) Magnitudes of average principal stretches ($\lambda_{1}$,
$\lambda_{2}$) computed from $\mathbf{F}_{\mathrm{avg}}$ with respect to the
movie frame number; (d) Average normalized sarcomere length with respect to
the movie frame number. We note that the definitions of
$\mathbf{F}_{\mathrm{avg}}$, $\lambda_{1}$, and $\lambda_{2}$ are introduced
in this text.

Fig. 1a shows a hiPSC-CM with flourescently labeled z-discs. Movies of these
flourescently labeled hiPSC-CMs contracting are incredibly information rich
[10]. Critically, with dynamic data it is possible to not only measure
structure but also measure aspects of how the structural components
functionally interact [11, 12, 13, 14, 4]. The focus of this work is on
extracting structural and functional information from movies of beating hiPSC-
CMs with flourescently labeled z-discs. To extract biologically relevant
quantitative information from these movies, we focus primarily on segmenting,
tracking, and analyzing individual sarcomeres. For context, each movie will
typically have over $100$ frames, and each frame will contain $100$s of
sarcomeres. Therefore, it is infeasible to perform segmentation and tracking
by hand at scale. In the first component of our computational framework, we
introduce a procedure for automated sarcomere segmentation and tracking.
Example results from segmentation and tracking are show in Fig. 1b.

Building on effective segmentation and tracking, the second contribution of
our computational framework is tools for automated analysis of hiPSC-CM
contractile behavior. To accomplish this, we first draw inspiration from
recent literature on defining summary metrics for describing cardiomyocyte
structure and cytoskeletal structure [15, 16] and mechanical deformation in
multiple cell types [17, 18]. Specifically, we treat the tracked sarcomeres as
fiducial markers and quantify the principal directions and magnitudes of
average cell contraction for each movie frame. An example of applying this
approach is shown in Fig. 1c, with Fig. 1d for comparison. Then, we move
beyond average metrics and extend our framework to readily perform
spatiotemporal analysis of sub-cellular sarcomere contraction. Specifically,
we treat the structurally disordered and complex hiPSC-CMs as spatial graphs
where z-discs are represented as nodes and sarcomeres are represented as
edges. With a spatial graph defined, performing spatial statistics based
analysis is straightforward. For example, it is possible to examine the
correlation between individual sarcomere contraction time series using both
euclidean and network distances.

In conjunction with segmentation, tracking, and analysis, our computational
framework – Sarc-Graph – also contains multiple visualization tools. The
remainder of this paper is focused on describing key components of the
framework. We note that all code is available free and open source on GitHub:
https://github.com/elejeune11/Sarc-Graph. The code is designed in a modular
way such that making major alterations to one component or performing only a
subset of the segmentation, tracking, and analysis steps should be
straightforward when appropriate. Looking forward, we anticipate that this
work is a starting point for a major effort in quantification and statistical
analysis of hiPSC-CM contractile behavior.

##  2 Methods

###  2.1 Data

In general, it is infeasible to perform segmentation and tracking by hand at
scale for movies of hiPSC-CMs with flourescently labeled z-discs. To address
the lack of “ground truth” data we invest significant effort in realistic
synthetic data generation to validate the segmentation and tracking components
of our computational framework. In addition, we show the results of applying
our framework to a selection of real experimental movies.

####  2.1.1 Synthetic data for software validation

In Fig. 2, we show the critical components of the synthetic data generation
pipeline. The first step is illustrated in Fig. 2a where we show an example of
generating the three-dimensional (3D) skeleton of z-disc and sarcomere
geometry for each movie frame. First, a baseline geometry is defined. Then,
the baseline geometry is deformed in space such that each sarcomere is
associated with a ground truth function for contraction (shown here as
fractional change in length) with respect to time. By starting with this 3D
skeleton geometry, it is straightforward to include examples with
inhomogeneous contraction, out of plane orientation and deformation, and
sarcomere chains that appear to intersect when projected into two-dimensional
(2D) space.

Figure 2:  Creating synthetic data to test the image analysis code: (a) Define
the geometry of z-discs (points) and sarcomeres (line segments) in three-
dimensional (3D) space and prescribe the ground truth deformation for each
sarcomere; (b) Approximate each z-disc as an oriented cylinder in 3D space;
(c) Convert the 3D geometry into a voxel array with higher resolution in the x
and y directions than the z direction; (d) Slice and project the voxel array
to obtain a two-dimensional (2D) image. Image noise and blur can be added to
the voxel array, 2D image, or both.

Given the 3D skeleton geometry in Fig. 2a, we then render each frame as
illustrated in Fig. 2b-d. First, each z-disc is treated as an oriented
cylinder with an associated height and radius. We note that with this approach
it is straightforward to include multiple z-disc sizes in the same image.
Then, the domain is converted into a voxel array such that the resolution is
representative of the experimental images of interest. The voxel array, a 3D
matrix, is then sliced around the simulated focal plane and projected into 2D.
Both image noise and blur can be added at the voxel array step, 2D image step,
or both. In the representative examples shown here, we add Gaussian blur to
the voxel array [19] and Perlin noise to the 2D image [20]. We note that the
entire synthetic data generation pipeline is available on GitHub. All code is
written in Python with heavy use of the numpy and matplotlib packages [21,
22].

####  2.1.2 Experimental data for software demonstration

|  key information on the experimental data |  key Sarc-Graph analysis challenges
---|---|---
E1 |  hiPSC-CMs cultured on a fibronectin coated glass substrate, vertically aligned fibers. |  large deformation
E2 |  hiPSC-CMs cultured on a fibronectin coated glass substrate, tangentially aligned fibers. |  large deformation
E3 |  Single paced hiPSC-CM constrained on a micropatterned island ($2000\,\mu m^{2}$) with fibronectin coating on top of soft polyacrylamide hydrogel ($7.9$ kPa). The method for making the hydrogel is described in [23]. |  low spatial resolution, large deformation
E4 |  Monolayer of paced hiPSC-CMs cultured on a fibronectin-coated glass substrate. |  substantial background signal
E5 |  Single paced hiPSC-CMs cultured on a fibronectin-coated glass substrate. |  low spatial resolution, small deformation
Table 1:  Summary of experimental examples included in this paper (E1, E2, E3,
E4, and E5). We note that Examples E1 and E2 have already been published and
made publicly available [10].

The first two examples of experimental data included with this paper (E1 and
E2) were previously made available in a publication by Toepfer et al. [10].
The protocol to introduce green fluorescent protein (GFP) onto z-disc protein
titin (TTN-GFP) in hiPSC-CMs is reported in [24]. Video imaging was conducted
on small clusters of hiPSC-CMs using a $100$X objective of a fluorescent
microscope with a minimum acquisition rate of 30 frames per second. The three
additional examples included with this publication were selected to showcase
the general applicability of Sarc-Graph. Cells in E3, E4, and E5 were
fluorescently labeled using lentiviral constructs lenti-GFP-actinin-2 or
lenti-mApple-actinin-2 as described in [23], and were electrically stimulated
at 1Hz using the IonOptix C-Pace EP Culture Pacer (IonOptix) during imaging.
Time-lapse videos of the cells contraction were acquired at 30 frames per
second using a 40x objective on a Nikon Eclipse Ti (Nikon Instruments, Inc.)
with an Evolve EMCCD Camera (Photometrics) or on a Zeiss Axiovert 200M
inverted spinning disk microscope with an ORCA-100 Camera (Hamamatsu),
equipped with a temperature and CO2 equilibrated environmental chamber. Key
details are summarized in Table 1. We note that, though not included as
examples with this publication, our computational framework has also been
tested on additional beating hiPSC-CM movies obtained by different researchers
with subtly different experimental set ups. We also note that for all examples
shown here, both experimental and synthetic, our framework does not require
any parameter tuning. All input parameters are identical across all cases
shown in this paper.

###  2.2 Code

Sarc-Graph is divided into seven modular components: file_pre_processing,
segmentation, tracking, time_series, spatial_graph, analysis_tools, and
functional_metrics. The first, file_pre_processing, is simply used to convert
each movie into a series of 2D numpy arrays [21], one for each frame of the
movie converted to grayscale. The code is designed such that the only new step
required to support a previously unsupported movie file format is a function
to convert each frame into a 2D numpy array. The remainder of the code can
loosely be grouped as segmentation which is necessary for extracting spatial
data, tracking and time_series which are necessary for extracting temporal
data, and spatial_graph, analysis_tools, and functional_metrics which are
necessary for synthesizing and interpreting the spatial and temporal data.

####  2.2.1 Segmentation

Figure 3:  Segmenting and tracking the z-discs and sarcomeres: (a) An example
of a raw two-dimensional (2D) input image; (b) The Laplacian of the input
image is used to detect the high gradients present at the edge of every z-disc
[25]; (c) z-discs are identified as closed contours where the value of the
Laplacian exceeds the threshold computed with Otsu’s method [26]; (d)
Sarcomeres are procedurally identified from the segmented z-discs; (e)
Sarcomeres and z-discs are tracked independently between frames with the
Python trackpy package [27]; (f) The algorithm to segment sarcomeres is based
on linking the approximately parallel z-discs to their closest neighbors in
the direction perpendicular to the z-disc; (g) Sarcomere properties are
computed from the pair of associated z-discs; (h) Tracking each sarcomere
leads to multiple spatially resolved time series.

Within the segmentation script, the function segmentation_all(folder_name,
gaussian_filter_size) will segment all z-discs and all sarcomeres from each
frame. The first parameter, “folder_name” specifies the folder that the data
is in. The second parameter, “gaussian_filter_size” specifies the standard
deviation for the Gaussian kernel applied to the image during z-disc
segmentation [28]. The default value for gaussian_filter_size $=1$ and that is
the value for all data shown in this paper. In our experience, a value of $2$
may be more appropriate for movies of hiPSC-CMs that have high levels of
irregular fluctuatig background signal, and thus it is left as a potential
tunable parameter with a default value of $1$. The key steps of segmentation
are illustrated in Fig. 3. z-disc segmentation is conducted by computing the
Laplacian of the image [25], applying a Gaussian filter to the Laplacian [28],
and then using the scikit-image “measure.find_contours()” function [29] to
identify contours with a level matching a scalar threshold computed via Otsu’s
method [26]. We note briefly that a global threshold determined with Otsu’s
method is sufficient (i.e., an adaptive threshold is not required) because
thresholding is performed on the Laplacian rather than the original image. As
illustrated in Fig. 3c, each contour represents a segmented z-disc and the
properties of each z-disc are computed algorithmically from the contours. Of
note, contour position is computed as the mean position of each pixel in the
contour, contour length is computed as the maximum possible distance between
two pixels in the contour, and contour endpoints are the locations of the
pixels that correspond to the maximum distance.

Once z-discs are segmented, sarcomeres are procedurally identified from the
segmented z-discs with a custom algorithm schematically illustrated in Fig.
3f. In brief, the steps of the algorithm are as follows:

  * 1.

Compute the distance between the center of each z-discs and the center of it’s
nearest neighbor.

  * 2.

Define median_neigh as the median distance between a z-disc and its nearest
neighbor.

  * 3.

For each z-disc, define a line that goes through the center of the z-disc and
is perpendicular to the axis defined by the contour endpoints. Then, define
the two points on the line that are median_neigh$/2$ distance from the z-disc
center. These points are referred to as “ghost_points.” This is illustrated in
Fig. 3f.

  * 4.

Find the nearest neighbor for each point in ghost_points.

  * 5.

If there is a two-way nearest neighbor match between a pair in ghost_points,
the pair of z-discs will correspond to a segmented sarcomere.

As illustrated in Fig. 3f, this algorithm links approximately parallel
z-discs. In order to flexibly accommodate potential image artifacts, the
criteria for “approximately parallel” is not strictly defined. However, in
most cases, there will not be a persistent two way nearest neighbor match
unless z-discs bound a convex quadrilateral, ideally a trapezoid. If
necessary, additional match rejection criteria could be integrated into the
algorithm. Once a sarcomere is identified, its properties are computed from
the properties of the corresponding paired z-discs. This is illustrated in
Fig. 3g. Of note, sarcomere width is computed as the mean length of the two
z-discs, sarcomere length is computed as the distance between the two z-disc
centers, and sarcomere angle is computed as the angle of the vector connecting
the two z-disc centers.

For the examples shown in this paper, the segmentation step takes order of
$10$s of seconds to a few minutes to run on a single laptop. We also note that
users working with still images rather than movies can still segment z-discs
and sarcomeres from a single frame, and can still create a spatial graph with
their data and perform several aspects of the spatial analysis described in
this paper. For generality, all length data is reported in units of pixels.

####  2.2.2 Tracking and processing time series data

The tracking component of Sarc-Graph is for independently tracking all z-discs
and all sarcomeres. Tracking the sarcomeres is necessary for reporting the
relative length change from frame to frame. Without tracking, it is only
possible to report population average relative length change which is
substantially less information rich because the sarcomeres do not necessarily
contract in perfect synchrony, and are not necessarily of identical lengths.
We track the z-discs and sarcomeres independently because we use this
information to construct a spatial graph, described later in the text.
Tracking is performed by calling the “run_all_tracking(folder_name,tp_depth)”
function. As stated previously, “folder_name” specifies the folder that the
data is in. The second parameter, “tp_depth” specifies the farthest distance
in pixels that a particle can travel between frames [27]. The default value
for tp_depth $=4$ and that is the value for all data shown in this paper. In
some cases, particularly for movies with lower spatial resolution, it is
necessary to change to tp_depth $=3$. Particle tracking is performed with the
Python package trackpy which is based on the Crocker–Grier algorithm [27, 30].
We note that we do not use trackpy for feature detection, instead we convert
each segmented z-disc and sarcomere to a pandas DataFrame for compatibility
with the trackpy framework [31]. The outcome of running run_all_tracking() is
a unique ID for each particle that is tracked for over $10\%$ of the movie.
With this information, it is possible to locate each tracked particle in space
across multiple frames, as illustrated in Fig. 3e where the position of each
sarcomere across all frames is plotted.

The time_series component of our computational framework processes the results
from tracking so that there is time series data describing the length,
normalized length $(L-\bar{L})/\bar{L}$, width, angle, and position for each
tracked sarcomere. This is accomplished by running the function
“timeseries_all(folder_name, keep_thresh)” where keep_thresh represents the
fraction of movie frames that a tracked sarcomere must be present in to get
processed. We recommend keep_thresh$=0.75$ for quantitative analysis of time
series data and keep_thresh$>1/k$ for visualization, where $k$ is the
approximate number of times the cell contracts during one movie. To process
the results from tracking, Gaussian process regression is used to interpolate
data across frames which temporarily lose signal [32]. Gaussian process
regression is implemented through python scikit-learn with a radial basis
function (RBF) and white noise kernel [33]. In addition to interpolating
between lost frames, Gaussian process regression has the added benefit of
adaptively smoothing the data without a need for additional dataset-specific
input parameters. The main result of running the timeseries_all function is
illustrated in Fig. 3h where normalized sarcomere length is plotted with
respect to frame number. Even with this synthetic data example, the individual
recovered time series deviate from the ground truth. This is for two main
reasons. First, for the ground truth geometry that we specified, some
sarcomeres are angled out of the image z-plane and are thus not perfectly
captured by a 2D image. Second, the limited resolution of the images
introduces artifacts where lengths less than 1 pixel cannot be perfectly
captured. However, the mean of all time series curves illustrated in Fig. 3h
is a near perfect fit to the ground truth.

For the examples shown in this paper, the tracking step takes order of $10$s
of seconds to a minute to run on a single laptop. The time_series step takes
on the order of $10$s of seconds to a few minutes to run. The time_series step
will output several user friendly text files describing sarcomere properties
with respect to time where each row corresponds to a tracked sarcomere and
each column corresponds to a movie frame. The normalized length data can be
automatically plotted with the function plot_normalized_tracked_timeseries()
from the analysis_tools file. For generality, all time data is reported in
units of frames.

Figure 4:  Data analysis from segmentation and tracking: (a) Tracked elements
are treated as fiducial markers and used to construct an average deformation
gradient $\mathbf{F}$ (see eqns. 1-4); (b) z-discs and sarcomeres are tracked
independently and used to construct a spatial graph (z-discs are “nodes” and
sarcomeres are “edges”, edge color corresponds to sarcomere orientation, node
color corresponds to correlation in sarcomere orientation) to measure
correlations based on network distance; (c) Automated analysis of time series
data is used to extract constants such as mean contraction time and number of
peaks.

####  2.2.3 Computing average deformation ($\mathbf{F}_{\mathrm{avg}}$)

Developing methods to describe data rich images and movies of cells with a
tractable set of output parameters is an active area of research [34, 16]. In
this paper specifically, we draw inspiration from recent work on summarizing
data from traction force microscopy experiments [18] and agent based model
simulations [35, 36] and propose a method to compute the mean deformation
gradient of each domain with respect to time. First, we define the standard
continuum mechanics deformation gradient $\mathbf{F}$ as follows:

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

where $d\mathbf{X}$ is a vector in the initial configuration, $d\mathbf{x}$ is
the vector in the deformed configuration, and $\mathbf{F}$ is the mapping
between them [37]. In the context of analyzing movies of contracting hiPSC-
CMs, we first define a set of $n$ vectors $\mathbf{v}$ that connect each
potential pair of tracked fiducial markers, as illustrated in Fig. 4a. The
direction and magnitude of each vector $\mathbf{v}$ will change as the
fiducial markers move between movie frames. With this definition, we can set
up the over-determined system of equations:

$\
\mathbf{F}_{\mathrm{avg}}\,\mathbf{\Lambda}_{0}=\mathbf{\Lambda}\qquad\mathrm{where}\qquad\mathbf{\Lambda}_{0}=[\mathbf{v}_{01},\mathbf{v}_{02},...,\mathbf{v}_{0n}]\qquad\mathrm{and}\qquad\mathbf{\Lambda}=[\mathbf{v}_{1},\mathbf{v}_{2},...,\mathbf{v}_{n}]$
(2)

where $\mathbf{F}_{\mathrm{avg}}$ is a $2\times 2$ matrix,
$\mathbf{\Lambda}_{0}$ is a $2\times n$ matrix of vectors in the initial
(reference) configuration movie frame, and $\mathbf{\Lambda}$ is a $2\times n$
matrix of vectors in the current (deformed) configuration movie frame. Note
that when the initial frame and the current frame are identical,
$\mathbf{\Lambda}_{0}=\mathbf{\Lambda}$ and
$\mathbf{F}_{\mathrm{avg}}=\mathbf{I}$. The initial configuration can be
defined as either the first frame of the movie, or a selected movie frame
where the cell is known to be in a relaxed state. Then, we can use the normal
equation to solve for the best fit average deformation gradient as:

$\
\mathbf{F}_{\mathrm{avg}}=\mathbf{\Lambda}\mathbf{\Lambda}_{0}^{\mathrm{T}}\,[\,\mathbf{\Lambda}_{0}\mathbf{\Lambda}_{0}^{\mathrm{T}}\,]^{-1}\,.$
(3)

Once computed, the components of $\mathbf{F}_{\mathrm{avg}}$ can be plotted
directly, or $\mathbf{F}_{\mathrm{avg}}$ can be manipulated to provide
information that is more convenient to compare and interpret. For example, we
can write

$\
\mathbf{F}_{\mathrm{avg}}=\mathbf{R}_{\mathrm{avg}}\mathbf{U}_{\mathrm{avg}}$
(4)

where $\mathbf{R}_{\mathrm{avg}}$ is a rotation tensor and
$\mathbf{U}_{\mathrm{avg}}$ is a symmetric positive definite tensor that
represents stretch. In the numerical setting, $\mathbf{R}_{\mathrm{avg}}$ and
$\mathbf{U}_{\mathrm{avg}}$ are computed with the scipy polar decomposition
function [28]. We can also compute the eigenvalues and eigenvectors of
$\mathbf{U}_{\mathrm{avg}}$, $\mathbf{\lambda}_{1}^{\mathrm{avg}}$,
$\mathbf{\lambda}_{2}^{\mathrm{avg}}$, $\mathbf{u}_{1}^{\mathrm{avg}}$, and
$\mathbf{u}_{2}^{\mathrm{avg}}$ respectively. In the numerical setting, this
is done with the numpy linalg package, and we always define
$\mathbf{\lambda}_{1}^{\mathrm{avg}}<\mathbf{\lambda}_{2}^{\mathrm{avg}}$
[21]. From a mechanical data interpretation perspective, the average stretch
values $\mathbf{\lambda}_{1}^{\mathrm{avg}}$ and
$\mathbf{\lambda}_{2}^{\mathrm{avg}}$ and the associated eigevectors are
particularly convenient because they contain information about both the
magnitude and direction of contraction. This is illustrated in Fig. 1b and 1c.

Because sarcomeres data is already processed through the time_series component
of our framework, we use the sarcomeres as our fiducial markers to define
$\mathbf{\Lambda}$ and $\mathbf{\Lambda}_{0}$. We note that this approach
would be just as effective with the z-discs chosen as fiducial markers. In
addition, if 3D data was available, this approach would still work with
$\mathbf{F}_{\mathrm{avg}}$ as a $3\times 3$ matrix and $\mathbf{\Lambda}$ and
$\mathbf{\Lambda}_{0}$ as $3\times n$ matrices. Within the analysis_tools
file, the function compute_F_whole_movie() will compute
$\mathbf{F}_{\mathrm{avg}}$ for each movie frame. The function
visualize_F_full_movie() will create a visualization of
$\mathbf{\lambda}_{1}^{\mathrm{avg}}$ and
$\mathbf{\lambda}_{2}^{\mathrm{avg}}$ as a function of the movie frame and
schematically illustrate $\mathbf{F}_{\mathrm{avg}}$ overlaid on the original
image data, as shown in Fig. 1. The computational cost of computing and
processing $\mathbf{F}_{\mathrm{avg}}$ for all frames is on the order of
seconds to $10$s of seconds on a single laptop. The computational cost of data
visualization is on the order of $10$s of seconds to minutes.

####  2.2.4 Comparing average deformation ($\mathbf{F}_{\mathrm{avg}}$) to
average sarcomere shortening ($\tilde{s}$, $s_{\mathrm{avg}}$) and
Orientational Order Parameter (OOP)

In addition to defining tensor $\mathbf{F}_{\mathrm{avg}}$ as a dynamic metric
that evolves over the course of the entire beating hiPSC-CM movie, it is
possible to formulate scalar metrics from $\mathbf{F}_{\mathrm{avg}}$ that are
more straightforward to directly compare to well established metrics for
describing hiPSC-CMs [7]. Here we focus on defining scalar metrics from
$\mathbf{F}_{\mathrm{avg}}$ that are directly comparable to average sarcomere
shortening ($\tilde{s}$, $s_{\mathrm{avg}}$) [10], and the Orientational Order
Parameter (OOP) [11, 16, 38].

We compute average sarcomere shortening in two different ways. In both cases,
we work with the normalized sarcomere length obtained from time_series
$y=(L-\bar{L})/\bar{L}$ where $L$ is sarcomere length in pixels. Given
normalized sarcomere length $y$ for each movie frame, we define shortening as:

$\ s=\frac{y_{\mathrm{max}}-y_{\mathrm{min}}}{y_{\mathrm{max}}+1}$ (5)

where $y_{\mathrm{max}}$ is the maximum normalized length and
$y_{\mathrm{min}}$ is the minimum normalized length. Then, we can define
$\tilde{s}$ as the median value of $s$ for all sarcomeres tracked in the
movie. We choose to use median rather than mean in this case because the mean
is more susceptible to outliers and thus $\tilde{s}$ is a better reflection of
the ground truth sarcomere deformation. In addition, we can compute the mean
normalized length time series $y_{\mathrm{avg}}$ and then compute
$s_{\mathrm{avg}}$ based on the single mean time series. If all sarcomeres are
contracting identically, $\tilde{s}$ and $s_{\mathrm{avg}}$ will be identical.
However, if this is not the case, the two values will likely differ and
$\tilde{s}$ will reflect average individual sarcomere behavior while
$s_{\mathrm{avg}}$ will reflect both individual sarcomeres and the (lack of)
synchrony between them. Functions to compute $\tilde{s}$ and
$s_{\mathrm{avg}}$ are given in the functional_metrics component of Sarc-
Graph.

In addition to computing mean sarcomere shortening, we specify the method for
computing OOP from hiPSC-CM movies that we implement in the functional_metrics
component of Sarc-Graph. When hiPSC-CMs are fully relaxed, individual
sarcomere chains are no longer under systolic tension and can thus become wavy
and lose local alignment. Therefore, the first step to computing OOP is to
define a frame of the movie where the lowest fraction of sarcomere chains will
be in their relaxed and potentially wavy configuration. To select this frame,
we compute $\det(\mathbf{F}_{\mathrm{avg}})$ for every movie frame with the
first frame of the movie used as the reference configuration. Then, we
identify the frame where $\det(\mathbf{F}_{\mathrm{avg}})$ is the greatest
(i.e. the sarcomeres are most relaxed) and select that frame as the reference
frame and recompute $\mathbf{F}_{\mathrm{avg}}$. With the appropriately
selected reference frame, we re-compute $\det(\mathbf{F}_{\mathrm{avg}})$ and
define the most contracted frame as the frame with the lowest value of
$\det(\mathbf{F}_{\mathrm{avg}})$.

We compute OOP for the static image corresponding to the most contracted frame
of the movie. To do this, we define structural tensor $\mathbb{T}$ following
the literature:

$\
\mathbb{T}\,=\,\biggl{<}\,2\begin{bmatrix}r_{x}^{i}r_{x}^{i}&r_{x}^{i}r_{y}^{i}\\\
r_{y}^{i}r_{x}^{i}&r_{y}^{i}r_{y}^{i}\end{bmatrix}\,-\,\begin{bmatrix}1&0\\\
0&1\end{bmatrix}\,\biggr{>}$ (6)

where $\mathbf{r}^{i}=[r_{x}^{i},r_{y}^{i}]$ is the unit vector describing the
orientation of the $i^{\mathrm{th}}$ sarcomere [11]. Given $\mathbb{T}$, with
eigenvalues $u_{\mathrm{max}}$ and $u_{\mathrm{min}}$ and corresponding unit
eigenvectors $\mathbf{v}_{\mathrm{max}}$ and $\mathbf{v}_{\mathrm{min}}$,
OOP$=u_{\mathrm{max}}$. When OOP$=0.0$, sarcomeres are oriented randomly. When
OOP$=1.0$, sarcomeres are perfectly aligned. Corresponding eigenvector
$\mathbf{v}_{\mathrm{max}}$ corresponds to the dominant direction of sarcomere
orientation.

In order to best relate $\mathbf{F}_{\mathrm{avg}}$ to OOP, we define two
scalar metrics $C_{\mathrm{iso}}$ and $C_{||}$ that describe average radial
contraction, and contraction in the direction of dominant sarcomere
orientation $\mathbf{v}_{\mathrm{max}}$ respectively. Unlike individual
sarcomere contraction or the average of individual sarcomere contraction,
$C_{\mathrm{iso}}$ and $C_{||}$ represent deformation throughout the whole
domain and not necessarily along the sarcomere axis. Metric $C_{\mathrm{iso}}$
is defined as

$\ C_{\mathrm{iso}}=1.0-\sqrt{\det(\mathbf{F}_{\mathrm{avg}})}$ (7)

where $\mathbf{F}_{\mathrm{avg}}$ is computed in the most contracted frame.
Practically, $C_{\mathrm{iso}}$ is the line shortening during contraction of
an equivalent system where contraction is not directionally dependent. When
$C_{\mathrm{iso}}=0.0$, no shortening occurs. As an example, a value
$C_{\mathrm{iso}}=0.10$ corresponds to an equivalent isotropic line shortening
of $10\%$ in all directions. Essentially, higher $C_{\mathrm{iso}}$ indicates
higher levels of hiPSC-CM contraction.

To quantify shortening in the dominant direction specified by
$\mathbf{v}_{\mathrm{max}}$ computed from $\mathbb{T}$, we consider the
relation $\mathbf{F}_{\mathrm{avg}}\mathbf{v}_{0}=\mathbf{v}_{\mathrm{max}}$
where $\mathbf{v}_{0}$ is the corresponding vector defined in the reference
frame. We can manipulate this equation to compute
$\mathbf{v}_{0}=\mathbf{F}_{\mathrm{avg}}^{-1}\mathbf{v}_{\mathrm{max}}$ and
then define $C_{||}$ as

$\
C_{||}=\frac{|\mathbf{v}_{0}|-|\mathbf{v}_{\mathrm{max}}|}{|\mathbf{v}_{0}|}\,.$
(8)

With this definition, $C_{||}$ is the fractional shortening of a line in the
direction of $\mathbf{v}_{\mathrm{max}}$ where direction is defined in the
contracted configuration. Conveniently, lower $\tilde{s}$ and
$s_{\mathrm{avg}}$, $C_{\mathrm{iso}}$, and $C_{||}$ all correspond to lower
levels of contraction, and higher $\tilde{s}$ and $s_{\mathrm{avg}}$,
$C_{\mathrm{iso}}$, and $C_{||}$ all correspond to higher levels of
contraction.

Functions to compute and visualize $\tilde{s}$, $s_{\mathrm{avg}}$, OOP,
$C_{\mathrm{iso}}$, and $C_{||}$ are all defined in the functional_metrics
script. The function “compute_metrics” will compute $\tilde{s}$,
$s_{\mathrm{avg}}$, OOP, $C_{\mathrm{iso}}$, and $C_{||}$ and create a
visualization of OOP and $\mathbf{F}_{\mathrm{avg}}$. The function
“visualize_lambda_as_functional_metric” will create a movie of
$\mathbf{\lambda}_{1}^{\mathrm{avg}}$ and
$\mathbf{\lambda}_{2}^{\mathrm{avg}}$ changing over the course of the movie.
Examples of this are provided as Supplementary Information. The computational
cost of computing $\mathbf{F}_{\mathrm{avg}}$ (including computing
$C_{\mathrm{iso}}$, and $C_{||}$) and $\tilde{s}$, $s_{\mathrm{avg}}$, and OOP
for all frames is on the order of seconds to $10$s of seconds on a single
laptop.

####  2.2.5 Preparing spatial graphs

In addition to analyzing average cell and sarcomere behavior (ex: average
z-disc and sarcomere morphological properties, $\mathbf{F}_{\mathrm{avg}}$),
we are interested in analyzing the spatiotemporal patterns that arise on the
sub-cellular level. We implemented the spatial_graph component of Sarc-Graph
to conveniently synthesize the outputs of the segmentation, tracking, and
time_series steps. Specifically, running the function create_spatial_graph()
will create a spatial graph where z-discs are represented as nodes and
sarcomeres are represented as edges. With a spatial graph structure, it is
then possible to better quantify spatiotemporal patterns in hiPSC-CM behavior.
For example, access to a spatial graph can help delineate sarcomere synchrony
that depends on Euclidean distance, and sarcomere synchrony that depends on
network distance. In addition, access to a spatial graph can help us better
quantify both global and local sarcomere alignment. This functionality is
illustrated in Fig. 4b.

In the run_all_tracking() step, z-discs and sarcomeres are tracked
independently. Every z-disc that is at least partially tracked (defined as
present in over $10\%$ of the movie frames) will have a unique local (frame-
specific) ID and a unique global ID. Each node of the spatial graph
corresponds to a unique global ID, and the local IDs specific to each frame
are linked to the corresponding global ID node. Each tracked sarcomere is
associated with two local z-disc IDs in each frame. To add sarcomeres to the
spatial graph, the global ID of the associated z-disc is first determined.
Then, either a new edge is created linking the two global IDs, or, if the edge
already exists, the weight of the edge is increased. The spatial graph is
created with the NetworkX python package, and all operations, such as
computing the shortest path between two nodes, rely on built-in NetworkX
functionality [39]. The function create_spatial_graph() runs in seconds on a
single laptop. Running create_spatial_graph() will also produce a schematic
drawing of the spatial graph.

####  2.2.6 Notable analysis and visualization tools

In addition to the functions to compute and visualize
$\mathbf{F}_{\mathrm{avg}}$ and other functional metrics
(compute_F_whole_movie(), visualize_F_full_movie(), compute_metrics,
visualize_lambda_as_functional_metric) and plot normalized sarcomere length
with respect to time (plot_normalized_tracked_timeseries()), the
analysis_tools component of Sarc-Graph contains multiple functions to further
visualize and analyze the data acquired from the hiPSC-CM movies. Example
outputs from these functions are included onn the Sarc-Graph GitHub page.
Additional key analysis and visualization functions are as follows:

  * •

visualize_segmentation(folder_name, gaussian_filter_size, frame_num):
Visualization of the segmented z-discs and sarcomeres overlaid on the original
image. The parameter gaussian_filter_size should match the value chosen in
segmentation. The default for gaussian_filter_size$\,=1$ and the default for
frame_num$\,=0$.

  * •

visualize_contract_anim_movie(folder_name, re_run_timeseries,
use_re_run_timeseries, keep_thresh=0.75): Visualization of the results of
tracking overlaid on the original movie of the contracting hiPSC-CMs. If
re_run_timeseries $=$ True and use_re_run_timeseries $=$ True this function
will re-run the timeseries_all() function with the new specified keep_thresh.
The saved tracking data will be designated for visualization purposes only.

  * •

cluster_timeseries_plot_dendrogram(folder_name,compute_dist_DTW,
compute_dist_euclidean): Visualize hierarchical clustering of all normalized
sarcomere length time series data and plot a dendrogram that illustrates the
clustering [28]. If compute_dist_DTW$\,=$True the similarity between each
sarcomere time series will be measured with a Dynamic Time Warping algorithm
[40]. Because this can be time consuming (order of minutes to $10$s of
minutes), we also provide an option to set compute_dist_euclidean$\,=$True and
instead base clustering on the Euclidean distance between time series.

  * •

plot_normalized_tracked_timeseries(folder_name): This function will create a
plot of the normalized length of the tracked sarcomeres with respect to frame
number.

  * •

plot_untracked_absolute_timeseries(folder_name): This function will create a
plot of the absolute sarcomere lengths in pixels with respect to frame number
for all segmented sarcomeres. The number of segmented sarcomeres will exceed
the number of tracked sarcomeres.

  * •

compute_timeseries_individual_parameters(folder_name): This function will
compute and save scalar values describing the normalized length time series
data for the tracked sarcomeres (contraction time, relaxation time, flat time,
period, offset, etc.) The results of this analysis will be saved in a
spreadsheet titled “timeseries_parameters_info.xlsx.”

  * •

analyze_J_full_movie(folder_name): This function will compute and save scalar
values describing the $\mathbf{F}_{\mathrm{avg}}$ time series data.
Specifically, it will compute time constants describing the plot of the scalar
Jacobian ($J=|\mathbf{F}_{\mathrm{avg}}|$) with respect to frame number. This
is illustrated in Fig. 3c.

  * •

compare_tracked_untracked(folder_name): This function will compare the tracked
and untracked populations through random sampling of the untracked population.
These plots are important for understanding if the tracked data is
significantly biased compared to the segmented data.

  * •

preliminary_spatial_temporal_correlation_info(folder_name): This function will
perform a preliminary analysis of spatial/temporal correlation within the
hiPSC-CMs. The main outcome is a plot of normalized cross-correlation score
between normalized sarcomere length time series data with respect to both
Euclidean and network distances.

##  3 Results

###  3.1 Synthetic data examples

Our first investigation into the performance of Sarc-Graph on synthetic data
is illustrated in Fig. 5. In this investigation, we define a sinusoidal chain
of $20$ sarcomeres that deform homogeneously following the time series
illustrated in Fig. 2a. Under baseline conditions (modest amplitude, little
angling out of plane, and low to moderate noise), Sarc-Graph is able to
successfully track all $20$ out of $20$ sarcomeres. This is illustrated in 5a.
Then, we adversely alter these baseline conditions to demonstrate the
limitations of our approach. In the first study, illustrated in Fig. 5b, we
alter the underlying 3D geometry of the sarcomere chain. As anticipated,
performance degrades when we increase the amplitude of the sine wave, and when
we increase the out of plane tilt of the sarcomere chain. That being said, we
are able to segment and track at least half of the $20$ sarcomeres in all but
four of the examples shown. In the second study, illustrated in Fig. 5c, we
adversely alter baseline conditions by adding noise to the rendered image.
Specifically, we add Perlin noise at varying levels of magnitude and with a
varying number of octaves [20]. From these results, it is clear that our
computational framework is quite robust to all tested types of low magnitude
noise, and Perlin noise with a low number of octaves (low frequency). However,
performance degrades in the presence of Perlin noise that has high gradients
and is similar in size and brightness to the simulated z-discs themselves. We
note that if a specific experimental dataset has a characteristic noise
pattern that can be removed with a known operation, implementing additional
steps in the segmentation component of our computational framework to address
it would be straightforward. In some cases, simply increasing the parameter
gaussian_filter_size will resolve the issue. In the experimental data that we
have tested the computational framework on thus far, our current segmentation
process is effective.

Figure 5:  For each image, the number of sarcomeres successfully segmented and
tracked (out of $20$ possible) is reported. Running the algorithm on synthetic
data shows that the code is robust to geometry and noise up to a limit: (a)
Illustration of successful z-disc and sarcomere segmentation and tracking for
a curved geometry in the presence of noise; (b) Performance degrades when the
baseline geometry moves partially out of plane (here slanted in the
z-direction), and when the sarcomere chain is too distorted (here high
amplitude to period ratio for sinusoidal geometry); (c) Performance degrades
in the presence of Perlin noise, in particular noise that is similar in size
and brightness to the z-discs themselves. Figure 6:  Example performance on
synthetic data of increasing complexity with a known ground truth:
Segmentation and tracking (marker color corresponds to sarcomere contraction
in the illustrated frame); Constructing a spatial graph (line color
corresponds to sarcomere orientation); Recovering individual time series data
(both Sarc-Graph output and ground truth shown); Recovering average
deformation behavior (both Sarc-Graph output and ground truth shown).

Our second investigation into the performance of Sarc-Graph on synthetic data
is illustrated in Fig. 6. In these examples, our aim is to design synthetic
data that contains many of challenges associated with analyzing real hiPSC-
CMs: out of plane deformation, sarcomere chain curvature, sarcomere chain
overlap, inhomogeneous deformation, variable z-disc size, and irregular
contraction. From S1 to S5 the examples are loosely organized from “least” to
“most” challenging for the computational framework to capture. The
corresponding movies are included in our GitHub repository. Briefly, example
S1 shows four sarcomere chains of variable curvature (decreasing from left to
right) and variable out of plane deformation (increasing from left to right)
deforming homogeneously. Our framework is able to recover $77$ out of $80$
sarcomeres, match the ground truth mean time series data, and match the ground
truth on $\mathbf{\lambda}_{1}^{\mathrm{avg}}$ and
$\mathbf{\lambda}_{2}^{\mathrm{avg}}$. This is consistent with the results
shown in Fig. 5. Example S2 shows two families of sarcomere chains, an
external ellipse and four internal chains that overlap in the center. Our
framework is able to recover $90$ out of $99$ sarcomeres (error occurs
primarily near the overlap region), match the ground truth mean time series
data, and match the ground truth on $\mathbf{\lambda}_{1}^{\mathrm{avg}}$ and
$\mathbf{\lambda}_{2}^{\mathrm{avg}}$.

Example S3 shows three closely positioned elliptical sarcomere chains that
deform inhomogeneously. Our framework is able to recover $34$ out of $40$
sarcomeres, however, because there is bias in which sarcomeres are recovered
(recovery is worse towards the bottom of the image frame), both the mean time
series data and the $\mathbf{\lambda}_{1}^{\mathrm{avg}}$ and
$\mathbf{\lambda}_{2}^{\mathrm{avg}}$ data are unable to perfectly capture the
ground truth. We note that the results for
$\mathbf{\lambda}_{1}^{\mathrm{avg}}$ and
$\mathbf{\lambda}_{2}^{\mathrm{avg}}$ are a better reflection of the ground
truth than the mean time series data. We also note that the amount of
sarcomere contraction simulated in S3 exceeds what is typically observed in
the experimental setting (over $30\%$ contraction).

Example S4 shows two families of sarcomere chains: an external ellipse, and
multiple partially overlapping and curved internal chains in different
z-planes. In Example S4, the sarcomeres deform inhomogeneously with a brief
abrupt jump in sarcomere deformation around frame $70$. Our framework is able
to recover $75$ out of $90$ sarcomeres, and is able to nearly recover the mean
time series data, and the $\mathbf{\lambda}_{1}^{\mathrm{avg}}$ and
$\mathbf{\lambda}_{2}^{\mathrm{avg}}$ data. We note that for both approaches,
the abrupt jump in sarcomere behavior seen in the ground truth is smoothed out
in the automated time series processing because during the abrupt jump
tracking temporarily fails for the entire outer ring of sarcomeres. There are
two main implications of this result. First, that large abrupt motion may
cause Sarc-Graph to fail. Second, that Sarc-Graph is able to track sarcomeres
on either end of a time period where signal is temporarily lost, and
interpolate data across the missing frames. Similar to example S3, the plot of
$\mathbf{\lambda}_{1}^{\mathrm{avg}}$ and
$\mathbf{\lambda}_{2}^{\mathrm{avg}}$ is a better reflection of the ground
truth than the mean time series curve.

Example S5 shows multiple sarcomere chains that are all angled out of plane,
appear to overlap at multiple points, and undergo inhomogeneous deformation.
In this case, we are only able to recover time series data for $27$ out of
$70$ sarcomeres. This results in mean time series data that is a poor
reflection of the ground truth. However, the
$\mathbf{\lambda}_{1}^{\mathrm{avg}}$ and
$\mathbf{\lambda}_{2}^{\mathrm{avg}}$ data is a much better match to the
ground truth. This is a positive sign towards the efficacy of our proposed
metric even when the number of sarcomeres recovered is relatively low. Though
only $27$ out of $70$ sarcomeres were fully tracked, segmented but only
partially tracked sarcomeres appear on the spatial graph shown in Fig. 6 S5.
Looking forward, additional criteria to further refine the spatial graph and
further delineate “probable” and “improbable” links can readily be added to
Sarc-Graph within the spatial_graph module.

| $\tilde{s}$ | $\tilde{s}^{\mathrm{GT}}$ | $s_{\mathrm{avg}}$ | $s_{\mathrm{avg}}^{\mathrm{GT}}$ | OOP | OOPGT | $C_{\mathrm{iso}}$ | $C_{\mathrm{iso}}^{\mathrm{GT}}$ | $C_{\mathrm{||}}$ | $C_{\mathrm{||}}^{\mathrm{GT}}$
---|---|---|---|---|---|---|---|---|---|---
S1 | $0.15$ | $0.15$ | $0.15$ | $0.15$ | $0.62$ | $0.63$ | $0.15$ | $0.15$ | $0.15$ | $0.15$
S2 | $0.085$ | $0.074$ | $0.074$ | $0.075$ | $0.080$ | $0.066$ | $0.076$ | $0.075$ | $0.076$ | $0.075$
S3 | $0.38$ | $0.40$ | $0.079$ | $0.24$ | $0.93$ | $0.88$ | $0.11$ | $0.16$ | $0.039$ | $0.070$
S4 | $0.058$ | $0.043$ | $0.042$ | $0.024$ | $0.16$ | $0.15$ | $0.04$ | $0.04$ | $0.039$ | $0.038$
S5 | $0.22$ | $0.22$ | $0.085$ | $0.042$ | $0.35$ | $0.21$ | $0.13$ | $0.16$ | $0.17$ | $0.20$
Table 2:  Comparison of Sarc-Graph computed $\tilde{s}$, $s_{\mathrm{avg}}$,
OOP, $C_{\mathrm{iso}}$, and $C_{\mathrm{||}}$ with the ground truth values of
these parameters for each synthetic data example.

In Table 2, we report median individual sarcomere shortening ($\tilde{s}$),
mean time series sarcomere shortening ($s_{\mathrm{avg}}$), Orientational
Order Parameter (OOP), equivalent isotropic contraction ($C_{\mathrm{iso}}$),
and contraction in the dominant direction of sarcomere orientation
($C_{\mathrm{||}}$) for each synthetic data example. For all of the synthetic
data cases, we are able to compute a ground truth for each of these parameters
from the known ground truth contraction, displacement, and orientation of the
sarcomere geometry. For examples S1 and S2, the parameters recovered by Sarc-
Graph are a near perfect match to the ground truth. However, for all other
samples minor discrepancies arise in $\tilde{s}$, OOP, $C_{\mathrm{iso}}$, and
$C_{\mathrm{||}}$. For example S4, $\tilde{s}$ is over estimated. For examples
S3 and S5, OOP is over estimated, $C_{\mathrm{iso}}$ is under estimated, and
$C_{\mathrm{||}}$ is under estimated. Consistent with the results shown in
Fig. 6, we believe that this error is due to bias in the region where
sarcomeres are recovered in S3, and failure of the algorithm for overlapping
and severely out of plane sarcomere chains in S4 and S5. We note that for all
examples but S1 and S2 – where sarcomere contraction is synchronous –
recovered $s_{\mathrm{avg}}$ is not a good match to the ground truth. This is
consistent with the normalized sarcomere length time series plot shown in Fig.
6. These errors indicate that care should be taken when comparing results
across beating hiPSC-CM movies that may have systemic differences. The
accessibility of our code for generating additional synthetic data can aid in
this preliminary check.

In all examples, most glaringly example S5, more sarcomeres are segmentated
and tracked for a small number of frames (and thus appear on the spatial
graph) then are tracked for enough frames to reconstruct the full time series.
The recovered time series ($\mathbf{\lambda}_{1}^{\mathrm{avg}}$ and
$\mathbf{\lambda}_{2}^{\mathrm{avg}}$) and parameters ($C_{\mathrm{iso}}$,
$C_{\mathrm{||}}$) associated with $\mathbf{F}_{\mathrm{avg}}$ are typically a
closer match to the synthetic data ground truth behavior than the mean curve
of the recovered time series data. Example S5 – with out of plane
inhomogeneous deformation and overlapping sarcomere chains – is perhaps the
best reflection of the challenges that we have seen in our experience thus far
with the experimental data. We note that the extensible code used to generate
the synthetic data is published on GitHub so implementing and testing
additional examples is straightforward.

Figure 7:  Example performance on experimental data E1 and E2: Raw image
visualization; z-disc and sarcomere segmentation; Sarcomeres that are tracked
for $1/3$ or more of the movie, red corresponds to a contracted state, blue
corresponds to a relaxed state (frame 20); Sarcomeres that are tracked for
$3/4$ or more of the movie, red corresponds to a contracted state, blue
corresponds to a relaxed state (frame 20), note that these sarcomeres are
included in the time series analysis; Individual time series data for
normalized sarcomere length; Average deformation behavior; Illustration of the
data represented as a spatial graph, color corresponds to sarcomere
orientation; Normalized sarcomere time series cross-correlation score plotted
as a function of the distance along the network. Supplementary Movies SI1 and
SI2 are included for further visualization.

###  3.2 Experimental data examples

Figure 8:  Example performance on experimental data E3 and E4: Raw image
visualization; z-disc and sarcomere segmentation; Sarcomeres that are tracked
for $1/3$ or more of the movie, red corresponds to a contracted state, blue
corresponds to a relaxed state (frame 50 for E3, frame 40 for E4); Sarcomeres
that are tracked for $3/4$ or more of the movie, red corresponds to a
contracted state, blue corresponds to a relaxed state (frame 50 for E3, frame
40 for E4), note that these sarcomeres are included in the time series
analysis; Individual time series data for normalized sarcomere length; Average
deformation behavior; Illustration of the data represented as a spatial graph,
color corresponds to sarcomere orientation; Normalized sarcomere time series
cross-correlation score plotted as a function of the distance along the
network. Supplementary Movies SI3 and SI4 are included for further
visualization. Figure 9:  Example performance on experimental data E5: Raw
image visualization; z-disc and sarcomere segmentation; Sarcomeres that are
tracked for $1/3$ or more of the movie, red corresponds to a contracted state,
blue corresponds to a relaxed state (frame 50); Sarcomeres that are tracked
for $3/4$ or more of the movie, red corresponds to a contracted state, blue
corresponds to a relaxed state (frame 50), note that these sarcomeres are
included in the time series analysis; Individual time series data for
normalized sarcomere length; Average deformation behavior; Illustration of the
data represented as a spatial graph, color corresponds to sarcomere
orientation; Normalized sarcomere time series cross-correlation score plotted
as a function of the distance along the network. Supplementary Movie SI5 is
included for further visualization.

Here, we show the performance of Sarc-Graph on the five experimental data
examples listed in Table 1. In Fig. 7 \- 9, we illustrate the key Sarc-Graph
outputs. In Table 3, we report $\tilde{s}$, $s_{\mathrm{avg}}$, OOP,
$C_{\mathrm{avg}}$, and $C_{\mathrm{||}}$ for each example. In the
Supplementary Information section, we include a movie for each sample that
shows both the change in length of the successfully tracked sarcomeres, and
how our computed metric $\mathbf{F}_{\mathrm{avg}}$ changes over the course of
the movie.

In Fig. 7, we show the performance of Sarc-Graph on two similar examples of
beating hiPSC-CM movies, examples E1 and E2. For example E1, we are able to
segment $\approx 500$ sarcomeres per frame and successfully track $74$
sarcomeres. The time series plot of absolute sarcomere length change with
respect to frame number shows three distinct contraction events during the
movie. Though the individual sarcomeres are clearly not perfectly in sync,
there is enough of a unifying pattern for these three peaks to emerge. These
three distinct contraction events are also reflected in the plot of
$\mathbf{\lambda}_{1}^{\mathrm{avg}}$ and
$\mathbf{\lambda}_{2}^{\mathrm{avg}}$. For example E2, we are able to segment
$\approx 500$ sarcomeres per frame and successfully track $60$ sarcomeres. The
time series plots of absolute sarcomere length change and
$\mathbf{\lambda}_{1}^{\mathrm{avg}}$ and
$\mathbf{\lambda}_{2}^{\mathrm{avg}}$ with respect to frame number also show
three distinct peaks. We also show the spatial graph representation of each
example and a plot of normalized cross-correlation between pairs of individual
sarcomere time series curves with respect to network distance.

Qualitatively, examples E1 and E2 are different in that the sarcomeres in
example E1 appear vertically aligned throughout while the sarcomeres in
example E2 appear tangentially aligned around the edge of the cell and
unaligned in the center. This qualitative observation is consistent with the
values of OOP reported in Table 3. In addition, we can use our computed metric
$\mathbf{F}_{\mathrm{avg}}$ to go beyond morphology alone and quantify the
function of the observed hiPSC-CMs. By looking at the time series plots of
$\mathbf{\lambda}_{1}^{\mathrm{avg}}$ and
$\mathbf{\lambda}_{2}^{\mathrm{avg}}$ with respect to frame number, we can see
that example E1 is deforming quite anisotropically
$\mathbf{\lambda}_{1}^{\mathrm{avg}}<\mathbf{\lambda}_{2}^{\mathrm{avg}}$
while the deformation of example E2 is much closer to isotropic where
$\mathbf{\lambda}_{1}^{\mathrm{avg}}\approx\mathbf{\lambda}_{2}^{\mathrm{avg}}$.
In example E1, $C_{\mathrm{||}}$ is meaningfully greater than
$C_{\mathrm{iso}}$ while in example E2 the two parameters are much closer in
value. Essentially, example E1 is experiencing more substantial oriented
contraction than example E2. We also note that for these two examples, where
$\tilde{s}$ is similar but OOP is different, metrics derived from
$\mathbf{F}_{\mathrm{avg}}$ are able to capture functional differences in
contraction.

| $\tilde{s}$ | $s_{\mathrm{avg}}$ | OOP | $C_{\mathrm{iso}}$ | $C_{\mathrm{||}}$
---|---|---|---|---|---
E1 | $0.17$ | $0.049$ | $0.48$ | $0.015$ | $0.026$
E2 | $0.15$ | $0.07$ | $0.29$ | $0.025$ | $0.028$
E3 | $0.18$ | $0.061$ | $0.68$ | $0.069$ | $0.041$
E4 | $0.13$ | $0.041$ | $0.16$ | $0.021$ | $0.027$
E5 | $0.13$ | $0.02$ | $0.38$ | $0.0060$ | $0.0066$
Table 3:  Sarc-Graph computed $\tilde{s}$, $s_{\mathrm{avg}}$, OOP,
$C_{\mathrm{iso}}$, and $C_{\mathrm{||}}$ for each experimental data example.

In Fig. 8, we show the performance of Sarc-Graph on two contrasting examples
of beating hiPSC-CM movies, examples E3 and E4. Example E3 is of a single
hiPSC-CM adhered to a patterned soft hydrogel substrate with highly aligned
sarcomeres and large deformation while Example E4 is of hiPSC-CMs in a
monolayer culture without a clear direction of sarcomere alignment. With Sarc-
Graph, we are able to segment approximately $100$, and track $44$ sarcomeres
from E3. From the plots of $\mathbf{\lambda}_{1}^{\mathrm{avg}}$ and
$\mathbf{\lambda}_{2}^{\mathrm{avg}}$, and scalar metrics OOP,
$C_{\mathrm{iso}}$, and $C_{\mathrm{||}}$, we observe high sarcomere alignment
and corresponding highly aligned contraction. Notably, the cell contracts
substantially in the direction perpendicular to fiber alignment, thus
exhibiting auxetic deformation during contraction. For example E4 we are able
to segment approximately $700$, and track $228$ sarcomeres. This is notable
because the conditions in E4 are far from ideal for image analysis. In both
cases, five distinct contraction events are observed over the course of the
movie and are visible from both the average normalized sarcomere length time
series and the $\mathbf{\lambda}_{1}^{\mathrm{avg}}$ and
$\mathbf{\lambda}_{2}^{\mathrm{avg}}$ time series.

In Fig. 9, we show the performance of Sarc-Graph on example E5. Example E5 is
of a hiPSC-CM on a glass slide with poor fiber alignment, curved fibers, and
relatively low deformation over the course of the movie. With Sarc-Graph, we
are able to segment approximately $800$ and track $539$ sarcomeres for example
E5. Though individual sarcomeres in this movie contract substantially
($\tilde{s}=0.13$), there is low overall synchrony, and overall average
deformation is quite low ($C_{\mathrm{iso}}=0.0060$). Despite low synchrony,
Sarc-Graph still detects $5$ distinct contractions in both the average time
series plot and the $\mathbf{\lambda}_{1}^{\mathrm{avg}}$ and
$\mathbf{\lambda}_{2}^{\mathrm{avg}}$ time series plot. Notably, overall
contraction is also slightly higher ($C_{\mathrm{||}}>C_{\mathrm{iso}}$) in
the dominant direction of sarcomere orientation. This example is notable in
because Sarc-Graph is able to segment and track a high number of sarcomeres
from an irregular geometry and detect subtle sarcomere behavior beyond the
limits of what has been accomplished with the previous state of the art [10].

##  4 Conclusion and outlook

The objective of this work is to provide an open source computational
framework to quantitatively analyze movies of beating hiPSC-CMs. To accomplish
this, we introduce tools to segment and track z-discs and sarcomeres, and
analyze their spatiotemporal behavior. Notably, we are able to automatically
segment and track a high number of sarcomeres across multiple experimental
conditions without any input parameter tuning, and we introduce two important
new approaches for the analysis of beating hiPSC-CM: a method for computing
average deformation gradient, and a method for treating the hiPSC-CMs as
spatial graphs. Looking forward, we see this work as an important tool for
substantial future study of hiPSC-CM behavior. Finally, our proposed
computational framework is heavily documented and has a modular extensible
design. We anticipate that many components of our open source software will be
directly useful to other researchers.

##  5 Additional information

Sarc-Graph, the code to generate the synthetic data, the example data, and
examples of the full analysis output can all be accessed through GitHub:
https://github.com/elejeune11/Sarc-Graph.

##  6 Supplementary Information

SI-1. A movie of E1 with a schematic illustration of deformation gradient
$\mathbf{F}_{\mathrm{avg}}$ and tracked sarcomeres overlayed (sarcomere color
corresponds to contraction level with red indicating the highest level of
contraction), magnitudes of average principal stretches ($\lambda_{1}$,
$\lambda_{2}$ computed from $\mathbf{F}_{\mathrm{avg}}$), and average
normalized sarcomere length, both with respect to the movie frame number. This
movie is a supplement to Fig. 7, sample E1.

SI-2. A movie of E2 with a schematic illustration of deformation gradient
$\mathbf{F}_{\mathrm{avg}}$ and tracked sarcomeres overlayed (sarcomere color
corresponds to contraction level with red indicating the highest level of
contraction), magnitudes of average principal stretches ($\lambda_{1}$,
$\lambda_{2}$ computed from $\mathbf{F}_{\mathrm{avg}}$), and average
normalized sarcomere length, both with respect to the movie frame number. This
movie is a supplement to Fig. 7, sample E2.

SI-3. A movie of E3 with a schematic illustration of deformation gradient
$\mathbf{F}_{\mathrm{avg}}$ and tracked sarcomeres overlayed (sarcomere color
corresponds to contraction level with red indicating the highest level of
contraction), magnitudes of average principal stretches ($\lambda_{1}$,
$\lambda_{2}$ computed from $\mathbf{F}_{\mathrm{avg}}$), and average
normalized sarcomere length, both with respect to the movie frame number. This
movie is a supplement to Fig. 8, sample E3.

SI-4. A movie of E4 with a schematic illustration of deformation gradient
$\mathbf{F}_{\mathrm{avg}}$ and tracked sarcomeres overlayed (sarcomere color
corresponds to contraction level with red indicating the highest level of
contraction), magnitudes of average principal stretches ($\lambda_{1}$,
$\lambda_{2}$ computed from $\mathbf{F}_{\mathrm{avg}}$), and average
normalized sarcomere length, both with respect to the movie frame number. This
movie is a supplement to Fig. 8, sample E4.

SI-5. A movie of E5 with a schematic illustration of deformation gradient
$\mathbf{F}_{\mathrm{avg}}$ and tracked sarcomeres overlayed (sarcomere color
corresponds to contraction level with red indicating the highest level of
contraction), magnitudes of average principal stretches ($\lambda_{1}$,
$\lambda_{2}$ computed from $\mathbf{F}_{\mathrm{avg}}$), and average
normalized sarcomere length, both with respect to the movie frame number. This
movie is a supplement to Fig. 9, sample E5. A still frame from this movie is
also shown in Fig. 1.

##  7 Acknowledgements

This work was supported by the CELL-MET Engineering Research Center NSF
ECC-1647837, the Boston University Hariri Junior Faculty Fellowship (EL), and
the David R. Dalton Career Development Professorship (EL). KZ acknowledges
fellowship support from the American Heart Association, award number
17PRE33660967.

## References

  * [1] Przemek A Gorski, Changwon Kho, and Jae Gyun Oh.  Measuring cardiomyocyte contractility and calcium handling in vitro.  In Experimental Models of Cardiovascular Diseases, pages 93–104. Springer, 2018.
  * [2] Côme Pasqualin, François Gannier, Angèle Yu, Claire O Malécot, Pierre Bredeloux, and Véronique Maupoil.  Sarcoptim for imagej: high-frequency online sarcomere length computing on stimulated cardiomyocytes.  American Journal of Physiology-Cell Physiology, 311(2):C277–C283, 2016.
  * [3] Kevin M Beussman, Marita L Rodriguez, Andrea Leonard, Nikita Taparia, Curtis R Thompson, and Nathan J Sniadecki.  Micropost arrays for measuring stem cell-derived cardiomyocyte contractility.  Methods, 94:43–50, 2016.
  * [4] Matthew Wheelwright, Zaw Win, Jennifer L Mikkila, Kamilah Y Amen, Patrick W Alford, and Joseph M Metzger.  Investigation of human ipsc-derived cardiac myocyte functional maturation by single cell traction force microscopy.  PLoS One, 13(4):e0194909, 2018.
  * [5] Kathleen A Clark, Abigail S McElhinny, Mary C Beckerle, and Carol C Gregorio.  Striated muscle cytoarchitecture: an intricate web of form and function.  Annual review of cell and developmental biology, 18(1):637–706, 2002.
  * [6] Francesco Silvio Pasqualini, Sean Paul Sheehy, Ashutosh Agarwal, Yvonne Aratyn-Schaus, and Kevin Kit Parker.  Structural phenotyping of stem cell-derived cardiomyocytes.  Stem cell reports, 4(3):340–347, 2015.
  * [7] Sean P Sheehy, Francesco Pasqualini, Anna Grosberg, Sung Jin Park, Yvonne Aratyn-Schaus, and Kevin Kit Parker.  Quality metrics for stem cell-derived cardiac myocytes.  Stem cell reports, 2(3):282–294, 2014.
  * [8] Alexandre JS Ribeiro, Olivier Schwab, Mohammad A Mandegar, Yen-Sin Ang, Bruce R Conklin, Deepak Srivastava, and Beth L Pruitt.  Multi-imaging method to assay the contractile mechanical output of micropatterned human ipsc-derived cardiac myocytes.  Circulation research, 120(10):1572–1583, 2017.
  * [9] Thomas Eschenhagen and Lucie Carrier.  Cardiomyopathy phenotypes in human-induced pluripotent stem cell-derived cardiomyocytes—a systematic review.  Pflügers Archiv-European Journal of Physiology, 471(5):755–768, 2019.
  * [10] Christopher N Toepfer, Arun Sharma, Marcelo Cicconet, Amanda C Garfinkel, Michael Mücke, Meraj Neyazi, Jon AL Willcox, Radhika Agarwal, Manuel Schmid, Jyoti Rao, et al.  Sarctrack: an adaptable software tool for efficient large-scale analysis of sarcomere function in hipsc-cardiomyocytes.  Circulation research, 124(8):1172–1183, 2019.
  * [11] Anna Grosberg, Po-Ling Kuo, Chin-Lin Guo, Nicholas A Geisse, Mark-Anthony Bray, William J Adams, Sean P Sheehy, and Kevin Kit Parker.  Self-organization of muscle cell structure and function.  PLoS computational biology, 7(2):e1001088, 2011.
  * [12] Jan David Kijlstra, Dongjian Hu, Nikhil Mittal, Eduardo Kausel, Peter van der Meer, Arman Garakani, and Ibrahim J Domian.  Integrated analysis of contractile kinetics, force generation, and electrical activity in single human stem cell-derived cardiomyocytes.  Stem cell reports, 5(6):1226–1238, 2015.
  * [13] Scott D Lundy, Wei-Zhong Zhu, Michael Regnier, and Michael A Laflamme.  Structural and functional maturation of cardiomyocytes derived from human pluripotent stem cells.  Stem cells and development, 22(14):1991–2002, 2013.
  * [14] Alexandre JS Ribeiro, Yen-Sin Ang, Ji-Dong Fu, Renee N Rivas, Tamer MA Mohamed, Gadryn C Higgs, Deepak Srivastava, and Beth L Pruitt.  Contractility of single cardiomyocytes differentiated from pluripotent stem cells depends on physiological shape and substrate stiffness.  Proceedings of the National Academy of Sciences, 112(41):12705–12710, 2015.
  * [15] Nancy K Drew, Mackenzie A Eagleson, Danny B Baldo Jr, Kevin Kit Parker, and Anna Grosberg.  Metrics for assessing cytoskeletal orientational correlations and consistency.  PLoS computational biology, 11(4):e1004190, 2015.
  * [16] Tessa Altair Morris, Jasmine Naik, Kirby Sinclair Fibben, Xiangduo Kong, Tohru Kiyono, Kyoko Yokomori, and Anna Grosberg.  Striated myocyte structural integrity: Automated analysis of sarcomeric z-discs.  PLoS computational biology, 16(3):e1007676, 2020.
  * [17] Susan E Leggett, Mohak Patel, Thomas M Valentin, Lena Gamboa, Amanda S Khoo, Evelyn Kendall Williams, Christian Franck, and Ian Y Wong.  Mechanophenotyping of 3d multicellular clusters using displacement arrays of rendered tractions.  Proceedings of the National Academy of Sciences, 117(11):5655–5663, 2020.
  * [18] David A Stout, Eyal Bar-Kochba, Jonathan B Estrada, Jennet Toyjanova, Haneesh Kesari, Jonathan S Reichner, and Christian Franck.  Mean deformation metrics for quantifying 3d cell–matrix interactions without requiring information about matrix material properties.  Proceedings of the National Academy of Sciences, 113(11):2898–2903, 2016.
  * [19] Richard A Haddad, Ali N Akansu, et al.  A class of fast gaussian binomial filters for speech and image processing.  IEEE Transactions on Signal Processing, 39(3):723–727, 1991.
  * [20] Ken Perlin.  An image synthesizer.  ACM Siggraph Computer Graphics, 19(3):287–296, 1985.
  * [21] Charles R. Harris, K. Jarrod Millman, Stéfan J. van der Walt, Ralf Gommers, Pauli Virtanen, David Cournapeau, Eric Wieser, Julian Taylor, Sebastian Berg, Nathaniel J. Smith, Robert Kern, Matti Picus, Stephan Hoyer, Marten H. van Kerkwijk, Matthew Brett, Allan Haldane, Jaime Fernández del Río, Mark Wiebe, Pearu Peterson, Pierre Gérard-Marchant, Kevin Sheppard, Tyler Reddy, Warren Weckesser, Hameer Abbasi, Christoph Gohlke, and Travis E. Oliphant.  Array programming with NumPy.  Nature, 585(7825):357–362, September 2020.
  * [22] J. D. Hunter.  Matplotlib: A 2d graphics environment.  Computing in Science & Engineering, 9(3):90–95, 2007.
  * [23] Anant Chopra, Matthew L Kutys, Kehan Zhang, William J Polacheck, Calvin C Sheng, Rebeccah J Luu, Jeroen Eyckmans, J Travis Hinson, Jonathan G Seidman, Christine E Seidman, et al.  Force generation via $\beta$-cardiac myosin, titin, and $\alpha$-actinin drives cardiac sarcomere assembly from cell-matrix adhesions.  Developmental cell, 44(1):87–96, 2018.
  * [24] Arun Sharma, Christopher N Toepfer, Tarsha Ward, Lauren Wasson, Radhika Agarwal, David A Conner, Johnny H Hu, and Christine E Seidman.  Crispr/cas9-mediated fluorescent tagging of endogenous proteins in human pluripotent stem cells.  Current protocols in human genetics, 96(1):21–11, 2018.
  * [25] G. Bradski.  The OpenCV Library.  Dr. Dobb’s Journal of Software Tools, 2000.
  * [26] Nobuyuki Otsu.  A threshold selection method from gray-level histograms.  IEEE transactions on systems, man, and cybernetics, 9(1):62–66, 1979.
  * [27] D Allan, C van der Wel, N Keim, TA Caswell, D Wieker, R Verweij, C Reid, T Bottaro, L Grueter, K Ramos, apiszcz, zoeith, RW Perry, F Boulogne, P Sinha, pfigliozzi, N Bruot, L Uieda, J Katins, H Mary, and A Ahmadia.  soft-matter/trackpy: Trackpy v0. 4.2 (2019).  URL https://doi. org/10.5281/zenodo, 3492186.
  * [28] Pauli Virtanen, Ralf Gommers, Travis E. Oliphant, Matt Haberland, Tyler Reddy, David Cournapeau, Evgeni Burovski, Pearu Peterson, Warren Weckesser, Jonathan Bright, Stéfan J. van der Walt, Matthew Brett, Joshua Wilson, K. Jarrod Millman, Nikolay Mayorov, Andrew R. J. Nelson, Eric Jones, Robert Kern, Eric Larson, C J Carey, İlhan Polat, Yu Feng, Eric W. Moore, Jake VanderPlas, Denis Laxalde, Josef Perktold, Robert Cimrman, Ian Henriksen, E. A. Quintero, Charles R. Harris, Anne M. Archibald, Antônio H. Ribeiro, Fabian Pedregosa, Paul van Mulbregt, and SciPy 1.0 Contributors.  SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python.  Nature Methods, 17:261–272, 2020.
  * [29] Stefan Van der Walt, Johannes L Schönberger, Juan Nunez-Iglesias, François Boulogne, Joshua D Warner, Neil Yager, Emmanuelle Gouillart, and Tony Yu.  scikit-image: image processing in python.  PeerJ, 2:e453, 2014.
  * [30] John C Crocker and David G Grier.  Methods of digital video microscopy for colloidal studies.  Journal of colloid and interface science, 179(1):298–310, 1996\.
  * [31] The pandas development team.  pandas-dev/pandas: Pandas, February 2020.
  * [32] Carl Edward Rasmussen and Christopher KI Williams.  Gaussian process regression for machine learning, 2006.
  * [33] Fabian Pedregosa, Gaël Varoquaux, Alexandre Gramfort, Vincent Michel, Bertrand Thirion, Olivier Grisel, Mathieu Blondel, Peter Prettenhofer, Ron Weiss, Vincent Dubourg, et al.  Scikit-learn: Machine learning in python.  the Journal of machine Learning research, 12:2825–2830, 2011.
  * [34] Anne E Carpenter, Thouis R Jones, Michael R Lamprecht, Colin Clarke, In Han Kang, Ola Friman, David A Guertin, Joo Han Chang, Robert A Lindquist, Jason Moffat, et al.  Cellprofiler: image analysis software for identifying and quantifying cell phenotypes.  Genome biology, 7(10):R100, 2006.
  * [35] Emma Lejeune, Berkin Dortdivanlioglu, Ellen Kuhl, and Christian Linder.  Understanding the mechanical link between oriented cell division and cerebellar morphogenesis.  Soft matter, 15(10):2204–2215, 2019.
  * [36] Emma Lejeune and Christian Linder.  Quantifying the relationship between cell division angle and morphogenesis through computational modeling.  Journal of theoretical biology, 418:1–7, 2017.
  * [37] Gerhard A Holzapfel.  Nonlinear solid mechanics a continuum approach for engineering, 2000.
  * [38] Akinori Umeno, Hiroko Kotani, Masakazu Iwasaka, and Shoogo Ueno.  Quantification of adherent cell orientation and morphology under strong magnetic fields.  IEEE transactions on magnetics, 37(4):2909–2911, 2001.
  * [39] Aric A. Hagberg, Daniel A. Schult, and Pieter J. Swart.  Exploring network structure, dynamics, and function using networkx.  In Gaël Varoquaux, Travis Vaught, and Jarrod Millman, editors, Proceedings of the 7th Python in Science Conference, pages 11 – 15, Pasadena, CA USA, 2008.
  * [40] Pavel Senin.  Dynamic time warping algorithm review.  Information and Computer Science Department University of Hawaii at Manoa Honolulu, USA, 855(1-23):40, 2008.