# A Data-Driven Computational Model for Engineered Cardiac Microtissues

Javiera Jilberto<EMAIL_ADDRESS>Samuel J. DePalma<EMAIL_ADDRESS>Jason
Lo<EMAIL_ADDRESS>Hiba Kobeissi<EMAIL_ADDRESS>Lani Quach
<EMAIL_ADDRESS>Emma Lejeune<EMAIL_ADDRESS>Brendon M. Baker
<EMAIL_ADDRESS>David Nordsletten<EMAIL_ADDRESS>Department of
Biomedical Engineering, University of Michigan, MI, USA  Department of Cardiac
Surgery, University of Michigan, MI, USA  Department of Biomedical
Engineering, School of Imaging Sciences and Biomedical Engineering, King’s
College London, London, UK  Department of Mechanical Engineering, Boston
University, MA, USA

###### Abstract

Engineered heart tissues (EHTs) present a potential solution to some of the
current challenges in the treatment of heart disease; however, the development
of mature, adult-like cardiac tissues remains elusive. Mechanical stimuli have
been observed to improve whole-tissue function and cardiomyocyte (CM)
maturation, although our ability to fully utilize these mechanisms is
hampered, in part, by our incomplete understanding of the mechanobiology of
EHTs. In this work, we leverage the experimental data produced by a
mechanically tunable experimental setup to generate tissue-specific
computational models of EHTs. Using imaging and functional data, our modeling
pipeline generates models with tissue-specific ECM and myofibril structure,
allowing us to estimate CM active stress. We use this experimental and
modeling pipeline to study different mechanical environments, where we
contrast the force output of the tissue with the computed active stress of
CMs. We show that the significant differences in measured experimental forces
can largely be explained by the levels of myofibril formation achieved by the
CMs in the distinct mechanical environments, with active stress showing more
muted variations across conditions. The presented model also enables us to
dissect the relative contributions of myofibrils and extracellular matrix to
tissue force output, a task difficult to address experimentally. These results
highlight the importance of tissue-specific modeling to augment EHT
experiments, providing deeper insights into the mechanobiology driving EHT
function.

###### keywords:

computational modeling , engineered heart tissues , cardiac biomechanics ,
mechanobiology

††journal: Acta Biomaterialia

##  1 Introduction

The development of engineered heart tissues (EHTs) for use in regenerative
therapies, drug testing, and disease modeling has the potential to improve the
life expectancy of millions of people that suffer from cardiac disease [1].
Current state-of-the-art EHTs are manufactured using cardiomyocytes (CMs)
derived from induced pluripotent stem cells (iPSCs) and scaffolds that mimic
the structure and mechanics of native cardiac tissue [2]. However, iPSC-CM
maturation in current EHT platforms remains a challenge [3]. This is evidenced
by the lack of hallmark attributes of mature CMs, such as myofibril alignment,
protein expression, calcium handling, and electrophysiological response, among
others [4]. Biophysical stimuli, such as electrical pacing [5, 6] and
mechanical loading [7, 8, 9], have been shown to enhance iPSC-CM maturation in
different EHT platforms. However, the underlying mechanisms driving these
observations remain incompletely understood [10], preventing scientists from
efficiently optimizing the application of these techniques.

Several in-vitro studies have attempted to elucidate the impact of different
mechanical microenvironmental perturbations on the maturity of iPSC-CMs. For
example, Leonard et al. [8] showed that increasing the resistance against
which the EHTs contract increases force generation. Similarly, Bliley et al.
[9] showed that EHTs grown under passive stretch evolve into more mature
tissues. Furthermore, other studies have shown that culturing CMs on hydrogels
with a modulus similar to that of the healthy adult myocardium enhances
electro-mechanical activity [11, 12]. DePalma et al. expanded upon this,
showing that iPSC-CMs adapt their contractile behavior in response to ECM
mechanics on synthetic, fibrous matrices that better mimic the anisotropic
mechanics of the cardiac ECM [7]. This and other work by Allen et al. [13]
also showed that increasing the anisotropy of the fibrous matrix results in
more aligned myofibrils. These studies highlight the impact of different
mechanical cues on EHT formation and function. However, the inherent
variability of EHTs – from their formation, maturation, ECM properties and
structure, myofibril formation, etc. – cloud our interpretation of each
mechanical parameter’s relevance and the underlying mechanobiological drivers.

A way of tackling this complex task and providing insight into the
multifaceted differences in EHTs is through biomechanical modeling.
Biomechanical models have been used to understand the mechanics of in-vitro
systems, helping decipher the mechanobiology behind the alignment of cells
[14, 15], the biomechanics of microtissue failure due to necking [16], the
cell force transmission through fibrous substrates [17, 18], and mechanics at
cell-to-cell junctions [19]. These models enable the examination of
experimental conditions – as well as exploration through in silico testing
whereby scenarios that would be virtually impossible to construct
experimentally can be studied. Biomechanical models – particularly at the
whole-organ level – have further enabled the integration and assimilation of
imaging and patient data [20, 21, 22, 23], providing pipelines for generating
_patient-specific models_ describing cardiac function. These pipelines enable
the integration of realistic structure/function and provide a platform for
understanding the significance of these factors on localized mechanics, such
as strain or stress.

To bypass many of the uncertainties associated with EHTs and delve into the
underlying impact of structure and function, we propose the integration of a
novel EHT platform with _tissue-specific_ computational models. In this study,
we leverage the experimental data obtained in the fibroTUG platform [24],
which uses electrospinning techniques to generate synthetic, fibrous matrices
with defined alignment and material characteristics. Imaging provides detailed
information on the ECM and myofiber architecture on a tissue-by-tissue basis.
Using image-based modeling and data assimilation, we create in-silico twins of
individual fibroTUGs and explore how different factors of the ECM and
myofibril structure impact localized cellular function and the resultant force
measures commonly reported in the literature. Through this integrated
experimental-computational approach, we observe that the resultant forces of
EHTs can be substantially biased by ECM alignment and mechanobiological
factors driving myofibril formation and function.

The paper is structured as follows: in the methods section, we detail the
process of creating the computational model from the experimental data. In the
results section, we validate our model by comparing the simulation results to
image-based algorithms and use experimental and non-experimental conditions to
analyze the role of the different mechanical variables in the iPSC-CMs stress
and force output relationship in EHTs. We follow with a discussion of the
model results and how we can use computational modeling approaches to explore
the mechanobiology of these tissues.

##  2 Materials and Methods

###  2.1 Experimental Setup

FibroTUG microtissues were fabricated as described in DePalma et. al (2023)
[24] (see also Fig. 1A). Briefly, dextran vinyl sulfone (DVS) fibers were
electrospun onto an array of polydimethylsiloxane (PDMS) posts attached to a
rotating mandrel. The bending stiffness of the posts $k_{\rm p}$ can be tuned
by altering the geometry of the posts and is measured experimentally, as
described in DePalma et al. (2023) [24]. By varying the rotational velocity of
the mandrel, the alignment of the fibers can also be controlled, resulting in
aligned fibrous matrices that exhibit high anisotropy or random matrices that
are more locally isotropic [7]. The DVS fibers present between posts were
stabilized by exposing them to UV light. Upon hydration, the unstabilized
fibers dissolve, leaving only the fibers suspended between two posts.
Secondary crosslinking in solutions with varying concentrations of
photoinitiator (LAP) resulted in matrices of varying stiffness (higher LAP
results in stiffer matrices). The stiffness of the matrix was characterized by
indenting the matrices incrementally, measuring the resulting post-deflection
(and thus the post-force $F_{\rm p}$), and then calculating the global strain
of the matrix (Fig. 1B). Further, images of the indented DVS fibers were taken
to record their organization and pair it with their force response. This setup
allows us to tune and control the post stiffness (soft, $k_{\rm p}=0.41$ N/m;
stiff, $k_{\rm p}=1.2$ N/m), fiber alignment (aligned or random), and fibrous
matrix stiffness (soft, LAP$=0.1$ mg/mL; stiff, LAP$=5.0$ mg/mL), parameters
that determine the mechanical environment where the iPSC-CMs develop [24].
While many permutations are possible, in this paper, we studied the following
permutations: soft fibers/soft posts, stiff fibers/soft posts, and soft
fibers/stiff posts for both aligned and random matrices, leading to a total of
six conditions.

