Modeling cardiomyocyte mechanics and autoregulation of contractility by mechano-chemo-transduction feedback.

The heart pumps blood into circulation against vascular resistance and actively regulates the contractile force to compensate for mechanical load changes. Our experimental data show that cardiomyocytes have a mechano-chemo-transduction (MCT) mechanism that increases intracellular Ca 2 + transient to enhance contractility in response to increased mechanical load. This study advances the cardiac excitation- Ca 2 + signaling-contraction (E-C) coupling model on conceptual and technical fronts. First, we developed analytical and computational models to perform 3-dimensional mechanical analysis of cardiomyocytes contracting in a viscoelastic medium under mechanical load. Next, we proposed an MCT feedback loop in the E-C coupling dynamic system to shift the feedforward paradigm of cardiac E-C coupling to an autoregulation model. Our combined modeling and experimental studies reveal that MCT enables autoregulation of E-C coupling and contractility in single cardiomyocytes, which underlies the heart's intrinsic autoregulation in compensatory response to load changes in order to maintain the stroke volume and cardiac output.

Introduction

The heart is a remarkably “smart” pump. In each heartbeat, cardiomyocytes generate contractile forces to pump blood into circulation against the mechanical load from vascular resistance. The blood pressure may change momentarily under different physiological conditions (i.e., posture change Munoz et al. (2021), lifting weights (Benn et al., 1996), exercise (Mitchell et al., 2005), etc.), and excessive loading may occur under various disease conditions (i.e., hypertension (Tello et al., 2019), pressure-overload (Jashari et al., 2015), dilated cardiomyopathy (Chandar et al., 2010), etc.). To maintain cardiac output, the heart has systemic and intrinsic regulations to adapt cardiomyocyte contraction to compensate for load changes. At the systemic level, the heart has autonomic neurohormonal regulations. At the organ level, the heart has intrinsic autoregulation as described by the Frank-Starling law and the Anrep effect (Von Anrep, 1912; Patterson et al., 1914; Starling, 1918). The Frank-Starling law describes the increase of contractile force in response to increased preload, i.e. stretching of the cardiomyocytes at the end of diastolic filling (Patterson et al., 1914; Starling, 1918). The Anrep effect describes the upregulation of contractility in response to increased afterload, i.e., increased systolic pressure against vascular resistance (Von Anrep, 1912). The current understanding is that the Frank-Starling mechanism is myofilament-based and does not involve changes in the transient; however, the Anrep effect is associated with changes in the cytosolic transient (Cingolani et al., 2012; Izu et al., 2020).

This adaptation to afterload at the cardiomyocyte level is the focus of this paper. Given that signaling is central to cardiac excitation- signaling-contraction (E-C) coupling (Bers, 2002), it is expected that mechanical load induced changes in signaling should impact many -dependent biochemical processes in cardiomyocytes. Despite the load adaptations at the organ and tissue levels, the underlying cellular and molecular mechanisms that transduce mechanical load to biochemical signals to regulate the contractility of cardiomyocytes remain largely unresolved (Moss and Fitzsimons, 2002; de Tombe et al., 2010; Bollensdorff et al., 2011; Cingolani et al., 2012).

Cardiac E-C coupling is governed by (at least) three nonlinear dynamic systems: (1) the electrical system that controls action potentials, which excites (2) the signaling system to cause intracellular signaling, which activates (3) the contractile system to cause cardiomyocyte contraction (Bers, 2001, 2002). The classical paradigm of cardiac E-C coupling emphasizes feedforward control from the action potential to the signaling and muscle contraction. How mechanical loading during contraction affects cardiomyocytes has been relatively understudied.

To investigate the effects of mechanical load on cardiomyocytes, we developed the Cell-in-Gel technology, in which live cardiomyocytes freshly isolated from the heart of an animal model are embedded in a viscoelastic hydrogel and electrically stimulated to undergo E-C coupling (Jian et al., 2014). We systematically adjusted the hydrogel stiffness to apply various levels of mechanical resistance to cardiomyocyte contraction, thereby simulating afterload at the single cell scale in a controlled way. Experimental studies showed that cardiomyocytes possess mechano-chemo-transduction (MCT) mechanisms that sense and transduce mechanical load to biochemical signals (Cingolani et al., 2012; Chen-Izu and Izu, 2017). MCT increases the cytosolic transient to enhance contractility in response to load increases, which provides a feedback loop in the E-C coupling. Our experimental study and conceptual model suggest that the cardiomyocyte can actively regulate contractility in compensation to load changes at the single-cell level (Shimkunas et al., 2021; Izu et al., 2021). However, quantitative modeling and rigorous mechanical analysis are needed to precisely understand how the MCT feedback in the E-C coupling dynamic system might enable autoregulation.

In this study, we aim to develop an autoregulation model by advancing cardiac modeling at technical and conceptual fronts. The technical innovation is to develop three-dimensional (3D) mechanical models of a single cardiomyocyte contracting under mechanical load in a viscoelastic medium. First, we use an analytical model using the viscoelastic Eshelby (VE Eshelby) inclusion analysis (Kazemi-Lari et al., 2021) of an idealized ellipsoid-shaped cell contracting under afterload in a viscoelastic medium, which provides exact, efficient, and fast calculations of the 3D strains and stresses in the cardiomyocyte during steady beat-to-beat contraction. Second, we develop a more refined numerical model using finite-element analysis (FEA) to provide detailed calculations of the 3D mechanical fields in realistic shaped cardiomyocytes. Third, we verify that the analytical model results are in good agreement with the spatially averaged values of the FEA model. These studies set a firm foundation for us to develop the autoregulation model by using the analytical model as a fast approximation to the average mechanical behavior of cardiomyocytes and also allow us to extract the parameters of the autoregulation model from the data measured in the Cell-in-Gel experiments. The goal is to create a realistic autoregulation model that can describe the autoregulation behavior of real cardiomyocytes.

A major conceptual development is to shift the classic paradigm of the feedforward E-C coupling model to a closed-loop feedback MCT model of autoregulation. We use the autoregulation model to propose a new mathematical framework to investigate how the cardiac E-C coupling dynamics might be regulated by MCT under changing afterload. Model simulations predict testable autoregulation behaviors in cardiomyocytes which are compared to experimental data.

Here, we present our findings that autoregulation of contractility naturally arises from the E-C coupling dynamic system with MCT feedback in cardiomyocytes at the single cardiomyocyte level. Our new autoregulatory E-C coupling model describes a fundamental MCT mechanism underlying the Anrep effect. It contributes to the heart’s intrinsic regulation of contractility in response to mechanical load changes under physiological and pathological conditions.

Model development

Experimental data