After defining matrix conditions, purified cultures of iPSC-CMs were seeded
onto the fibroTUG and cultured in this environment for seven days. On day 7,
time-lapse videos of the microtissue’s spontaneous contractions were acquired
(see Video S1). These videos were processed to obtain post-displacement curves
as seen in Fig. 1C. Finally, immunofluorescence staining is used to image cell
nuclei, titin, and the DVS fibers (see Fig. 1D).

###  2.2 Image Processing

This subsection details the process of extracting information from the DVS
fiber and titin images (Fig. 2) to obtain quantities that characterize the
structure of the matrix and the myofibril network that can then be projected
into a 2D fibroTUG model.

DVS fibers. The processing starts by creating a mask of the fibers that is
used to compute the local fiber density $\rho_{\rm f}$, alignment
$\mathbf{f}_{0}$, and dispersion $\kappa_{\rm{f}}$ as shown in Fig. 2A. The
methods used to compute these quantities are detailed in Supplementary
Information Section S1. The fiber density $\rho_{f}$ takes values between 0
(no fibers) and 1 (fiber), allowing us to define the mechanical presence of
fibers. The fiber alignment vector, $\mathbf{f}_{0}$, provides the direction
of the local stiffness anisotropy in the direction of the fibers. Finally, the
dispersion $\kappa_{\rm f}$ is a parameter used in continuum models [25] to
represent regions where fibers follow a distribution around a mean vector. In
our case, it allows the local stiffness to move from anisotropic (no
dispersion, $\kappa_{\rm{f}}=0$) to isotropic ($\kappa_{\rm{f}}=0.5$)
depending on the local fibers.

Titin. Myofibrils were identified via immunofluorescent imaging of iPSC-CMs
containing a titin-GFP reporter, which allows for the visualization of the
sarcomeres’ Z-discs [26]. When images of titin were available for the whole
domain, tissue-specific fields describing the structure of the myofibril
network were computed. This was done following a similar strategy to the
processing of the DVS fibers, where the titin images were used to find the
myofibril density $\rho_{\rm m}$, alignment $\mathbf{m}_{0}$, and dispersion
$\mathbf{\kappa}_{\rm m}$ (Fig. 2B). The steps to compute these quantities
from the images are presented in Supplementary Information Section S2. The
myofibril density $\rho_{\rm m}$ (taking values from 0 to 1) allows us to
define the contractile regions. The alignment vector $\mathbf{m}_{0}$ defines
the direction of the contraction and the dispersion $\kappa_{\rm m}$ activates
an isotropic contraction in regions of unorganized myofibrils.

To obtain the visualization of the titin in the whole tissue shown in Fig. 1D,
several smaller images were stitched together since a high magnification is
needed to have a clear view of the Z-discs. This is a slow process, and since
the main objective was to quantify myofibril structural characteristics, we
decided to accelerate the process by only imaging the center of the tissues
and statistically characterize the myofibril alignment and density. Given the
six different experimental conditions, the myofibril orientation was analyzed
and compiled for all the images available ($N\geq 7$ per condition). A Von
Mises distribution [25] was fit to the histogram of angles (measured relative
to the post-to-post direction, see Supplementary Information Section S2.1 for
more details). The distribution is characterized by a parameter $\xi$
(assuming the mean is zero), with high values of $\xi$ indicating myofibrils
oriented preferentially in the post-to-post direction. The resulting
probability density functions (PDFs) and the data histograms are shown in Fig.
3A. The myofibril density was computed for these partial images and then a
mean value was computed per condition (see Fig. 3B). The mean density values
were then normalized by the aligned, soft fibers/soft post condition value.
Table 1 shows the final density values used. Given a fiber network and
geometry, we used the probability fitting and the density measurements
specific for each condition to generate computational myofibril fields that
followed these two parameters on top of the fiber fields. The procedure to
create these fields is shown in Supplementary Information Section S3. One
important thing to notice is that whenever this approach was taken, no
myofibril dispersion was considered ($\kappa_{\rm m}=0$), as it is difficult
to generate computationally in a meaningful way. The impact of these
considerations was studied in Section S3.1 and S3.2 of the Supplementary
Information.

###  2.3 Biomechanical Model

To model microtissue mechanics, we used a constrained mixture continuum
mechanics framework [27, 28, 29]. Due to the thinness of the fibroTUG tissues
($\sim 12\mu m$ compared to the $\sim 400\mu m$ in length), we consider the
tissue domain $\Omega\subset\mbox{$\mathbb{R}$}^{2}$. The boundary at the post
is denoted by $\Gamma_{p}$. The reference coordinates of a given point in
$\Omega$ are denoted by $\mathbf{X}$. Under internal and external loads, this
point moves to a deformed position $\mathbf{x}=\mathbf{X}+\mathbf{u}$, where
$\mathbf{u}$ is the displacement field. The deformation gradient tensor
$\mathbf{F}=\nabla\mathbf{u}+\mathbf{I}$ describes the deformation of the
material with respect to the reference coordinates and $J=\det{\mathbf{F}}$
the relative volume change [28]. We further define
$\mathbf{C}=\mathbf{F}^{T}\mathbf{F}$ as the right Green Cauchy deformation
tensors. Constitutive relations are often expressed in terms of the invariants
of $\mathbf{C}$ [30]. In this work, we considered the following,

$I_{1}=\mathbf{C}:\mathbf{I},\qquad\bar{I}_{1}=J^{-2/3}I_{1},$

where $\bar{I}_{1}$ is the isochoric version of $I_{1}$ [31]. Furthermore,
invariants describing the deformation in the fiber and myofibril direction
are,

$\qquad I_{\rm 4f}=\mathbf{C}:\mathbf{f}_{0}\otimes\mathbf{f}_{0},\qquad
I_{\rm 4m}=\mathbf{C}:\mathbf{m}_{0}\otimes\mathbf{m}_{0},$

and these are further modified to account for local dispersion $\kappa_{\rm
f}$, $\kappa_{\rm m}$ [25, 32],

$\displaystyle I_{\rm 4f}^{*}=\kappa_{\rm f}I_{1}+(1-2\kappa_{\rm f})I_{\rm
4f},\qquad I_{\rm 4m}^{*}=\kappa_{\rm m}I_{1}+(1-2\kappa_{\rm m})I_{\rm 4m}.$

When the fiber dispersion $\kappa_{\rm f}=0$, $I_{4\rm f}^{*}=I_{4\rm f}$,
i.e., the local response of $I_{4\rm f}^{*}$ is fully anisotropic, whereas
when $\kappa_{\rm f}=0.5$, $I_{4\rm f}^{*}=I_{1}/2$, and the material behaves
locally as a fully isotropic material. The same is true for the myofibril
invariant $I_{4\rm m}^{*}$.