The Cell-in-Gel system, shown schematically in Figure 1A, was invented by our group to control afterload at the single-cell level and to study how cardiomyocyte E-C coupling is affected by afterload. It is composed of live cardiomyocytes embedded in a biocompatible constraining viscoelastic hydrogel matrix that has similar mechanical properties as the extracellular matrix of the heart tissue (Pislaru et al., 2014). The hydrogel consists of a polyvinyl alcohol (PVA) polymer crosslinked with 4-boronate-polyethylene glycol (4B-PEG) crosslinker. The 4B-PEG also attaches to the glycans on the cell surface, effectively tethering the cell surface to the hydrogel matrix. The hydrogel is constantly perfused with Tyrode’s solution that provides a constant physiological milieu to keep the cell alive during the experiments.

Figure 1

Cell-in-Gel experiments to study cardiomyocyte contraction under mechanical load

(A–I) (A) Schematic of the Cell-in-Gel system showing a single-cell embedded in a 3D hydrogel and electrically excited to beat, (B) simultaneous measurements of contraction and transients using the gel-dissolve protocol, (C) approximation of the Cell-in-Gel system using the Eshelby inclusion BVP and the 3D schematic of the equivalent ellipsoidal cell (inclusion), (D-E) a representative experimentally measured load-free (LF) strain history cycle of a cell and its power spectrum, (F) the proposed autoregulation model with Mechano-Chemo-Transduction (MCT) feedback, (G–I) comparison of LF (and MCT-free) and in-gel responses from the gel-dissolve (self-control) experiment, and (J) simulated histories (dashed lines) of the MCT-free and MCT eigenstrains (β) calculated from Equations 8, 9a, and 9b.

The Cell-in-Gel system provides a 3D viscoelastic medium with several advantages. By varying the crosslinker density (C) between experiments, the stiffness of the hydrogel can be adjusted to apply different levels of afterload on the contracting cell (Jian et al., 2014). The hydrogel can also be dissolved away, releasing the cell from the gel into the aqueous solution to beat freely (called the “load-free” condition).

We performed experiments using rabbit ventricular myocytes to investigate the mechanical load effects on mammalian cardiomyocytes by using freshly isolated left ventricular myocytes from mouse and rabbit hearts (Jian et al., 2014; Shimkunas et al., 2021; Hegyi et al., 2021). As shown in one such experiment of Figure 1B, the cardiomyocyte was electrically stimulated to undergo E-C coupling, first in the hydrogel under mechanical afterload, and then after the hydrogel was dissolved to release the cell into the Tyrode’s solution (load-free). Such self-control experiments allow us to compare the transient and contraction of the same cell under mechanical load versus load-free, which avoids cell to cell variability and provides self-consistent information for model calibration.

The viscoelastic behavior of the hydrogel in the Cell-in-Gel system was obtained from the rheology data as described previously (Kazemi-Lari et al., 2021). For the purpose of mechanical modeling, the viscoelastic behavior is approximated by the generalized Maxwell viscoelastic model, where the relaxation modulus, , is mathematically represented by a Prony series according towhere and are the respective instantaneous short-term and long-term shear moduli, () are relaxation time constants, and represent phase fractions that apportion the amount of material that relaxes according to each time constant (Wineman and Rajagopal, 2000). Fitting to the experimental rheology data was achieved by selecting time constants at equal increments along a logarithmic scale, s, and then fitting the parameters () and (Kazemi-Lari et al., 2021). For most of the hydrogel samples, the distribution of ξ versus (relaxation spectrum) resembled a Gaussian distribution, therefore the relaxation spectra for all hydrogels were constrained to lie alongwhere is the mean time constant, is the SD(characteristic width) of the distribution, and is a normalizing constant to make . In the current study, the presented mechanical analysis is performed in a hydrogel with crosslink density of approximately 8% corresponding to the coefficients of the viscoelastic model as listed in Table 1.

Table 1

Parameter values of the generalized Maxwell viscoelastic model used for mathematical modeling of the hydrogel at 8% crosslink density

G_0 (kPa)	G_∞ (kPa)	τ_0 (s)	σ_0	ξ_0
4.723	≈ 0	0.0783	0.571	3.36224

In addition, the elastic modulus of the cardiomyocyte was estimated using the available data in the literature from the uniaxial force-stretch measurements on a single ventricular myocyte stimulated at 1.0 Hz frequency (Prosser et al., 2013). According to the extracted stress-strain response curve, the cell had an approximate Young’s modulus of kPa when electrically stimulated (Kazemi-Lari et al., 2021). Treating the cell as an incompressible linear elastic solid, the effective shear modulus was kPa. To our knowledge, no viscoelastic data on a single cardiomyocyte currently exists in the literature.

Analytical model: Viscoelastic Eshelby inclusion problem

The mechanical modeling is first analytically explored in the context of the viscoelastic extension of the classical Eshelby inclusion boundary-value problem (BVP) that has a known closed-form solution. For this purpose, the cell is approximated as an ellipsoidal inclusion that is embedded in and perfectly bonded to an infinite gel matrix. The inclusion undergoes an eigenstrain, i.e., a spontaneous inelastic strain that simulates cell contraction and creates residual stresses in both the inclusion and matrix. Importantly, an ellipsoidal inclusion is the only known shape that has a uniform interior strain and stress field as an exact solution at a given time instant. The analysis treats the cardiac muscle cell as an isotropic continuum to calculate the average stress and strain of the cell in a homogenized sense (Kazemi-Lari et al., 2021).

To address the cardiomyocyte’s autoregulation mechanisms in response to changes in mechanical afterload, we next propose a feedback model to capture the up-regulation of cytosolic transients to elevated mechanical stress in the cell. The increased leads to enhanced cell contractility in the mathematical form of an amplified eigenstrain (contractility).

To verify the reasonableness of the proposed VE Eshelby approach, a 3D finite-element mechanical analysis is also performed on the actual geometry of a typical rabbit left ventricular cardiomyocyte. The outer envelope of the 3D geometry of the cell is obtained using surface reconstruction of a stack of two-dimensional (2D) images captured through the thickness of the cell. The cell is again assumed to be an isotropic solid with homogenized properties.

Our analytical model is based on the well-known Eshelby inclusion analysis (Eshelby, 1957, 1959), where the exact solution for the stress and strain fields of an elastic matrix with an ellipsoidal-shape inclusion was developed. The original formulation addressed the problem where a subregion (the inclusion occupying a subvolume ) of an infinite, linear elastic solid undergoes a spontaneous change of shape (transformation strain or eigenstrain). Because the inclusion is constrained by the surrounding matrix which occupies the complementary volume , a residual state of stress is created inside and outside the inclusion and the inclusion is restrained from achieving its stress-free configuration. In this study, we use a viscoelastic extension of the original model to better quantify the residual stress fields in the Cell-in-Gel system (Shaw et al., 2013; Kazemi-Lari et al., 2021). We envision the cardiac cell as an ellipsoidal inclusion attempting to contract (in response to electrical excitation), but its deformation is constrained by the surrounding hydrogel matrix which only permits a partial contraction to be achieved. Consistent with the cardiomyocyte in the Cell-in-Gel experiment (Figure 1C), the ellipsoidal inclusion has dimensions and , resulting in an aspect ratio of . We assume the cell membrane is perfectly bonded to the gel and the remote boundaries of the gel are load-free.

The cell is assumed to undergo a 3D inelastic strain , called the eigenstrain. This simulates the strain of the cell when it contracts load-free (unconstrained by any matrix). In other words, the eigenstrain () of the cell is interpreted here as the cardiac muscle cell’s intrinsic contractility or the cell’s potential for contraction if there was no constraining gel (load-free). The mechanical resistance of the gel, however, causes the magnitude of the cell strain () to be somewhat less. This tug-of-war causes self-equilibrated residual stress fields in the inclusion and the matrix. As was obtained from a self-control experiment, a typical experimentally measured axial load-free strain (β), which is magnitude of the fractional shortening of the cell but with a negative sign, is shown in Figure 1D.

With respect to the Cartesian coordinate system (Scalar quantities are denoted in normal type, whereas bold-face symbols denote vector or higher-order tensor quantities. The Einstein summation convention is used for repeated lower case latin subscripts within a term, as in .) shown in Figure 1C(bottom), the displacements , strains , and stresses for the elasto-static problem arewhere is the position vector with Cartesian components and orthonormal base vectors (). are the constant components of the 4th-order elastic stiffness tensor . The quantities and are position-dependent components of 3rd-order and 4th-order tensors, respectively, as originally defined in (Eshelby, 1957, 1959).

Moreover, according to the theory of linear viscoelasticity, the constitutive (stress-strain) response of the material can be expressed as a convolution integral in terms of a relaxation modulus and a prescribed strain history. That is, for a given shear strain history , the shear stress history is expressed in terms of the relaxation shear modulus aswhere is the strain rate. To account for the viscosity of the cell and gel in the framework of the Eshelby inclusion analysis, all mechanical field quantities become time-dependent. Because we are only interested in the steady state (long-term) mechanical response of the cell contraction, we can exploit the correspondence principle of linear viscoelasticity (Wineman and Rajagopal, 2000) to solve the original time-dependent equations in the frequency domain and substantially reduce the computational cost. The solutions are then transformed back into the time-domain to get the cyclic time-history of the field quantities. Taking the Fourier transform of the field equations in (Equation 3a), (Equation 3b), (Equation 3c), (Equation 3d), the displacements, strains, and stresses can be obtained in frequency domain by

An arbitrary periodic eigenstrain history can be expressed as the complex Fourier series.where the coefficients are complex-valued . Decomposing the load-free strain history into n = 128 discrete frequencies, the power spectrum of is shown in Figure 1E on a logarithmic scale. In addition, once the stress in the inclusion is known, the mechanical power output of the cell (inclusion) can be calculated. This is related to the product of the instantaneous cell stress and strain rate , asand the work done per cycle (period T) is

The volume of the ellipsoidal inclusion is (see Figure 1C(bottom)), and is the start time of a cycle.

Autoregulation model

The developed VE Eshelby model assumes the cardiomyocyte to be a contractile but otherwise passive (MCT-free) inclusion, and it analyzes the in-gel cell response based on the load-free cell strain history (eigenstrain). Based on experimental data, we hypothesize that the MCT mechanism effectively alters the eigenstrain history, , in response to the mechanical afterload on the cell. That is, the real cell senses the increased load and increases the transient to enhance contractility. The MCT- gain is clearly observed in our experiments, but the exact MCT mechanism is incompletely understood. A well-calibrated mechanical analysis should provide clues to inform the ongoing investigation of MCT mechanisms and functional consequences in cardiomyocytes.

Our approach to mathematically capture the cardiomyocyte’s autoregulatory response is depicted in the block diagram of Figure 1F. The complete feedback loop requires two additional mathematical “links” (bold green boxes) to the VE Eshelby analysis. The first link is the known enhancement of contractility because of the amplified transient ( -Contractility Coupling), and the second is the gain because of the mechano-chemo-transduction (Stress- MCT). Each is described below.

Ca2+-Contractility Coupling

This first link is described by the following first-order kinetic model that connects the cytosolic concentration to the eigenstrain β (contractility),where is a time constant, is a functional of the transient , and the initial condition is at the start of the cycle. The right-hand side of Equation (8) is a “forcing function” (f), chosen to capture the system dynamics. Figures 1G and 1H show measured time histories of the strain and respectively, comparing load-free and in-gel responses. Note that the experiment in Figure 1A began with the cell contracting in the cross-linked gel and then the gel was dissolved, allowing the cell to contract in load-free. It should be emphasized that the plots show representative cycles for each case (in-gel and load-free) measured on the same cell (self-control experiment), and this provides a direct comparison that avoids complications from cell-to-cell variability.

Figure 1G shows that the cell actually achieves a slightly larger contraction in-gel than when load-free. Owing to the presence of the viscoelastic gel, however, the cell is unable to fully recover its deformation and a residual strain of almost –4 % remains at the end of each cycle. The in-gel transient exhibits an MCT- gain () compared to the load-free case in Figure 1H, indicating a significant upregulation of the intracellular transient, as measured using the standard Fura-2 fluorescence ratiometric method (Jian et al., 2014). Cross-plots of strain versus for the LF and in-gel cases are provided in Figure 1I, showing a much larger hysteresis for the in-gel case.

From the observed behaviors in Figures 1G–1I, the following nonlinear function is proposed for f,where is the average of the transient during one actuation cycle and are positive constants. The magnitude of f in Equation (9a) is essentially proportional to when , yet when because of the exponential in the denominator becoming large. The resulting response is very nearly periodic, , because we choose a small value for () and during the last half of the time interval.

Overall, the coupled eigenstrain- model in Equations 8, 9a, and 9b has four parameters . Identification of these parameters was done recursively using a numerical nonlinear least-squares minimization algorithm. The model parameters were first calibrated to the measured load-free axial strain and histories, giving the values in Table 2. The resulting model was then tested using the measured in-gel history to predict the upregulated eigenstrain.

Table 2

Eigenstrain- parameter values

a_1 (s)	b_0	b_1	b_2
0.086	0.044	3.226	19.735

Figure 1J demonstrates an excellent fit of the model’s MCT-free eigenstrain, (dashed blue line), against the experimentally measured load-free axial strain, Exp  (solid blue line). The load-free eigenstrain is the calculated response of the model, Equations 8, 9a, and 9b, to the load-free transient in Figure 1H. The upregulated eigenstrain of the model, , in Figure 1J is the calculated response of the model to the in-gel transient of Figure 1H. This accounts for mechano-chemo-transduction, which produces a peak MCT eigenstrain of compared to the peak load-free (MCT-free) eigenstrain of . It also produces faster upstroke of the transient, with the peak occurring about 0.2 s sooner. The contractility of the cell is enhanced in both its amplitude and rate.