As mentioned in Section 2.1, a fibroTUG tissue consists of two components, DVS
fibers and iPSC-CMs. We assumed the two components work in parallel and,
therefore, the strain energy density of the tissue $\Psi$ is given by the sum
of the energy of the fibers $\Psi_{ECM}$ and the iPSC-CMs $\Psi_{CM}$,

$\Psi=\Psi_{ECM}+\Psi_{CM},$ (1)

The strain energy function $\Psi_{\rm ECM}$ is given by the strain energy
density describing the fiber mechanics $\Psi_{\rm f}$ and a term that delivers
numerical stability on the areas where there are no fibers $\Psi_{\rm st}$,

$\displaystyle\Psi_{ECM}=\Psi_{\rm f}(\mathbf{C};\rho_{\rm
f},\mathbf{f}_{0},\kappa_{\rm f})+\Psi_{\rm st}(\mathbf{C}).$

The first term corresponds to a modified neofiber material law [33] and
integrates the structural information obtained from the DVS fiber images,

$\displaystyle\Psi_{\rm f}(\mathbf{C};\rho_{\rm f},\mathbf{f}_{0},\kappa_{\rm
f})=\rho_{\rm f}\left(\frac{C_{1}}{4}(I_{4\rm
f}^{*}-1)^{2}+\frac{C_{2}}{2}(\bar{I}_{1}-2)\right),$ (2)

where $C_{1}$ is the stiffness in the fiber direction and $C_{2}$ is the
isotropic stiffness. This formulation was chosen since the experimental data
from the passive stretching of the fibers (Fig. 1B) showed a close-to-linear
behavior that is well captured by this material law. The stabilization term
is,

$\Psi_{\rm st}=\frac{K}{2}\left[(J-1)^{2}+(\ln
J)^{2}\right]+\mu\operatorname{tr}(\mathbf{E}^{2}).$ (3)

The first term inside the parenthesis penalizes volumetric changes, while the
second corresponds to the deviatoric term of a Saint-Venant Kirchhoff
material. The parameters $K=10^{-3}$ kPa and $\mu=10^{-2}$ kPa are chosen to
be $\ll C_{1}$, so the mechanical response of the ECM is dominated by
$\Psi_{\rm f}$. A sensitivity analysis confirmed that these parameters are
minimally important (see Section S6 of the Supplementary Information for more
details).

The contractile iPSC-CMs were modeled by a passive component representing the
bulk stiffness of the cells ($\Psi_{\rm c}$) and an active contraction
component ($\Psi_{\rm a}$),

$\Psi_{\rm CM}=\Psi_{\rm c}+\Psi_{\rm a}.$ (4)

The passive component is modeled using the 2D version of the compressible
Neohookean material law presented in Pascon et al. (2019) [34], which has the
following strain energy density function,

$\Psi_{\rm c}=\frac{K_{c}}{2}[\ln J]^{2}+\frac{\mu_{\rm
c}}{2}\left(I_{1}-2-2\ln J\right).$ (5)

We use $\mu_{\rm c}=2$ kPa, which is close to values derived from stress-
strain curves obtained from isolated iPSC-CMs [35]. Due to the 2D nature of
the EHTs, the compressible modulus was assumed negligible (i.e. $K_{c}=0$) due
to the ease for cells to deform out-of-plane under in-plane compression.

The active component is given by,

$\Psi_{a}=\int_{0}^{I_{\rm 4m}^{*}}\eta\rho_{\rm m}\phi(s)\;{\rm d}s.$ (6)

Here, $\phi$ is a function that models the length-dependent behavior of
cardiomyocytes, and it is taken from [36] (see Supplementary Information
Section S4 for more details). The parameter $\eta$ controls the magnitude of
the activation in time of the iPSC-CMs. Note that we assume the activation
only occurs where there are myofibrils (hence the multiplication by $\rho_{\rm
m}$), but the passive response of the iPSC-CMs, is present everywhere in the
tissue. This is because iPSC-CMs nuclei appear evenly distributed in the
images of the tissues, but only some develop myofibrils.

The total Cauchy stress $\bm{\sigma}$ is computed from $\Psi$ using,

$\bm{\sigma}=J^{-1}\frac{\partial\Psi}{\partial\mathbf{F}}\mathbf{F}^{T},$ (7)

and, under a quasi-static regime, the stress balance is given by,

$\nabla\cdot\bm{\sigma}=0.$ (8)

Throughout the paper, we assess the mechanics of the tissues using the active
stress and the strain in the myofibril direction,

$\sigma_{\rm a,m}=\bm{\sigma}_{\rm
a}:\mathbf{m}\otimes\mathbf{m},\qquad\varepsilon_{\rm
m}=\frac{1}{2}\left(\mathbf{C}-\mathbf{I}\right):\mathbf{m}\otimes\mathbf{m}.$
(9)

where $\bm{\sigma}_{\rm a}=J^{-1}\frac{\partial\Psi_{\rm
a}}{\partial\mathbf{F}}\mathbf{F}^{T}$.

###  2.4 Data Assimilation

Once the structure of each tissue is defined, the only unknowns remaining in
the system are the material parameters of the constitutive law of the fibers,
$C_{1}$ and $C_{2}$, and of the active component of the iPSC-CMs, $\eta$. For
simplicity, we considered these parameters to be constant within a fibroTUG.
To find these tissue-specific values, we assimilate the functional data shown
in Fig. 1C-D using a parameter identification strategy introduced by [20].
Briefly, this technique integrates material parameters into the set of state
variables, enabling the addition of constraints to match measured data. This
method allows solving for both the displacement field and additional material
parameters in the same forward simulation. In our case, we divide the
assimilation into two steps. First, we use the data from the fiber indentation
test (Fig. 1C) to find the stiffness parameters of the fibers, and then we use
the post-displacement trace (Fig. 1D) to find the active parameter, $\eta$.
The corresponding boundary value problem equations are described in Section S5
of the Supplementary Information.

##  3 Results

###  3.1 Model validation

The objective of this first set of simulations was to assess the ability of
the proposed pipeline and model to capture fibroTUG mechanics. To do so, three
complete datasets of the soft fibers/soft post condition were analyzed. The
fiber stiffness was set to the mean value found for the fibrous samples of
this condition, $C_{1}=3.84$ kPa (see Section 3.2). By design, the post-
displacement (and, therefore, the post-force) data were exactly matched by the
simulations as shown in Fig. 4A. The data assimilation enabled us to identify
the parameter $\eta$ that reproduces the post-displacement and force curves.
Fig. 4B shows the mean $\sigma_{\rm a,m}$ at each time point for each tissue.
To assess the ability of the model to capture local deformations, we validated
the results of the simulations against the results from MicroBundleCompute, an
image tracking software developed specifically for measuring the internal
displacements of contracting EHTs [37]. For the three tissues, correlation
plots were calculated between the predicted displacements of the simulation
and the measured displacements from the video. Fig. 4C and D show the measured
and predicted displacement fields in the post-to-post direction (Y) and the
perpendicular direction (X) for one of the three datasets. Fig. 4E and F show
the displacement’s correlation in Y and X, respectively. The $R^{2}$ parameter
for the correlation in the Y direction is always greater than 0.89, while in
the X direction, a positive correlation is observed, with $R^{2}$ values close
to 0.4.

###  3.2 Active stress prediction for experimental conditions