As a baseline study on the VE Eshelby inclusion model integrated with the proposed -Contractility coupling dynamics, a single cardiomyocyte was exercised in a self-control experiment and the developed VE Eshelby inclusion model was simulated using both as-measured load-free eigenstrain () and predicted MCT eigenstrain () to capture the influence of upregulation on the cell deformation and stress. The MCT-free simulated results (obtained via ) are shown as gray lines for reference. Figures 2A and 2B show the simulated upregulated () axial strain and stress histories during a periodic cycle. Peak strains and stresses are marked by unfilled circles.

Figure 2

Modeling results for the baseline case using the analytical VE Eshelby simulation, with and without MCT autoregulation feedback

(A–H) (A) strain history, (B) stress history, (C) mechanical power output by the cell, (D) mechanical work done by the cell on the surrounding, (E and F) comparison of the measured gain to the in-gel MCT simulation calculated from (Equation 10a), (Equation 10b), and (G and H) simulated MCT cell strain and stress responses (in-gel) including cell viscoelasticity, showing good agreement with the measured cell strain response.

In the MCT-free analysis, the simulated and measured strain timing and amplitude in-gel show considerable discrepancies. The agreement is much better between the measured and simulated strain histories in Figure 2A when MCT effects are taken into account using . The magnitude of the simulated peak contraction is much closer and the time lag discrepancy is only 0.04 s. Importantly, the maximum predicted axial stress resulting from the MCT eigenstrain is almost twice the MCT-free case in Figure 2B.

The effect of MCT upregulation is even more significant on mechanical energy. A comparison of the mechanical power and work histories between the MCT and MCT-free simulations are plotted in Figures 2C and 2D. The increases in mechanical power and work are quite large, because the strain magnitude, strain rate, and stress magnitudes are all enhanced in the MCT model. The calculated mechanical work delivered by the upregulated cell at the end of a cycle is pJ, which is about larger than the MCT-free cell.

Stress- Ca2+ MCT

To complete the feedback loop, this second link relates the sensed afterload to the increase in Ca2+ transient. As mentioned earlier, we do not yet know the molecular identity or location of cell mechanosensors. For simplicity, our phenomenological model assumes that MCT- gain is a function of the cell’s internal axial stress , because this is the dominant stress component and the stress is uniform throughout the ellipsoid. Accordingly, the input of the proposed model is the time-dependent MCT-upregulated axial stress, , and the output of the model is the increase in Ca2+ transient, according to,where is the difference (MCT- gain) between the in-gel and load-free cytosolic in Figure 1G at every moment, is the minimum value of (maximum compressive axial stress) within each beating cycle, and , , , , , are five constant parameters. The time ranges between zero and the period of each cycle ( 2.096 s here). The calibrated parameters for the current cell, identified recursively are listed in Table 3. The quantity in brackets on the right-hand side is always positive because of the addition of to the cell stress. Note also the explicit time dependence () in the denominator of g grows large as time progresses and quickly drives the right-hand forcing function to zero.

Table 3

gain-stress parameter values

c_1 (s)	d_0 (kPa^−1)	d_1 (kPa)	d_2 (kPa)	d_3 (kPa^−3s^−4))
0.199	0.103	8.584	1.063	0.126

Figures 2E and 2F compares the Ca2+ gain of the model with experimentally measured data. The predicted time-dependent Ca2+ gain ()—by a numerical solution of Equation (Equation 10a), (Equation 10b)—matches well with the measured history (Figure 2E). It is also worth mentioning that the MCT-upregulated axial stress inside the cell (Figure 2B) leads the Ca2+ gain time-response by about 0.06 s. The resulting Ca2+ gain versus the (simulated) cell axial stress in Figure 2E exhibits an oddly hysteretic response, but the important segment in the response appears to be the rising, nearly linear, portion (see the arrow).

With these two proposed links established, the MCT process forms a closed-loop where the mechanical force is sensed by the cell’s mechanosensors and the cytosolic Ca2+ transient is increased by . In turn, the upregulated Ca2+ transient magnifies the eigenstrain (β-MCT) and enhances the contractility of the cell.

Cell viscoelasticity

Until now, the cell’s passive mechanical properties have been modeled as purely elastic, but even better agreement with the experimental data can be achieved by incorporating the viscosity of the cell itself (Figures 2G and 2H). This final step is of secondary importance (compared to the MCT effect) and is somewhat speculative, because the viscosity of the cell is unknown. In the absence of better information, we chose the simplest possible viscoelastic constitutive model, a three-parameter viscoelastic solid (Equation (1) with ) with kPa, kPa, and s. The same from Figure 2A is used as the input to the simulation. As seen in Figure 2G, this produces a remarkably good match with the measured strain response of the cell (see ), such that the calculated peak strain of is quite close to the measured value of and the time lag discrepancy is almost eliminated. The fact that the simplest possible viscoelastic model generates such good agreement provides us confidence that the assumed cell’s viscosity parameter values are reasonable. Accordingly, Figure 2H provides our best estimate of the MCT cell’s axial stress history during the in-gel portion of the Cell-in-Gel experiment. The peak stress of kPa occurs at s, and the minimum stress (maximum compression) of kPa occurs at s.

Numerical model: Finite-element analysis

The 3D mechanical field is also investigated numerically for the actual irregular geometries of left ventricular cardiomyocytes. The cardiomyocyte’s 3D geometry is obtained by surface reconstruction of the cell’s outer envelope using a stack of 2D images captured through the thickness of the cell. The 2D confocal images were acquired using standard optical sectioning techniques on the Olympus FluoView FV1000 confocal microscope with 0.5 μm resolution in the Z-direction (through thickness). The outlines of the cell exterior were segmented using 3D-Doctor (Able Software Corp.) and the 3D rendered surface was smoothed at multiple layers to avoid unnecessary geometrical and computational complexities.

Because an analytical solution is not available to such a complicated shape, an FEA was performed in COMSOL Multiphysics (COMSOL Multiphysics, 2020) using tetrahedral elements with quadratic shape functions. A fine unstructured computational mesh was constructed for the cell domain and a gradually growing mesh (toward the outer boundaries) was developed for the gel domain. Because we are only interested in the steady-state cyclic response of the cell, the FEA was performed in the frequency domain and the results were transformed back into the time-domain using the inverse DFT. A mesh convergence study was also performed to measure the sensitivity of the response with respect to the size of the computational grid and gel’s domain. In the finite-element model, the cardiomyocyte geometry was encompassed by a cuboid gel domain that extended about 5 times the dimensions of the cell in each direction and had viscoelastic properties similar to the hydrogel. Even with the improvement from the DFT, the FEA computations were intensive compared to the VE Eshelby calculations. Each FEA simulation took about 5 hours on an Intel Xenon quad-core desktop computer, whereas a VE Eshelby simulation took only a few seconds in Wolfram Mathematica (Wolfram Research, Inc., 0000) on a laptop computer.

Modeling results

FEA results of a baseline case

Here the finite-element (FE) mechanical modeling results of the single-cell contraction inside a viscoelastic hydrogel are presented. FE simulations are performed on a few actual geometries of typical young rabbit heart’s left ventricular cardiomyocytes and a representative cell geometry is selected for demonstration and detailed analysis. Although not shown here, the FEA was also performed on a few other irregular muscle cell geometries from different regions of the left ventricular wall and very similar conclusions and findings were drawn. The FE modeling is performed in the frequency domain with a frequency range of 0–30.057 Hz and frequency increments of 0.477 Hz. At each frequency, a quasi-static linear viscoelastic solution is stored and the augmented results are then transformed into the time-domain using inverse DFT. For comparison, the VE Eshelby model was also used to calculate the multiaxial stress and strain time responses of the corresponding ellipsoid-shape approximation of the cardiomyocyte with identical length, surface area, and volume as the actual geometry. The cell contraction was calculated using the upregulated eigenstrain () in both modeling approaches. The cell was treated as a contractile, elastic object with a shear modulus of kPa contracting inside a viscoelastic hydrogel with kPa and a crosslink density of %.

As a baseline case of actual cardiomyocyte geometry, the mechanical analysis of a typical rabbit ventricular myocyte as shown in Figures 3A and 3B, is studied here. The FE simulation was carried out on a 64-bit operating system with Intel Xenon E−2136 processor (@ 3.30 GHz) and 32 GB memory. For the baseline FEA model, about 290k full-integration 3D continuum elements were used and the solution time was about 5.5 h.

Figure 3

Simulation results of the numerical model obtained by FEA of the baseline case

For a Figure360 author presentation of this figure, see https://doi.org/10.1016/j.isci.2022.104667.

(A–I) (A) 3D geometry of the cardiomyocyte used in FEA, (B) confocal micrograph of the cardiomyocyte, (C and D) axial stress and strain fields at different cross sections along the cell’s axis from FEA of the cardiomyocyte contraction in 8% gel, (E and F) average axial strain and stress histories inside the cell and their comparison against the VE Eshelby simulation results for an equivalent ellipsoid with 8.75 and 4.65, (G) traction vector on cardiyomycyte’s surface, (H) distribution of normal surface traction, (I) distribution of shear surface traction, and (J) time history of the magnitude of force components summed over the cell surface. See also Video S1.

Figures 3C and 3D show the contour plots of the axial stress and strain at several cross-sectional slices along the longitudinal axis of the cell. As expected, the strain and stress fields within the inclusion are no longer uniform, now gently varying mostly along the longitudinal (x) axis. The contour plots are shown at the time of maximum contraction ( s) where the average axial strain is about % and the average axial stress is 1.49 kPa. It can be seen that the axial strain and stress are closer to their corresponding average values near the mid-span of the cell where the axial cross section is almost constant with x. Moving toward the two ends of the cell, a higher gradient is observed in values of both axial strain and stress within each cross section. This is because of the fact that a typical cardiomyocyte has jagged narrow ends.

The time history of the average axial strain and stress values are also shown in Figures 3E and 3F using the red dashed lines. The green shaded bands in Figures 3E and 3F represent the 90% cumulative distribution of each corresponding quantity. It is observed that as the contraction increases, the distribution of strain and stress values become broader. In addition, the black solid lines depict the analytically obtained histories from the VE Eshelby inclusion problem using the ellipsoidal approximation. For this cell, the ellipsoidal inclusion dimensions are μm, μm, and μm. The cell’s aspect ratios in the equivalent ellipsoidal approximation are 8.75 and 4.65. It can be seen that the results of the VE Eshelby inclusion problem agree quite well with the FEA average values of axial strain and stress at each instant of time. The predicted axial strain and stress histories from the VE Eshelby inclusion model agree within 8% and 4%, respectively, of the corresponding average histories from the FEA results. Specifically, the maximum contraction occurs at s at % strain from the FEA and it occurs at s at % strain in the VE Eshelby analysis. Moreover, the peak cell’s average axial stress of 3.13 kPa occurs at s in both simulations.

To give a sense of the magnitude and direction of the traction on the surface of the cardiomyocyte, a scaled schematic of the traction vector distribution on the entire cell’s surface is provided in Figure 3G. The traction distribution is taken at the moment of cell’s peak average axial stress ( s). Figures 3H and 3I show the spatial distribution of the magnitude of the normal and shear components of the traction vector, respectively, on the surface of the cell. Figure 3H shows that the magnitude of the normal traction is near zero for most of the surface except for a small region near the end of the cell where it is maximum. Figure 3I shows that the shear traction is small near the mid-span of the cell and gradually grows and reaches its maximum value near the two ends of the cell. It can be concluded that surface tractions are mostly governed by the normal stress near the two ends and by shear stress almost anywhere else. Such a nonuniform distribution of surface tractions may provide clues regarding the spatial distribution and sensitivity of mechanosensors along the surface of cardiac muscle cells (Video S1). Moreover, at each instant of time, the sum of the magnitudes of the surface shear component of traction is larger than the normal component as shown in Figure 3J. This might suggest a higher sensitivity of surface mechanosensors to shear stresses and a potential role shear tractions play on MCT.

Figure 4 shows the snapshots, within the plane in the cell and gel, of several mechanical field quantities, including the total displacement, longitudinal strain, longitudinal stress, shear strain, and shear stress at the time of maximum contraction. For comparison, the same field quantities are also evaluated analytically using the VE Eshelby inclusion model of the equivalent ellipsoid. The mechanical field quantities in the equivalent ellipsoid model, although always symmetric with respect to the plane (mid-span), agree reasonably well with the corresponding FEA results of the actual geometry.

Figure 4

Comparing the simulation results using numerical and analytic models for the baseline case

(A–H) Snapshots of the mechanical field quantities at the time of maximum contraction on plane and their comparison with the VE Eshelby inclusion results for the equivalent ellipsoidal approximation of the cardiomyocyte. (A and B) total displacement, (C and D) axial strain, (E and F) axial stress, (G and H) shear strain, and (I-J) shear stress field. See also Video S2 and S3.

Figures 4A and 4B present the contour plot of the magnitude of displacement field normalized by the cell’s half-length, a. Streamlines are overlaid on the equivalent ellipsoidal (Figure 4B) to show the directions of displacements during contraction, showing how the gel is drawn inward toward the ends of the cell at yet is pushed away from the cell at the mid-span . The contour plots also show that nonzero displacements in the gel are localized to the vicinity of the cell, and rapidly approach zero by or so. Thus, the map provides useful information about the expected displacements in the gel, which could be compared to the displacements of the micro-beads in the Cell-in-Gel experiments (Figure 1).