In this section, we focus on understanding the effect of different mechanical
environments on the iPSC-CMs ability to exert contractile stress. Images of
five fibrous matrices for each condition were processed as described in 2.2
(with the exception of the random, soft fibers/soft post, where only four
matrices with the paired indentation test were available). The fiber stiffness
for each available matrix was identified using the indentation experimental
data as detailed in 2.4. As expected, the mean fiber stiffness for the stiff
matrix cases (mean $C_{1}=19.47$/$C_{1}=24.66$ for aligned/random matrices) is
about four times stiffer than the soft matrix case (mean
$C_{1}=3.84$/$C_{1}=4.26$ for aligned/random matrices). More detailed results
of this process can be found in Supplementary Information Section S8. We
computationally generated five myofibril fields on top of each matrix
following the probabilistic characterization shown in Fig. 3. Further, five
representative samples (i.e., that have a similar mean and standard deviation
than the whole set) of post-displacement curves were used as input in the
active data assimilation step. This gives us $N=125$ in-silico models per
experimental condition (except for the random, soft fibers/soft post condition
where $N=100$).

Fig. 5A-C shows the results for aligned matrices, including bar plots for the
post-force and the mean active stress $\sigma_{\rm a,m}$, and the $\sigma_{\rm
a,m}$ field at maximum contraction for aligned matrices. Fig. 5D-F shows the
same for random matrices. Both $F_{\rm p}$ and $\sigma_{\rm a,m}$ follow the
same trends when comparing condition to condition. For instance, in the
aligned case, the magnitude of both $F_{\rm p}$ and $\sigma_{\rm a,m}$ are
highest in the soft fibers/soft post-condition and lowest in the case with
stiff fibers. However, the relative differences between conditions in
$\sigma_{\rm a,m}$ are lower compared to $F_{\rm p}$. For example, the force
output drops from $3.01\,\mu$N in the soft fibers/soft post case (grey bars in
5A) to only $1.19\,\mu$N in the case with stiffer fibers (red bars in 5A),
which corresponds to a -60.5% relative change. Contrarily, the maximum
$\sigma_{\rm a,m}$ falls from 2.45 to 1.55 kPa between these two cases, only a
-36.8% variation.

###  3.3 Isolating the effect of mechanical variations

To investigate the origin of the $F_{\rm p}/\sigma_{\rm a,m}$ relative
differences, we performed several simulations outside the parameter range of
current experimental conditions. Specifically, we studied variations in post
stiffness, fiber stiffness, fiber alignment, and myofibril density. To do so,
we took a single post-force curve (Fig. 6A) and computed the active stress
needed to generate that force output given a certain parameter set. In this
experiment, higher values of $\sigma_{\rm a,m}$ indicate that the tissue is
less efficient in transmitting cell stress to force output as this was held
constant across all conditions. We performed this test using five random
matrices and five aligned matrices with five computationally generated
myofibril fields each (Fig. 6B). The alignment of myofibrils was taken from
the soft post/soft fibers case (gray curves in Fig. 3A). The stiffness in the
fiber direction was set to be $C_{1}=3.84$ and $C_{1}=19.47$ kPa, for the soft
fiber and the stiff fiber case, respectively.

First, we considered a constant and uniform myofibril density $\rho_{\rm
m}=1$. Results of this experiment are shown in Fig. 6C, for aligned matrices
and Fig. 6E, for random matrices. On aligned matrices, the active stress
necessary to generate the force output represented in Fig. 6A was
$2.12,2.60,1.95$ kPa for the soft fibers/soft post, stiff fibers, and stiff
post conditions, respectively. Similar relative trends are observed for the
random case, although with higher values compared to their aligned
counterparts. Second, to study the effect of myofibril density, we considered
the case where $\rho_{\rm m}=\rho_{\rm data}$, which is the specific density
computed from the data (shown in Table 1). Fig. 6D shows the results of these
simulations. The magnitude of $\sigma_{\rm a,m}$ is elevated in those cases
with lower $\rho_{\rm m}$, when compared with the $\rho_{\rm m}=1$ simulations
for both aligned and random matrices.

Finally, to understand the reason behind the performance of random matrices
compared to aligned matrices, we investigated the relative importance of
myofibril alignment and fiber alignment in the transmission of iPSC-CM stress
to post-force. To do so, we simulated mixed, synthetic scenarios where aligned
myofibril fields ($\xi=3.07$) were generated on top of random fiber fields
(Rf-Am) and vice-versa, random myofibrils ($\xi=1.67$) on aligned fibers (Af-
Rm) as shown in Fig. 7A. These cases were compared with the experimental
scenarios, aligned fibers with aligned myofibrils (Af-Am) and random fibers
with random myofibrils (Rf-Rm), which correspond to the bar plots shown in
Fig. 6C,E. We performed the same experiment as in the previous section with
these two mixed cases. Results for the soft fibers/soft post case are shown in
Fig. 7B. The results for the other conditions (with stiffer fibers or stiffer
posts) are very similar and are shown in Supplementary Information Section S9.
The difference in $\sigma_{\rm a,m}$ between the experimental conditions Af-Am
and Rf-Rm is 0.846 kPa. When random myofibrils are used on top of aligned
fibers (Af-Rm), that difference is reduced to 0.731 kPa (-13.6%). When
myofibrils are aligned in random matrices (Rf-Am) $\sigma_{\rm a,m}$ the
reduction is similar (-16.7%). The highest variation is produced by changing
the matrix alignment. For example, the difference is reduced by 86.4% when we
fixed the myofibril alignment to be aligned and changed the matrix alignment
(i.e., from Af-Am to Rf-Am).

##  4 Discussion

In this work, we developed a pipeline to generate data-driven computational
models of EHTs. By combining comprehensive experimental data with
biomechanical computational models, we are able to model EHT mechanics and
compute metrics to estimate iPSC-CM function. The pipeline developed allows us
to model the explicit fiber structure and myofibril network and to infer the
mechanical properties of these components from functional experimental data.
This is facilitated by the use of the fibroTUG platform, which provides great
control over the mechanical environment and enables obtaining detailed imaging
and functional data [24]. The biomechanical computational model augments the
analysis and conclusions obtained from the experiments. It also enables the
test of conditions that are experimentally not feasible to achieve, allowing
us to decipher the effect of the different variables of the system in the
tissue mechanics. Below, we discuss the results of the different simulations
performed using the model.

###  4.1 Model Validation

For the three validation datasets, the correlation between the displacements
predicted by the simulations and the measurements of the MicroBundleCompute
software in the post-to-post direction is strong, with high $R^{2}$ values
(mean 0.938) and slopes close to 1 (0.967, 1.047, 0.970 for the three tissues,
see Fig. 4E). The correlations in the perpendicular direction are also
positive but with slopes further from 1 (1.495, 0.694, 0.623 for the three
tissues) and a mean $R^{2}=0.372$ (see Fig. 4F). The decrease in correlation
in the X directions is explained by several reasons. First and foremost, the
videos are taken using brightfield imaging [24], which mainly shows the
fiber’s deformations. Our modeling approach uses a constrained mixture
framework, where the kinematics of fibers and cells are homogenized, meaning
that the computed displacements represent the average displacement of both
these components. Furthermore, the displacements in the X direction are
smaller and harder to track as demonstrated by a simple exercise where two
people track different features across the contraction cycle. This test showed
a decrease in the measured displacement correlation, from $R^{2}=0.9$ in the Y
direction to $R^{2}=0.72$ in the X direction (see Supplementary Information
Section S7 for more details). For these reasons, we believe that imaging-to-
modeling comparison is not direct, but it does allow us to confirm that our
model is meaningfully capturing the main features of fibroTUG kinematics.

###  4.2 Active stress prediction for experimental conditions

In Fig. 5, the computational pipeline was used to assess the mechanics of
fibroTUGs under different mechanical environments. Both $F_{\rm p}$ and
$\sigma_{\rm a,m}$ show significant differences between the conditions (Fig.
5A-B,D-E). Since $\sigma_{\rm a,m}$ is a value that considers the structure of
all the tissue components, the differences in this value indicate that
modifying the mechanical environment where the cells develop will influence
their contractile maturity. These differences are also captured by $F_{\rm
p}$, but the relative differences are influenced by a multitude of factors,
including fiber mechanics, myofibrillar density and alignment, and contractile
maturity. This proves that the force output is the reflection of different
variables in the system, not only the iPSC-CMs active stress. This highlights
the importance of giving proper context to force output, which is often
treated as a direct surrogate of iPSC-CM stress [6, 8, 38].

###  4.3 Isolating the effect of mechanical variations

One key advantage of using computational models to study EHT mechanics is
their flexibility to study scenarios that are difficult, if not impossible, to
create in-vitro. These simulations can shed light on the influence of
different variables of a system by allowing us to change their values
individually and measure the variation in an output of interest. We performed
controlled simulations where one mechanical variable was modified at a time to
assess the changes in the $F_{\rm p}$/$\sigma_{\rm a,m}$ relationship. In Fig.
6C and E, we can see that when $\rho_{\rm m}=1$, changing the ECM stiffness
and the boundary conditions will have a slight impact on the capacity of the
tissue to translate cell active stress into force output. For example, for
stiffer fibers, a higher amount of the work done by the iPSC-CMs is lost in
tugging the less deformable matrix. Conversely, when the post is stiffer, the
magnitude of $\sigma_{\rm a,m}$ is slightly lower, corresponding to a more
efficient transduction of myofibrillar stress to the total output force.

When $\rho_{\rm m}=\rho_{\rm data}$, the conditions with lower myofibril
formation need to generate higher values of $\sigma_{\rm a,m}$. This is not
surprising because there will be less myofibril area to generate the same
force, which means that the existing sarcomeres need to generate more stress
to compensate. However, our model allows us to measure that effect and compare
it with the effect of other mechanical variations. The observed changes due to
poor myofibril formation are much higher than those observed due to post or
fiber stiffness (observed by the change of the bar plot magnitude between Fig.
6C-E and Fig. 6D-F). Furthermore, the cases that have lower $\rho_{\rm m}$
(and that need higher $\sigma_{\rm a,m}$ to generate the same force) correlate
with the cases that have lower force output experimentally (Fig. 5A-B). This
indicates that myofibril density is one of the most important mechanical
parameters explaining EHT force output. The importance of myofibril formation
has also been observed in single-cell models [39].

When evaluating the effect of fiber alignment, we can see that the magnitude
of $\sigma_{\rm a,m}$ is higher in random than aligned matrices. This makes
sense, as in these cases, both fibers and myofibrils are less aligned in the
post-to-post direction (see, for example, the PDFs in Fig. 3A). Here, a lot of
the work performed by the iPSC-CMs gets lost in pulling the transverse
direction. In Fig. 7, we decouple the alignment of fibers and myofibrils with
the intention of understanding the relative importance of myofibril alignment
versus fiber alignment. The results show that aligning the myofibrils on top
of a random matrix (Am-Rf) or vice-versa (Rm-Af) only explains a small part of
the difference between the all-aligned (Af-Am) and the all-random condition
(Rf-Rm). For this reason, we conclude that, _for the range of observed
variations in myofibril alignment_ , the matrix structure has a higher effect
on the force output. It could be that other mechanical environments generate a
higher disarray of myofibrils making this parameter more dominant.

The results of the model complement nicely with other biomarkers studied in
this platform [24]. For example, the aligned soft fibers/soft post condition
showed a higher proportion of connexin 43, more mature forms of myosin
(MLC-2v), and lower beats per minute, to mention a few. Our model shows that
this condition and the aligned soft fibers/stiff post are the most efficient
in transmitting force. However, the case with soft posts presents much higher
myofibril strains (see Fig. 6C), which has shown to be important in iPSC-CM
maturation [40]. These observations show the benefits of using a combined
computational-experimental approach to understand different aspects that are
involved in the EHT function and how they relate to the iPSC-CM maturation.

###  4.4 Limitations

The computational generation of myofibril fields enabled the flexibility to
explore non-experimental scenarios to understand the importance of myofibril
alignment. However, generating myofibrils networks that are representative of
real sarcomere organization is challenging. We performed several tests to
assess our approach using the validation datasets. We tested the importance of
including myofibril dispersion $\kappa_{\rm m}$ and studied the prediction
error induced when using computationally generated myofibril fields instead of
imaged-based ones. Detailed results of these tests are shown in Supplementary
Information Section S3. Briefly, we showed that not including $\kappa_{\rm m}$
will mainly impact the deformations in X, with these simulations showing less
necking than expected. However, the active stress prediction remains very
similar. The test of the artificial myofibril fields showed that the sarcomere
strain values are the most affected quantity, but the active stress prediction
is, again, very similar to the one computed using image-generated myofibril
fields. The results of these tests show that we can use our proposed method to
correctly assess changes in $\sigma_{\rm a,m}$.

A few assumptions were made in the model not based on the experimental data.
The passive stiffness of the iPSC-CMs, isolated from the fibers, was not
measured in our experiments. This parameter is tricky to measure, as it is
known that the cell will change its stiffness during development [41] and that
this value also depends on the temperature and the state of its contractile
apparatus [42]. Therefore, for simplicity, we considered the passive stiffness
of the cells to be 2 kPa, a value similar to the one measured by Ballan et al.
[35] in iPSC-CMs, and lower than what is usually measured in mature CMs [43].
We further assume that the passive response of the cell was isotropic. This is
because no studies were found on the anisotropy of iPSC-CMs and since the
iPSC-CMs considered in this study are still immature (thus, the cytoskeleton
is less organized), we believe that our assumptions of an isotropic, less
stiff passive response than in adult CMs are reasonable.

Another modeling assumption is that all iPSC-CMs have a synchronous activation
and that all their contraction can be modeled using a single $\eta$ parameter.
This means that the spatial heterogeneity observed in the $\sigma_{\rm a,m}$
fields in Fig. 5C,F are the product of only the length-dependent function
$\phi$ in Eq. (6). The simultaneous activation assumption is backed up by
these tissues being very small, and the calcium transients measured
experimentally did not show spatial differences [24]. However, it is possible
(and probably expected) that not all the iPSC-CMs in one tissue have the same
maturity, which could be modeled by having an $\eta$ field instead of a single
scalar value. For example, instead of only using the post-displacement/force
as a constraint, we could use the displacement measurements from the
MicroBundleCompute software to find a local $\eta$. Different strategies to do
so will be explored in the future. Nevertheless, we believe that using a
single $\eta$ parameter allows us to compute an average activation and that
the overall regional differences will not affect the conclusions of this
paper.

Finally, we modeled the fibroTUG in 2D, which also forces us to consider a
material with low compressibility. This decision is based on the low thickness
to length ratio that is $\sim 3\%$, and because only 2D projected images were
available for the fibers, myofibrils, and displacement tracking. Future
studies will aim to generate a pipeline to reconstruct the full tissue
geometry from 3D stacks obtained from high-magnification imaging techniques
and assess the validity of the 2D model.