Contour plots of the axial strain () and axial stress () are provided in Figures 4C–4F, respectively. A clear difference between the two analyses is that the stress and strain fields in the ellipsoidal geometry approximation are uniform at any given time instant but are less uniform in the actual geometry with gradients near the corners of the cell (Videos S2 and S3). In both approaches, the fields in the gel (matrix) are nonuniform with a local strain and stress concentration just outside the cell’s ends, but these decay quickly to zero away from the cell. In the gel, both plots exhibit “hot spots” near the cell ends and at the mid-span (although much less severe). The maximum stresses and strains in the gel are quite local, less than about to an extent near the cell ends.

Figures 4G and 4H show the contour plots of the shear strain in the (). Owing to the diagonal structure of (only normal eigenstrains) and the axisymmetry of the approximated ellipsoidal geometry, the shear strain inside the inclusion is equal to zero. Moreover, despite the nonuniformity of the shear strain inside the actual cell geometry, its value is quite small everywhere inside the cell. The shear strain is, however, much larger just outside the cell (in the gel) with values growing with increasing . Near the two ends, the shear strain tends to zero again right outside the cell. In addition, the shear strain field is more pronounced and extends further outside the actual cell geometry. Similar behavior is observed in the shear stress field, as shown in Figures 4I and 4J. The shear stress in the ellipsoidal approximation is identically zero, whereas it is nonzero in the actual geometry with high gradients near the corners of the cell. The shear stress field inside the gel is quite similar for both actual and ellipsoidal cardiomyocyte geometries.

The aforementioned studies show that the analytical model solutions (assuming an ellipsoid-shaped cell) can accurately capture the cell-averaged values of the mechanical quantities predicted by the FEA model of real cardiomyocytes (rod-shaped with irregular geometry). This allowed us to develop the autoregulation model in Equations (Equation 8), (Equation 9a), (Equation 9b), (Equation 10b), (Equation 10a) by using the analytical model as a fast approximation of the average mechanical behavior of cardiomyocytes in the Cell-in-Gel experiments.

Autoregulation of contractility by MCT feedback

We conducted experiments and performed model simulations to test the hypothesis that MCT feedback enables cardiomyocytes to autoregulate transient and contractility to maintain contraction amplitude despite the increase in the afterload (Figure 5). Systematic experiments were performed to study the increasing mechanical load (by mixing 10% PVA with a series of crosslinker concentrations, CL%, to increase hydrogel density) in relation to transient and contraction (Shimkunas et al., 2021). Figure 5B shows that the transient amplitude increases with increasing load until reaching a peak (in 0% 10% range), and then declines under higher load (10% 15%). We can see in Figure 5C that with increased transient, cardiomyocytes are able to maintain a relatively stable contraction amplitude despite the increasing load (0% 10%) but fail under higher loads (10% 15%). These experimental data are recapitulated in our model simulations. Our autoregulation model predicts the peak transient, contraction, and stress by the autoregulation dynamics of Equations (Equation 8), (Equation 9a), (Equation 9b), (Equation 10a), for crosslinker concentrations, C, ranging between 0% 15% (Figures 5D–5G). For this purpose, typical values of aspect ratio ( 5), load-free peak (), and peak load-free contraction ( %) were chosen.

Figure 5

Autoregulation of contractility: model predictions and experimental verification

(A) Autoregulation concept for excitation-contraction coupling by MCT feedback. Experimental data from rabbit ventricular myocytes are shown in (B) and (C). Cells were embedded in the gels of different stiffness by mixing 10% PVA with a crosslinker of indicated concentrations (CL%).

(B) shows the transient peak (indexed by Fura-2 ratio) in relation to gel stiffness.

(C) represents the contraction amplitude (fractional shortening %) in relation to gel stiffness. Statistical tests: The bars show the mean and SEM of each group with indicated number of cells/animals. All groups passed the D’Agostino-Pearson normality test. One-way ANOVA test is used for multiple groups comparison, and Student’s t test for pairwise comparison of neighboring groups.

(D–G) Model prediction of (D) transient peak, (E) contraction amplitude, (F) axial stress history, and (G) peak axial stress as a function of crosslinker concentration, C %.

A comparison of model predictions (Figures 5D and 5E) and experimental data (Figures 5B and 5C) shows that the simulated peak and peak contraction amplitude agree reasonably well with experimentally measured data. Specifically, a similar bell-shaped relationship between the peak transient and hydrogel stiffness is observed in model predictions (Figure 5D) with increasing crosslink concentration (C). From analytical modeling, the peak of is predicted to occur at % where the peak contraction amplitude is %. In addition, the model prediction falls mostly within one SD of the mean measured data. Modeling predictions in Figure 5E suggest that autoregulation of contractility occurs in cardiomyocytes contracting in hydrogels with crosslink concentration up to % where cells are able to maintain contraction amplitudes. Beyond this, the MCT- gain and cell contraction decrease progressively as the gel resistance overwhelms the cell contraction under excessive afterload.

Moreover, in our previous work (Kazemi-Lari et al., 2021), we performed a parametric study to correlate the hydrogel stiffness (crosslinker concentration, C%) with the axial stress but without the MCT effect. Now, we have added MCT- gain in the autoregulation model and correlated the hydrogel stiffness with the peak axial stress during the contraction cycle (Figure 5F). Our autoregulation model simulation predicts higher internal stress in the cardiomyocytes in-gel as a result of autoregulation of contractility. Figure 5G shows that the peak axial stress increases as the gel stiffness increases and reaches a maximum of 3.14 kPa at around % beyond which it asymptotically decreases to an approximate value of 2.5 kPa. This trend is consistent with the trends observed in the model predictions of the peak transient and contraction amplitude. Therefore, model simulation empowers the calculation of the 3D mechanical stress in cardiomyocytes under various load conditions, which cannot yet be measured experimentally.

Discussion

The classical paradigm of cardiac E-C coupling emphasizes feedforward control: the action potential activates the signaling which in turn activates the muscle contraction. Based on recent experiment data revealing MCT mechanisms (Shimkunas et al., 2021; Hegyi et al., 2021), the mathematical models developed herein delineate how closing the MCT-feedback loop in E-C coupling enables autoregulation of contractility in cardiomyocytes.

To put our studies into the context of literature, most previous experiments in the cardiac E-C coupling field used cardiomyocytes bathed in aqueous solutions under load-free condition. Those experimental studies significantly contributed to the current understanding of E-C coupling, especially in the areas of cardiac electrophysiology and signaling. Nevertheless, the mechanical load effects on cardiomyocytes were missing in those load-free experiments.

In efforts to apply mechanical load on cardiomyocytes, several techniques have been developed in the field. The usage and limitations of various techniques for applying mechanical loading at the single-cell level have been analyzed in a review paper (Chen-Izu and Izu, 2017). One technique is to use a pair of micro-cantilevers (i.e., carbon fibers, glass fibers, etc.) attached to the cell surface to stretch the rod-shaped cardiomyocyte along the long axis (Le Guennec et al., 1990; Alvarez et al., 1999; Calaghan and White, 2004; Iribe et al., 2009; Prosser et al., 2011). Such a one-dimensional stretching method applies longitudinal tension on the cell but no transverse compression or surface traction on lateral surfaces. Another method is to attach cells to a membrane and apply 2D stretch. However, the 2D-stretching method is limited to neonatal cardiomyocytes or stem-cells that can be cultured and attached to the membrane, because the rod-shaped adult cardiomyocytes do not attach well. It is important to know that different experimental settings apply different types of mechanical stress on the cell, and these may activate different mechanosensors and MCT pathways.