##  5 Conclusions

This work introduces a data-driven biomechanical computational model of EHTs.
The implemented pipeline leverages the experimental data of the fibroTUG
platform, which was introduced specifically to study iPSC-CMs function due to
mechanics. With the model, we are able to measure the iPSC-CM active stress
that replicates the experimental observations. Unlike the tissue force output,
this value is a direct measure of iPSC-CM contractile function. We used the
model to assess the effects of different mechanical environments on iPSC-CMs
stress generation. The results followed the same trends as the force output,
indicating maturation differences, but with lower relative differences. This
highlights that EHT force generation is the product of different variables
besides iPSC-CM contraction. We explore this fact with our model and show that
myofibril density is one of the key factors explaining the experimental
differences observed. The developed framework opens the door to many other
applications in the cardiac tissue engineering field, and it shows how a
combination of in-silico with in-vitro approaches can help us better
understand the mechanobiology of EHTs.

##  6 Acknowledgements

JJ acknowledges the support of ANID and Fulbright Chile through the Becas BIO
program. DN acknowledges funding from the Engineering and Physical Sciences
Research Council Healthcare Technology Challenge Award (EP/R003866/1). BMB
acknowledges the support from the National Science Foundation through
CBET-2033654. All authors acknowledge the support from the National Science
Foundation through the Nanosystems Engineering Research Center for Directed
Multiscale Assembly of Cellular Metamaterials with Nanoscale Precision (CELL-
MET, EEC-1647837).

Table 1: Results of the myofibril density characterization (mean $\pm$ standard deviation) for the different experimental conditions. The normalized density is the mean density of each case divided by the aligned soft fibers/soft post condition mean density (i.e., divided by 0.643). | Aligned Fibers | Random Fibers
---|---|---
Condition | Density | Normalized Density | Density | Normalized Density
Soft Fibers / Soft Post | 0.643 $\pm$ 0.104 | 1.000 | 0.395 $\pm$ 0.180 | 0.613
Stiff Fibers / Soft Post | 0.444 $\pm$ 0.126 | 0.690 | 0.391 $\pm$ 0.060 | 0.608
Soft Fibers / Stiff Post | 0.593 $\pm$ 0.128 | 0.922 |  0.563 $\pm$ 0.119 | 0.875
Figure 1: Experimental setup and data output. (A) Diagram of FibroTUG
manufacturing process showing the electrospun fibers, the selective
crosslinking, and the cell seeding process. (B) Force vs. strain curves
obtained by performing an indentation test in the fibers only. (C) Post
displacement traces obtained by tracking the post position in videos of the
contractile EHTs. (D) Immunofluorescence images of the DVS fibers and the
titin protein that indicates the location of Z-discs.  Figure 2: Image
processing results for (A) the DVS fibers showing (from left to right) fiber
density ($\rho_{\rm f}$), fiber alignment (shown as the absolute angle respect
the post-to-post direction, $|\theta_{\rm f}|$), and fiber dispersion
$\kappa_{\rm f}$, and for (B) For the titin images showing myofibril density
$\rho_{\rm m}$, myofibril alignment (as the absolute angle respect the post-
to-post direction, $|\theta_{\rm m}|$), and myofibril dispersion $\kappa_{\rm
m}$.  Figure 3: Myofibril characterization. (A) Histograms and their
respective Von Mises fit for each condition studied. (B) Example image of the
density computation for the different cases. The bar plots at the right show
the density value computed for all the images available. All data presented as
mean ± standard deviation. Figure 4: Results of the pipeline for three
complete datasets. Each color represents a different tissue. (A) Post
displacement data (dots) and model (line). (B) Simulated active stress. (C)
MicroBundleCompute displacement outputs (left) and simulation results (right)
in the post-to-post direction Y. (D) Measured displacements (left) and
simulation results in the X direction. (E) Correlation between the tracking
and the simulation displacements in Y and (F) in X. (C) and (D) correspond to
the case pictured in cyan. Figure 5: Results at maximum contraction for the
experimental conditions of interest. (A) Force output, (B) mean $\sigma_{\rm
a,m}$, and (C) $\sigma_{\rm a,m}$ field at maximum contraction for aligned
fibrous matrix. (D) Force output, (E) mean $\sigma_{\rm a,m}$, and (F)
$\sigma_{\rm a,m}$ field at maximum contraction for random fibrous matrix. All
data presented as mean $\pm$ standard deviation. * $p<0.0001$ by unpaired
t-tests. Figure 6: Impact of different mechanical conditions on the post-
force/active stress relationship. (A) Post-force and post-displacement for
each condition. Note that the force is kept the same for all conditions, while
the post-displacement varies in the case of stiff posts. (B) Procedure to
generate in-silico samples (see supplementary Information S3 for more
details). An aligned and a random matrix were built from the images, and
computationally myofibril fields were created using the corresponding PDF.
Results for an aligned condition (C) with $\rho_{\rm m}=1$, and (D) with
$\rho_{\rm m}$ equal to the specific density obtained from data for each case
($\rho_{\rm data}$). (E) and (F) show the results for random matrices with
$\rho_{\rm m}=1$ and $\rho_{\rm m}=\rho_{\rm data}$. The values below the bar
plots in (D) and (F) correspond to the mean density used in these cases (see
Table 1). All data presented as mean $\pm$ standard deviation. * $p<0.0001$ by
unpaired t-tests. Figure 7: Procedure for the myofibril vs fiber alignment
test. Simulated conditions that have aligned fibers with random myofibrils
(Af-Rm) and vice-versa (Rf-Am) in-silico tissues are generated by using
fibrous matrices generated from images and the PDFs describing myofibril
alignment. These simulations are compared with the aligned (Af-Am) and random
(Rf-Rm) experimental conditions. (B) Mean $\sigma_{\rm a,m}$ at maximum
contraction for the Af-Am, Af-Rm, Rf-Am, Rf-Rm, simulations using soft
fibers/soft posts. All data presented as mean $\pm$ standard deviation. *
$p<0.0001$ by unpaired t-tests.

## References

  * [1] C. W. Tsao, A. W. Aday, Z. I. Almarzooq, A. Alonso, A. Z. Beaton, M. S. Bittencourt, A. K. Boehme, A. E. Buxton, A. P. Carson, Y. Commodore-Mensah, M. S. Elkind, K. R. Evenson, C. Eze-Nliam, J. F. Ferguson, G. Generoso, J. E. Ho, R. Kalani, S. S. Khan, B. M. Kissela, K. L. Knutson, D. A. Levine, T. T. Lewis, J. Liu, M. S. Loop, J. Ma, M. E. Mussolino, S. D. Navaneethan, A. M. Perak, R. Poudel, M. Rezk-Hanna, G. A. Roth, E. B. Schroeder, S. H. Shah, E. L. Thacker, L. B. VanWagner, S. S. Virani, J. H. Voecks, N.-Y. Wang, K. Yaffe, S. S. M. and, Heart disease and stroke statistics—2022 update: A report from the american heart association, Circulation 145 (8) (feb 2022).  doi:10.1161/cir.0000000000001052.
  * [2] Y. Zhao, G. Eng, B. W. Lee, M. Radisic, G. Vunjak-Novakovic, Chapter 32 - cardiac tissue engineering, in: R. Lanza, R. Langer, J. P. Vacanti, A. Atala (Eds.), Principles of Tissue Engineering (Fifth Edition), fifth edition Edition, Academic Press, 2020, pp. 593–616.  doi:https://doi.org/10.1016/B978-0-12-818422-6.00033-2.
URL https://www.sciencedirect.com/science/article/pii/B9780128184226000332

  * [3] M. F. Tenreiro, A. F. Louro, P. M. Alves, M. Serra, Next generation of heart regenerative therapies: progress and promise of cardiac tissue engineering, npj Regenerative Medicine 6 (1) (jun 2021).  doi:10.1038/s41536-021-00140-4.
  * [4] E. Karbassi, A. Fenix, S. Marchiano, N. Muraoka, K. Nakamura, X. Yang, C. E. Murry, Cardiomyocyte maturation: advances in knowledge and implications for regenerative medicine, Nature Reviews Cardiology 17 (6) (2020) 341–359.  doi:10.1038/s41569-019-0331-x.
  * [5] J.-L. Ruan, N. L. Tulloch, M. V. Razumova, M. Saiget, V. Muskheli, L. Pabon, H. Reinecke, M. Regnier, C. E. Murry, Mechanical stress conditioning and electrical stimulation promote contractility and force maturation of induced pluripotent stem cell-derived human cardiac tissue, Circulation 134 (20) (2016) 1557–1567.  doi:10.1161/circulationaha.114.014998.
  * [6] K. Ronaldson-Bouchard, S. P. Ma, K. Yeager, T. Chen, L. Song, D. Sirabella, K. Morikawa, D. Teles, M. Yazawa, G. Vunjak-Novakovic, Advanced maturation of human cardiac tissue grown from pluripotent stem cells, Nature 556 (7700) (2018) 239–243.  doi:10.1038/s41586-018-0016-3.
  * [7] S. J. DePalma, C. D. Davidson, A. E. Stis, A. S. Helms, B. M. Baker, Microenvironmental determinants of organized iPSC-cardiomyocyte tissues on synthetic fibrous matrices, Biomaterials Science 9 (1) (2021) 93–107.  doi:10.1039/d0bm01247e.
  * [8] A. Leonard, A. Bertero, J. D. Powers, K. M. Beussman, S. Bhandari, M. Regnier, C. E. Murry, N. J. Sniadecki, Afterload promotes maturation of human induced pluripotent stem cell derived cardiomyocytes in engineered heart tissues, Journal of Molecular and Cellular Cardiology 118 (2018) 147–158.  doi:10.1016/j.yjmcc.2018.03.016.
  * [9] J. M. Bliley, M. C. S. C. Vermeer, R. M. Duffy, I. Batalov, D. Kramer, J. W. Tashman, D. J. Shiwarski, A. Lee, A. S. Teplenin, L. Volkers, B. Coffin, M. F. Hoes, A. Kalmykov, R. N. Palchesko, Y. Sun, J. D. H. Jongbloed, N. Bomer, R. A. de Boer, A. J. H. Suurmeijer, D. A. Pijnappels, M. C. Bolling, P. van der Meer, A. W. Feinberg, Dynamic loading of human engineered heart tissue enhances contractile function and drives a desmosome-linked disease phenotype, Science Translational Medicine 13 (603) (jul 2021).  doi:10.1126/scitranslmed.abd1817.
  * [10] C. Tu, B. S. Chao, J. C. Wu, Strategies for improving the maturity of human induced pluripotent stem cell-derived cardiomyocytes, Circulation Research 123 (5) (2018) 512–514.  doi:10.1161/circresaha.118.313472.
  * [11] A. Körner, M. Mosqueira, M. Hecker, N. D. Ullrich, Substrate stiffness influences structural and functional remodeling in induced pluripotent stem cell-derived cardiomyocytes, Frontiers in Physiology 12 (aug 2021).  doi:10.3389/fphys.2021.710619.
  * [12] A. J. S. Ribeiro, Y.-S. Ang, J.-D. Fu, R. N. Rivas, T. M. A. Mohamed, G. C. Higgs, D. Srivastava, B. 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) (2015) 12705–12710.  doi:10.1073/pnas.1508073112.
  * [13] A. C. Allen, E. Barone, N. Momtahan, C. O. Crosby, C. Tu, W. Deng, K. Polansky, J. Zoldan, Temporal impact of substrate anisotropy on differentiating cardiomyocyte alignment and functionality, Tissue Engineering Part A 25 (19-20) (2019) 1426–1437.  doi:10.1089/ten.tea.2018.0258.
  * [14] N. Thavandiran, N. Dubois, A. Mikryukov, S. Massé, B. Beca, C. A. Simmons, V. S. Deshpande, J. P. McGarry, C. S. Chen, K. Nanthakumar, G. M. Keller, M. Radisic, P. W. Zandstra, Design and formulation of functional pluripotent stem cell-derived cardiac microtissues, Proceedings of the National Academy of Sciences 110 (49) (nov 2013).  doi:10.1073/pnas.1311120110.
  * [15] S. Das, A. Ippolito, P. McGarry, V. S. Deshpande, Cell reorientation on a cyclically strained substrate, PNAS Nexus 1 (5) (sep 2022).  doi:10.1093/pnasnexus/pgac199.
  * [16] H. Wang, A. A. Svoronos, T. Boudou, M. S. Sakar, J. Y. Schell, J. R. Morgan, C. S. Chen, V. B. Shenoy, Necking and failure of constrained 3d microtissues induced by cellular tension, Proceedings of the National Academy of Sciences 110 (52) (2013) 20923–20928.  doi:10.1073/pnas.1313662110.
  * [17] P. Robison, M. A. Caporizzo, H. Ahmadzadeh, A. I. Bogush, C. Y. Chen, K. B. Margulies, V. B. Shenoy, B. L. Prosser, Detyrosinated microtubules buckle and bear load in contracting cardiomyocytes, Science 352 (6284) (apr 2016).  doi:10.1126/science.aaf0659.
  * [18] H. Wang, A. Abhilash, C. S. Chen, R. G. Wells, V. B. Shenoy, Long-range force transmission in fibrous matrices enabled by tension-driven alignment of fibers, Biophysical Journal 107 (11) (2014) 2592–2603.  doi:10.1016/j.bpj.2014.09.044.
  * [19] W. Ronan, R. M. McMeeking, C. S. Chen, J. P. McGarry, V. S. Deshpande, Cooperative contractility: The role of stress fibres in the regulation of cell-cell junctions, Journal of Biomechanics 48 (3) (2015) 520–528.  doi:10.1016/j.jbiomech.2014.11.025.
  * [20] R. Miller, E. Kerfoot, C. Mauger, T. F. Ismail, A. A. Young, D. A. Nordsletten, An implementation of patient-specific biventricular mechanics simulations with a deep learning and computational pipeline, Frontiers in Physiology 12 (sep 2021).  doi:10.3389/fphys.2021.716597.
  * [21] L. Asner, M. Hadjicharalambous, R. Chabiniok, D. Peressutti, E. Sammut, J. Wong, G. Carr-White, R. Razavi, A. King, N. Smith, J. Lee, D. Nordsletten, Patient-specific modeling for left ventricular mechanics using data-driven boundary energies, Computer Methods in Applied Mechanics and Engineering 314 (2017) 269–295.  doi:10.1016/j.cma.2016.08.002.
  * [22] J. H. Bracamonte, S. K. Saunders, J. S. Wilson, U. T. Truong, J. S. Soares, Patient-specific inverse modeling of in vivo cardiovascular mechanics with medical image-derived kinematics as input data: Concepts, methods, and applications, Applied Sciences 12 (8) (2022) 3954.  doi:10.3390/app12083954.
  * [23] G. Balaban, H. Finsberg, S. Funke, T. F. Håland, E. Hopp, J. Sundnes, S. Wall, M. E. Rognes, In vivo estimation of elastic heterogeneity in an infarcted human heart, Biomechanics and Modeling in Mechanobiology 17 (5) (2018) 1317–1329.  doi:10.1007/s10237-018-1028-5.
  * [24] S. J. DePalma, J. Jillberto, A. E. Stis, D. D. Huang, J. Lo, C. D. Davidson, M. E. Jewett, H. Kobeissi, A. Chowdhury, C. S. Chen, E. Lejeune, A. S. Helms, D. A. Nordsletten, B. M. Baker, Microenvironmental mechanical inputs dictate the assembly and function of ipsc-derived cardiomyocyte microtissues, in preparation (2023).
  * [25] A. J. Schriefl, A. J. Reinisch, S. Sankaran, D. M. Pierce, G. A. Holzapfel, Quantitative assessment of collagen fibre orientations from two-dimensional images of soft biological tissues, Journal of The Royal Society Interface 9 (76) (2012) 3081–3093.  doi:10.1098/rsif.2012.0339.
  * [26] A. Sharma, C. N. Toepfer, T. Ward, L. Wasson, R. Agarwal, D. A. Conner, J. H. Hu, C. E. Seidman, CRISPR/cas9-mediated fluorescent tagging of endogenous proteins in human pluripotent stem cells, Current Protocols in Human Genetics 96 (1) (jan 2018).  doi:10.1002/cphg.52.
  * [27] J. D. Humphfrey, K. R. Rajagopal, A Constrained Mixture Model for Growth and Remodeling of Soft Tissues, Mathematical Models and Methods in Applied Sciences 12 (03) (2002) 407–430.  doi:10.1142/s0218202502001714.
  * [28] G. A. Holzapfel, Nonlinear Solid Mechanics, John Wiley & Sons, 2000.
  * [29] J. Bonet, R. D. Wood, Nonlinear Continuum Mechanics for Finite Element Analysis, Cambridge University Press, 2008.  doi:10.1017/cbo9780511755446.
  * [30] G. A. Holzapfel, R. W. Ogden, Constitutive modelling of passive myocardium: a structurally based framework for material characterization, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 367 (1902) (2009) 3445–3475.  doi:10.1098/rsta.2009.0091.
  * [31] D. Nolan, A. Gower, M. Destrade, R. Ogden, J. McGarry, A robust anisotropic hyperelastic formulation for the modelling of soft tissue, Journal of the Mechanical Behavior of Biomedical Materials 39 (2014) 48–60.  doi:10.1016/j.jmbbm.2014.06.016.
  * [32] G. A. Holzapfel, R. W. Ogden, S. Sherifova, On fibre dispersion modelling of soft biological tissues: a review, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 475 (2224) (2019) 20180736.  doi:10.1098/rspa.2018.0736.
  * [33] M. Hadjicharalambous, R. Chabiniok, L. Asner, E. Sammut, J. Wong, G. Carr-White, J. Lee, R. Razavi, N. Smith, D. Nordsletten, Analysis of passive cardiac constitutive laws for parameter estimation using 3d tagged MRI, Biomechanics and Modeling in Mechanobiology 14 (4) (2014) 807–828.  doi:10.1007/s10237-014-0638-9.
  * [34] J. P. Pascon, Large deformation analysis of plane-stress hyperelastic problems via triangular membrane finite elements, International Journal of Advanced Structural Engineering 11 (3) (2019) 331–350.  doi:10.1007/s40091-019-00234-w.
  * [35] N. Ballan, N. Shaheen, G. M. Keller, L. Gepstein, Single-cell mechanical analysis of human pluripotent stem cell-derived cardiomyocytes for drug testing and pathophysiological studies, Stem Cell Reports 15 (3) (2020) 587–596.  doi:10.1016/j.stemcr.2020.07.006.
  * [36] R. Kerckhoffs, P. Bovendeerd, F. Prinzen, K. Smits, T. Arts, Intra- and interventricular asynchrony of electromechanics in the ventricularly paced heart, Journal of Engineering Mathematics 47 (3/4) (2003) 201–216.  doi:10.1023/b:engi.0000007972.73874.da.
  * [37] H. Kobeissi, J. Jilberto, M. Çağatay Karakan, X. Gao, S. J. DePalma, S. L. Das, L. Quach, J. Urquia, B. M. Baker, C. S. Chen, D. Nordsletten, E. Lejeune, Microbundlecompute: Automated segmentation, tracking, and analysis of sub-domain deformation in cardiac microbundles, in preparation (2023).
  * [38] T. Boudou, W. R. Legant, A. Mu, M. A. Borochin, N. Thavandiran, M. Radisic, P. W. Zandstra, J. A. Epstein, K. B. Margulies, C. S. Chen, A microfabricated platform to measure and manipulate the mechanics of engineered cardiac microtissues, Tissue Engineering Part A 18 (9-10) (2012) 910–919.  doi:10.1089/ten.tea.2011.0341.
  * [39] A. S. V. Roest, C. Liu, M. M. Morck, K. B. Kooiker, G. Jung, D. Song, A. Dawood, A. Jhingran, G. Pardon, S. Ranjbarvaziri, G. Fajardo, M. Zhao, K. S. Campbell, B. L. Pruitt, J. A. Spudich, K. M. Ruppel, D. Bernstein, Hypertrophic cardiomyopathy $\beta$-cardiac myosin mutation (p710r) leads to hypercontractility by disrupting super relaxed state, Proceedings of the National Academy of Sciences 118 (24) (jun 2021).  doi:10.1073/pnas.2025030118.
  * [40] Y.-C. Tsan, S. J. DePalma, Y.-T. Zhao, A. Capilnasiu, Y.-W. Wu, B. Elder, I. Panse, K. Ufford, D. L. Matera, S. Friedline, T. S. O’Leary, N. Wubshet, K. K. Y. Ho, M. J. Previs, D. Nordsletten, L. L. Isom, B. M. Baker, A. P. Liu, A. S. Helms, Physiologic biomechanics enhance reproducible contractile development in a stem cell derived cardiac muscle platform, Nature Communications 12 (1) (oct 2021).  doi:10.1038/s41467-021-26496-1.
  * [41] S. Lahmers, Y. Wu, D. R. Call, S. Labeit, H. Granzier, Developmental control of titin isoform expression and passive stiffness in fetal and neonatal myocardium, Circulation Research 94 (4) (2004) 505–513.  doi:10.1161/01.res.0000115522.52554.86.
  * [42] N. M. King, M. Methawasin, J. Nedrud, N. Harrell, C. S. Chung, M. Helmes, H. Granzier, Mouse intact cardiac myocyte mechanics: cross-bridge and titin-based stress in unactivated cells, Journal of General Physiology 137 (1) (2010) 81–91.  doi:10.1085/jgp.201010499.
  * [43] M. A. Caporizzo, C. Y. Chen, K. Bedi, K. B. Margulies, B. L. Prosser, Microtubules increase diastolic stiffness in failing human cardiomyocytes and myocardium, Circulation 141 (11) (2020) 902–915.  doi:10.1161/circulationaha.119.043930.