We invented a 3D Cell-in-Gel system where freshly isolated ventricular myocytes were embedded in a viscoelastic hydrogel polymer matrix (Jian et al., 2014). We carried out experiments to measure the intracellular transient and contraction while the cardiomyocytes were paced to undergo E-C coupling under mechanical load in the hydrogel. Next, the hydrogel was dissolved to release the cell into Tyrode’s solution and the same cell was contracting load-free in the solution. Such self-control experiments allow us to compare the transient and contraction of the same cell under mechanical load versus load-free, which provide data for model parameter extraction. To further investigate autoregulation phenomena, we systematically tuned the hydrogel stiffness to apply different afterloads on cells. Experimental data show that cardiomyocytes were able to increase the transient in accordance with increased gel resistance to maintain the contraction amplitude, suggesting autoregulation of contractility at the single-cell level.

To develop mathematical models for analyzing the force and deformation fields in cardiomyocytes during contraction in a 3D viscoelastic environment, we started from the elastic Eshelby inclusion model (Shaw et al., 2013) by approximating the cell as an elastic ellipsoidal inclusion. Later, we developed a viscoelastic generalization of the Eshelby inclusion analysis (Kazemi-Lari et al., 2021), where the model produces an exact solution of the time-dependent mechanical fields during beat-to-beat contraction of the cardiomyocyte in the viscoelastic environment of the surrounding hydrogel matrix.

Here, we further exploited the analytical viscoelastic model to develop a numerical model using FEA to account for the actual geometries of real cardiomyocytes. The goal is to provide accurate mechanical analyses of the cell contraction in a 3D viscoelastic environment. We performed modeling using analytical and numerical methods and the two approaches showed good agreements. In developing the model, we initially assumed the cardiomyocyte to be a contractile but otherwise passive (MCT-free) object. However, the experimental data from cell-in-gel experiments show that the real cardiomyocytes can adapt to load increases by upregulating the cytosolic transient to enhance contractility. This is because of the fact that a living cardiomyocyte is not a passive object but has an MCT mechanism to actively regulate its contraction in response to load changes.

We conducted mechanical analyses of different load conditions and used the autoregulation model to predict the mechanical stresses inside the cell and on the cell surface (see Videos S1, S2, and S3). Model simulations show that when the cardiomyocyte is contracting in an aqueous solution, the axial strain magnitude (same as the fractional shortening) is large, but the stress inside the cell is minimal, and the surface traction is also negligible. In contrast, when the cardiomyocyte is contracting in a hydrogel, mechanical analyses show a knockdown of contraction amplitude, if the cell is a passive viscoelastic object without any MCT regulation. This hypothetical MCT-free case serves as a comparison to understand discrepancies between this assumption and the experimental measurements from real cardiomyocytes. The experimental data clearly show upregulations of transient and contraction, revealing the necessity to add MCT into the model to describe the living cardiomyocytes that can actively regulate cell function. Based on experimental data, we developed an autoregulation model by incorporating a load-dependent MCT gain that leads to enhanced MCT-eigenstrain. The eigenstrain history was made dependent on the measured transients. The model parameters were obtained by fitting the experimental data of transients and load-free strain. Then, the model was used to predict the MCT-eigenstrain history based on in-gel transients. Using this upregulated eigenstrain, the VE Eshelby model accurately predicted the enhanced contractility of the cell in response to afterload. Furthermore, we used the model to simulate the transient and contraction amplitude in cardiomyocytes embedded in hydrogels of different stiffness. The simulation results predicted that the transient should increase in accordance with increasing gel stiffness, which should enhance contractility to maintain the contraction amplitude in a range of mechanical loads, but exceedingly high loads would knock down the contraction amplitude and decrease the transient. These non-trivial model predictions have been confirmed by experimental data.

Our mathematical models provide new computational tools to calculate the 3D mechanical fields in cells and study how cardiomyocytes autoregulate E-C coupling and contraction under mechanical loads. Simulation results show that the cardiomyocyte contracting in a viscoelastic environment should experience 3D strain and stress including (to various amounts) longitudinal tension, transverse compression, and shear and normal surface stresses. Parametric analyses show how the mechanical stresses are related to physiological and pathological conditions such as the load level, the cell aspect ratio and shape, the stiffness of the cell and of the environment.

The mechanical stress increases when the cell is contracting under higher resistance or after overload induced damage/remodeling, which would be associated with pathological conditions such as hypertension (Tello et al., 2019), pressure-overload (Jashari et al., 2015), and heart failure with reduced ejection fraction (HFrEF) (Chandar et al., 2010). The stress also increases when the cell is in a stiffer environment, which would be associated with fibrosis (Naser et al., 2021; Neff and Bradshaw, 2021; Frangogiannis, 2021), infarction (Villalba-Orero et al., 2021), and heart failure with preserved ejection fraction (HFpEF) (Zile et al., 2015). Interestingly, the stress increases with the cell’s aspect ratio (length to width ratio). High mechanical stress, among many other factors, may contribute to the high prevalence of atrial fibrillations (Virani et al., 2021). Furthermore, model analyses show that cardiomyocytes under higher stress also need more energy to generate the work output during contraction, which would increase the demand for ATP synthesis in mitochondria. However, mitochondria function is impaired in heart failure (Garbincius and Elrod, 2021), which would further exacerbate cardiac dysfunction. Therefore, our models provide useful tools to quantify mechanical stresses and energetic demands in cardiomyocytes and study how impaired autoregulation leads to cardiac dysfunction and heart diseases.

Limitations of the study

The current autoregulation model also has limitations. Here we first focus on closing the MCT feedback loop in E-C coupling (Figure 1F) to start the effort to build an autoregulation model. Based on existing experimental data, we first develop a phenomenological model to capture the autoregulation phenomena. However, this model does not yet include detailed MCT mechanisms which need to be addressed in future model development. For example, a mechanistic model should include MCT signaling pathways from mechano-sensors (i.e., dystroglycans, integrins, etc.) to chemo-transducers (i.e., NOSs, Noxs, CaMKII, etc.) to effector molecules (i.e., ion channels, handling proteins, contractile proteins, etc.), as well as MCT induced post-translational modification of E-C coupling proteins that underlie the autoregulation phenomena. A more complete model should also incorporate other feedback loops including the sarcomere level feedback and mechano-electro-transduction. However, experimental studies on these feedback mechanisms and MCT pathways still need to be done to provide necessary data for model development. In fact, our Cell-in-Gel system is designed to embed cardiomyocytes in a hydrogel that enables applying mechanical load on the internal and cell-surface mechanosensors (through surface tractions). We will also use the FEA model developed here to quantify the 3D stains and stresses in cardiomyocyte’s architecture, which is essential for identifying putative mechanosensors inside the cell and on the cell surface. Our FEA model at the single-cell level will also pave the foundation for developing tissue-level and organ-level FEA models in the future. Currently, both experimental and modeling studies of mechanical load effects on cardiomyocytes and mechano-transduction mechanisms are still at a pioneering early stage. Our studies aim to fill critical knowledge gaps and build foundations for understanding how the heart autoregulates contractility in response to mechanical load changes to maintain cardiac output in health and diseases.

Perspective

Biological cells, as the basic unit of life and function, are controlled by dynamic systems that integrate a myriad of biochemical processes and molecular interactions. Deciphering the control mechanisms at the single-cell level is essential for understanding how the cell can sense its environment and regulate the biochemical processes to achieve cell function. The cardiomyocyte presents an important and informative case study for understanding how nonlinear dynamic systems integrate to control cardiac excitation-contraction and autoregulate the contractile function in connection with the mechanical load changes to maintain cardiac output. Given that the bioelectrical system, the signaling system, and the contractile/motility system are ubiquitous in biological cells, our work in understanding how mechanical stresses affect these systems in cardiomyocytes may have broader implications for understanding the mechano-transduction in other cell types.

STAR★Methods

Key resources table

Resource availability

Lead contact

Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contact, Ye Chen-Izu (ychenizu@ucdavis.edu).

Materials availability

This study focuses on mathematical modeling and did not generate new physical materials. All software listed in the key resources table can be obtained commercially.

Method details

All computational model development and simulation details are described in this paper. Further details of the previous development of VE Eshelby inclusion model can be found in Shaw et al. (2013); Kazemi-Lari et al. (2021).

This modeling paper uses our previously published experimental data which are cited in the paper. We summarize the essential methods used in the experiments.

Animal and cells

New Zealand White rabbits, male, 4-6 months of age, were purchased from Charles River Laboratories (Wilmington, MA). Standard enzymatic techniques were used to isolate the left ventricular myocytes (Horvath et al., 2013; Hegyi et al., 2018). All experiments were performed at room temperature and within 10 hours after cell isolation.

Cell-in-Gel system

We used our Cell-in-Gel system to embed freshly isolated rabbit ventricular myocytes in a 3-dimensional (3-D) hydrogel made of polyvinyl alcohol (PVA, 89-98 kDa, hydrolyzed, Sigma Aldrich) and 4-armed boronate-polyethylene glycol (4B-PEG) crosslinker (Jian et al., 2014).

Gel-forming protocol: The cell suspension was first mixed with PVA solution (10% wt.) at 3:7 ratio. The PVA-cell suspension was mixed with the crosslinker, which embeds cells in the viscoelastic hydrogel within minutes. The 4B-PEG crosslinker also spontaneously binds to the cis-diols of glycosylated molecules on the cell surface, thereby tethering the cell to the gel polymer matrix (Jian et al., 2014).

Gel-dissolve protocol: Sorbitol (1% w/v) was added into the BTy solution to compete off the crosslinker binding to PVA, and thereby dissolving the gel. Note that it takes about 15-20 minutes to dissolve the gel and refresh the bath solution.

Normalization protocol: We performed paired measurements in the same cell to normalize measurements obtained in gel (after load) to measurements obtained in solution (load free). These self-control experiments using the gel-forming and gel-dissolve protocols helped to minimize potential errors due to cell-to-cell variability in dye loading and use of Flou-5N which is a single wavelength fluorescent indicator.

Tuning the gel stiffness: The gel stiffness is tunable by adjusting the crosslinker to the PVA ratio. We used the cell suspension containing 10% PVA to mix with the solution containing 10% crosslinker at equal volume (1:1 ratio) for experiments unless noted otherwise. To apply different levels of mechanical load on the cell, we systematically tuned the hydrogel stiffness ranging from soft to hard gel using a series of crosslinker concentrations of 5%, 7.5%, 9%, 10%, 11%, 12.5% and 15%, each mixing with 10% PVA solution at 1:1 ratio.

Measurement of cytosolic transient

Cytosolic concentration, , was measured using the Fura-2 ratiometric method as previously described (Jian et al., 2014). The IonOptix Calcium & Contractility System (IonOptix, USA) was mounted on an Olympus X71 inverted microscope. Cells were perfused continuously with BTy solution and electrically stimulated at 0.5 Hz. The cells were imaged using wide-field epi-fluorescence microscopy with 40x objective (UPlanSApo 40X, water-immersion, NA 1.15, with correction for the No.1 glass coverslip thickness). Excitation light was generated with a xenon arc lamp and passed through a HyperSwitch with a galvanometer-driven mirror switching at 500 Hz. Fura-2 was excited by dual excitation lights at 340 nm and 380 nm wave lengths; the fluorescence emission passed through a D510/40 nm bandpass filter and was captured in a photomultiplier tube (PMT). The whole-cell fluorescence signal was acquired from the ROI encompassing a single cell in its entirety. Photobleaching of the dye during experiments was monitored, best-fit to a linear function, and corrected in the data analyses. Background fluorescence was measured from an area adjacent to the cell, which was then subtracted from the total fluorescence signal measured in each experiment. The ratio of Fura-2 emission signals from 340 nm and 380 nm was calculated as a measure of .

Measurement of cardiomyocyte contraction

Cardiomyocyte contraction was measured using the IonOptix Calcium & Contractility System with a high-speed camera (MyoCam-S, 1000 frames/sec), concurrent with the fluorescence measurements of signals. Cells were continuously perfused with BTy solution and electrically stimulated at 0.5 Hz. The cardiomyocyte’s sarcomere length (SL) was measured from the striation pattern of the cell using the real-time fast Fourier transform (FFT). The cell contraction was measured by the fractional shortening of SL calculated as the percentage of the SL shortening to the diastolic SL.

Quantification and statistical analysis

Data sets from each experimental group were analyzed to calculate the mean value (Ave.), standard deviation (SD), sample number (n), and standard error of mean (SEM). Statistical comparison between data sets of two different groups used Student’s t-test for either unpaired or for paired measurements in the case of self-control experiment. One-way ANOVA with Bonferroni post hoc test was used to compare multiple groups. Statistical significance of the difference between experimental groups was tested using the standard 95% confidence level.

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Software and algorithms
WOLFRAM MATHEMATICA, Version: 11.1.1.0	Wolfram	https://www.wolfram.com/mathematica
COMSOL Multiphysics®, Version: 5.5	Comsol	https://www.comsol.com
3D-DOCTOR, Version: Trial	Able Software Corp.	https://www.ablesw.com/3d-doctor
Deposited data
Measurements of Ca^2+ transient and contraction	This paper	DRYAD: https://dx.doi.org/10.25338/B8T935